Chaste  Release::3.4
VertexMeshWriter.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "VertexMeshWriter.hpp"
37 #include "Version.hpp"
38 #include "Cylindrical2dVertexMesh.hpp"
39 #include "Toroidal2dVertexMesh.hpp"
40 
44 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
46 {
51 };
52 
54 // Implementation
56 
57 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
59  const std::string& rBaseName,
60  const bool clearOutputDir)
61  : AbstractMeshWriter<ELEMENT_DIM, SPACE_DIM>(rDirectory, rBaseName, clearOutputDir),
62  mpMesh(NULL),
63  mpIters(new MeshWriterIterators<ELEMENT_DIM, SPACE_DIM>),
64  mpNodeMap(NULL),
65  mNodeMapCurrentIndex(0)
66 {
67  mpIters->pNodeIter = NULL;
68  mpIters->pElemIter = NULL;
69 
70 #ifdef CHASTE_VTK
71  // Dubious, since we shouldn't yet know what any details of the mesh are.
72  mpVtkUnstructedMesh = vtkUnstructuredGrid::New();
73 #endif //CHASTE_VTK
74 }
75 
76 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
78 {
79  if (mpIters->pNodeIter)
80  {
81  delete mpIters->pNodeIter;
82  delete mpIters->pElemIter;
83  }
84 
85  delete mpIters;
86 
87  if (mpNodeMap)
88  {
89  delete mpNodeMap;
90  }
91 
92 #ifdef CHASTE_VTK
93  // Dubious, since we shouldn't yet know what any details of the mesh are.
94  mpVtkUnstructedMesh->Delete(); // Reference counted
95 #endif //CHASTE_VTK
96 }
97 
98 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
100 {
101  if (mpMesh)
102  {
103  // Sanity check
104  assert(this->mNumNodes == mpMesh->GetNumNodes());
105 
106  std::vector<double> coordinates(SPACE_DIM+1);
107 
108  // Get the node coordinates using the node iterator (thus skipping deleted nodes)
109  for (unsigned j=0; j<SPACE_DIM; j++)
110  {
111  coordinates[j] = (*(mpIters->pNodeIter))->GetPoint()[j];
112  }
113  coordinates[SPACE_DIM] = (*(mpIters->pNodeIter))->IsBoundaryNode();
114 
115  ++(*(mpIters->pNodeIter));
116 
117  return coordinates;
118  }
119  else
120  {
122  }
123 }
124 
125 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
127 {
128  // This method should only be called in 3D
129  assert(SPACE_DIM == 3);
130  assert(mpMesh);
131  assert(this->mNumElements == mpMesh->GetNumElements());
132 
133  // Create data structure for this element
134  VertexElementData elem_data;
135 
136  // Store node indices owned by this element
137  elem_data.NodeIndices.resize((*(mpIters->pElemIter))->GetNumNodes());
138  for (unsigned j=0; j<elem_data.NodeIndices.size(); j++)
139  {
140  unsigned old_index = (*(mpIters->pElemIter))->GetNodeGlobalIndex(j);
141  elem_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
142  }
143 
144  // Store faces owned by this element
145  elem_data.Faces.resize((*(mpIters->pElemIter))->GetNumFaces());
146  for (unsigned i=0; i<elem_data.Faces.size(); i++)
147  {
148  // Get pointer to this face
149  VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = (*(mpIters->pElemIter))->GetFace(i);
150 
151  // Create data structure for this face
152  ElementData face_data;
153 
154  // Store this face's index
155  face_data.AttributeValue = p_face->GetIndex();
156 
157  // Store node indices owned by this face
158  face_data.NodeIndices.resize(p_face->GetNumNodes());
159  for (unsigned j=0; j<face_data.NodeIndices.size(); j++)
160  {
161  unsigned old_index = p_face->GetNodeGlobalIndex(j);
162  face_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
163  }
164 
165  // Store this face
166  elem_data.Faces[i] = face_data;
167 
169  }
170 
171  // Set attribute
172  elem_data.AttributeValue = (unsigned)(*(mpIters->pElemIter))->GetAttribute();
173  ++(*(mpIters->pElemIter));
174 
175  return elem_data;
176 }
177 
178 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
180 {
182 
183  if (mpMesh)
184  {
185  assert(this->mNumElements == mpMesh->GetNumElements());
186 
187  ElementData elem_data;
188  elem_data.NodeIndices.resize((*(mpIters->pElemIter))->GetNumNodes());
189  for (unsigned j=0; j<elem_data.NodeIndices.size(); j++)
190  {
191  unsigned old_index = (*(mpIters->pElemIter))->GetNodeGlobalIndex(j);
192  elem_data.NodeIndices[j] = mpMesh->IsMeshChanging() ? mpNodeMap->GetNewIndex(old_index) : old_index;
193  }
194 
195  // Set attribute
196  elem_data.AttributeValue = (*(mpIters->pElemIter))->GetAttribute();
197  ++(*(mpIters->pElemIter));
198 
199  return elem_data;
200  }
201  else
202  {
204  }
205 }
206 
207 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
209 {
210 #ifdef CHASTE_VTK
211  assert(SPACE_DIM==3 || SPACE_DIM == 2);
212 
213  // Create VTK mesh
214  MakeVtkMesh(rMesh);
215 
216  // Now write VTK mesh to file
217  assert(mpVtkUnstructedMesh->CheckAttributes() == 0);
218  vtkXMLUnstructuredGridWriter* p_writer = vtkXMLUnstructuredGridWriter::New();
219 #if VTK_MAJOR_VERSION >= 6
220  p_writer->SetInputData(mpVtkUnstructedMesh);
221 #else
222  p_writer->SetInput(mpVtkUnstructedMesh);
223 #endif
224  // Uninitialised stuff arises (see #1079), but you can remove valgrind problems by removing compression:
225  // **** REMOVE WITH CAUTION *****
226  p_writer->SetCompressor(NULL);
227  // **** REMOVE WITH CAUTION *****
228 
229  std::string vtk_file_name = this->mpOutputFileHandler->GetOutputDirectoryFullPath() + this->mBaseName;
230  if (stamp != "")
231  {
232  vtk_file_name += "_" + stamp;
233  }
234  vtk_file_name += ".vtu";
235 
236  p_writer->SetFileName(vtk_file_name.c_str());
237  //p_writer->PrintSelf(std::cout, vtkIndent());
238  p_writer->Write();
239  p_writer->Delete(); // Reference counted
240 #endif //CHASTE_VTK
241 }
242 
249 template<>
251 {
252 #ifdef CHASTE_VTK
253  // Create VTK mesh
254  if (dynamic_cast<Toroidal2dVertexMesh*>(&rMesh))
255  {
256  MutableVertexMesh<2, 2>* p_mesh_for_vtk = static_cast<Toroidal2dVertexMesh*>(&rMesh)->GetMeshForVtk();
257  MakeVtkMesh(*p_mesh_for_vtk);
258 
259  // Avoid memory leak
260  delete p_mesh_for_vtk;
261  }
262  else if (dynamic_cast<Cylindrical2dVertexMesh*>(&rMesh))
263  {
264  MutableVertexMesh<2, 2>* p_mesh_for_vtk = static_cast<Cylindrical2dVertexMesh*>(&rMesh)->GetMeshForVtk();
265  MakeVtkMesh(*p_mesh_for_vtk);
266 
267  // Avoid memory leak
268  delete p_mesh_for_vtk;
269  }
270  else
271  {
272  MakeVtkMesh(rMesh);
273  }
274 
275  // Now write VTK mesh to file
276  assert(mpVtkUnstructedMesh->CheckAttributes() == 0);
277  vtkXMLUnstructuredGridWriter* p_writer = vtkXMLUnstructuredGridWriter::New();
278 #if VTK_MAJOR_VERSION >= 6
279  p_writer->SetInputData(mpVtkUnstructedMesh);
280 #else
281  p_writer->SetInput(mpVtkUnstructedMesh);
282 #endif
283  // Uninitialised stuff arises (see #1079), but you can remove valgrind problems by removing compression:
284  // **** REMOVE WITH CAUTION *****
285  p_writer->SetCompressor(NULL);
286  // **** REMOVE WITH CAUTION *****
287 
288  std::string vtk_file_name = this->mpOutputFileHandler->GetOutputDirectoryFullPath() + this->mBaseName;
289  if (stamp != "")
290  {
291  vtk_file_name += "_" + stamp;
292  }
293  vtk_file_name += ".vtu";
294 
295  p_writer->SetFileName(vtk_file_name.c_str());
296  //p_writer->PrintSelf(std::cout, vtkIndent());
297  p_writer->Write();
298  p_writer->Delete(); // Reference counted
299 #endif //CHASTE_VTK
300 }
301 
302 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
304 {
305 #ifdef CHASTE_VTK
306  // Make the Vtk mesh
307  vtkPoints* p_pts = vtkPoints::New(VTK_DOUBLE);
308  p_pts->GetData()->SetName("Vertex positions");
309  for (unsigned node_num=0; node_num<rMesh.GetNumNodes(); node_num++)
310  {
311  c_vector<double, SPACE_DIM> position = rMesh.GetNode(node_num)->rGetLocation();
312  if (SPACE_DIM==2)
313  {
314  p_pts->InsertPoint(node_num, position[0], position[1], 0.0);
315  }
316  else
317  {
318  p_pts->InsertPoint(node_num, position[0], position[1], position[2]);
319  }
320  }
321 
322  mpVtkUnstructedMesh->SetPoints(p_pts);
323  p_pts->Delete(); // Reference counted
325  iter != rMesh.GetElementIteratorEnd();
326  ++iter)
327  {
328  vtkCell* p_cell;
329  if (SPACE_DIM == 2)
330  {
331  p_cell = vtkPolygon::New();
332  }
333  else
334  {
335  p_cell = vtkConvexPointSet::New();
336  }
337  vtkIdList* p_cell_id_list = p_cell->GetPointIds();
338  p_cell_id_list->SetNumberOfIds(iter->GetNumNodes());
339  for (unsigned j=0; j<iter->GetNumNodes(); ++j)
340  {
341  p_cell_id_list->SetId(j, iter->GetNodeGlobalIndex(j));
342  }
343  mpVtkUnstructedMesh->InsertNextCell(p_cell->GetCellType(), p_cell_id_list);
344  p_cell->Delete(); // Reference counted
345  }
346 #endif //CHASTE_VTK
347 }
348 
349 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
350 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddCellData(std::string dataName, std::vector<double> dataPayload)
351 {
352 #ifdef CHASTE_VTK
353  vtkDoubleArray* p_scalars = vtkDoubleArray::New();
354  p_scalars->SetName(dataName.c_str());
355  for (unsigned i=0; i<dataPayload.size(); i++)
356  {
357  p_scalars->InsertNextValue(dataPayload[i]);
358  }
359 
360  vtkCellData* p_cell_data = mpVtkUnstructedMesh->GetCellData();
361  p_cell_data->AddArray(p_scalars);
362  p_scalars->Delete(); // Reference counted
363 #endif //CHASTE_VTK
364 }
365 
366 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
367 void VertexMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddPointData(std::string dataName, std::vector<double> dataPayload)
368 {
369 #ifdef CHASTE_VTK
370  vtkDoubleArray* p_scalars = vtkDoubleArray::New();
371  p_scalars->SetName(dataName.c_str());
372  for (unsigned i=0; i<dataPayload.size(); i++)
373  {
374  p_scalars->InsertNextValue(dataPayload[i]);
375  }
376 
377  vtkPointData* p_point_data = mpVtkUnstructedMesh->GetPointData();
378  p_point_data->AddArray(p_scalars);
379  p_scalars->Delete(); // Reference counted
380 #endif //CHASTE_VTK
381 }
382 
384 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
386 {
387  this->mpMeshReader = NULL;
388  mpMesh = &rMesh;
389 
390  this->mNumNodes = mpMesh->GetNumNodes();
391  this->mNumElements = mpMesh->GetNumElements();
392 
393  typedef typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator NodeIterType;
394  mpIters->pNodeIter = new NodeIterType(mpMesh->GetNodeIteratorBegin());
395 
396  typedef typename VertexMesh<ELEMENT_DIM,SPACE_DIM>::VertexElementIterator ElemIterType;
397  mpIters->pElemIter = new ElemIterType(mpMesh->GetElementIteratorBegin());
398 
399  // Set up node map if we might have deleted nodes
400  mNodeMapCurrentIndex = 0;
401  if (mpMesh->IsMeshChanging())
402  {
403  mpNodeMap = new NodeMap(mpMesh->GetNumAllNodes());
404  for (NodeIterType it = mpMesh->GetNodeIteratorBegin(); it != mpMesh->GetNodeIteratorEnd(); ++it)
405  {
406  mpNodeMap->SetNewIndex(it->GetIndex(), mNodeMapCurrentIndex++);
407  }
408  }
409  WriteFiles();
410 }
411 
412 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
414 {
415  std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
416 
417  // Write node file
418  std::string node_file_name = this->mBaseName + ".node";
419  out_stream p_node_file = this->mpOutputFileHandler->OpenOutputFile(node_file_name);
420 
421  // Write the node header
422  unsigned num_attr = 0;
423  unsigned max_bdy_marker = 1; // as we include boundary node information in the node file
424  unsigned num_nodes = this->GetNumNodes();
425 
426  *p_node_file << num_nodes << "\t";
427  *p_node_file << SPACE_DIM << "\t";
428  *p_node_file << num_attr << "\t";
429  *p_node_file << max_bdy_marker << "\n";
430  *p_node_file << std::setprecision(6);
431 
432  // Write each node's data
433  for (unsigned item_num=0; item_num<num_nodes; item_num++)
434  {
435  std::vector<double> current_item = this->GetNextNode();
436  *p_node_file << item_num;
437  for (unsigned i=0; i<SPACE_DIM+1; i++)
438  {
439  *p_node_file << "\t" << current_item[i];
440  }
441  *p_node_file << "\n";
442  }
443  *p_node_file << comment << "\n";
444  p_node_file->close();
445 
446  // Write element file
447  std::string element_file_name = this->mBaseName + ".cell";
448  out_stream p_element_file = this->mpOutputFileHandler->OpenOutputFile(element_file_name);
449 
450  // Write the element header
451  num_attr = 1; //Always write element attributes
452  unsigned num_elements = this->GetNumElements();
453  *p_element_file << num_elements << "\t" << num_attr << "\n";
454 
455  // Write each element's data
456  for (unsigned item_num=0; item_num<num_elements; item_num++)
457  {
458  if (SPACE_DIM == 2) // In 2D, write the node indices owned by this element
459  {
460  // Get data for this element
461  ElementData elem_data = this->GetNextElement();
462 
463  // Get the node indices owned by this element
464  std::vector<unsigned> node_indices = elem_data.NodeIndices;
465 
466  // Write this element's index and the number of nodes owned by it to file
467  *p_element_file << item_num << "\t" << node_indices.size();
468 
469  // Write the node indices owned by this element to file
470  for (unsigned i=0; i<node_indices.size(); i++)
471  {
472  *p_element_file << "\t" << node_indices[i];
473  }
474 
475  *p_element_file << "\t" << elem_data.AttributeValue;
476 
477  // New line
478  *p_element_file << "\n";
479  }
480  else // 3D
481  {
482  assert(SPACE_DIM == 3);
483 
484  // Get data for this element
485  VertexElementData elem_data_with_faces = this->GetNextElementWithFaces();
486 
487  // Get the node indices owned by this element
488  std::vector<unsigned> node_indices = elem_data_with_faces.NodeIndices;
489 
490  // Write this element's index and the number of nodes owned by it to file
491  *p_element_file << item_num << "\t" << node_indices.size();
492 
493  // Write the node indices owned by this element to file
494  for (unsigned i=0; i<node_indices.size(); i++)
495  {
496  *p_element_file << "\t" << node_indices[i];
497  }
498 
499  // Get the faces owned by this element (if any)
500  std::vector<ElementData> faces = elem_data_with_faces.Faces;
501 
502  // Write the number of faces owned by this element to file
503  *p_element_file << "\t" << faces.size();
504 
505  for (unsigned j=0; j<faces.size(); j++)
506  {
507  // Get the node indices owned by this face
508  std::vector<unsigned> face_node_indices = faces[j].NodeIndices;
509 
510  // Write this face's index and the number of nodes owned by it to file
511  *p_element_file << "\t" << faces[j].AttributeValue << "\t" << face_node_indices.size();
512 
513  // Write the node indices owned by this element to file
514  for (unsigned i=0; i<face_node_indices.size(); i++)
515  {
516  *p_element_file << "\t" << face_node_indices[i];
517  }
518 
520  }
521 
522  *p_element_file << "\t" << elem_data_with_faces.AttributeValue;
523 
524  // New line
525  *p_element_file << "\n";
526  }
527  }
528 
529  *p_element_file << comment << "\n";
530  p_element_file->close();
531 }
532 
534 // Explicit instantiation
536 
537 template class VertexMeshWriter<1,1>;
538 template class VertexMeshWriter<1,2>;
539 template class VertexMeshWriter<1,3>;
540 template class VertexMeshWriter<2,2>;
541 template class VertexMeshWriter<2,3>;
542 template class VertexMeshWriter<3,3>;
void AddCellData(std::string dataName, std::vector< double > dataPayload)
MeshWriterIterators< ELEMENT_DIM, SPACE_DIM > * mpIters
void WriteFilesUsingMesh(VertexMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
Node< SPACE_DIM > * GetNode(unsigned index) const
virtual unsigned GetNumNodes() const
Definition: VertexMesh.cpp:611
AbstractMesh< ELEMENT_DIM, SPACE_DIM >::NodeIterator * pNodeIter
virtual std::vector< double > GetNextNode()
vtkUnstructuredGrid * mpVtkUnstructedMesh
std::vector< ElementData > Faces
std::vector< unsigned > NodeIndices
VertexMeshWriter(const std::string &rDirectory, const std::string &rBaseName, const bool clearOutputDir=true)
void SetNewIndex(unsigned oldIndex, unsigned newIndex)
Definition: NodeMap.cpp:66
VertexElement< ELEMENT_DIM-1, SPACE_DIM > * GetFace(unsigned index) const
VertexElementData GetNextElementWithFaces()
VertexElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
Definition: VertexMesh.hpp:668
void WriteVtkUsingMesh(VertexMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, std::string stamp="")
VertexMesh< ELEMENT_DIM, SPACE_DIM >::VertexElementIterator * pElemIter
void MakeVtkMesh(VertexMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
std::vector< double > GetNextNode()
ElementData GetNextElement()
static std::string GetProvenanceString()
virtual ElementData GetNextElement()
VertexElementIterator GetElementIteratorEnd()
Definition: VertexMesh.hpp:675
std::vector< unsigned > NodeIndices
void AddPointData(std::string dataName, std::vector< double > dataPayload)