Chaste  Release::3.4
QuadraticMeshHelper.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 "QuadraticMeshHelper.hpp"
37 
38 #define SEEK_TO_CONTENT(methNameDirect, methNameIncrement, index) \
39  if (index > 0u) { \
40  if (rMeshReader.IsFileFormatBinary()) { \
41  rMeshReader.methNameDirect(index - 1u); \
42  } else { \
43  for (unsigned i=0; i<index-1u; ++i) { \
44  rMeshReader.methNameIncrement(); \
45  } } }
46 
47 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
48 void SeekToBoundaryElement(AbstractMeshReader<ELEMENT_DIM, SPACE_DIM>& rMeshReader,
49  unsigned boundaryElementIndex)
50 {
51  SEEK_TO_CONTENT(GetFaceData, GetNextFaceData, boundaryElementIndex);
52 }
53 
54 template<unsigned DIM>
56  AbstractMeshReader<DIM,DIM>* pMeshReader)
57 {
58  assert(pMesh);
59  assert(pMeshReader);
60 
61  if (pMesh->GetNumLocalElements() > 0u)
62  {
63  pMeshReader->Reset();
64 
65  // Create a set of element indices we own
66  std::set<unsigned> owned_element_indices;
68  iter != pMesh->GetElementIteratorEnd();
69  ++iter)
70  {
71  owned_element_indices.insert(iter->GetIndex());
72  }
73 
74  const std::vector<unsigned>& r_node_perm = pMesh->rGetNodePermutation();
75 
76  // Add the extra nodes (1 extra node in 1D, 3 in 2D, 6 in 3D) to the element data
77  for (typename AbstractMeshReader<DIM,DIM>::ElementIterator iter = pMeshReader->GetElementIteratorBegin(owned_element_indices);
78  iter != pMeshReader->GetElementIteratorEnd();
79  ++iter)
80  {
81  std::vector<unsigned> nodes = iter->NodeIndices;
82  assert(nodes.size()==(DIM+1)*(DIM+2)/2);
83  Element<DIM,DIM>* p_element = pMesh->GetElement(iter.GetIndex());
84  assert(p_element->GetNumNodes()==DIM+1); // Element is initially linear
85 
86  // Add extra nodes to make it a quad element
87  for (unsigned j=DIM+1; j<(DIM+1)*(DIM+2)/2; j++)
88  {
89  unsigned node_index = nodes[j];
90  if (!r_node_perm.empty())
91  {
92  node_index = r_node_perm[node_index];
93  }
94  Node<DIM>* p_node = pMesh->GetNodeOrHaloNode(node_index);
95  p_element->AddNode(p_node);
96  p_node->AddElement(p_element->GetIndex());
97  p_node->MarkAsInternal();
98  }
99  }
100  }
101 }
102 
103 template<unsigned DIM>
105  AbstractMeshReader<DIM,DIM>* pMeshReader)
106 {
107  assert(pMesh);
108  assert(pMeshReader);
109  // We only have boundary elements in 2d or 3d
110  if (DIM > 1 && pMesh->GetNumLocalBoundaryElements() > 0u)
111  {
112  // If the data is on disk our job is easy
113  if (pMeshReader->GetOrderOfBoundaryElements() == 2u)
114  {
115  // The work should have been done in the linear constructor, but let's check
116  // that the first face has more than DIM nodes.
117  assert((*pMesh->GetBoundaryElementIteratorBegin())->GetNumNodes()==DIM*(DIM+1)/2);
118  return;
119  }
120  else
121  {
122  AddNodesToBoundaryElements(pMesh, pMeshReader);
123  }
124  }
125 }
126 
127 template<unsigned DIM>
129  AbstractMeshReader<DIM,DIM>* pMeshReader)
130  {
131  // Loop over all boundary elements, find the equivalent face from all
132  // the elements, and add the extra nodes to the boundary element
133  bool boundary_element_file_has_containing_element_info = false;
134 
135  if (pMeshReader)
136  {
137  boundary_element_file_has_containing_element_info = pMeshReader->GetReadContainingElementOfBoundaryElement();
138  }
139 
140  if (DIM > 1)
141  {
142  if (boundary_element_file_has_containing_element_info)
143  {
144  pMeshReader->Reset();
145  }
146 
148  // we may need to skip through the boundary element file searching for containing elements hints.
149  // This counter keeps track of our position in the file.
150  unsigned next_face_on_file = 0u;
151 
154  iter != pMesh->GetBoundaryElementIteratorEnd();
155  ++iter)
156  {
157 
158  // collect the nodes of this boundary element in a set
159  std::set<unsigned> boundary_element_node_indices;
160  for (unsigned i=0; i<DIM; i++)
161  {
162  boundary_element_node_indices.insert( (*iter)->GetNodeGlobalIndex(i) );
163  }
164 
165  bool found_this_boundary_element = false;
166 
167  // Loop over elements surrounding this face, then loop over each face of that element, and see if it matches
168  // this boundary element.
169  // Note, if we know what elem it should be in (boundary_element_file_has_containing_element_info==true)
170  // we will reset elem_index immediately (below)
171  Node<DIM>* p_representative_node = (*iter)->GetNode(0);
172  for (typename Node<DIM>::ContainingElementIterator element_iter = p_representative_node->ContainingElementsBegin();
173  element_iter != p_representative_node->ContainingElementsEnd();
174  ++element_iter)
175  {
176  unsigned elem_index = *element_iter;
177 
178  // We know what elem it should be in (but we'll still check the node indices match in case)
179  if (boundary_element_file_has_containing_element_info)
180  {
181  unsigned face_index = (*iter)->GetIndex();
183  do
184  {
185  elem_index = pMeshReader->GetNextFaceData().ContainingElement;
186  next_face_on_file++;
187  }
188  while (face_index >= next_face_on_file);
189  }
190 
191  Element<DIM,DIM>* p_element = pMesh->GetElement(elem_index);
192 
193  // For each element, loop over faces (the opposites to a node)
194  for (unsigned face=0; face<DIM+1; face++)
195  {
196  // Collect the node indices for this face
197  std::set<unsigned> node_indices;
198  for (unsigned local_node_index=0; local_node_index<DIM+1; local_node_index++)
199  {
200  if (local_node_index != face)
201  {
202  node_indices.insert( p_element->GetNodeGlobalIndex(local_node_index) );
203  }
204  }
205 
206  assert(node_indices.size()==DIM);
207 
208  // See if this face matches the boundary element, and add internal nodes if so
209  if (node_indices == boundary_element_node_indices)
210  {
211  QuadraticMeshHelper<DIM>::AddExtraBoundaryNodes(pMesh, *iter, p_element, face);
212 
213  found_this_boundary_element = true;
214  break;
215  }
216  }
217 
218  // If the containing element info was given, we must have found the face first time
219  if (boundary_element_file_has_containing_element_info && !found_this_boundary_element)
220  {
221  #define COVERAGE_IGNORE
222  //std::cout << (*iter)->GetIndex() << " " << pMeshReader->GetNextFaceData().ContainingElement << "\n";
223  EXCEPTION("Boundary element " << (*iter)->GetIndex()
224  << "wasn't found in the containing element given for it "
225  << elem_index);
226  #undef COVERAGE_IGNORE
227  }
228 
229  if (found_this_boundary_element)
230  {
231  break;
232  }
233  }
234 
235  if (!found_this_boundary_element)
236  {
237  #define COVERAGE_IGNORE
238  EXCEPTION("Unable to find a face of an element which matches one of the boundary elements");
239  #undef COVERAGE_IGNORE
240  }
241  }
242  }
243 }
244 
245 template<unsigned DIM>
247 {
248 #ifndef NDEBUG
249  unsigned expected_num_nodes = DIM*(DIM+1)/2;
252  iter != pMesh->GetBoundaryElementIteratorEnd();
253  ++iter)
254  {
255  assert((*iter)->GetNumNodes() == expected_num_nodes);
256  }
257 #endif
258 }
259 
260 template<unsigned DIM>
262  BoundaryElement<DIM-1,DIM>* pBoundaryElement,
263  Node<DIM>* pNode)
264 {
265  assert(DIM > 1);
266 
267  // Add node to the boundary node list
268  if (!pNode->IsBoundaryNode())
269  {
270  pNode->SetAsBoundaryNode();
271  pMesh->mBoundaryNodes.push_back(pNode);
272  }
273  // Add it to the boundary element
274  pBoundaryElement->AddNode(pNode);
275 }
276 
277 template<unsigned DIM>
279  BoundaryElement<DIM-1,DIM>* pBoundaryElement,
280  Element<DIM,DIM>* pElement,
281  unsigned internalNode)
282 {
283  assert(DIM > 1);
284  assert(internalNode >= DIM+1);
285  assert(internalNode < (DIM+1)*(DIM+2)/2);
286  Node<DIM>* p_internal_node = pElement->GetNode(internalNode);
287  AddNodeToBoundaryElement(pMesh, pBoundaryElement, p_internal_node);
288 }
289 
290 template<unsigned DIM>
292  BoundaryElement<DIM-1,DIM>* pBoundaryElement,
293  Element<DIM,DIM>* pElement,
294  unsigned nodeIndexOppositeToFace)
295 {
296  assert(DIM!=1);
297  if (DIM==2)
298  {
299  assert(nodeIndexOppositeToFace<3);
300  // the single internal node of the element's face will be numbered 'face+3'
301  AddNodeToBoundaryElement(pMesh, pBoundaryElement, pElement, nodeIndexOppositeToFace+3);
302  }
303  else
304  {
305  assert(DIM==3);
306 
307  unsigned b_elem_n0 = pBoundaryElement->GetNodeGlobalIndex(0);
308  unsigned b_elem_n1 = pBoundaryElement->GetNodeGlobalIndex(1);
309 
310  unsigned offset;
311  bool reverse;
312 
313  if (nodeIndexOppositeToFace==0)
314  {
315  // face opposite to node 0 = {1,2,3}, with corresponding internals {9,8,5}
316  HelperMethod1(b_elem_n0, b_elem_n1, pElement, 1, 2, 3, offset, reverse);
317  HelperMethod2(pMesh, pBoundaryElement, pElement, 9, 8, 5, offset, reverse);
318  }
319  else if (nodeIndexOppositeToFace==1)
320  {
321  // face opposite to node 1 = {2,0,3}, with corresponding internals {7,9,6}
322  HelperMethod1(b_elem_n0, b_elem_n1, pElement, 2, 0, 3, offset, reverse);
323  HelperMethod2(pMesh, pBoundaryElement, pElement, 7, 9, 6, offset, reverse);
324  }
325  else if (nodeIndexOppositeToFace==2)
326  {
327  // face opposite to node 2 = {0,1,3}, with corresponding internals {8,7,4}
328  HelperMethod1(b_elem_n0, b_elem_n1, pElement, 0, 1, 3, offset, reverse);
329  HelperMethod2(pMesh, pBoundaryElement, pElement, 8, 7, 4, offset, reverse);
330  }
331  else
332  {
333  assert(nodeIndexOppositeToFace==3);
334  // face opposite to node 3 = {0,1,2}, with corresponding internals {5,6,4}
335  HelperMethod1(b_elem_n0, b_elem_n1, pElement, 0, 1, 2, offset, reverse);
336  HelperMethod2(pMesh, pBoundaryElement, pElement, 5, 6, 4, offset, reverse);
337  }
338  }
339 }
340 
342 // two unpleasant helper methods for AddExtraBoundaryNodes()
344 
345 #define COVERAGE_IGNORE
346 template<unsigned DIM>
347 void QuadraticMeshHelper<DIM>::HelperMethod1(unsigned boundaryElemNode0, unsigned boundaryElemNode1,
348  Element<DIM,DIM>* pElement,
349  unsigned node0, unsigned node1, unsigned node2,
350  unsigned& rOffset,
351  bool& rReverse)
352 {
353  if (pElement->GetNodeGlobalIndex(node0)==boundaryElemNode0)
354  {
355  rOffset = 0;
356  if (pElement->GetNodeGlobalIndex(node1)==boundaryElemNode1)
357  {
358  rReverse = false;
359  }
360  else
361  {
362  assert(pElement->GetNodeGlobalIndex(node2)==boundaryElemNode1);
363  rReverse = true;
364  }
365  }
366  else if (pElement->GetNodeGlobalIndex(node1)==boundaryElemNode0)
367  {
368  rOffset = 1;
369  if (pElement->GetNodeGlobalIndex(node2)==boundaryElemNode1)
370  {
371  rReverse = false;
372  }
373  else
374  {
375  assert(pElement->GetNodeGlobalIndex(node0)==boundaryElemNode1);
376  rReverse = true;
377  }
378  }
379  else
380  {
381  assert(pElement->GetNodeGlobalIndex(node2)==boundaryElemNode0);
382  rOffset = 2;
383  if (pElement->GetNodeGlobalIndex(node0)==boundaryElemNode1)
384  {
385  rReverse = false;
386  }
387  else
388  {
389  assert(pElement->GetNodeGlobalIndex(node1)==boundaryElemNode1);
390  rReverse = true;
391  }
392  }
393 }
394 #undef COVERAGE_IGNORE
395 
396 
397 #define COVERAGE_IGNORE
398 template<unsigned DIM>
400  BoundaryElement<DIM-1,DIM>* pBoundaryElement,
401  Element<DIM,DIM>* pElement,
402  unsigned internalNode0, unsigned internalNode1, unsigned internalNode2,
403  unsigned offset,
404  bool reverse)
405 {
406  if (offset==1)
407  {
408  unsigned temp = internalNode0;
409  internalNode0 = internalNode1;
410  internalNode1 = internalNode2;
411  internalNode2 = temp;
412  }
413  else if (offset == 2)
414  {
415  unsigned temp = internalNode0;
416  internalNode0 = internalNode2;
417  internalNode2 = internalNode1;
418  internalNode1 = temp;
419  }
420 
421  if (reverse)
422  {
423  unsigned temp = internalNode1;
424  internalNode1 = internalNode2;
425  internalNode2 = temp;
426  }
427 
428  AddNodeToBoundaryElement(pMesh, pBoundaryElement, pElement, internalNode0);
429  AddNodeToBoundaryElement(pMesh, pBoundaryElement, pElement, internalNode1);
430  AddNodeToBoundaryElement(pMesh, pBoundaryElement, pElement, internalNode2);
431 }
432 #undef COVERAGE_IGNORE
433 
434 // Explicit instantiation
437 
438 template class QuadraticMeshHelper<1>;
439 template class QuadraticMeshHelper<2>;
440 template class QuadraticMeshHelper<3>;
static void AddNodesToBoundaryElements(AbstractTetrahedralMesh< DIM, DIM > *pMesh, AbstractMeshReader< DIM, DIM > *pMeshReader)
ElementIterator GetElementIteratorBegin()
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
std::vector< Node< SPACE_DIM > * > mBoundaryNodes
virtual unsigned GetNumLocalBoundaryElements() const
static void AddInternalNodesToElements(AbstractTetrahedralMesh< DIM, DIM > *pMesh, AbstractMeshReader< DIM, DIM > *pMeshReader)
Definition: Node.hpp:58
void SetAsBoundaryNode(bool value=true)
Definition: Node.cpp:127
void MarkAsInternal()
Definition: Node.cpp:359
BoundaryElementIterator GetBoundaryElementIteratorBegin() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
virtual Node< SPACE_DIM > * GetNodeOrHaloNode(unsigned index) const
unsigned ContainingElement
const std::vector< unsigned > & rGetNodePermutation() const
static void AddInternalNodesToBoundaryElements(AbstractTetrahedralMesh< DIM, DIM > *pMesh, AbstractMeshReader< DIM, DIM > *pMeshReader)
bool IsBoundaryNode() const
Definition: Node.cpp:165
#define EXCEPTION(message)
Definition: Exception.hpp:143
static void AddExtraBoundaryNodes(AbstractTetrahedralMesh< DIM, DIM > *pMesh, BoundaryElement< DIM-1, DIM > *pBoundaryElement, Element< DIM, DIM > *pElement, unsigned nodeIndexOppositeToFace)
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
static void HelperMethod2(AbstractTetrahedralMesh< DIM, DIM > *pMesh, BoundaryElement< DIM-1, DIM > *pBoundaryElement, Element< DIM, DIM > *pElement, unsigned internalNode0, unsigned internalNode1, unsigned internalNode2, unsigned offset, bool reverse)
static void HelperMethod1(unsigned boundaryElemNode0, unsigned boundaryElemNode1, Element< DIM, DIM > *pElement, unsigned node0, unsigned node1, unsigned node2, unsigned &rOffset, bool &rReverse)
virtual unsigned GetNumLocalElements() const
virtual ElementData GetNextFaceData()=0
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
ElementIterator GetElementIteratorEnd()
ContainingElementIterator ContainingElementsBegin() const
Definition: Node.hpp:444
virtual void Reset()=0
virtual bool GetReadContainingElementOfBoundaryElement()
void AddElement(unsigned index)
Definition: Node.cpp:269
virtual unsigned GetOrderOfBoundaryElements()
void AddNode(Node< SPACE_DIM > *pNode)
unsigned GetIndex() const
Definition: Node.cpp:159
static void CheckBoundaryElements(AbstractTetrahedralMesh< DIM, DIM > *pMesh)
ContainingElementIterator ContainingElementsEnd() const
Definition: Node.hpp:452
static void AddNodeToBoundaryElement(AbstractTetrahedralMesh< DIM, DIM > *pMesh, BoundaryElement< DIM-1, DIM > *pBoundaryElement, Element< DIM, DIM > *pElement, unsigned internalNode)
BoundaryElementIterator GetBoundaryElementIteratorEnd() const