Chaste  Release::2018.1
MixedDimensionMesh.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 
36 #include "MixedDimensionMesh.hpp"
37 #include "Exception.hpp"
38 
39 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
41  : DistributedTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::DistributedTetrahedralMesh(partitioningMethod)
42 {
43 }
44 
45 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
47 {
48  for (unsigned i=0; i<mCableElements.size(); i++)
49  {
50  delete mCableElements[i];
51  }
52 }
53 
54 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
56 {
58  //Note that the above method may permute the node (for a parMETIS partition) after construction
59  //We have to convert to the permuted form first
60 
61  // Add cable elements
62  mNumCableElements = rMeshReader.GetNumCableElements();
63  //this->mCableElements.reserve(mNumCableElements);
64 
65  for (unsigned element_index=0; element_index < mNumCableElements; element_index++)
66  {
67  ElementData element_data = rMeshReader.GetNextCableElementData();
68  //Convert the node indices from the original to the permuted
69  if (!this->mNodePermutation.empty())
70  {
71  for (unsigned j=0; j<2; j++) // cables are always 1d
72  {
73  element_data.NodeIndices[j] = this->mNodePermutation[ element_data.NodeIndices[j] ];
74  }
75  }
76 
77  //Determine if we own any nodes on this cable element
78  bool node_owned = false;
79  for (unsigned j=0; j<2; j++) // cables are always 1d
80  {
81  try
82  {
83  this->SolveNodeMapping(element_data.NodeIndices[j]);
84  node_owned = true;
85  break;
86  }
87  catch (Exception &)
88  {
89  //We deal with non-owned nodes in the next part
90  }
91  }
92 
93  //If we don't locally own either node, then we don't construct the cable
94  if (node_owned)
95  {
96  std::vector<Node<SPACE_DIM>*> nodes;
97  nodes.reserve(2u);
98 
99  for (unsigned j=0; j<2; j++) // cables are always 1d
100  {
101  //Note that if we own one node on a cable element then we are likely to own the other.
102  //If not, we are likely to have a halo.
103  //If not, (free-running Purkinje with monodomain mesh?), then this will terminate.
104  try
105  {
106  nodes.push_back(this->GetNodeOrHaloNode(element_data.NodeIndices[j]) );
107  }
108  // LCOV_EXCL_START
109  catch (Exception&)
110  {
112  }
113  // LCOV_EXCL_STOP
114  }
115 
116  Element<1u, SPACE_DIM>* p_element = new Element<1u,SPACE_DIM>(element_index, nodes, false);
117  RegisterCableElement(element_index);
118  this->mCableElements.push_back(p_element);
119  for (unsigned node_index=0; node_index<p_element->GetNumNodes(); ++node_index)
120  {
121  mNodeToCablesMapping.insert(std::pair<Node<SPACE_DIM>*, Element<1u, SPACE_DIM>*>(
122  p_element->GetNode(node_index), p_element));
123  }
124 
125  if (rMeshReader.GetNumCableElementAttributes() > 0)
126  {
127  assert(rMeshReader.GetNumCableElementAttributes() == 1);
128  p_element->SetAttribute(element_data.AttributeValue);
129  }
130  }
131  }
132 
133  rMeshReader.Reset();
134 }
135 
136 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
138 {
139  mCableElementsMapping[index] = this->mCableElements.size();
140 }
141 
142 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
144 {
145  return mNumCableElements;
146 }
147 
148 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
150 {
151  return mCableElements.size();
152 }
153 
154 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
156 {
157  std::map<unsigned, unsigned>::const_iterator element_position = mCableElementsMapping.find(globalElementIndex);
158 
159  if (element_position == mCableElementsMapping.end())
160  {
161  EXCEPTION("Requested cable element " << globalElementIndex << " does not belong to processor " << PetscTools::GetMyRank());
162  }
163 
164  unsigned index = element_position->second;
165 
166  return mCableElements[index];
167 }
168 
169 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
171 {
172  //This should not throw in the distributed parallel case
173  try
174  {
175  unsigned tie_break_index = this->GetCableElement(elementIndex)->GetNodeGlobalIndex(0);
176 
177  //if it is in my range
178  if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(tie_break_index))
179  {
180  return true;
181  }
182  else
183  {
184  return false;
185  }
186  }
187  catch (Exception &)
188  {
189  //We don't own this cable element
190  return false;
191  }
192 }
193 
194 
195 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
197 {
198  return mNodeToCablesMapping.equal_range(pNode);
199 }
200 
201 
202 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
204 {
205  return mCableElements.begin();
206 }
207 
208 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
210 {
211  return mCableElements.end();
212 }
213 
214 // Explicit instantiation
215 template class MixedDimensionMesh<1,1>;
216 template class MixedDimensionMesh<1,2>;
217 template class MixedDimensionMesh<1,3>;
218 template class MixedDimensionMesh<2,2>;
219 template class MixedDimensionMesh<2,3>;
220 template class MixedDimensionMesh<3,3>;
221 
222 // Serialization for Boost >= 1.36
CableElementIterator GetCableElementIteratorEnd() const
void SetAttribute(double attribute)
unsigned GetNumLocalCableElements() const
Definition: Node.hpp:58
void RegisterCableElement(unsigned index)
std::pair< NodeCableIterator, NodeCableIterator > CableRangeAtNode
virtual unsigned GetNumCableElementAttributes() const
#define EXCEPTION(message)
Definition: Exception.hpp:143
virtual void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
#define NEVER_REACHED
Definition: Exception.hpp:206
std::vector< unsigned > NodeIndices
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
virtual ElementData GetNextCableElementData()
CableElementIterator GetCableElementIteratorBegin() const
virtual void Reset()=0
unsigned GetNumNodes() const
unsigned GetNumCableElements() const
bool CalculateDesignatedOwnershipOfCableElement(unsigned globalElementIndex)
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
virtual unsigned GetNumCableElements() const
MixedDimensionMesh(DistributedTetrahedralMeshPartitionType::type partitioningMethod=DistributedTetrahedralMeshPartitionType::METIS_LIBRARY)
CableRangeAtNode GetCablesAtNode(const Node< SPACE_DIM > *pNode)
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
Element< 1u, SPACE_DIM > * GetCableElement(unsigned globalElementIndex) const
std::vector< Element< 1, SPACE_DIM > * >::const_iterator CableElementIterator