Chaste  Release::3.4
XdmfMeshWriter.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 <sstream>
37 #include <map>
38 
39 #include "XdmfMeshWriter.hpp"
40 #include "DistributedTetrahedralMesh.hpp"
41 #include "Version.hpp"
42 
43 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
45  const std::string& rBaseName,
46  const bool clearOutputDir)
47  : AbstractTetrahedralMeshWriter<ELEMENT_DIM, SPACE_DIM>(rDirectory, rBaseName, clearOutputDir),
48  mNumberOfTimePoints(1u),
49  mTimeStep(1.0)
50 {
51 }
52 
53 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
55  bool keepOriginalElementIndexing)
56 {
57 #ifdef _MSC_VER
58  EXCEPTION("XDMF is not supported under Windows at present.");
59 #else
60  assert(keepOriginalElementIndexing);
61  this->mpDistributedMesh = dynamic_cast<DistributedTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* >(&rMesh);
62  bool mesh_is_distributed = (this->mpDistributedMesh != NULL) && PetscTools::IsParallel();
63 
65  {
66  // Write main test Grid collection (to be later replaced by temporal collection)
67  // Write references to geometry and topology chunk(s)
68  unsigned num_chunks = 1;
69  if (mesh_is_distributed)
70  {
71  num_chunks = PetscTools::GetNumProcs();
72  }
73  WriteXdmfMasterFile(num_chunks);
74  }
75  if (!mesh_is_distributed && !PetscTools::AmMaster())
76  {
77  //If the mesh is not distributed then the master knows everything and will write the geometry/topology as a single chunk
78  PetscTools::Barrier("XdmfMeshWriter wait for chunks to be written");
79  return;
80  }
81 
82  // Geometry
83  std::stringstream local_geometry_file_name;
84  local_geometry_file_name << this->mBaseName << "_geometry_"<< PetscTools::GetMyRank() <<".xml";
85  out_stream geometry_file = this->mpOutputFileHandler->OpenOutputFile(local_geometry_file_name.str());
86  std::string geom_type = "XYZ";
87  if (SPACE_DIM == 2)
88  {
89  geom_type = "XY";
90  }
91  (*geometry_file) << "<Geometry GeometryType=\""<< geom_type <<"\">\n";
92  unsigned num_nodes = rMesh.GetNumNodes();
93  if (this->mpDistributedMesh)
94  {
95  num_nodes = this->mpDistributedMesh->GetNumLocalNodes() + this->mpDistributedMesh->GetNumHaloNodes();
96  }
97 
98  (*geometry_file) << "\t<DataItem Format=\"XML\" Dimensions=\""<< num_nodes <<" "<< SPACE_DIM <<"\" DataType=\"Float\">";
99 
100  // Map a global node index into a local index (into mNodes and mHaloNodes as if they were concatenated)
101  std::map<unsigned, unsigned> global_to_node_index_map;
102 
103  //Node index that we are writing to the chunk (index into mNodes and mHaloNodes as if they were concatenated)
104  unsigned index = 0;
105 
106  // Owned nodes come first
108  iter != rMesh.GetNodeIteratorEnd();
109  ++iter)
110  {
111  global_to_node_index_map[iter->GetIndex()] = index;
112  index++;
113  (*geometry_file) << "\n\t\t";
114  c_vector<double, SPACE_DIM> current_item = (iter)->rGetLocation();
115  for (unsigned j=0; j<SPACE_DIM; j++)
116  {
117  (*geometry_file) << current_item[j] << "\t";
118  }
119  }
120 
121  // Halo nodes
122  if (this->mpDistributedMesh)
123  {
124  for (typename DistributedTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::HaloNodeIterator halo_iter=this->mpDistributedMesh->GetHaloNodeIteratorBegin();
125  halo_iter != this->mpDistributedMesh->GetHaloNodeIteratorEnd();
126  ++halo_iter)
127  {
128  global_to_node_index_map[(*halo_iter)->GetIndex()] = index;
129  index++;
130  (*geometry_file) << "\n\t\t";
131  c_vector<double, SPACE_DIM> current_item = (*halo_iter)->rGetLocation();
132  for (unsigned j=0; j<SPACE_DIM; j++)
133  {
134  (*geometry_file) << current_item[j] << "\t";
135  }
136  }
137  }
138  (*geometry_file) << "\n";
139 
140  (*geometry_file) << "\t</DataItem>\n";
141  (*geometry_file) << "</Geometry>\n";
142  (*geometry_file) << "<!-- " + ChasteBuildInfo::GetProvenanceString() + "-->\n";
143  geometry_file->close();
144 
145  // Topology
146  std::stringstream local_topology_file_name;
147  local_topology_file_name << this->mBaseName << "_topology_"<< PetscTools::GetMyRank() <<".xml";
148  out_stream topology_file = this->mpOutputFileHandler->OpenOutputFile(local_topology_file_name.str());
149  std::string top_type = "Tetrahedron";
150  if (SPACE_DIM == 2)
151  {
152  top_type = "Triangle";
153  }
154  unsigned num_elems = rMesh.GetNumElements();
155  if (this->mpDistributedMesh)
156  {
157  num_elems = this->mpDistributedMesh->GetNumLocalElements();
158  }
159  (*topology_file) << "<Topology TopologyType=\""<< top_type <<"\" NumberOfElements=\""<< num_elems <<"\">\n";
160  (*topology_file) << "\t<DataItem Format=\"XML\" Dimensions=\""<< num_elems <<" "<< ELEMENT_DIM+1 <<"\">";
162  elem_iter != rMesh.GetElementIteratorEnd();
163  ++elem_iter)
164  {
165  (*topology_file) << "\n\t\t";
166  for (unsigned j=0; j<ELEMENT_DIM+1; j++)
167  {
168  unsigned local_index = global_to_node_index_map[ elem_iter->GetNodeGlobalIndex(j) ];
169  (*topology_file) << local_index <<"\t";
170  }
171  }
172  (*topology_file) << "\n";
173 
174  (*topology_file) << "\t</DataItem>\n";
175  (*topology_file) << "</Topology>\n";
176  (*topology_file) << "<!-- " + ChasteBuildInfo::GetProvenanceString() + "-->\n";
177  topology_file->close();
178  PetscTools::Barrier("XdmfMeshWriter wait for chunks to be written");
179 #endif
180 }
181 
182 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
184 {
185 #ifdef _MSC_VER
186  EXCEPTION("XDMF is not supported under Windows at present.");
187 #else
188  // This method is only called when there is no mesh. We are writing from a reader.
189  if (PetscTools::AmMaster())
190  {
191  WriteXdmfMasterFile();
192 
193  // Geometry
194  out_stream geometry_file = this->mpOutputFileHandler->OpenOutputFile(this->mBaseName + "_geometry_0.xml");
195  std::string geom_type = "XYZ";
196  if (SPACE_DIM == 2)
197  {
198  geom_type = "XY";
199  }
200  (*geometry_file) << "<Geometry GeometryType=\""<< geom_type <<"\">\n";
201  (*geometry_file) << "\t<DataItem Format=\"XML\" Dimensions=\""<< this->GetNumNodes() <<" "<< SPACE_DIM <<"\" DataType=\"Float\">";
202  for (unsigned item_num=0; item_num<this->GetNumNodes(); item_num++)
203  {
204  (*geometry_file) << "\n\t\t";
205  std::vector<double> current_item = this->GetNextNode();
206  for (unsigned j=0; j<SPACE_DIM; j++)
207  {
208  (*geometry_file) << current_item[j]<<"\t";
209  }
210  }
211  (*geometry_file) << "\n";
212 
213  (*geometry_file) << "\t</DataItem>\n";
214  (*geometry_file) << "</Geometry>\n";
215  (*geometry_file) << "<!-- " + ChasteBuildInfo::GetProvenanceString() + "-->\n";
216  geometry_file->close();
217 
218  // Topology
219  out_stream topology_file = this->mpOutputFileHandler->OpenOutputFile(this->mBaseName + "_topology_0.xml");
220  std::string top_type = "Tetrahedron";
221  if (SPACE_DIM == 2)
222  {
223  top_type = "Triangle";
224  }
225  (*topology_file) << "<Topology TopologyType=\""<< top_type <<"\" NumberOfElements=\""<< this->GetNumElements() <<"\">\n";
226  (*topology_file) << "\t<DataItem Format=\"XML\" Dimensions=\""<< this->GetNumElements() <<" "<< ELEMENT_DIM+1 <<"\">";
227  for (unsigned item_num=0; item_num<this->GetNumElements(); item_num++)
228  {
229  (*topology_file) << "\n\t\t";
230  std::vector<unsigned> current_item = this->GetNextElement().NodeIndices;
231  for (unsigned j=0; j<ELEMENT_DIM+1; j++)
232  {
233  (*topology_file) << current_item[j]<<"\t";
234  }
235  }
236  (*topology_file) << "\n";
237 
238  (*topology_file) << "\t</DataItem>\n";
239  (*topology_file) << "</Topology>\n";
240  (*topology_file) << "<!-- " + ChasteBuildInfo::GetProvenanceString() + "-->\n";
241  topology_file->close();
242  }
243 #endif
244 }
245 
246 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
248 {
249 #ifndef _MSC_VER
250  assert(PetscTools::AmMaster());
251  // Define namespace symbols
252  XERCES_CPP_NAMESPACE_USE
253 
254  // Initialize Xerces.
255  XMLPlatformUtils::Initialize();
256 
257  DOMImplementation* p_DOM_implementation = DOMImplementationRegistry::getDOMImplementation(X("core"));
258 
259  DOMDocumentType* p_DOM_document_type = p_DOM_implementation->createDocumentType(X("Xdmf"),0,X("Xdmf.dtd"));
260  DOMDocument* p_DOM_document = p_DOM_implementation->createDocument(0, X("Xdmf"), p_DOM_document_type);
261  DOMElement* p_root_element = p_DOM_document->getDocumentElement();
262  p_root_element->setAttribute(X("Version"), X("2.0"));
263  p_root_element->setAttribute(X("xmlns:xi"), X("http://www.w3.org/2001/XInclude"));
264 
265  DOMElement* p_domain_element = p_DOM_document->createElement(X("Domain"));
266  p_root_element->appendChild(p_domain_element);
267 
268  // Temporal collection
269  DOMElement* p_grid_temp_collection_element = p_DOM_document->createElement(X("Grid"));
270  p_grid_temp_collection_element->setAttribute(X("CollectionType"), X("Temporal"));
271  p_grid_temp_collection_element->setAttribute(X("GridType"), X("Collection"));
272  p_domain_element->appendChild(p_grid_temp_collection_element);
273 
274  // Time values
275  DOMElement* p_time_element = p_DOM_document->createElement(X("Time"));
276  p_time_element->setAttribute(X("TimeType"), X("HyperSlab"));
277  p_grid_temp_collection_element->appendChild(p_time_element);
278 
279  DOMElement* p_time_dataitem_element = p_DOM_document->createElement(X("DataItem"));
280  p_time_dataitem_element->setAttribute(X("Format"),X("XML"));
281  p_time_dataitem_element->setAttribute(X("NumberType"),X("Float"));
282  p_time_dataitem_element->setAttribute(X("Dimensions"),X("3"));
283  p_time_element->appendChild(p_time_dataitem_element);
284 
285  std::stringstream time_stream;
286  time_stream << "0.0 " << mTimeStep << " " << mNumberOfTimePoints;
287  DOMText* p_time_text = p_DOM_document->createTextNode(X(time_stream.str()));
288  p_time_dataitem_element->appendChild(p_time_text);
289 
290  for(unsigned t=0; t<mNumberOfTimePoints; ++t)
291  {
292  DOMElement* p_grid_collection_element = p_DOM_document->createElement(X("Grid"));
293  p_grid_collection_element->setAttribute(X("CollectionType"), X("Spatial"));
294  p_grid_collection_element->setAttribute(X("GridType"), X("Collection"));
295  //p_grid_collection_element->setAttribute(X("Name"), X("spatial_collection"));
296  p_grid_temp_collection_element->appendChild(p_grid_collection_element);
297 
298  if (t==0)
299  {
300  for (unsigned chunk=0; chunk<numberOfChunks; chunk++)
301  {
302  std::stringstream chunk_stream;
303  chunk_stream << chunk;
304 
305  DOMElement* p_grid_element = p_DOM_document->createElement(X("Grid"));
306  p_grid_element->setAttribute(X("GridType"), X("Uniform"));
307  p_grid_element->setAttribute(X("Name"), X("Chunk_" + chunk_stream.str()));
308  p_grid_collection_element->appendChild(p_grid_element);
309 
310  //DOMElement* p_geom_element = p_DOM_document->createElement(X("Geometry"));
311  //p_geom_element->setAttribute(X("Reference"),X("/Xdmf/Domain/Geometry[1]"));
312  DOMElement* p_geom_element = p_DOM_document->createElement(X("xi:include"));
313  p_geom_element->setAttribute(X("href"), X(this->mBaseName+"_geometry_"+chunk_stream.str()+".xml"));
314  p_grid_element->appendChild(p_geom_element);
315  //DOMElement* p_topo_element = p_DOM_document->createElement(X("Topology"));
316  //p_topo_element->setAttribute(X("Reference"),X("/Xdmf/Domain/Topology[1]"));
317  DOMElement* p_topo_element = p_DOM_document->createElement(X("xi:include"));
318  p_topo_element->setAttribute(X("href"), X(this->mBaseName+"_topology_"+chunk_stream.str()+".xml"));
319  p_grid_element->appendChild(p_topo_element);
320 
321  /*
322  * p_grid_element may now need an Attribute (node data). Call Annotate,
323  * which here does nothing, but in pde can be overloaded to print variables
324  */
325  AddDataOnNodes(p_grid_element, p_DOM_document, t);
326  }
327  }
328  else // t>0
329  {
330  for (unsigned chunk=0; chunk<numberOfChunks; chunk++)
331  {
332  std::stringstream chunk_stream;
333  chunk_stream << chunk;
334 
335  DOMElement* p_grid_element = p_DOM_document->createElement(X("Grid"));
336  p_grid_element->setAttribute(X("GridType"), X("Subset"));
337  p_grid_element->setAttribute(X("Section"), X("All"));
338  p_grid_collection_element->appendChild(p_grid_element);
339 
340  /*
341  * p_grid_element may now need an Attribute (node data). Call Annotate,
342  * which here does nothing, but in pde can be overloaded to print variables
343  */
344  AddDataOnNodes(p_grid_element, p_DOM_document, t);
345  DOMElement* p_grid_ref_element = p_DOM_document->createElement(X("Grid"));
346  p_grid_ref_element->setAttribute(X("GridType"), X("Uniform"));
347  p_grid_ref_element->setAttribute(X("Reference"), X("XML"));
348  //p_grid_ref_element->setAttribute(X("Name"), X("Chunk_" + chunk_stream.str()));
349 
350  DOMText* p_ref_text = p_DOM_document->createTextNode(X("/Xdmf/Domain/Grid/Grid/Grid[@Name=\"Chunk_"+chunk_stream.str()+"\"]"));
351  p_grid_ref_element->appendChild(p_ref_text);
352  p_grid_element->appendChild(p_grid_ref_element);
353  }
354  }
355  }
356  // Create a Comment node, and then append this to the root element.
357  DOMComment* p_provenance_comment = p_DOM_document->createComment(X(" "+ChasteBuildInfo::GetProvenanceString()));
358  p_DOM_document->appendChild(p_provenance_comment);
359 
360  XMLFormatTarget* p_target = new LocalFileFormatTarget(X(this->mpOutputFileHandler->GetOutputDirectoryFullPath() + this->mBaseName+".xdmf"));
361 
362 #if _XERCES_VERSION >= 30000
363  DOMLSSerializer *p_serializer = ((DOMImplementationLS*)p_DOM_implementation)->createLSSerializer();
364  p_serializer->getDomConfig()->setParameter(XMLUni::fgDOMWRTFormatPrettyPrint, true);
365  DOMLSOutput *p_output = ((DOMImplementationLS*)p_DOM_implementation)->createLSOutput();
366  p_output->setByteStream(p_target);
367  p_serializer->write(p_DOM_document, p_output);
368 #else
369  DOMWriter* p_serializer = ((DOMImplementationLS*)p_DOM_implementation)->createDOMWriter();
370  p_serializer->setFeature(XMLUni::fgDOMWRTFormatPrettyPrint, true);
371  p_serializer->writeNode(p_target, *p_DOM_document);
372 #endif
373 
374  // Cleanup
375  p_serializer->release();
376  p_DOM_document->release();
377  delete p_target;
378  XMLPlatformUtils::Terminate();
379 #endif // _MSC_VER
380 }
381 
383 // Explicit instantiation
385 
386 template class XdmfMeshWriter<1,1>;
387 template class XdmfMeshWriter<1,2>;
388 template class XdmfMeshWriter<1,3>;
389 template class XdmfMeshWriter<2,2>; // Actually used
390 template class XdmfMeshWriter<2,3>;
391 template class XdmfMeshWriter<3,3>; // Actually used
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
static void Barrier(const std::string callerId="")
Definition: PetscTools.cpp:134
XdmfMeshWriter(const std::string &rDirectory, const std::string &rBaseName, const bool clearOutputDir=true)
virtual unsigned GetNumElements() const
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::vector< Node< SPACE_DIM > * >::const_iterator HaloNodeIterator
static bool AmMaster()
Definition: PetscTools.cpp:120
NodeIterator GetNodeIteratorEnd()
virtual unsigned GetNumNodes() const
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
NodeIterator GetNodeIteratorBegin(bool skipDeletedNodes=true)
static bool IsParallel()
Definition: PetscTools.cpp:97
static std::string GetProvenanceString()
static unsigned GetNumProcs()
Definition: PetscTools.cpp:108
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
void WriteXdmfMasterFile(unsigned numberOfChunks=1u)