Chaste  Release::2018.1
AbstractTetrahedralMesh.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 <limits>
37 #include "AbstractTetrahedralMesh.hpp"
38 
40 // Implementation
42 
43 
44 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
46 {
47  unsigned lo=this->GetDistributedVectorFactory()->GetLow();
48  unsigned hi=this->GetDistributedVectorFactory()->GetHigh();
49  for (unsigned element_index=0; element_index<mElements.size(); element_index++)
50  {
51  Element<ELEMENT_DIM, SPACE_DIM>* p_element = mElements[element_index];
52  p_element->SetOwnership(false);
53  for (unsigned local_node_index=0; local_node_index< p_element->GetNumNodes(); local_node_index++)
54  {
55  unsigned global_node_index = p_element->GetNodeGlobalIndex(local_node_index);
56  if (lo<=global_node_index && global_node_index<hi)
57  {
58  p_element->SetOwnership(true);
59  break;
60  }
61  }
62  }
63 }
64 
65 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
67  : mMeshIsLinear(true)
68 {
69 }
70 
71 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
73 {
74  // Iterate over elements and free the memory
75  for (unsigned i=0; i<mElements.size(); i++)
76  {
77  delete mElements[i];
78  }
79  // Iterate over boundary elements and free the memory
80  for (unsigned i=0; i<mBoundaryElements.size(); i++)
81  {
82  delete mBoundaryElements[i];
83  }
84 }
85 
86 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
88 {
89  return mElements.size();
90 }
91 
92 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
94 {
95  return GetNumElements();
96 }
97 
98 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
100 {
101  return mElements.size();
102 }
103 
104 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
106 {
107  return mBoundaryElements.size();
108 }
109 
110 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
112 {
113  return GetNumBoundaryElements();
114 }
115 
116 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
118 {
119  return mBoundaryElements.size();
120 }
121 
122 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
124 {
125  return 0u;
126 }
127 
128 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
130 {
131  return this->GetNumNodes();
132 }
133 
134 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
136 {
137  return this->GetNumAllNodes();
138 }
139 
140 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
142 {
143  unsigned local_index = SolveElementMapping(index);
144  return mElements[local_index];
145 }
146 
147 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
149 {
150  unsigned local_index = SolveBoundaryElementMapping(index);
151  return mBoundaryElements[local_index];
152 }
153 
154 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
156 {
157  return mBoundaryElements.begin();
158 }
159 
160 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
162 {
163  return mBoundaryElements.end();
164 }
165 
166 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
168  unsigned elementIndex,
169  c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian,
170  double& rJacobianDeterminant,
171  c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian) const
172 {
173  mElements[SolveElementMapping(elementIndex)]->CalculateInverseJacobian(rJacobian, rJacobianDeterminant, rInverseJacobian);
174 }
175 
176 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
178  unsigned elementIndex,
179  c_vector<double, SPACE_DIM>& rWeightedDirection,
180  double& rJacobianDeterminant) const
181 {
182  mBoundaryElements[SolveBoundaryElementMapping(elementIndex)]->CalculateWeightedDirection(rWeightedDirection, rJacobianDeterminant );
183 }
184 
185 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
187 {
188  if (ELEMENT_DIM <= 1)
189  {
190  //If the ELEMENT_DIM of the mesh is 1 then the boundary will have ELEMENT_DIM = 0
191  EXCEPTION("1-D mesh has no boundary normals");
192  }
193  for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryElementIterator face_iter = this->GetBoundaryElementIteratorBegin();
194  face_iter != this->GetBoundaryElementIteratorEnd();
195  ++face_iter)
196  {
197  //Form a set for the boundary element node indices
198  std::set<unsigned> boundary_element_node_indices;
199  for (unsigned i=0; i<ELEMENT_DIM; i++)
200  {
201  boundary_element_node_indices.insert( (*face_iter)->GetNodeGlobalIndex(i) );
202  }
203 
204  Node<SPACE_DIM>* p_opposite_node = nullptr;
205  Node<SPACE_DIM>* p_representative_node = (*face_iter)->GetNode(0);
206  for (typename Node<SPACE_DIM>::ContainingElementIterator element_iter = p_representative_node->ContainingElementsBegin();
207  element_iter != p_representative_node->ContainingElementsEnd();
208  ++element_iter)
209  {
210  Element<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*element_iter);
211  //Form a set for the element node indices
212  std::set<unsigned> element_node_indices;
213  for (unsigned i=0; i<=ELEMENT_DIM; i++)
214  {
215  element_node_indices.insert( p_element->GetNodeGlobalIndex(i) );
216  }
217 
218  std::vector<unsigned> difference(ELEMENT_DIM);
219 
220  std::vector<unsigned>::iterator set_iter = std::set_difference(
221  element_node_indices.begin(),element_node_indices.end(),
222  boundary_element_node_indices.begin(), boundary_element_node_indices.end(),
223  difference.begin());
224  if (set_iter - difference.begin() == 1)
225  {
226  p_opposite_node = this -> GetNodeOrHaloNode(difference[0]);
227  break;
228  }
229  }
230  assert(p_opposite_node != nullptr);
231 
232  // Vector from centroid of face to opposite node
233  c_vector<double, SPACE_DIM> into_mesh = p_opposite_node->rGetLocation() - (*face_iter)->CalculateCentroid();
234  c_vector<double, SPACE_DIM> normal = (*face_iter)->CalculateNormal();
235 
236  if (inner_prod(into_mesh, normal) > 0.0)
237  {
238  EXCEPTION("Inward facing normal in boundary element index "<<(*face_iter)->GetIndex());
239  }
240  }
241 }
242 
243 
244 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
246 {
247  assert(ELEMENT_DIM == 1); // LCOV_EXCL_LINE
248 
249  for (unsigned node_index=0; node_index<=width; node_index++)
250  {
251  Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, node_index==0 || node_index==width, node_index);
252  this->mNodes.push_back(p_node); // create node
253  if (node_index==0) // create left boundary node and boundary element
254  {
255  this->mBoundaryNodes.push_back(p_node);
256  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(0, p_node) );
257  }
258  if (node_index==width) // create right boundary node and boundary element
259  {
260  this->mBoundaryNodes.push_back(p_node);
261  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(1, p_node) );
262  }
263  if (node_index > 0) // create element
264  {
265  std::vector<Node<SPACE_DIM>*> nodes;
266  nodes.push_back(this->mNodes[node_index-1]);
267  nodes.push_back(this->mNodes[node_index]);
268  this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(node_index-1, nodes) );
269  }
270  }
271 
272  this->RefreshMesh();
273 }
274 
275 
276 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
277 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructRectangularMesh(unsigned width, unsigned height, bool stagger)
278 {
279  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
280  assert(ELEMENT_DIM == 2); // LCOV_EXCL_LINE
281 
282  //Construct the nodes
283  unsigned node_index=0;
284  for (unsigned j=0; j<height+1; j++)
285  {
286  for (unsigned i=0; i<width+1; i++)
287  {
288  bool is_boundary=false;
289  if (i==0 || j==0 || i==width || j==height)
290  {
291  is_boundary=true;
292  }
293  //Check in place for parallel
294  assert(node_index==(width+1)*(j) + i);
295  Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index++, is_boundary, i, j);
296  this->mNodes.push_back(p_node);
297  if (is_boundary)
298  {
299  this->mBoundaryNodes.push_back(p_node);
300  }
301  }
302  }
303 
304  //Construct the boundary elements
305  unsigned belem_index=0;
306  //Top
307  for (unsigned i=0; i<width; i++)
308  {
309  std::vector<Node<SPACE_DIM>*> nodes;
310  nodes.push_back(this->mNodes[height*(width+1)+i+1]);
311  nodes.push_back(this->mNodes[height*(width+1)+i]);
312  assert(belem_index==i);
313  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
314  }
315  //Right
316  for (unsigned j=1; j<=height; j++)
317  {
318  std::vector<Node<SPACE_DIM>*> nodes;
319  nodes.push_back(this->mNodes[(width+1)*j-1]);
320  nodes.push_back(this->mNodes[(width+1)*(j+1)-1]);
321  assert(belem_index==width+j-1);
322  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
323  }
324  //Bottom
325  for (unsigned i=0; i<width; i++)
326  {
327  std::vector<Node<SPACE_DIM>*> nodes;
328  nodes.push_back(this->mNodes[i]);
329  nodes.push_back(this->mNodes[i+1]);
330  assert(belem_index==width+height+i);
331  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
332  }
333  //Left
334  for (unsigned j=0; j<height; j++)
335  {
336  std::vector<Node<SPACE_DIM>*> nodes;
337  nodes.push_back(this->mNodes[(width+1)*(j+1)]);
338  nodes.push_back(this->mNodes[(width+1)*(j)]);
339  assert(belem_index==2*width+height+j);
340  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,nodes));
341  }
342 
343  //Construct the elements
344  unsigned elem_index = 0;
345  for (unsigned j=0; j<height; j++)
346  {
347  for (unsigned i=0; i<width; i++)
348  {
349  unsigned parity=(i+(height-j))%2;//Note that parity is measured from the top-left (not bottom left) for historical reasons
350  unsigned nw=(j+1)*(width+1)+i; //ne=nw+1
351  unsigned sw=(j)*(width+1)+i; //se=sw+1
352  std::vector<Node<SPACE_DIM>*> upper_nodes;
353  upper_nodes.push_back(this->mNodes[nw]);
354  upper_nodes.push_back(this->mNodes[nw+1]);
355  if (stagger==false || parity == 1)
356  {
357  upper_nodes.push_back(this->mNodes[sw+1]);
358  }
359  else
360  {
361  upper_nodes.push_back(this->mNodes[sw]);
362  }
363  assert(elem_index==2*(j*width+i));
364  this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++,upper_nodes));
365  std::vector<Node<SPACE_DIM>*> lower_nodes;
366  lower_nodes.push_back(this->mNodes[sw+1]);
367  lower_nodes.push_back(this->mNodes[sw]);
368  if (stagger==false ||parity == 1)
369  {
370  lower_nodes.push_back(this->mNodes[nw]);
371  }
372  else
373  {
374  lower_nodes.push_back(this->mNodes[nw+1]);
375  }
376  this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++,lower_nodes));
377  }
378  }
379 
380  this->RefreshMesh();
381 }
382 
383 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
385  unsigned height,
386  unsigned depth)
387 {
388  assert(SPACE_DIM == 3); // LCOV_EXCL_LINE
389  assert(ELEMENT_DIM == 3); // LCOV_EXCL_LINE
390  //Construct the nodes
391 
392  unsigned node_index = 0;
393  for (unsigned k=0; k<depth+1; k++)
394  {
395  for (unsigned j=0; j<height+1; j++)
396  {
397  for (unsigned i=0; i<width+1; i++)
398  {
399  bool is_boundary = false;
400  if (i==0 || j==0 || k==0 || i==width || j==height || k==depth)
401  {
402  is_boundary = true;
403  }
404  assert(node_index == (k*(height+1)+j)*(width+1)+i);
405  Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index++, is_boundary, i, j, k);
406 
407  this->mNodes.push_back(p_node);
408  if (is_boundary)
409  {
410  this->mBoundaryNodes.push_back(p_node);
411  }
412  }
413  }
414  }
415 
416  // Construct the elements
417 
418  unsigned elem_index = 0;
419  unsigned belem_index = 0;
420  unsigned element_nodes[6][4] = {{0, 1, 5, 7}, {0, 1, 3, 7},
421  {0, 2, 3, 7}, {0, 2, 6, 7},
422  {0, 4, 6, 7}, {0, 4, 5, 7}};
423 /* Alternative tessellation - (gerardus)
424  * Note that our method (above) has a bias that all tetrahedra share a
425  * common edge (the diagonal 0 - 7). In the following method the cube is
426  * split along the "face diagonal" 1-2-5-6 into two prisms. This also has a bias.
427  *
428  unsigned element_nodes[6][4] = {{ 0, 6, 5, 4},
429  { 0, 2, 6, 1},
430  { 0, 1, 6, 5},
431  { 1, 2, 3, 7},
432  { 1, 2, 6, 7},
433  { 1, 6, 7, 5 }};
434 */
435  std::vector<Node<SPACE_DIM>*> tetrahedra_nodes;
436 
437  for (unsigned k=0; k<depth; k++)
438  {
439  if (k!=0)
440  {
441  // height*width squares on upper face, k layers of 2*height+2*width square aroun
442  assert(belem_index == 2*(height*width+k*2*(height+width)) );
443  }
444  for (unsigned j=0; j<height; j++)
445  {
446  for (unsigned i=0; i<width; i++)
447  {
448  // Compute the nodes' index
449  unsigned global_node_indices[8];
450  unsigned local_node_index = 0;
451 
452  for (unsigned z = 0; z < 2; z++)
453  {
454  for (unsigned y = 0; y < 2; y++)
455  {
456  for (unsigned x = 0; x < 2; x++)
457  {
458  global_node_indices[local_node_index] = i+x+(width+1)*(j+y+(height+1)*(k+z));
459 
460  local_node_index++;
461  }
462  }
463  }
464 
465  for (unsigned m = 0; m < 6; m++)
466  {
467  // Tetrahedra #m
468 
469  tetrahedra_nodes.clear();
470 
471  for (unsigned n = 0; n < 4; n++)
472  {
473  tetrahedra_nodes.push_back(this->mNodes[global_node_indices[element_nodes[m][n]]]);
474  }
475 
476  assert(elem_index == 6 * ((k*height+j)*width+i)+m );
477  this->mElements.push_back(new Element<ELEMENT_DIM,SPACE_DIM>(elem_index++, tetrahedra_nodes));
478  }
479 
480  // Are we at a boundary?
481  std::vector<Node<SPACE_DIM>*> triangle_nodes;
482 
483  if (i == 0) //low face at x==0
484  {
485  triangle_nodes.clear();
486  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
487  triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
488  triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
489  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
490  triangle_nodes.clear();
491  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
492  triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
493  triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
494  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
495  }
496  if (i == width-1) //high face at x=width
497  {
498  triangle_nodes.clear();
499  triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
500  triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
501  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
502  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
503  triangle_nodes.clear();
504  triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
505  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
506  triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
507  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
508  }
509  if (j == 0) //low face at y==0
510  {
511  triangle_nodes.clear();
512  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
513  triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
514  triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
515  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
516  triangle_nodes.clear();
517  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
518  triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
519  triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
520  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
521  }
522  if (j == height-1) //high face at y=height
523  {
524  triangle_nodes.clear();
525  triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
526  triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
527  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
528  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
529  triangle_nodes.clear();
530  triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
531  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
532  triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
533  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
534  }
535  if (k == 0) //low face at z==0
536  {
537  triangle_nodes.clear();
538  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
539  triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
540  triangle_nodes.push_back(this->mNodes[global_node_indices[2]]);
541  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
542  triangle_nodes.clear();
543  triangle_nodes.push_back(this->mNodes[global_node_indices[0]]);
544  triangle_nodes.push_back(this->mNodes[global_node_indices[1]]);
545  triangle_nodes.push_back(this->mNodes[global_node_indices[3]]);
546  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
547  }
548  if (k == depth-1) //high face at z=depth
549  {
550  triangle_nodes.clear();
551  triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
552  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
553  triangle_nodes.push_back(this->mNodes[global_node_indices[5]]);
554  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
555  triangle_nodes.clear();
556  triangle_nodes.push_back(this->mNodes[global_node_indices[4]]);
557  triangle_nodes.push_back(this->mNodes[global_node_indices[6]]);
558  triangle_nodes.push_back(this->mNodes[global_node_indices[7]]);
559  this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>(belem_index++,triangle_nodes));
560  }
561  }//i
562  }//j
563  }//k
564 
565  this->RefreshMesh();
566 }
567 
568 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
569 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ConstructRegularSlabMesh(double spaceStep, double width, double height, double depth)
570 {
571  assert(spaceStep>0.0);
572  assert(width>0.0);
573  if (ELEMENT_DIM > 1)
574  {
575  assert(height>0.0);
576  }
577  if (ELEMENT_DIM > 2)
578  {
579  assert(depth>0.0);
580  }
581  unsigned num_elem_x=(unsigned)((width+0.5*spaceStep)/spaceStep); //0.5*spaceStep is to ensure that rounding down snaps to correct number
582  unsigned num_elem_y=(unsigned)((height+0.5*spaceStep)/spaceStep);
583  unsigned num_elem_z=(unsigned)((depth+0.5*spaceStep)/spaceStep);
584 
585  //Make it obvious that actual_width_x etc. are temporaries used in spotting for exception
586  {
587  double actual_width_x=num_elem_x*spaceStep;
588  double actual_width_y=num_elem_y*spaceStep;
589  double actual_width_z=num_elem_z*spaceStep;
590 
591  //Note here that in ELEMENT_DIM > 1 cases there may be a zero height or depth - in which case we don't need to use relative comparisons
592  // Doing relative comparisons with zero is okay - if we avoid division by zero.
593  // However, it's best not to test whether " fabs( 0.0 - 0.0) > DBL_EPSILON*0.0 "
594  if (fabs (actual_width_x - width) > DBL_EPSILON*width
595  ||( height!= 0.0 && fabs (actual_width_y - height) > DBL_EPSILON*height)
596  ||( depth != 0.0 && fabs (actual_width_z - depth) > DBL_EPSILON*depth ))
597  {
598  EXCEPTION("Space step does not divide the size of the mesh");
599  }
600  }
601  switch (ELEMENT_DIM)
602  {
603  case 1:
604  this->ConstructLinearMesh(num_elem_x);
605  break;
606  case 2:
607  this->ConstructRectangularMesh(num_elem_x, num_elem_y); // Stagger=default value
608  break;
609  default:
610  case 3:
611  this->ConstructCuboid(num_elem_x, num_elem_y, num_elem_z);
612  }
613  this->Scale(spaceStep, spaceStep, spaceStep);
614 }
615 
616 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
618  unsigned dimension, double spaceStep,
619  double width, double height, double depth)
620 {
621  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
622  if (dimension >= SPACE_DIM)
623  {
624  EXCEPTION("Cannot split on non-existent dimension");
625  }
626  // Rotate the width -> height -> depth around (if dimension is non-default)
627  if (SPACE_DIM == 2 && dimension == 0)
628  {
629  double temp = height ;
630  height = width;
631  width = temp;
632  }
633  else if (SPACE_DIM == 3)
634  {
635  unsigned rotate_perm = SPACE_DIM - 1u - dimension; // How many shuffles to get the user's axis to the top
636  for (unsigned i=0; i<rotate_perm; i++)
637  {
638  double temp = depth;
639  depth = height;
640  height = width;
641  width = temp;
642  }
643  }
644  this->ConstructRegularSlabMesh(spaceStep, width, height, depth);
645  if (SPACE_DIM == 2 && dimension == 0)
646  {
647  // Rotate the positions back again x -> y -> x
648  // this->Rotate(M_PI_2);
649  c_matrix<double, 2, 2> axis_rotation = zero_matrix<double>(2, 2);
650  axis_rotation(0,1)=1.0;
651  axis_rotation(1,0)=-1.0;
652  this->Rotate(axis_rotation);
653  this->Translate(0.0, width); // Formerly known as height, but we rotated it
654  }
655  else if (SPACE_DIM == 3 && dimension == 0)
656  {
657  // this->RotateZ(M_PI_2);
658  // this->RotateY(M_PI_2);
659  // RotY * RotZ = [0 0 1; 1 0 0; 0 1 0] x->y->z->x
660  //this->Translate(depth /*old width*/, width /*old height*/, 0.0);
661  c_matrix<double, 3, 3> axis_permutation = zero_matrix<double>(3, 3);
662  axis_permutation(0, 2)=1.0;
663  axis_permutation(1, 0)=1.0;
664  axis_permutation(2, 1)=1.0;
665  this->Rotate(axis_permutation);
666  }
667  else if (SPACE_DIM == 3 && dimension == 1)
668  {
669  // this->RotateY(-M_PI_2);
670  // this->RotateZ(-M_PI_2);
671  // // RotZ' after RotY' = RotZ' * RotY' = [0 1 0; 0 0 1; 1 0 0] x->z->y->x
672  // this->Translate(height /*old width*/, 0.0, width /*old depth*/);
673  c_matrix<double, 3, 3> axis_permutation = zero_matrix<double>(3, 3);
674  axis_permutation(0, 1)=1.0;
675  axis_permutation(1, 2)=1.0;
676  axis_permutation(2, 0)=1.0;
677  this->Rotate(axis_permutation);
678  }
679 
680 }
681 
682 
683 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
685 {
686  // This may throw in the distributed parallel case
687  unsigned tie_break_index = this->GetBoundaryElement(faceIndex)->GetNodeGlobalIndex(0);
688 
689  // if it is in my range
690  if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(tie_break_index))
691  {
692  return true;
693  }
694  else
695  {
696  return false;
697  }
698 }
699 
700 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
702 {
703  // This may throw in the distributed parallel case
704  unsigned tie_break_index = this->GetElement(elementIndex)->GetNodeGlobalIndex(0);
705 
706  // if it is in my range
707  if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(tie_break_index))
708  {
709  return true;
710  }
711  else
712  {
713  return false;
714  }
715 }
716 
717 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
719 {
720 
721  if (this->mNodes.size() == 0u)
722  {
723 // LCOV_EXCL_START
724  /*
725  * Coverage of this block requires a mesh regular slab mesh with the number of
726  * elements in the primary dimension less than (num_procs - 1), e.g. a 1D mesh
727  * one element wide with num_procs >=3.
728  * Note that if a process owns no nodes, then it will still need to enter the
729  * collective call to MatMPIAIJSetPreallocation. In PetscTools::SetupMat, the
730  * rowPreallocation parameter uses the special value zero to indicate no preallocation.
731  */
732 
733  // This process owns no nodes and thus owns none of the mesh
734  return (1u);
735 // LCOV_EXCL_STOP
736  }
737 
738  unsigned nodes_per_element = this->mElements[0]->GetNumNodes(); //Usually ELEMENT_DIM+1, except in Quadratic case
739  if (ELEMENT_DIM <= 2u)
740  {
741  /*
742  * Note that we start assuming that each internal node is connected to 1 element.
743  * This is to avoid the trivial situation in which there are no internal nodes (a
744  * single line or a single triangle. We need the minimum connectivity to be 2 (in 1D) or
745  * 3 (in 2D) even if there is only one element.
746  */
747  unsigned max_num = 1u; // See note above.
748  unsigned boundary_max_num = 0u;
749  for (unsigned local_node_index=0; local_node_index<this->mNodes.size(); local_node_index++)
750  {
751  unsigned num = this->mNodes[local_node_index]->GetNumContainingElements();
752  if (this->mNodes[local_node_index]->IsBoundaryNode()==false && num>max_num)
753  {
754  max_num = num;
755  }
756  if (this->mNodes[local_node_index]->IsBoundaryNode() && num>boundary_max_num)
757  {
758  boundary_max_num = num;
759  }
760  }
761  bool linear = (nodes_per_element == ELEMENT_DIM + 1);
762  /*
763  * In 1d each containing element is connected to one node (or 2 if quadratic), add to this the node itself
764  * and the connectivity is GetNumContainingElements() + 1 or 2*GetNumContainingElements() + 1
765  */
766  if (ELEMENT_DIM == 1)
767  {
768  if (linear)
769  {
770  return max_num+1;
771  }
772  else
773  {
774  return 2*max_num+1;
775  }
776  }
777  // Not returned ...else if (ELEMENT_DIM == 2)
778  /*
779  * In 2d each containing element adds one connected node (since one node will be shared by a previous element)
780  * this leads to a connectivity of GetNumContainingElements() + 1 (or 3*GetNumContainingElements() + 1) in the quadratic case
781  *
782  * If the node is on a boundary then the one elements will have an unpaired node and the connectivity is
783  * GetNumContainingElements() + 1 (or 3*(GetNumContainingElements() + 3 for quadratic)
784  */
785  if (linear)
786  {
787  return std::max(max_num+1, boundary_max_num+2);
788  }
789  else
790  {
791  return std::max(3*max_num+1, 3*boundary_max_num+3);
792  }
793  }
794 
795  /*
796  * In 3d there are many more cases. In general a non-boundary node has fewer connecting nodes than it has elements.
797  * A node on the boundary may have even fewer, unless it is on a corner and has more faces than it has elements.
798  * We can, in the linear case estimate an upper bound as max(elements, faces)+2.
799  * However for the sake of accuracy we are going for a brute force solution.
800  *
801  * This may prove to be a bottle-neck...
802  */
803 
804  std::set<unsigned> forward_star_nodes; // Used to collect each node's neighbours
805  unsigned max_connectivity = 0u;
806  for (unsigned local_node_index=0; local_node_index<this->mNodes.size(); local_node_index++)
807  {
808  forward_star_nodes.clear();
809 
810  for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[local_node_index]->ContainingElementsBegin();
811  it != this->mNodes[local_node_index]->ContainingElementsEnd();
812  ++it)
813  {
814  Element<ELEMENT_DIM, SPACE_DIM>* p_elem = this->GetElement(*it);
815  for (unsigned i=0; i<nodes_per_element; i++)
816  {
817  forward_star_nodes.insert(p_elem->GetNodeGlobalIndex(i));
818  }
819  }
820  if (forward_star_nodes.size() > max_connectivity)
821  {
822  max_connectivity = forward_star_nodes.size();
823  }
824  }
825  return max_connectivity;
826 }
827 
828 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
829 void AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetHaloNodeIndices(std::vector<unsigned>& rHaloIndices) const
830 {
831  // Make sure the output vector is empty
832  rHaloIndices.clear();
833 }
834 
835 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
837 {
838  for (unsigned i=0; i<rOtherMesh.GetNumNodes(); i++)
839  {
840  Node<SPACE_DIM>* p_node = rOtherMesh.GetNode(i);
841  assert(!p_node->IsDeleted());
842  const c_vector<double, SPACE_DIM>& location=p_node->rGetLocation();
843  bool is_boundary=p_node->IsBoundaryNode();
844 
845  Node<SPACE_DIM>* p_node_copy = new Node<SPACE_DIM>(i, location, is_boundary);
846  this->mNodes.push_back( p_node_copy );
847  if (is_boundary)
848  {
849  this->mBoundaryNodes.push_back( p_node_copy );
850  }
851  }
852 
853  for (unsigned i=0; i<rOtherMesh.GetNumElements(); i++)
854  {
855  Element<ELEMENT_DIM, SPACE_DIM>* p_elem = rOtherMesh.GetElement(i);
856  assert(!p_elem->IsDeleted());
857  std::vector<Node<SPACE_DIM>*> nodes_for_element;
858  for (unsigned j=0; j<p_elem->GetNumNodes(); j++)
859  {
860  nodes_for_element.push_back(this->mNodes[ p_elem->GetNodeGlobalIndex(j) ]);
861  }
862  Element<ELEMENT_DIM, SPACE_DIM>* p_elem_copy = new Element<ELEMENT_DIM, SPACE_DIM>(i, nodes_for_element);
863  p_elem_copy->RegisterWithNodes();
864  this->mElements.push_back(p_elem_copy);
865  }
866 
867  for (unsigned i=0; i<rOtherMesh.GetNumBoundaryElements(); i++)
868  {
869  BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_b_elem = rOtherMesh.GetBoundaryElement(i);
870  assert(!p_b_elem->IsDeleted());
871  std::vector<Node<SPACE_DIM>*> nodes_for_element;
872  for (unsigned j=0; j<p_b_elem->GetNumNodes(); j++)
873  {
874  nodes_for_element.push_back(this->mNodes[ p_b_elem->GetNodeGlobalIndex(j) ]);
875  }
876  BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_b_elem_copy = new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(i, nodes_for_element);
877  p_b_elem_copy->RegisterWithNodes();
878  this->mBoundaryElements.push_back(p_b_elem_copy);
879  }
880  this->RefreshMesh();
881 
882 }
883 
884 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
886  std::vector<std::vector<unsigned> >& rNodesToSendPerProcess,
887  std::vector<std::vector<unsigned> >& rNodesToReceivePerProcess)
888 {
889  assert( rNodesToSendPerProcess.empty() );
890  assert( rNodesToReceivePerProcess.empty() );
891 
892  //Initialise vectors of sets for the exchange data
893  std::vector<std::set<unsigned> > node_sets_to_send_per_process;
894  std::vector<std::set<unsigned> > node_sets_to_receive_per_process;
895 
896  node_sets_to_send_per_process.resize(PetscTools::GetNumProcs());
897  node_sets_to_receive_per_process.resize(PetscTools::GetNumProcs());
898  std::vector<unsigned> global_lows = this->GetDistributedVectorFactory()->rGetGlobalLows();
899 
900  for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = this->GetElementIteratorBegin();
901  iter != this->GetElementIteratorEnd();
902  ++iter)
903  {
904  std::vector <unsigned> nodes_on_this_process;
905  std::vector <unsigned> nodes_not_on_this_process;
906  //Calculate local and non-local node indices
907  for (unsigned i=0; i<ELEMENT_DIM+1; i++)
908  {
909  unsigned node_index=iter->GetNodeGlobalIndex(i);
910  if (this->GetDistributedVectorFactory()->IsGlobalIndexLocal(node_index))
911  {
912  nodes_on_this_process.push_back(node_index);
913  }
914  else
915  {
916  nodes_not_on_this_process.push_back(node_index);
917  }
918  }
919  /*
920  * If this is a TetrahedralMesh (not distributed) then it's possible that we own none
921  * of the nodes in this element. In that case we must skip the element.
922  */
923  if (nodes_on_this_process.empty())
924  {
925  continue; //Move on to the next element.
927  }
928  // If there are any non-local nodes on this element then we need to add to the data exchange
929  if (!nodes_not_on_this_process.empty())
930  {
931  for (unsigned i=0; i<nodes_not_on_this_process.size(); i++)
932  {
933  // Calculate who owns this remote node
934  unsigned remote_process=global_lows.size()-1;
935  for (; global_lows[remote_process] > nodes_not_on_this_process[i]; remote_process--)
936  {
937  }
938 
939  // Add this node to the correct receive set
940  node_sets_to_receive_per_process[remote_process].insert(nodes_not_on_this_process[i]);
941 
942  // Add all local nodes to the send set
943  for (unsigned j=0; j<nodes_on_this_process.size(); j++)
944  {
945  node_sets_to_send_per_process[remote_process].insert(nodes_on_this_process[j]);
946  }
947  }
948  }
949  }
950 
951  for (unsigned process_number = 0; process_number < PetscTools::GetNumProcs(); process_number++)
952  {
953  std::vector<unsigned> process_send_vector( node_sets_to_send_per_process[process_number].begin(),
954  node_sets_to_send_per_process[process_number].end() );
955  std::sort(process_send_vector.begin(), process_send_vector.end());
956 
957  rNodesToSendPerProcess.push_back(process_send_vector);
958 
959  std::vector<unsigned> process_receive_vector( node_sets_to_receive_per_process[process_number].begin(),
960  node_sets_to_receive_per_process[process_number].end() );
961  std::sort(process_receive_vector.begin(), process_receive_vector.end());
962 
963  rNodesToReceivePerProcess.push_back(process_receive_vector);
964  }
965 
966 }
967 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
969 {
970  c_vector<double, 2> min_max;
971  min_max[0] = DBL_MAX;
972  min_max[1] = 0.0;
973  for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator ele_iter = GetElementIteratorBegin();
974  ele_iter != GetElementIteratorEnd();
975  ++ele_iter)
976  {
977  c_vector<double, 2> ele_min_max = ele_iter->CalculateMinMaxEdgeLengths();
978  if (ele_min_max[0] < min_max[0])
979  {
980  min_max[0] = ele_min_max[0];
981  }
982  if (ele_min_max[1] > min_max[1])
983  {
984  min_max[1] = ele_min_max[1];
985  }
986  }
987  return min_max;
988 }
989 
990 
991 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
993  bool strict,
994  std::set<unsigned> testElements,
995  bool onlyTryWithTestElements)
996 {
997  for (std::set<unsigned>::iterator iter=testElements.begin(); iter!=testElements.end(); iter++)
998  {
999  assert(*iter<this->GetNumElements());
1000  if (this->mElements[*iter]->IncludesPoint(rTestPoint, strict))
1001  {
1002  assert(!this->mElements[*iter]->IsDeleted());
1003  return *iter;
1004  }
1005  }
1006 
1007  if (!onlyTryWithTestElements)
1008  {
1009  for (unsigned i=0; i<this->mElements.size(); i++)
1010  {
1011  if (this->mElements[i]->IncludesPoint(rTestPoint, strict))
1012  {
1013  assert(!this->mElements[i]->IsDeleted());
1014  return i;
1015  }
1016  }
1017  }
1018 
1019  // If it's in none of the elements, then throw
1020  std::stringstream ss;
1021  ss << "Point [";
1022  for (unsigned j=0; (int)j<(int)SPACE_DIM-1; j++)
1023  {
1024  ss << rTestPoint[j] << ",";
1025  }
1026  ss << rTestPoint[SPACE_DIM-1] << "] is not in ";
1027  if (!onlyTryWithTestElements)
1028  {
1029  ss << "mesh - all elements tested";
1030  }
1031  else
1032  {
1033  ss << "set of elements given";
1034  }
1035  EXCEPTION(ss.str());
1036 }
1037 
1038 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1040  std::set<unsigned> testElements)
1041 {
1042  assert(testElements.size() > 0);
1043  EXCEPT_IF_NOT(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE // CalculateInterpolationWeights hits an assertion otherwise
1044 
1045  double max_min_weight = -std::numeric_limits<double>::infinity();
1046  unsigned closest_index = 0;
1047  for (std::set<unsigned>::iterator iter = testElements.begin();
1048  iter != testElements.end();
1049  iter++)
1050  {
1051  c_vector<double, ELEMENT_DIM+1> weight = this->mElements[*iter]->CalculateInterpolationWeights(rTestPoint);
1052  double neg_weight_sum = 0.0;
1053  for (unsigned j=0; j<=ELEMENT_DIM; j++)
1054  {
1055  if (weight[j] < 0.0)
1056  {
1057  neg_weight_sum += weight[j];
1058  }
1059  }
1060  if (neg_weight_sum > max_min_weight)
1061  {
1062  max_min_weight = neg_weight_sum;
1063  closest_index = *iter;
1064  }
1065  }
1066  assert(!this->mElements[closest_index]->IsDeleted());
1067  return closest_index;
1068 }
1069 
1070 // Explicit instantiation
1071 template class AbstractTetrahedralMesh<1,1>;
1072 template class AbstractTetrahedralMesh<1,2>;
1073 template class AbstractTetrahedralMesh<1,3>;
1074 template class AbstractTetrahedralMesh<2,2>;
1075 template class AbstractTetrahedralMesh<2,3>;
1076 template class AbstractTetrahedralMesh<3,3>;
virtual unsigned GetNumBoundaryElements() const
virtual unsigned GetNumLocalBoundaryElements() const
Definition: Node.hpp:58
BoundaryElementIterator GetBoundaryElementIteratorBegin() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
virtual unsigned GetMaximumNodeIndex()
virtual bool CalculateDesignatedOwnershipOfBoundaryElement(unsigned faceIndex)
void ConstructRegularSlabMesh(double spaceStep, double width, double height=0, double depth=0)
bool IsBoundaryNode() const
Definition: Node.cpp:164
Node< SPACE_DIM > * GetNode(unsigned index) const
virtual unsigned GetNumElements() const
#define EXCEPTION(message)
Definition: Exception.hpp:143
virtual unsigned GetNumCableElements() const
unsigned CalculateMaximumNodeConnectivityPerProcess() const
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
void SetOwnership(bool ownership)
virtual unsigned GetNumNodes() const
virtual void GetHaloNodeIndices(std::vector< unsigned > &rHaloIndices) const
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
void RegisterWithNodes()
Definition: Element.cpp:65
bool IsDeleted() const
ContainingElementIterator ContainingElementsBegin() const
Definition: Node.hpp:485
virtual void GetWeightedDirectionForBoundaryElement(unsigned elementIndex, c_vector< double, SPACE_DIM > &rWeightedDirection, double &rJacobianDeterminant) const
unsigned GetNumNodes() const
virtual unsigned GetNumVertices() const
unsigned GetContainingElementIndex(const ChastePoint< SPACE_DIM > &rTestPoint, bool strict=false, std::set< unsigned > testElements=std::set< unsigned >(), bool onlyTryWithTestElements=false)
#define EXCEPT_IF_NOT(test)
Definition: Exception.hpp:158
virtual void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true)
BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * GetBoundaryElement(unsigned index) const
virtual void ConstructCuboid(unsigned width, unsigned height, unsigned depth)
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:139
void CalculateNodeExchange(std::vector< std::vector< unsigned > > &rNodesToSendPerProcess, std::vector< std::vector< unsigned > > &rNodesToReceivePerProcess)
unsigned GetNumAllBoundaryElements() const
unsigned GetNearestElementIndexFromTestElements(const ChastePoint< SPACE_DIM > &rTestPoint, std::set< unsigned > testElements)
virtual bool CalculateDesignatedOwnershipOfElement(unsigned elementIndex)
virtual void ConstructLinearMesh(unsigned width)
static unsigned GetNumProcs()
Definition: PetscTools.cpp:108
std::vector< BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * >::const_iterator BoundaryElementIterator
virtual c_vector< double, 2 > CalculateMinMaxEdgeLengths()
void ConstructFromMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rOtherMesh)
bool IsDeleted() const
Definition: Node.cpp:412
ContainingElementIterator ContainingElementsEnd() const
Definition: Node.hpp:493
void ConstructRegularSlabMeshWithDimensionSplit(unsigned dimension, double spaceStep, double width, double height=0, double depth=0)
BoundaryElementIterator GetBoundaryElementIteratorEnd() const