VtkWriter.hh 11.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13 14 15 16 17 18
#include <list>
#include <vector>
#include "DOFVector.h"

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
19 20
  using namespace std;

Thomas Witkowski's avatar
Thomas Witkowski committed
21
  template<typename T>
22
  void VtkWriter::writeFileToStream(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
23 24 25 26 27
  {
    int nVertices = (*dataCollector)[0]->getNumberVertices();
    int nElements = (*dataCollector)[0]->getNumberElements();
    int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);

28
    if (dim == 2 && degree == 2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
29 30
      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
      nElements *= 4;
31
    } else if (dim == 2 && degree == 3) {
Thomas Witkowski's avatar
Thomas Witkowski committed
32 33
      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
      nElements *= 9;
34
    } else if (dim == 2 && degree == 4) {
Thomas Witkowski's avatar
Thomas Witkowski committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
      nVertices += (*dataCollector)[0]->getNumberInterpPoints();
      nElements *= 16;
    }
    file << "<?xml version=\"1.0\"?>\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "  <UnstructuredGrid>\n";
    file << "    <Piece NumberOfPoints=\"" << nVertices 
	 << "\" NumberOfCells=\"" << nElements << "\">\n";
    file << "      <Points>\n";
    file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    
    writeVertexCoords(file);
    
    file << "        </DataArray>\n";
    file << "      </Points>\n";
    file << "      <Cells>\n";
    file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
    
Thomas Witkowski's avatar
Thomas Witkowski committed
53
    for (int i = 0; i < nElements; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
      file << " " << (i + 1) * vertices << "\n";

    file << "        </DataArray>\n";
    file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";

    for (int i = 0; i < nElements; i++) {
      switch (vertices) {
      case 2:
	file << " 3\n";
	break;
      case 3: 
	file << " 5\n";
	break;
      case 4:
	file << " 10\n";
	break;		    
      default:
	break;
      }	   
    }

    file << "        </DataArray>\n";
    file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
77

Thomas Witkowski's avatar
Thomas Witkowski committed
78 79 80 81 82 83 84
    writeConnectivity(file);
    
    file << "        </DataArray>\n";    
    file << "      </Cells>\n";
    file << "      <PointData>\n";
    
    for (int i = 0; i < static_cast<int>(dataCollector->size()); i++) {
85 86
      file << "        <DataArray type=\"Float32\" Name=\"" 
	   << (*dataCollector)[i]->getValues()->getName() 
Thomas Witkowski's avatar
Thomas Witkowski committed
87
	   << "\" format=\"ascii\">\n";
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89 90 91 92 93 94 95 96 97 98
      writeVertexValues(file, i);
      
      file << "        </DataArray>\n";
    }
    
    file << "      </PointData>\n";
    file << "    </Piece>\n";
    file << "  </UnstructuredGrid>\n";
    file << "</VTKFile>\n";
  }
99

Thomas Witkowski's avatar
Thomas Witkowski committed
100 101
  
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
102
  void VtkWriter::writeVertexCoords(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
103
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
104
    DOFVector<list<VertexInfo> > *vertexInfos = 
105
      (*dataCollector)[0]->getVertexInfos();
Thomas Witkowski's avatar
Thomas Witkowski committed
106
    DOFVector<list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
107 108 109 110 111
    int counter = 0;

    // For all DOFs of vertices, write the coordinates.
    for (it.reset(); !it.end(); ++it) {
      // for all vertex infos of this DOF
Thomas Witkowski's avatar
Thomas Witkowski committed
112
      for (list<VertexInfo>::iterator it2 = it->begin(); it2 != it->end(); ++it2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
113 114
	it2->outputIndex = counter++;
	writeCoord(file, it2->coords);
115
      }     
Thomas Witkowski's avatar
Thomas Witkowski committed
116 117 118 119
    }

    // For the second dim case, write also the interpolation points.
    if ((dim == 2) && (degree > 1)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
120 121
      DOFVector<list< WorldVector<double> > > *interpPointCoords = (*dataCollector)[0]->getInterpPointCoords();
      DOFVector<list< WorldVector<double> > >::Iterator pointIt(interpPointCoords, USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
122
      
123
      counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
124
      for (pointIt.reset(); !pointIt.end(); ++pointIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
125
	for (list< WorldVector<double> >::iterator it2 = pointIt->begin(); 
126 127
	     it2 != pointIt->end(); ++it2) {
	  counter++;
Thomas Witkowski's avatar
Thomas Witkowski committed
128
	  writeCoord(file, *it2);
129
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
130 131 132 133
      }
    }
  }

134

Thomas Witkowski's avatar
Thomas Witkowski committed
135
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
136
  void VtkWriter::writeVertexValues(T &file, int componentNo)
Thomas Witkowski's avatar
Thomas Witkowski committed
137 138 139
  {
    DOFVector<int> *interpPointInd = (*dataCollector)[componentNo]->getInterpPointInd();
    DOFVector<double> *values = (*dataCollector)[componentNo]->getValues();
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    DOFVector<list<WorldVector<double> > > *dofCoords = (*dataCollector)[componentNo]->getDofCoords();
Thomas Witkowski's avatar
Thomas Witkowski committed
141 142 143

    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
    DOFVector<double>::Iterator valueIt(values, USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
144
    DOFVector<list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS);
145

Thomas Witkowski's avatar
Thomas Witkowski committed
146
    file << scientific;
147
    file.precision(15);
148
         
Thomas Witkowski's avatar
Thomas Witkowski committed
149 150 151 152 153
    // Write the values for all vertex DOFs.
    for (intPointIt.reset(), valueIt.reset(), coordIt.reset();
	 !intPointIt.end(); 
	 ++intPointIt, ++valueIt, ++coordIt) {

154 155 156
      if (*intPointIt == -2) 
	for (int i = 0; i < static_cast<int>(coordIt->size()); i++) 
	  file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";      
Thomas Witkowski's avatar
Thomas Witkowski committed
157 158 159 160
    }        

    // For the second dim case, write also the values of the interpolation points.
    if ((dim == 2) && (degree > 1)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
161
      DOFVector<list<WorldVector<double> > >::Iterator 
Thomas Witkowski's avatar
Thomas Witkowski committed
162 163 164 165 166 167 168
	interpCoordIt((*dataCollector)[componentNo]->getInterpPointCoords(), USED_DOFS);
    
      for (intPointIt.reset(), valueIt.reset(), interpCoordIt.reset();
	   !intPointIt.end(); 
	   ++intPointIt, ++valueIt, ++interpCoordIt) {      
	
	if (*intPointIt >= 0) {
169
	  for (unsigned int i = 0; i < interpCoordIt->size(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
170 171 172 173 174
	    file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";
	}
      }
    }    
  }
175

Thomas Witkowski's avatar
Thomas Witkowski committed
176 177
  
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
178
  void VtkWriter::writeConnectivity(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
179 180 181 182 183 184 185 186 187 188
  {
    // For the second dim case, and if higher order Lagrange elements are used,
    // write the connectivity by extra functions.
    if ((dim == 2) && (degree == 2)) {
      writeConnectivity_dim2_degree2(file);
    } else if ((dim == 2) && (degree == 3)) {
      writeConnectivity_dim2_degree3(file);
    } else if ((dim == 2) && (degree == 4)) {
      writeConnectivity_dim2_degree4(file);   
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
189 190
      list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
      list<ElementInfo>::iterator elementIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
191 192 193 194 195 196 197 198 199 200 201 202
      int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);
      
      for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
	// for all vertices
	for (int i = 0; i < vertices; i++) {
	  file << " " << elementIt->vertexInfo[i]->outputIndex;
	}
	file << "\n";
      } 
    }
  }
  
203

Thomas Witkowski's avatar
Thomas Witkowski committed
204
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
205
  void VtkWriter::writeConnectivity_dim2_degree2(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
206
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
207 208
    list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
    list<ElementInfo>::iterator elementIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210 211
    vector<vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
    vector<vector<int> >::iterator pointIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236

    int nVertices = (*dataCollector)[0]->getNumberVertices();

    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
	 pointIt != interpPoints->end(); 
	 ++pointIt, ++elementIt) {
    
      file << " " << elementIt->vertexInfo[0]->outputIndex
	   << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[2] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[2]->outputIndex 
	   << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[1] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[1]->outputIndex
	   << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[2] + nVertices << "\n";

      file << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[1] + nVertices 
	   << " " << (*pointIt)[2] + nVertices << "\n";
    }
  }

237

Thomas Witkowski's avatar
Thomas Witkowski committed
238
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
239
  void VtkWriter::writeConnectivity_dim2_degree3(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
240
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
241 242
    list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
    list<ElementInfo>::iterator elementIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244 245
    vector<vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
    vector<vector<int> >::iterator pointIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292

    int nVertices = (*dataCollector)[0]->getNumberVertices();

    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
	 pointIt != interpPoints->end(); 
	 ++pointIt, ++elementIt) {
                          
      file << " " << elementIt->vertexInfo[0]->outputIndex
	   << " " << (*pointIt)[3] + nVertices
	   << " " << (*pointIt)[4] + nVertices << "\n";

      file << " " << (*pointIt)[4] + nVertices
	   << " " << (*pointIt)[5] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << (*pointIt)[3] + nVertices
	   << " " << (*pointIt)[4] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << (*pointIt)[2] + nVertices
	   << " " << (*pointIt)[3] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[1]->outputIndex
	   << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[5] + nVertices << "\n";

      file << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[6] + nVertices
	   << " " << (*pointIt)[5] + nVertices << "\n";

      file << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[2] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[2]->outputIndex 
	   << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[2] + nVertices << "\n";
     
      
    }
  }

293

Thomas Witkowski's avatar
Thomas Witkowski committed
294
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
295
  void VtkWriter::writeConnectivity_dim2_degree4(T &file)
Thomas Witkowski's avatar
Thomas Witkowski committed
296
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
297 298
    list<ElementInfo> *elements = (*dataCollector)[0]->getElementInfos();
    list<ElementInfo>::iterator elementIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
299

Thomas Witkowski's avatar
Thomas Witkowski committed
300 301
    vector<vector<int> > *interpPoints = (*dataCollector)[0]->getInterpPoints();
    vector<vector<int> >::iterator pointIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376

    int nVertices = (*dataCollector)[0]->getNumberVertices();

    for (pointIt = interpPoints->begin(), elementIt = elements->begin(); 
	 pointIt != interpPoints->end(); 
	 ++pointIt, ++elementIt) {

      file << " " << elementIt->vertexInfo[0]->outputIndex 
	   << " " << (*pointIt)[5] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << (*pointIt)[5] + nVertices
	   << " " << (*pointIt)[9] + nVertices
	   << " " << (*pointIt)[6] + nVertices << "\n";

      file << " " << (*pointIt)[6] + nVertices
	   << " " << (*pointIt)[7] + nVertices
	   << " " << (*pointIt)[9] + nVertices << "\n";

      file << " " << (*pointIt)[7] + nVertices
	   << " " << (*pointIt)[9] + nVertices
	   << " " << (*pointIt)[10] + nVertices << "\n";

      file << " " << (*pointIt)[7] + nVertices
	   << " " << (*pointIt)[8] + nVertices
	   << " " << (*pointIt)[10] + nVertices << "\n";

      file << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[8] + nVertices
	   << " " << (*pointIt)[10] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[1]->outputIndex 
	   << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[8] + nVertices << "\n";

      file << " " << (*pointIt)[4] + nVertices
	   << " " << (*pointIt)[5] + nVertices
	   << " " << (*pointIt)[9] + nVertices << "\n";

      file << " " << (*pointIt)[4] + nVertices
	   << " " << (*pointIt)[9] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << (*pointIt)[9] + nVertices
	   << " " << (*pointIt)[10] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[10] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << (*pointIt)[0] + nVertices
	   << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[10] + nVertices << "\n";

      file << " " << (*pointIt)[3] + nVertices
	   << " " << (*pointIt)[4] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << (*pointIt)[2] + nVertices
	   << " " << (*pointIt)[3] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << (*pointIt)[1] + nVertices
	   << " " << (*pointIt)[2] + nVertices
	   << " " << (*pointIt)[11] + nVertices << "\n";

      file << " " << elementIt->vertexInfo[2]->outputIndex 
	   << " " << (*pointIt)[2] + nVertices
	   << " " << (*pointIt)[3] + nVertices << "\n";

    }
  }

}