Chaste  Release::2018.1
AbstractTetrahedralMesh.hpp
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 
36 #ifndef ABSTRACTTETRAHEDRALMESH_HPP_
37 #define ABSTRACTTETRAHEDRALMESH_HPP_
38 
39 #include "ChasteSerialization.hpp"
40 #include "ClassIsAbstract.hpp"
42 #include <boost/serialization/vector.hpp>
43 #include <boost/serialization/base_object.hpp>
44 #include <boost/serialization/split_member.hpp>
45 
46 #include <vector>
47 #include <string>
48 #include <cassert>
49 #include <boost/foreach.hpp>
50 
51 #include "AbstractMesh.hpp"
52 #include "BoundaryElement.hpp"
53 #include "Element.hpp"
54 #include "GenericMeshReader.hpp"
55 #include "AbstractMeshReader.hpp"
56 #include "TrianglesMeshReader.hpp"
57 #include "TrianglesMeshWriter.hpp"
58 #include "ArchiveLocationInfo.hpp"
59 #include "FileFinder.hpp"
60 
61 
63 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
65 
70 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
71 class AbstractTetrahedralMesh : public AbstractMesh<ELEMENT_DIM, SPACE_DIM>
72 {
73  friend class AbstractConductivityTensors<ELEMENT_DIM, SPACE_DIM>; //A class which needs a global to local element mapping
74  friend class CentroidWriter; //A test class which needs access to mElements in order to check that local/global indices match
75 protected:
80 
81 private:
91  virtual unsigned SolveElementMapping(unsigned index) const = 0;
92 
101  virtual unsigned SolveBoundaryElementMapping(unsigned index) const = 0;
102 
104  friend class boost::serialization::access;
105 
117  template<class Archive>
118  void save(Archive & archive, const unsigned int version) const
119  {
120  archive & boost::serialization::base_object<AbstractMesh<ELEMENT_DIM,SPACE_DIM> >(*this);
121  archive & mMeshIsLinear;
122  // Create a mesh writer pointing to the correct file and directory
125  false);
126  // Binary meshes have similar content to the original Triangle/Tetgen format, but take up less space on disk
127  mesh_writer.SetWriteFilesAsBinary();
128 
129  // Archive the mesh permutation, so we can just copy the original mesh files whenever possible
130  bool permutation_available = (this->rGetNodePermutation().size() != 0);
131  archive & permutation_available;
132 
133  if (permutation_available)
134  {
135  const std::vector<unsigned>& rPermutation = this->rGetNodePermutation();
136  archive & rPermutation;
137  }
138 
139  if (!this->IsMeshOnDisk() || this->mMeshChangesDuringSimulation)
140  {
141  mesh_writer.WriteFilesUsingMesh(*(const_cast<AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>*>(this)));
142  }
143  else
144  {
145  unsigned order_of_element = (mMeshIsLinear?1:2);
146  unsigned& order_of_boundary_element = order_of_element;
147 
148  // Mesh in disc, copy it to the archiving folder
149  std::string original_file=this->GetMeshFileBaseName();
150  std::shared_ptr<AbstractMeshReader<ELEMENT_DIM, SPACE_DIM> > p_original_mesh_reader
151  = GenericMeshReader<ELEMENT_DIM, SPACE_DIM>(original_file, order_of_element, order_of_boundary_element);
152 
153  if (p_original_mesh_reader->IsFileFormatBinary())
154  {
155  // Mesh is in binary format, we can just copy the files across ignoring the mesh reader
156  if (PetscTools::AmMaster())
157  {
158  FileFinder mesh_base(this->GetMeshFileBaseName());
159  FileFinder mesh_folder = mesh_base.GetParent();
160  std::string mesh_leaf_name = mesh_base.GetLeafNameNoExtension();
161  std::vector<FileFinder> mesh_files = mesh_folder.FindMatches(mesh_leaf_name + ".*");
163  BOOST_FOREACH(const FileFinder& r_mesh_file, mesh_files)
164  {
166  dest_dir);
167  ABORT_IF_THROWS(r_mesh_file.CopyTo(dest_file));
168  }
169  }
170  }
171  else
172  {
173  // Mesh in text format, use the mesh writer to "binarise" it
174  mesh_writer.WriteFilesUsingMeshReaderAndMesh(*p_original_mesh_reader,
176  }
177  }
178 
179  // Make sure that the files are written before slave processes proceed
180  PetscTools::Barrier("AbstractTetrahedralMesh::save");
181  }
182 
189  template<class Archive>
190  void load(Archive & archive, const unsigned int version)
191  {
192  archive & boost::serialization::base_object<AbstractMesh<ELEMENT_DIM,SPACE_DIM> >(*this);
193  archive & mMeshIsLinear;
194 
195  bool permutation_available=false;
196  std::vector<unsigned> permutation;
197 
198  if (version > 0)
199  {
200  archive & permutation_available;
201 
202  if (permutation_available)
203  {
204  archive & permutation;
205  }
206  }
207 
208  // Store the DistributedVectorFactory loaded from the archive
210  this->mpDistributedVectorFactory = nullptr;
211 
212  // Check whether we're migrating, or if we can use the original partition for the mesh
213  DistributedVectorFactory* p_our_factory = nullptr;
214  if (p_factory)
215  {
216  p_our_factory = p_factory->GetOriginalFactory();
217  }
218  if (p_our_factory && p_our_factory->GetNumProcs() == p_factory->GetNumProcs())
219  {
220  // Specify the node distribution
221  this->SetDistributedVectorFactory(p_our_factory);
222  }
223  else
224  {
225  // Migrating; let the mesh re-partition if it likes
227  p_our_factory = nullptr;
228  }
229 
230  if (mMeshIsLinear)
231  {
232  // I am a linear mesh
234 
235  if (permutation_available)
236  {
237  mesh_reader.SetNodePermutation(permutation);
238  }
239 
240  this->ConstructFromMeshReader(mesh_reader);
241  }
242  else
243  {
244  // I am a quadratic mesh and need quadratic information from the reader
246  this->ConstructFromMeshReader(mesh_reader);
247  }
248 
249  // Make sure we're using the correct vector factory
250  if (p_factory)
251  {
252  if (!this->mpDistributedVectorFactory)
253  {
254  // If we're not using a DistributedTetrahedralMesh, ConstructFromMeshReader won't set
255  // this->mpDistributedVectorFactory.
256  this->mpDistributedVectorFactory = p_factory;
257  }
258  else
259  {
260  // We need to update p_factory to match this->mpDistributedVectorFactory, and then use
261  // p_factory, since the rest of the code (e.g. AbstractCardiacPde) will be using p_factory.
262  p_factory->SetFromFactory(this->mpDistributedVectorFactory);
263  if (p_our_factory != this->mpDistributedVectorFactory)
264  {
265  // Avoid memory leak
266  delete this->mpDistributedVectorFactory;
267  }
268  this->mpDistributedVectorFactory = p_factory;
269  }
270  }
271  }
272  BOOST_SERIALIZATION_SPLIT_MEMBER()
273 
274 protected: // Give access of these variables to subclasses
275 
277  std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> mElements;
278 
280  std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> mBoundaryElements;
281 
289  void SetElementOwnerships();
290 
291 public:
292 
294  // Iterators //
296 
298  typedef typename std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *>::const_iterator BoundaryElementIterator;
299 
301  class ElementIterator;
302 
308  inline ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true);
309 
313  inline ElementIterator GetElementIteratorEnd();
314 
316  // Methods //
318 
323 
327  virtual ~AbstractTetrahedralMesh();
328 
329 
333  virtual unsigned GetNumElements() const;
334 
338  virtual unsigned GetNumLocalElements() const;
339 
343  virtual unsigned GetNumBoundaryElements() const;
344 
348  virtual unsigned GetNumLocalBoundaryElements() const;
349 
353  unsigned GetNumAllElements() const;
354 
358  unsigned GetNumAllBoundaryElements() const;
359 
365  virtual unsigned GetNumCableElements() const;
366 
371  virtual unsigned GetNumVertices() const;
372 
379  virtual unsigned GetMaximumNodeIndex();
380 
387  Element<ELEMENT_DIM, SPACE_DIM>* GetElement(unsigned index) const;
388 
395  BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* GetBoundaryElement(unsigned index) const;
396 
397 
404  virtual void ConstructFromMeshReader(AbstractMeshReader<ELEMENT_DIM, SPACE_DIM>& rMeshReader)=0;
405 
414  void ConstructFromMesh(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rOtherMesh);
415 
416 
420  BoundaryElementIterator GetBoundaryElementIteratorBegin() const;
421 
426  BoundaryElementIterator GetBoundaryElementIteratorEnd() const;
427 
436  virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian,
437  double& rJacobianDeterminant,
438  c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian) const;
439 
448  virtual void GetWeightedDirectionForBoundaryElement(unsigned elementIndex,
449  c_vector<double, SPACE_DIM>& rWeightedDirection,
450  double& rJacobianDeterminant) const;
451 
459  void CheckOutwardNormals();
460 
473  virtual void ConstructLinearMesh(unsigned width);
474 
490  virtual void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true);
491 
503  virtual void ConstructCuboid(unsigned width, unsigned height, unsigned depth);
504 
515  void ConstructRegularSlabMesh(double spaceStep, double width, double height=0, double depth=0);
516 
533  void ConstructRegularSlabMeshWithDimensionSplit(unsigned dimension, double spaceStep, double width, double height=0, double depth=0);
534 
535 
543  virtual bool CalculateDesignatedOwnershipOfBoundaryElement( unsigned faceIndex );
544 
552  virtual bool CalculateDesignatedOwnershipOfElement( unsigned elementIndex );
553 
564 
571  virtual void GetHaloNodeIndices(std::vector<unsigned>& rHaloIndices) const;
572 
586  void CalculateNodeExchange( std::vector<std::vector<unsigned> >& rNodesToSendPerProcess,
587  std::vector<std::vector<unsigned> >& rNodesToReceivePerProcess);
588 
589 
598  virtual c_vector<double, 2> CalculateMinMaxEdgeLengths();
599 
613  unsigned GetContainingElementIndex(const ChastePoint<SPACE_DIM>& rTestPoint,
614  bool strict=false,
615  std::set<unsigned> testElements=std::set<unsigned>(),
616  bool onlyTryWithTestElements = false);
617 
623  unsigned GetNearestElementIndexFromTestElements(const ChastePoint<SPACE_DIM>& rTestPoint,
624  std::set<unsigned> testElements);
625 
627  // Nested classes //
629 
633  class ElementIterator
634  {
635  public:
641  inline Element<ELEMENT_DIM, SPACE_DIM>& operator*();
642 
647  inline Element<ELEMENT_DIM, SPACE_DIM>* operator->();
648 
654  inline bool operator!=(const typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator& rOther);
655 
660  inline ElementIterator& operator++();
661 
673  typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator elementIter,
674  bool skipDeletedElements=true);
675 
676  private:
678  AbstractTetrahedralMesh& mrMesh;
679 
681  typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator mElementIter;
682 
685 
690  inline bool IsAtEnd();
691 
696  inline bool IsAllowedElement();
697  };
698 };
699 
701 
702 namespace boost {
703 namespace serialization {
710 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
711 struct version<AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM> >
712 {
715 };
716 } // namespace serialization
717 } // namespace boost
718 
719 
721 // ElementIterator class implementation - most methods are inlined //
723 
724 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
726  bool skipDeletedElements)
727 {
728  return ElementIterator(*this, mElements.begin(), skipDeletedElements);
729 }
730 
731 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
733 {
734  return ElementIterator(*this, mElements.end());
735 }
736 
737 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
739 {
740  assert(!IsAtEnd());
741  return **mElementIter;
742 }
743 
744 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
746 {
747  assert(!IsAtEnd());
748  return *mElementIter;
749 }
750 
751 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
753 {
754  return mElementIter != rOther.mElementIter;
755 }
756 
757 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
759 {
760  do
761  {
762  ++mElementIter;
763  }
764  while (!IsAtEnd() && !IsAllowedElement());
765 
766  return (*this);
767 }
768 
769 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
772  typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator elementIter,
773  bool skipDeletedElements)
774  : mrMesh(rMesh),
775  mElementIter(elementIter),
776  mSkipDeletedElements(skipDeletedElements)
777 {
778  if (mrMesh.mElements.size() == 0)
779  {
780  // Cope with empty meshes
781  mElementIter = mrMesh.mElements.end();
782  }
783  else
784  {
785  // Make sure we start at an allowed element
786  if (mElementIter == mrMesh.mElements.begin() && !IsAllowedElement())
787  {
788  ++(*this);
789  }
790  }
791 }
792 
793 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
795 {
796  return mElementIter == mrMesh.mElements.end();
797 }
798 
799 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
801 {
802  return !(mSkipDeletedElements && (*this)->IsDeleted());
803 }
804 
805 #endif /*ABSTRACTTETRAHEDRALMESH_HPP_*/
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
virtual unsigned GetNumBoundaryElements() const
#define TEMPLATED_CLASS_IS_ABSTRACT_2_UNSIGNED(T)
std::vector< Element< ELEMENT_DIM, SPACE_DIM > * > mElements
bool operator!=(const typename AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM >::ElementIterator &rOther)
virtual unsigned GetNumLocalBoundaryElements() const
#define ABORT_IF_THROWS(block)
Definition: Exception.hpp:225
static void Barrier(const std::string callerId="")
Definition: PetscTools.cpp:134
BoundaryElementIterator GetBoundaryElementIteratorBegin() const
static std::string GetMeshFilename()
virtual unsigned GetMaximumNodeIndex()
virtual bool CalculateDesignatedOwnershipOfBoundaryElement(unsigned faceIndex)
const std::vector< unsigned > & rGetNodePermutation() const
void ConstructRegularSlabMesh(double spaceStep, double width, double height=0, double depth=0)
void load(Archive &archive, const unsigned int version)
virtual unsigned GetNumElements() const
virtual unsigned GetNumCableElements() const
unsigned CalculateMaximumNodeConnectivityPerProcess() const
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
static bool AmMaster()
Definition: PetscTools.cpp:120
FileFinder CopyTo(const FileFinder &rDest) const
Definition: FileFinder.cpp:308
std::vector< Element< ELEMENT_DIM, SPACE_DIM > * >::iterator mElementIter
Element< ELEMENT_DIM, SPACE_DIM > & operator*()
std::string GetMeshFileBaseName() const
virtual void GetHaloNodeIndices(std::vector< unsigned > &rHaloIndices) const
DistributedVectorFactory * GetOriginalFactory()
virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian) const
virtual unsigned GetNumLocalElements() const
virtual unsigned SolveElementMapping(unsigned index) const =0
Element< ELEMENT_DIM, SPACE_DIM > * operator->()
bool mMeshChangesDuringSimulation
DistributedVectorFactory * mpDistributedVectorFactory
virtual void GetWeightedDirectionForBoundaryElement(unsigned elementIndex, c_vector< double, SPACE_DIM > &rWeightedDirection, double &rJacobianDeterminant) const
ElementIterator(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, typename std::vector< Element< ELEMENT_DIM, SPACE_DIM > * >::iterator elementIter, bool skipDeletedElements=true)
std::vector< BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * > mBoundaryElements
void SetNodePermutation(std::vector< unsigned > &rPermutationVector)
std::string GetExtension() const
Definition: FileFinder.cpp:250
#define CHASTE_VERSION_CONTENT(N)
virtual unsigned GetNumVertices() const
void SetFromFactory(DistributedVectorFactory *pFactory)
unsigned GetContainingElementIndex(const ChastePoint< SPACE_DIM > &rTestPoint, bool strict=false, std::set< unsigned > testElements=std::set< unsigned >(), bool onlyTryWithTestElements=false)
virtual void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true)
std::vector< FileFinder > FindMatches(const std::string &rPattern) const
Definition: FileFinder.cpp:430
bool IsMeshOnDisk() const
BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * GetBoundaryElement(unsigned index) const
virtual void ConstructCuboid(unsigned width, unsigned height, unsigned depth)
virtual void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)=0
FileFinder GetParent() const
Definition: FileFinder.cpp:255
void CalculateNodeExchange(std::vector< std::vector< unsigned > > &rNodesToSendPerProcess, std::vector< std::vector< unsigned > > &rNodesToReceivePerProcess)
Forward declaration which is going to be used for friendship.
unsigned GetNumAllBoundaryElements() const
gcov doesn't like this file...
std::string GetLeafNameNoExtension() const
Definition: FileFinder.cpp:245
void save(Archive &archive, const unsigned int version) const
virtual unsigned SolveBoundaryElementMapping(unsigned index) const =0
unsigned GetNearestElementIndexFromTestElements(const ChastePoint< SPACE_DIM > &rTestPoint, std::set< unsigned > testElements)
virtual bool CalculateDesignatedOwnershipOfElement(unsigned elementIndex)
virtual void SetDistributedVectorFactory(DistributedVectorFactory *pFactory)
virtual void ConstructLinearMesh(unsigned width)
std::vector< BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * >::const_iterator BoundaryElementIterator
static std::string GetArchiveDirectory()
virtual void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
static std::string GetArchiveRelativePath()
virtual c_vector< double, 2 > CalculateMinMaxEdgeLengths()
void ConstructFromMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rOtherMesh)
void ConstructRegularSlabMeshWithDimensionSplit(unsigned dimension, double spaceStep, double width, double height=0, double depth=0)
void WriteFilesUsingMeshReaderAndMesh(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
BoundaryElementIterator GetBoundaryElementIteratorEnd() const