Chaste  Release::2018.1
TrianglesMeshReader.cpp
1 /*
2 
3 Copyright (c) 2005-2018, 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 #include <cassert>
36 #include <sstream>
37 #include <iostream>
38 
39 #include "TrianglesMeshReader.hpp"
40 #include "Exception.hpp"
41 
42 static const char* NODES_FILE_EXTENSION = ".node";
43 static const char* ELEMENTS_FILE_EXTENSION = ".ele";
44 static const char* FACES_FILE_EXTENSION = ".face";
45 static const char* EDGES_FILE_EXTENSION = ".edge";
46 static const char* NCL_FILE_EXTENSION = ".ncl";
47 static const char* CABLE_FILE_EXTENSION = ".cable";
48 
50 // Implementation
52 
53 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
55  unsigned orderOfElements,
56  unsigned orderOfBoundaryElements,
57  bool readContainingElementForBoundaryElements)
58  : mFilesBaseName(pathBaseName),
59  mNodeItemWidth(0),
60  mElementItemWidth(0),
61  mFaceItemWidth(0),
62  mNumNodes(0),
63  mNumElements(0),
64  mNumFaces(0),
65  mNumCableElements(0),
66  mNodesRead(0),
67  mElementsRead(0),
68  mCableElementsRead(0),
69  mFacesRead(0),
70  mBoundaryFacesRead(0),
71  mNclItemsRead(0),
72  mNumNodeAttributes(0),
73  mNumElementAttributes(0),
74  mNumFaceAttributes(0),
75  mNumCableElementAttributes(0),
76  mOrderOfElements(orderOfElements),
77  mOrderOfBoundaryElements(orderOfBoundaryElements),
78  mEofException(false),
79  mReadContainingElementOfBoundaryElement(readContainingElementForBoundaryElements),
80  mFilesAreBinary(false),
81  mMeshIsHexahedral(false),
82  mNodeFileReadBuffer(nullptr),
83  mElementFileReadBuffer(nullptr),
84  mFaceFileReadBuffer(nullptr),
85  mNodePermutationDefined(false)
86 {
87  // Only linear and quadratic elements
88  assert(orderOfElements==1 || orderOfElements==2);
90  {
91  EXCEPTION("Boundary element file should not have containing element info if it is quadratic");
92  }
93  if (mOrderOfElements==1)
94  {
95  mNodesPerElement = ELEMENT_DIM+1;
96  }
97  else
98  {
99  assert(SPACE_DIM==ELEMENT_DIM); // LCOV_EXCL_LINE
100  mNodesPerElement = (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2;
101  }
102 
104  {
105  mNodesPerBoundaryElement = ELEMENT_DIM;
106  }
107  else
108  {
109  assert(SPACE_DIM==ELEMENT_DIM); // LCOV_EXCL_LINE
110  mNodesPerBoundaryElement = ELEMENT_DIM*(ELEMENT_DIM+1)/2;
111  }
112 
113  mIndexFromZero = false; // Initially assume that nodes are not numbered from zero
114 
115  OpenFiles();
116  ReadHeaders();
117 }
118 
119 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
121 {
122  delete[] mNodeFileReadBuffer;
123  delete[] mElementFileReadBuffer;
124  delete[] mFaceFileReadBuffer;
125 }
126 
127 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
129 {
130  return mNumElements;
131 }
132 
133 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
135 {
136  return mNumNodes;
137 }
138 
139 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
141 {
142  return mNumFaces;
143 }
144 
145 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147 {
148  return mNumCableElements;
149 }
150 
151 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
153 {
154  return mNumElementAttributes;
155 }
156 
157 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
159 {
160  return mNumFaceAttributes;
161 }
162 
163 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
165 {
166  return mNumCableElementAttributes;
167 }
168 
169 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
171 {
172  CloseFiles();
173 
174  mNodesRead = 0;
175  mElementsRead = 0;
176  mFacesRead = 0;
177  mBoundaryFacesRead = 0;
178  mCableElementsRead = 0;
179  mNclItemsRead = 0;
180  mEofException = false;
181 
182  OpenFiles();
183  ReadHeaders();
184 }
185 
186 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
188 {
189  std::vector<double> ret_coords(SPACE_DIM);
190 
191  mNodeAttributes.clear(); // clear attributes for this node
192  GetNextItemFromStream(mNodesFile, mNodesRead, ret_coords, mNumNodeAttributes, mNodeAttributes);
193 
194  mNodesRead++;
195  return ret_coords;
196 }
197 
198 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
200 {
201  ElementData element_data;
202  element_data.NodeIndices.resize(mNodesPerElement);
203  element_data.AttributeValue = 0.0; // If an attribute is not read this stays as zero, otherwise overwritten.
204 
205  std::vector<double> element_attributes;
206  GetNextItemFromStream(mElementsFile, mElementsRead, element_data.NodeIndices, mNumElementAttributes, element_attributes);
207 
208  if (mNumElementAttributes > 0)
209  {
210  element_data.AttributeValue = element_attributes[0];
211  }
212 
213  EnsureIndexingFromZero(element_data.NodeIndices);
214 
215  mElementsRead++;
216 
217  if (mNodePermutationDefined)
218  {
219  for (std::vector<unsigned>::iterator node_it = element_data.NodeIndices.begin();
220  node_it != element_data.NodeIndices.end();
221  ++ node_it)
222  {
223  assert(*node_it < mPermutationVector.size());
224  *node_it = mPermutationVector[*node_it];
225  }
226  }
227 
228  return element_data;
229 }
230 
231 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
233 {
234  ElementData element_data;
235  element_data.NodeIndices.resize(2u);
236  element_data.AttributeValue = 0; // If an attribute is not read this stays as zero, otherwise overwritten.
237 
238  std::vector<double> cable_element_attributes;
239  GetNextItemFromStream(mCableElementsFile, mCableElementsRead, element_data.NodeIndices, mNumCableElementAttributes, cable_element_attributes);
240  if (mNumCableElementAttributes > 0)
241  {
242  element_data.AttributeValue = cable_element_attributes[0];
243  }
244 
245  EnsureIndexingFromZero(element_data.NodeIndices);
246 
247  mCableElementsRead++;
248 
249  // Node permutation can only be done with binary data...
250 // if (mNodePermutationDefined)
251 // {
252 // for (std::vector<unsigned>::iterator node_it = element_data.NodeIndices.begin();
253 // node_it != element_data.NodeIndices.end();
254 // ++ node_it)
255 // {
256 // assert(*node_it < mPermutationVector.size());
257 // *node_it = mPermutationVector[*node_it];
258 // }
259 // }
260 
261  return element_data;
262 }
263 
264 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
266 {
267  ElementData face_data;
268  std::vector<unsigned> ret_indices;
269 
270  // In the first case there's no file, all the nodes are set as faces
271  if (ELEMENT_DIM == 1)
272  {
273  ret_indices.push_back( mOneDimBoundary[mBoundaryFacesRead] );
274  }
275  else
276  {
277  ret_indices.resize(mNodesPerBoundaryElement);
278 
279  assert(ELEMENT_DIM != 0); // LCOV_EXCL_LINE //Covered in earlier exception, but needed in loop guard here.
280  do
281  {
282  face_data.AttributeValue = 1.0; // If an attribute is not read this stays as one, otherwise overwritten.
283 
284  std::vector<double> face_attributes; //will store face attributes, if any
285  if (mReadContainingElementOfBoundaryElement)
286  {
287  assert(mNumFaceAttributes == 0);
288  GetNextItemFromStream(mFacesFile, mFacesRead, ret_indices, 1, face_attributes);
289 
290  if (face_attributes.size() > 0)
291  {
292  face_data.ContainingElement = (unsigned) face_attributes[0];// only one face attribute registered for the moment
293  }
294 
295  }
296  else
297  {
298  GetNextItemFromStream(mFacesFile, mFacesRead, ret_indices, mNumFaceAttributes,
299  face_attributes);
300 
301  if (mNumFaceAttributes > 0)
302  {
303  face_data.AttributeValue = face_attributes[0]; //only one face attribute registered for the moment
304  }
305  }
306 
307  EnsureIndexingFromZero(ret_indices);
308 
309  mFacesRead++;
310  }
311  while (ELEMENT_DIM==2 && face_data.AttributeValue==0.0); //In triangles format we ignore internal edges (which are marked with attribute 0)
312  }
313 
314  mBoundaryFacesRead++;
315 
316  if (mNodePermutationDefined)
317  {
318  for (std::vector<unsigned>::iterator node_it = ret_indices.begin();
319  node_it != ret_indices.end();
320  ++ node_it)
321  {
322  assert(*node_it < mPermutationVector.size());
323  *node_it = mPermutationVector[*node_it];
324  }
325  }
326 
327  face_data.NodeIndices = ret_indices;
328 
329  return face_data;
330 }
331 
332 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
333 std::vector<double> TrianglesMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNode(unsigned index)
334 {
335  if (!mFilesAreBinary)
336  {
337  EXCEPTION("Random access is only implemented in mesh readers for binary mesh files.");
338  }
339  if (index >= mNumNodes)
340  {
341  EXCEPTION("Node does not exist - not enough nodes.");
342  }
343 
344  if (mNodePermutationDefined)
345  {
346  assert(index<mInversePermutationVector.size());
347  index = mInversePermutationVector[index];
348  }
349 
350  // Put the file stream pointer to the right location
351  if (index > mNodesRead)
352  {
353  // This is a monotonic (but non-contiguous) read. Let's assume that it's more efficient
354  // to seek from the current position rather than from the start of the file
355  mNodesFile.seekg( mNodeItemWidth*(index-mNodesRead), std::ios_base::cur);
356  }
357  else if (mNodesRead != index)
358  {
359  mNodesFile.seekg(mNodeFileDataStart + mNodeItemWidth*index, std::ios_base::beg);
360  }
361 
362  mNodesRead = index; // Allow GetNextNode() to note the position of the item after this one
363  // Read the next item.
364  return GetNextNode();
365 }
366 
367 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
369 {
370  if (!mFilesAreBinary)
371  {
372  EXCEPTION("Random access is only implemented in mesh readers for binary mesh files.");
373  }
374  if (index >= mNumElements)
375  {
376  EXCEPTION("Element " << index << " does not exist - not enough elements (only " << mNumElements << ").");
377  }
378 
379  // Put the file stream pointer to the right location
380  if (index > mElementsRead)
381  {
382  // This is a monotonic (but non-contiguous) read. Let's assume that it's more efficient
383  // to seek from the current position rather than from the start of the file
384  mElementsFile.seekg( mElementItemWidth*(index-mElementsRead), std::ios_base::cur);
385  }
386  else if (mElementsRead != index)
387  {
388  mElementsFile.seekg(mElementFileDataStart + mElementItemWidth*index, std::ios_base::beg);
389  }
390 
391  mElementsRead = index; // Allow GetNextElementData() to note the position of the item after this one
392  // Read the next item.
393  return GetNextElementData();
394 }
395 
396 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
398 {
399  if (!mFilesAreBinary)
400  {
401  EXCEPTION("Random access is only implemented in mesh readers for binary mesh files.");
402  }
403  if (index >=mNumFaces)
404  {
405  EXCEPTION("Face does not exist - not enough faces.");
406  }
407 
408  /*
409  *
410  if (index > mFacesRead)
411  {
412  // This would be a monotonic (but non-contiguous) read. But we don't actually read faces with this access pattern.
414  mFacesFile.seekg( mFaceItemWidth*(index-mFacesRead), std::ios_base::cur);
415  }
416  else
417  */
418  // Put the file stream pointer to the right location
419  if (mFacesRead != index)
420  {
421  mFacesFile.seekg(mFaceFileDataStart + mFaceItemWidth*index, std::ios_base::beg);
422  }
423  mFacesRead = index; // Allow next call to mark the position in the file stream
424 
425  // Read the next item
426  return GetNextFaceData();
427 }
428 
429 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
431 {
432  if (!mFilesAreBinary)
433  {
434  EXCEPTION("NCL file functionality is only implemented in mesh readers for binary mesh files.");
435  }
436 
437  if (!mNclFileAvailable)
438  {
439  EXCEPTION("No NCL file available for this mesh.");
440  }
441  if (index >= mNumNodes)
442  {
443  EXCEPTION("Connectivity list does not exist - not enough nodes.");
444  }
445 
446  if (mNodePermutationDefined)
447  {
448  assert(index < mInversePermutationVector.size());
449  index = mInversePermutationVector[index];
450  }
451 
452  // Put the file stream pointer to the right location
453  if (index > mNclItemsRead)
454  {
455  // This is a monotonic (but non-contiguous) read. Let's assume that it's more efficient
456  // to seek from the current position rather than from the start of the file
457  mNclFile.seekg( mNclItemWidth*(index-mNclItemsRead), std::ios_base::cur);
458  }
459  else if ( mNclItemsRead != index )
460  {
461  mNclFile.seekg(mNclFileDataStart + mNclItemWidth*index, std::ios_base::beg);
462  }
463 
464  // Read the next item
465  std::vector<unsigned> containing_element_indices;
466  containing_element_indices.resize(mMaxContainingElements);
467 
468  std::vector<double> dummy; // unused here
469  GetNextItemFromStream(mNclFile, index, containing_element_indices, 0, dummy);
470  mNclItemsRead = index + 1; //Ready for the next call
471 
472  EnsureIndexingFromZero(containing_element_indices);
473 
474  unsigned num_containing_elements = mMaxContainingElements;
475  while ( containing_element_indices[num_containing_elements-1] == UINT_MAX )
476  {
477  num_containing_elements--;
478  }
479 
480  containing_element_indices.resize(num_containing_elements);
481 
482  return containing_element_indices;
483 }
484 
485 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
487 {
488  OpenNodeFile();
489  OpenElementsFile();
490  OpenFacesFile();
491  OpenNclFile();
492  OpenCableElementsFile();
493 }
494 
495 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
497 {
498  // Nodes definition
499  std::string file_name = mFilesBaseName + NODES_FILE_EXTENSION;
500  mNodesFile.open(file_name.c_str(), std::ios::binary);
501  if (!mNodesFile.is_open())
502  {
503  EXCEPTION("Could not open data file: " + file_name);
504  }
505 }
506 
507 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
509 {
510  // Elements definition
511  std::string file_name;
512  if (ELEMENT_DIM == SPACE_DIM)
513  {
514  file_name = mFilesBaseName + ELEMENTS_FILE_EXTENSION;
515  }
516  else
517  {
518  if (ELEMENT_DIM == 1)
519  {
520  file_name = mFilesBaseName + EDGES_FILE_EXTENSION;
521  }
522  else if (ELEMENT_DIM == 2)
523  {
524  file_name = mFilesBaseName + FACES_FILE_EXTENSION;
525  }
526  else
527  {
528  EXCEPTION("Can't have a zero-dimensional mesh in a one-dimensional space");
529  }
530  }
531 
532  mElementsFile.open(file_name.c_str(), std::ios::binary);
533  if (!mElementsFile.is_open())
534  {
535  EXCEPTION("Could not open data file: " + file_name);
536  }
537 }
538 
539 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
541 {
542  // Faces/edges definition
543  std::string file_name;
544  if (ELEMENT_DIM == 3)
545  {
546  file_name = mFilesBaseName + FACES_FILE_EXTENSION;
547  }
548  else if (ELEMENT_DIM == 2)
549  {
550  file_name = mFilesBaseName + EDGES_FILE_EXTENSION;
551  }
552  else //if (ELEMENT_DIM == 1)
553  {
554  // There is no file, data will be read from the node file (with boundaries marked)
555  return;
556  }
557 
558  mFacesFile.open(file_name.c_str(), std::ios::binary);
559  if (!mFacesFile.is_open())
560  {
561  EXCEPTION("Could not open data file: " + file_name);
562  }
563 }
564 
565 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
567 {
568  std::string file_name = mFilesBaseName + NCL_FILE_EXTENSION;
569  mNclFile.open(file_name.c_str(), std::ios::binary);
570 
571  mNclFileAvailable = mNclFile.is_open();
572 }
573 
574 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
576 {
577  std::string file_name = mFilesBaseName + CABLE_FILE_EXTENSION;
578  mCableElementsFile.open(file_name.c_str(), std::ios::binary);
579  if (!mCableElementsFile.is_open())
580  {
581  mNumCableElements = 0u;
582  mNumCableElementAttributes = 0u;
583  }
584 }
585 
586 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
588 {
589  return mNodeAttributes;
590 }
591 
592 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
594 {
595  /*
596  * Reading node file header
597  */
598  std::string buffer;
599  GetNextLineFromStream(mNodesFile, buffer);
600  std::stringstream node_header_line(buffer);
601  unsigned dimension;
602  node_header_line >> mNumNodes >> dimension >> mNumNodeAttributes >> mMaxNodeBdyMarker;
603  if (SPACE_DIM != dimension)
604  {
605  EXCEPTION("SPACE_DIM != dimension read from file ");
606  }
607 
608  // Is there anything else on the header line?
609  std::string extras;
610  node_header_line >> extras;
611  if (extras == "BIN")
612  {
613  mFilesAreBinary = true;
614  mNodeFileDataStart = mNodesFile.tellg(); // Record the position of the first byte after the header.
615  mNodeItemWidth = SPACE_DIM * sizeof(double);
616 
617  // We enforce that all binary files (written by Chaste) are indexed from zero
618  mIndexFromZero = true;
619  }
620  else
621  {
622  // #1621 - ncl files are only supported in binary read mode.
623  assert(!mNclFileAvailable);
624 
625  // Get the next line to see if it is indexed from zero or not
626  GetNextLineFromStream(mNodesFile, buffer);
627  std::stringstream node_first_line(buffer);
628  unsigned first_index;
629  node_first_line >> first_index;
630  assert(first_index == 0 || first_index == 1);
631  mIndexFromZero = (first_index == 0);
632 
633  // Close, reopen, skip header
634  mNodesFile.close();
635  OpenNodeFile();
636  GetNextLineFromStream(mNodesFile, buffer);
637  }
638 
639  /*
640  * Reading element file header
641  */
642  GetNextLineFromStream(mElementsFile, buffer);
643  std::stringstream element_header_line(buffer);
644 
645  unsigned extra_attributes = 0;
646 
647  if (ELEMENT_DIM == SPACE_DIM)
648  {
649  element_header_line >> mNumElements >> mNumElementNodes >> mNumElementAttributes;
650 
651  extra_attributes = mNumElementAttributes;
652 
653  // Is there anything else on the header line?
654  std::string element_extras;
655  element_header_line >> element_extras;
656  if (element_extras == "BIN")
657  {
658  // Double check for binaryness
659  assert (mFilesAreBinary);
660  }
661  else if (element_extras == "HEX")
662  {
663  mMeshIsHexahedral = true;
664  if (ELEMENT_DIM == 2)
665  {
666  mNodesPerElement = 4;
667  mNodesPerBoundaryElement = 2;
668  }
669  if (ELEMENT_DIM == 3)
670  {
671  mNodesPerElement = 8;
672  mNodesPerBoundaryElement = 4;
673  }
674  }
675  else
676  {
677  assert (element_extras == "");
678  }
679 
680  // The condition below allows the writer to cope with a NodesOnlyMesh
681  if (mNumElements != 0)
682  {
683  if (mNumElementNodes != mNodesPerElement)
684  {
685  EXCEPTION("Number of nodes per elem, " << mNumElementNodes << ", does not match "
686  << "expected number, " << mNodesPerElement << " (which is calculated given "
687  << "the order of elements chosen, " << mOrderOfElements << " (1=linear, 2=quadratics)");
688  }
689  }
690  }
691  else
692  {
693  // Note that .face files don't have the number of nodes in a face element in the header (its element dim +1)
694  element_header_line >> mNumElements >> mNumFaceAttributes;
695 
696  extra_attributes = mNumFaceAttributes;
697 
698  if (ELEMENT_DIM == 1 || ELEMENT_DIM == 2)
699  {
700  mNumElementAttributes = mNumFaceAttributes;
701  }
702 
703  // Is there anything else on the header line?
704  std::string element_extras;
705  element_header_line >> element_extras;
706  if (element_extras == "BIN")
707  {
708  // Double check for binaryness
709  assert (mFilesAreBinary);
710  }
711 
712  mNodesPerElement = ELEMENT_DIM+1;
713  }
714 
715  if (mFilesAreBinary)
716  {
717  mElementFileDataStart = mElementsFile.tellg(); // Record the position of the first byte after the header.
718  mElementItemWidth = mNodesPerElement*sizeof(unsigned) + extra_attributes*sizeof(double);
719  }
720 
721  /*
722  * Reading face/edge file header.
723  * The condition below allows the writer to cope with a NodesOnlyMesh.
724  */
725  if (mNumElements != 0)
726  {
727  if (ELEMENT_DIM == 1)
728  {
729  GetOneDimBoundary();
730  mNumFaces = mOneDimBoundary.size();
731  }
732  else
733  {
734  GetNextLineFromStream(mFacesFile, buffer);
735  std::stringstream face_header_line(buffer);
736 
737  face_header_line >> mNumFaces >> mNumFaceAttributes;
738  assert(mNumFaceAttributes==0 || mNumFaceAttributes==1);
739 
740  /*
741  * If mNumFaceAttributes=1 then loop over and set mNumFaces to be
742  * the number of faces which are marked as boundary faces.
743  * Double check for binaryness.
744  */
745  std::string face_extras;
746  face_header_line >> face_extras;
747  assert (mFilesAreBinary == (face_extras == "BIN"));
748  if (mNumFaceAttributes==1)
749  {
750  unsigned num_boundary_faces = 0;
751  bool end_of_file=false;
752  while (!end_of_file)
753  {
754  try
755  {
756  GetNextFaceData();
757  num_boundary_faces++;
758  }
759  catch(Exception& e)
760  {
761  if (mEofException)
762  {
763  end_of_file = true;
764  }
765  else
766  {
767  throw e;
768  }
769  }
770  }
771  mNumFaces = num_boundary_faces;
772 
775  // if (mNumFaces==0)
776  // {
777  // EXCEPTION("No boundary elements found. NOTE: elements in face/edge file with an attribute value of 0 are considered to be internal (non-boundary) elements");
778  // }
779 
780  // close the file, reopen, and skip the header again
781  mFacesFile.close();
782  mFacesFile.clear(); // Older versions of gcc don't explicitly reset "fail" and "eof" flags in std::ifstream after calling close()
783  OpenFacesFile();
784  GetNextLineFromStream(mFacesFile, buffer);
785  mFacesRead = 0;
786  mBoundaryFacesRead = 0;
787  }
788  }
789  }
790 
791  if (mFilesAreBinary)
792  {
793  mFaceFileDataStart = mFacesFile.tellg(); // Record the position of the first byte after the header.
794  mFaceItemWidth = ELEMENT_DIM*sizeof(unsigned) + mNumFaceAttributes*sizeof(double);
795  }
796 
797  /*
798  * Read NCL file (if one is available)
799  */
800  if (mNclFileAvailable)
801  {
802  GetNextLineFromStream(mNclFile, buffer);
803  std::stringstream ncl_header_line(buffer);
804  unsigned num_nodes_in_file;
805  ncl_header_line >> num_nodes_in_file >> mMaxContainingElements;
806 
807  if (mNumNodes != num_nodes_in_file)
808  {
809  EXCEPTION("NCL file does not contain the correct number of nodes for mesh");
810  }
811 
812  mNclFileDataStart = mNclFile.tellg(); // Record the position of the first byte after the header
813  mNclItemWidth = mMaxContainingElements * sizeof(unsigned);
814  }
815 
816  /*
817  * Read cable file (if one is available)
818  */
819  if (mCableElementsFile.is_open())
820  {
821  GetNextLineFromStream(mCableElementsFile, buffer);
822  std::stringstream cable_header_line(buffer);
823  unsigned num_nodes_per_cable_element;
824  cable_header_line >> mNumCableElements >> num_nodes_per_cable_element >> mNumCableElementAttributes;
825  assert(num_nodes_per_cable_element == 2u);
826  mCableElementsRead = 0u;
827  }
828 }
829 
830 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
832 {
833  mNodesFile.close();
834  mElementsFile.close();
835  mFacesFile.close();
836  mNclFile.close();
837  mCableElementsFile.close();
838 }
839 
840 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
841 void TrianglesMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNextLineFromStream(std::ifstream& rFileStream, std::string& rRawLine)
842 {
843  bool line_is_blank;
844  mEofException = false;
845  do
846  {
847  getline(rFileStream, rRawLine, '\n');
848  if (rFileStream.eof())
849  {
850  mEofException = true;
851  EXCEPTION("File contains incomplete data: unexpected end of file.");
852  }
853 
854  // Get rid of any comment
855  rRawLine = rRawLine.substr(0, rRawLine.find('#',0));
856 
857  line_is_blank = (rRawLine.find_first_not_of(" \t",0) == std::string::npos);
858  }
859  while (line_is_blank);
860 }
861 
862 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
863 template<class T_DATA>
864 void TrianglesMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNextItemFromStream(std::ifstream& rFileStream, unsigned expectedItemNumber,
865  std::vector<T_DATA>& rDataPacket, const unsigned& rNumAttributes, std::vector<double>& rAttributes)
866 {
867  if (mFilesAreBinary)
868  {
869  if (!rDataPacket.empty()) // Avoid MSVC 10 assertion
870  {
871  rFileStream.read((char*)&rDataPacket[0], rDataPacket.size()*sizeof(T_DATA));
872  }
873  if (rNumAttributes > 0)
874  {
875  for (unsigned i = 0; i < rNumAttributes; i++)
876  {
877  double attribute;
878  rFileStream.read((char*) &attribute, sizeof(double));
879  rAttributes.push_back(attribute);
880  }
881  }
882  }
883  else
884  {
885  std::string buffer;
886  GetNextLineFromStream(rFileStream, buffer);
887  std::stringstream buffer_stream(buffer);
888 
889  unsigned item_index;
890  buffer_stream >> item_index;
891 
892  // If we are indexing from zero our expected item number is one larger
893  expectedItemNumber += mIndexFromZero ? 0 : 1;
894 
895  if (item_index != expectedItemNumber)
896  {
897  if (!mIndexFromZero)
898  {
899  // To fix the exception message to agree with file format
900  expectedItemNumber--;
901  }
902  EXCEPTION("Data for item " << expectedItemNumber << " missing");
903  }
904 
905  for (unsigned i=0; i<rDataPacket.size(); i++)
906  {
907  buffer_stream >> rDataPacket[i];
908  }
909 
910  if (rNumAttributes > 0)
911  {
912  for (unsigned i = 0; i < rNumAttributes; i++)
913  {
914  double attribute;
915  buffer_stream >> attribute;
916  if (buffer_stream.fail())
917  {
918  EXCEPTION("Error in reading attribute index " << i << " (out of " << rNumAttributes << ") in one of the files in " << mFilesBaseName);
919  }
920  rAttributes.push_back(attribute);
921  }
922  }
923  }
924 }
925 
926 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
928 {
929  return mFilesBaseName;
930 }
931 
932 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
934 {
935  assert(ELEMENT_DIM == 1); // LCOV_EXCL_LINE
936  mNumFaceAttributes = 0;
937  if (!mOneDimBoundary.empty())
938  {
939  // We have already read this and have reset the reader (probably from the mesh class)
940  return;
941  }
942  std::vector<unsigned> node_indices(2);
943  std::vector<double> dummy_attribute; // unused
944 
945  // Count how many times we see each node
946  std::vector<unsigned> node_count(mNumNodes); // Covers the case if it's indexed from 1
947  for (unsigned element_index=0; element_index<mNumElements;element_index++)
948  {
949  GetNextItemFromStream(mElementsFile, element_index, node_indices, mNumElementAttributes, dummy_attribute);
950  if (!mIndexFromZero)
951  {
952  // Adjust so we are indexing from zero
953  node_indices[0]--;
954  node_indices[1]--;
955  }
956  node_count[node_indices[0]]++;
957  node_count[node_indices[1]]++;
958  }
959 
960  // Find the ones which are terminals (only one mention)
961  for (unsigned node_index=0; node_index<mNumNodes;node_index++)
962  {
963  if (node_count[node_index] == 1u)
964  {
965  mOneDimBoundary.push_back(node_index);
966  }
967  }
968 
969  // Close the file, reopen, and skip the header again
970  mElementsFile.close();
971  mElementsFile.clear(); // Older versions of gcc don't explicitly reset "fail" and "eof" flags in std::ifstream after calling close()
972  OpenElementsFile();
973  std::string buffer;
974  GetNextLineFromStream(mElementsFile, buffer);
975 }
976 
977 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
979 {
980  if (!mIndexFromZero) // If node indices do not start at zero move them all down one so they do
981  {
982  for (unsigned i=0; i<rNodeIndices.size(); i++)
983  {
984  rNodeIndices[i]--;
985  }
986  }
987 }
988 
989 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
991 {
992  return mFilesAreBinary;
993 }
994 
995 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
997 {
998  return mNclFileAvailable;
999 }
1000 
1001 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1003 {
1004  mNodeFileReadBuffer = new char[bufferSize];
1005  mElementFileReadBuffer = new char[bufferSize];
1006  mFaceFileReadBuffer = new char[bufferSize];
1007 
1008  mNodesFile.rdbuf()->pubsetbuf(mNodeFileReadBuffer, bufferSize);
1009  mElementsFile.rdbuf()->pubsetbuf(mElementFileReadBuffer, bufferSize);
1010  mFacesFile.rdbuf()->pubsetbuf(mFaceFileReadBuffer, bufferSize);
1011 }
1012 
1013 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1014 void TrianglesMeshReader<ELEMENT_DIM, SPACE_DIM>::SetNodePermutation(std::vector<unsigned>& rPermutationVector)
1015 {
1016  if (!mFilesAreBinary)
1017  {
1018  // It would be too inefficient otherwise...
1019  EXCEPTION("Permuted read can only be used with binary files since it requires random access to the node file.");
1020  }
1021 
1022  mNodePermutationDefined = true;
1023  mPermutationVector = rPermutationVector;
1024  mInversePermutationVector.resize(mPermutationVector.size());
1025  for (unsigned index=0; index<mPermutationVector.size(); index++)
1026  {
1027  mInversePermutationVector[mPermutationVector[index]]=index;
1028  }
1029 }
1030 
1031 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1033 {
1034  return(mNodePermutationDefined);
1035 }
1036 
1040 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1042 {
1043  return mPermutationVector;
1044 }
1045 
1046 // Explicit instantiation
1047 template class TrianglesMeshReader<0,1>;
1048 template class TrianglesMeshReader<1,1>;
1049 template class TrianglesMeshReader<1,2>;
1050 template class TrianglesMeshReader<1,3>;
1051 template class TrianglesMeshReader<2,2>;
1052 template class TrianglesMeshReader<2,3>;
1053 template class TrianglesMeshReader<3,3>;
1054 
1055 
1060 template void TrianglesMeshReader<0,1>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1061 template void TrianglesMeshReader<0,1>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1062 template void TrianglesMeshReader<1,1>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1063 template void TrianglesMeshReader<1,1>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1064 template void TrianglesMeshReader<1,2>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1065 template void TrianglesMeshReader<1,2>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1066 template void TrianglesMeshReader<1,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1067 template void TrianglesMeshReader<1,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1068 template void TrianglesMeshReader<2,2>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1069 template void TrianglesMeshReader<2,2>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1070 template void TrianglesMeshReader<2,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1071 template void TrianglesMeshReader<2,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
1072 template void TrianglesMeshReader<3,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<unsigned>&, const unsigned&, std::vector<double>&);
1073 template void TrianglesMeshReader<3,3>::GetNextItemFromStream(std::ifstream&, unsigned, std::vector<double> &, const unsigned&, std::vector<double>&);
ElementData GetNextCableElementData()
unsigned GetNumFaceAttributes() const
unsigned GetNumNodes() const
unsigned ContainingElement
void GetNextLineFromStream(std::ifstream &rFileStream, std::string &rRawLine)
#define EXCEPTION(message)
Definition: Exception.hpp:143
TrianglesMeshReader(std::string pathBaseName, unsigned orderOfElements=1, unsigned orderOfBoundaryElements=1, bool readContainingElementsForBoundaryElements=false)
std::vector< unsigned > GetContainingElementIndices(unsigned index)
unsigned GetNumElementAttributes() const
unsigned GetNumCableElementAttributes() const
unsigned GetNumCableElements() const
unsigned GetNumElements() const
std::vector< unsigned > NodeIndices
ElementData GetNextElementData()
void SetNodePermutation(std::vector< unsigned > &rPermutationVector)
void SetReadBufferSize(unsigned bufferSize)
std::vector< double > GetNextNode()
void EnsureIndexingFromZero(std::vector< unsigned > &rNodeIndices)
ElementData GetFaceData(unsigned index)
unsigned GetNumFaces() const
const std::vector< unsigned > & rGetNodePermutation()
std::string GetMeshFileBaseName()
std::vector< double > GetNodeAttributes()
void GetNextItemFromStream(std::ifstream &rFileStream, unsigned expectedItemNumber, std::vector< T_DATA > &rDataPacket, const unsigned &rNumAttributes, std::vector< double > &rAttributes)
std::vector< double > GetNode(unsigned index)
ElementData GetElementData(unsigned index)