Chaste Commit::f2ff7ee04e70ac9d06c57344df8d017dbb12b97b
VertexMesh.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "VertexMesh.hpp"
37#include "MutableMesh.hpp"
38#include "RandomNumberGenerator.hpp"
40
41template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
43 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> vertexElements)
44 : mpDelaunayMesh(nullptr)
45{
46
47 // Reset member variables and clear mNodes and mElements
48 Clear();
49
50 // Populate mNodes and mElements
51 for (unsigned node_index = 0; node_index < nodes.size(); node_index++)
52 {
53 Node<SPACE_DIM>* p_temp_node = nodes[node_index];
54 this->mNodes.push_back(p_temp_node);
55 }
56 for (unsigned elem_index = 0; elem_index < vertexElements.size(); elem_index++)
57 {
58 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
59 mElements.push_back(p_temp_vertex_element);
60 }
61
62 // In 3D, populate mFaces
63 if (SPACE_DIM == 3)
64 {
65 // Use a std::set to keep track of which faces have been added to mFaces
66 std::set<unsigned> faces_counted;
67
68 // Loop over mElements
69 for (unsigned elem_index = 0; elem_index < mElements.size(); elem_index++)
70 {
71 // Loop over faces of this element
72 for (unsigned face_index = 0; face_index < mElements[elem_index]->GetNumFaces(); face_index++)
73 {
74 VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* p_face = mElements[elem_index]->GetFace(face_index);
75 unsigned global_index = p_face->GetIndex();
76
77 // If this face is not already contained in mFaces, add it, and update faces_counted
78 if (faces_counted.find(global_index) == faces_counted.end())
79 {
80 mFaces.push_back(p_face);
81 faces_counted.insert(global_index);
82 }
83 }
84 }
85 }
86
87 // Register elements with nodes
88 for (unsigned index = 0; index < mElements.size(); index++)
89 {
91
92 unsigned element_index = p_element->GetIndex();
93 unsigned num_nodes_in_element = p_element->GetNumNodes();
94
95 for (unsigned node_index = 0; node_index < num_nodes_in_element; node_index++)
96 {
97 p_element->GetNode(node_index)->AddElement(element_index);
98 }
99 }
100
102
103 this->mMeshChangesDuringSimulation = false;
104}
105
106template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
109 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> vertexElements)
110 : mpDelaunayMesh(nullptr)
111{
112 // Reset member variables and clear mNodes, mFaces and mElements
113 Clear();
114
115 // Populate mNodes mFaces and mElements
116 for (unsigned node_index = 0; node_index < nodes.size(); node_index++)
118 Node<SPACE_DIM>* p_temp_node = nodes[node_index];
119 this->mNodes.push_back(p_temp_node);
120 }
121
122 for (unsigned face_index = 0; face_index < faces.size(); face_index++)
123 {
124 VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* p_temp_face = faces[face_index];
125 mFaces.push_back(p_temp_face);
127
128 for (unsigned elem_index = 0; elem_index < vertexElements.size(); elem_index++)
129 {
130 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
131 mElements.push_back(p_temp_vertex_element);
132 }
133
134 // Register elements with nodes
135 for (unsigned index = 0; index < mElements.size(); index++)
136 {
137 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = mElements[index];
138 for (unsigned node_index = 0; node_index < p_temp_vertex_element->GetNumNodes(); node_index++)
139 {
140 Node<SPACE_DIM>* p_temp_node = p_temp_vertex_element->GetNode(node_index);
141 p_temp_node->AddElement(p_temp_vertex_element->GetIndex());
142 }
143 }
144
145 this->mMeshChangesDuringSimulation = false;
146}
147
152template <>
153VertexMesh<2, 2>::VertexMesh(TetrahedralMesh<2, 2>& rMesh, bool isPeriodic, bool isBounded)
154 : mpDelaunayMesh(&rMesh)
155{
156 //Note !isPeriodic is not used except through polymorphic calls in rMesh
157
158 // Reset member variables and clear mNodes, mFaces and mElements
159 Clear();
160
161 if (!isBounded)
162 {
163 unsigned num_elements = mpDelaunayMesh->GetNumAllNodes();
164 unsigned num_nodes = mpDelaunayMesh->GetNumAllElements();
165
166 // Allocate memory for mNodes and mElements
167 this->mNodes.reserve(num_nodes);
168
169 // Create as many elements as there are nodes in the mesh
170 mElements.reserve(num_elements);
171 for (unsigned elem_index = 0; elem_index < num_elements; elem_index++)
173 VertexElement<2, 2>* p_element = new VertexElement<2, 2>(elem_index);
174 mElements.push_back(p_element);
175 }
176
177 // Populate mNodes
179
180 // Loop over elements of the Delaunay mesh (which are nodes/vertices of this mesh)
181 for (unsigned i = 0; i < num_nodes; i++)
182 {
183 // Loop over nodes owned by this triangular element in the Delaunay mesh
184 // Add this node/vertex to each of the 3 vertex elements
185 for (unsigned local_index = 0; local_index < 3; local_index++)
186 {
187 unsigned elem_index = mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
188 unsigned num_nodes_in_elem = mElements[elem_index]->GetNumNodes();
189 unsigned end_index = num_nodes_in_elem > 0 ? num_nodes_in_elem - 1 : 0;
190
191 mElements[elem_index]->AddNode(this->mNodes[i], end_index);
192 }
193 }
194 }
195 else // Is Bounded
196 {
197 // First create an extended mesh to include points extended from the boundary
198 std::vector<Node<2> *> nodes;
199 for (typename TetrahedralMesh<2,2>::NodeIterator node_iter = mpDelaunayMesh->GetNodeIteratorBegin();
200 node_iter != mpDelaunayMesh->GetNodeIteratorEnd();
201 ++node_iter)
202 {
203 nodes.push_back(new Node<2>(node_iter->GetIndex(), node_iter->rGetLocation(),node_iter->IsBoundaryNode()));
204 }
205
206 // Add new nodes
207 unsigned new_node_index = mpDelaunayMesh->GetNumNodes();
208 for (TetrahedralMesh<2,2>::ElementIterator elem_iter = mpDelaunayMesh->GetElementIteratorBegin();
209 elem_iter != mpDelaunayMesh->GetElementIteratorEnd();
210 ++elem_iter)
211 {
212 for (unsigned j=0; j<3; j++)
213 {
214 Node<2>* p_node_a = mpDelaunayMesh->GetNode(elem_iter->GetNodeGlobalIndex(j));
215 Node<2>* p_node_b = mpDelaunayMesh->GetNode(elem_iter->GetNodeGlobalIndex((j+1)%3));
216
217 std::set<unsigned> node_a_element_indices = p_node_a->rGetContainingElementIndices();
218 std::set<unsigned> node_b_element_indices = p_node_b->rGetContainingElementIndices();
219
220 std::set<unsigned> shared_elements;
221 std::set_intersection(node_a_element_indices.begin(),
222 node_a_element_indices.end(),
223 node_b_element_indices.begin(),
224 node_b_element_indices.end(),
225 std::inserter(shared_elements, shared_elements.begin()));
226
227
228 /*
229 * Note using boundary nodes to identify the boundary edges won't work with
230 * triangles which have 3 boundary nodes
231 * if ((p_node_a->IsBoundaryNode() && p_node_b->IsBoundaryNode()))
232 */
233
234 if (shared_elements.size() == 1) // It's a boundary edge
235 {
236 c_vector<double,2> edge = p_node_b->rGetLocation() - p_node_a->rGetLocation();
237 double edge_length = norm_2(edge);
238 c_vector<double,2> normal_vector;
239
240 normal_vector[0]= edge[1];
241 normal_vector[1]= -edge[0];
242
243 double dij = norm_2(normal_vector);
244 assert(dij>1e-5); //Sanity check
245 normal_vector /= dij;
246
247 double extra_node_scaling = 1.0; // increase to add more points per external edge (makes rounder cells)
248
249 int num_sections = ceil(edge_length*extra_node_scaling);
250 for (int section=0; section<=num_sections; section++)
251 {
252 double ratio = (double)section/(double)num_sections;
253 c_vector<double,2> new_node_location = normal_vector + ratio*p_node_a->rGetLocation() + (1-ratio)*p_node_b->rGetLocation();
254 nodes.push_back(new Node<2>(new_node_index, new_node_location));
255 new_node_index++;
256 }
257 }
258 }
259 }
260 MutableMesh<2,2> extended_mesh(nodes);
261
262 unsigned num_elements = mpDelaunayMesh->GetNumAllNodes();
263 unsigned num_nodes = extended_mesh.GetNumAllElements();
264
265 // Allocate memory for mNodes and mElements
266 this->mNodes.reserve(num_nodes);
267
268 // Create as many elements as there are nodes in the mesh
269 mElements.reserve(num_elements);
270 for (unsigned elem_index = 0; elem_index < num_elements; elem_index++)
271 {
272 VertexElement<2, 2>* p_element = new VertexElement<2, 2>(elem_index);
273 mElements.push_back(p_element);
274 }
275
276 // Populate mNodes
278
279 // Loop over elements of the Delaunay mesh (which are nodes/vertices of this mesh)
280 for (unsigned i = 0; i < num_nodes; i++)
281 {
282 // Loop over nodes owned by this triangular element in the Delaunay mesh
283 // Add this node/vertex to each of the 3 vertex elements
284 for (unsigned local_index = 0; local_index < 3; local_index++)
285 {
286 unsigned elem_index = extended_mesh.GetElement(i)->GetNodeGlobalIndex(local_index);
288 if (elem_index < num_elements)
289 {
290 unsigned num_nodes_in_elem = mElements[elem_index]->GetNumNodes();
291 unsigned end_index = num_nodes_in_elem > 0 ? num_nodes_in_elem - 1 : 0;
293 mElements[elem_index]->AddNode(this->mNodes[i], end_index);
294 }
295 }
296 }
297 }
299 // Reorder mNodes anticlockwise
300 for (unsigned elem_index = 0; elem_index < mElements.size(); elem_index++)
301 {
307 std::vector<std::pair<double, unsigned> > index_angle_list;
308 for (unsigned local_index = 0; local_index < mElements[elem_index]->GetNumNodes(); local_index++)
309 {
310 c_vector<double, 2> vectorA = mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
311 c_vector<double, 2> vectorB = mElements[elem_index]->GetNodeLocation(local_index);
312 c_vector<double, 2> centre_to_vertex = mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
313
314 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
315 unsigned global_index = mElements[elem_index]->GetNodeGlobalIndex(local_index);
316
317 std::pair<double, unsigned> pair(angle, global_index);
318 index_angle_list.push_back(pair);
319 }
320
321 // Sort the list in order of increasing angle
322 sort(index_angle_list.begin(), index_angle_list.end());
323
324 // Create a new Voronoi element and pass in the appropriate Nodes, ordered anticlockwise
325 VertexElement<2, 2>* p_new_element = new VertexElement<2, 2>(elem_index);
326 for (unsigned count = 0; count < index_angle_list.size(); count++)
328 unsigned local_index = count > 1 ? count - 1 : 0;
329 p_new_element->AddNode(mNodes[index_angle_list[count].second], local_index);
330 }
331
332 // Replace the relevant member of mElements with this Voronoi element
333 delete mElements[elem_index];
334 mElements[elem_index] = p_new_element;
335 }
336
337 this->mMeshChangesDuringSimulation = false;
338}
348template <>
350 : mpDelaunayMesh(&rMesh)
351{
352 // Reset member variables and clear mNodes, mFaces and mElements
353 Clear();
354
355 unsigned num_nodes = mpDelaunayMesh->GetNumAllElements();
356
357 // Allocate memory for mNodes
358 this->mNodes.reserve(num_nodes);
359
360 // Populate mNodes
362
363 std::map<unsigned, VertexElement<3, 3>*> index_element_map;
364 unsigned face_index = 0;
365 unsigned element_index = 0;
366
367 // Loop over each edge of the Delaunay mesh and populate mFaces and mElements
368 for (TetrahedralMesh<3, 3>::EdgeIterator edge_iterator = mpDelaunayMesh->EdgesBegin();
369 edge_iterator != mpDelaunayMesh->EdgesEnd();
370 ++edge_iterator)
371 {
372 Node<3>* p_node_a = edge_iterator.GetNodeA();
373 Node<3>* p_node_b = edge_iterator.GetNodeB();
374
375 if (!(p_node_a->IsBoundaryNode() && p_node_b->IsBoundaryNode()))
376 {
377 std::set<unsigned>& node_a_element_indices = p_node_a->rGetContainingElementIndices();
378 std::set<unsigned>& node_b_element_indices = p_node_b->rGetContainingElementIndices();
379 std::set<unsigned> edge_element_indices;
380
381 std::set_intersection(node_a_element_indices.begin(),
382 node_a_element_indices.end(),
383 node_b_element_indices.begin(),
384 node_b_element_indices.end(),
385 std::inserter(edge_element_indices, edge_element_indices.begin()));
386
387 c_vector<double, 3> edge_vector;
388 edge_vector = p_node_b->rGetLocation() - p_node_a->rGetLocation();
389
390 c_vector<double, 3> mid_edge;
391 mid_edge = edge_vector * 0.5 + p_node_a->rGetLocation();
392
393 unsigned element0_index = *(edge_element_indices.begin());
394
395 c_vector<double, 3> basis_vector1;
396 basis_vector1 = mNodes[element0_index]->rGetLocation() - mid_edge;
397
398 c_vector<double, 3> basis_vector2;
399 basis_vector2 = VectorProduct(edge_vector, basis_vector1);
400
406 std::vector<std::pair<double, unsigned> > index_angle_list;
407
408 // Loop over each element containing this edge (i.e. those containing both nodes of the edge)
409 for (std::set<unsigned>::iterator index_iter = edge_element_indices.begin();
410 index_iter != edge_element_indices.end();
411 ++index_iter)
412 {
413 // Calculate angle
414 c_vector<double, 3> vertex_vector = mNodes[*index_iter]->rGetLocation() - mid_edge;
416 double local_vertex_dot_basis_vector1 = inner_prod(vertex_vector, basis_vector1);
417 double local_vertex_dot_basis_vector2 = inner_prod(vertex_vector, basis_vector2);
418
419 double angle = atan2(local_vertex_dot_basis_vector2, local_vertex_dot_basis_vector1);
420
421 std::pair<double, unsigned> pair(angle, *index_iter);
422 index_angle_list.push_back(pair);
423 }
424
425 // Sort the list in order of increasing angle
426 sort(index_angle_list.begin(), index_angle_list.end());
428 // Create face
429 VertexElement<2, 3>* p_face = new VertexElement<2, 3>(face_index);
430 face_index++;
431 for (unsigned count = 0; count < index_angle_list.size(); count++)
432 {
433 unsigned local_index = count > 1 ? count - 1 : 0;
434 p_face->AddNode(mNodes[index_angle_list[count].second], local_index);
435 }
436
437 // Add face to list of faces
438 mFaces.push_back(p_face);
439
440 // Add face to appropriate elements
441 if (!p_node_a->IsBoundaryNode())
442 {
443 unsigned node_a_index = p_node_a->GetIndex();
444 if (index_element_map[node_a_index])
445 {
446 // If there is already an element with this index, add the face to it...
447 index_element_map[node_a_index]->AddFace(p_face);
448 }
449 else
450 {
451 // ...otherwise create an element, add the face to it, and add to the map
452 mVoronoiElementIndexMap[node_a_index] = element_index;
453 VertexElement<3, 3>* p_element = new VertexElement<3, 3>(element_index);
454 element_index++;
455 p_element->AddFace(p_face);
456 index_element_map[node_a_index] = p_element;
457 }
458 }
459 if (!p_node_b->IsBoundaryNode())
460 {
461 unsigned node_b_index = p_node_b->GetIndex();
462 if (index_element_map[node_b_index])
463 {
464 // If there is already an element with this index, add the face to it...
465 index_element_map[node_b_index]->AddFace(p_face);
466 }
467 else
468 {
469 // ...otherwise create an element, add the face to it, and add to the map
470 mVoronoiElementIndexMap[node_b_index] = element_index;
471 VertexElement<3, 3>* p_element = new VertexElement<3, 3>(element_index);
472 element_index++;
473 p_element->AddFace(p_face);
474 index_element_map[node_b_index] = p_element;
475 }
476 }
478 }
479
480 // Populate mElements
481 unsigned elem_count = 0;
482 for (std::map<unsigned, VertexElement<3, 3>*>::iterator element_iter = index_element_map.begin();
483 element_iter != index_element_map.end();
484 ++element_iter)
485 {
486 mElements.push_back(element_iter->second);
487 mVoronoiElementIndexMap[element_iter->first] = elem_count;
488 elem_count++;
489 }
490
491 this->mMeshChangesDuringSimulation = false;
492}
498template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
500 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*>& rElements)
501{
502 // Build a list of unique edges from nodes within all the elements
503 for (auto elem : rElements)
504 {
505 elem->SetEdgeHelper(&mEdgeHelper);
506 elem->BuildEdges();
507 }
508}
509
510template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
512{
513 return mEdgeHelper.GetNumEdges();
514}
515
516template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
518{
519 return mEdgeHelper.GetEdge(index);
520}
521
522template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
524{
525 return mEdgeHelper;
526}
528template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
530{
531 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
532 c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
533 double jacobian_det;
534
535 // Loop over elements of the Delaunay mesh and populate mNodes
536 for (unsigned i = 0; i < rMesh.GetNumElements(); i++)
537 {
538 // Calculate the circumcentre of this element in the Delaunay mesh
539 rMesh.GetInverseJacobianForElement(i, jacobian, jacobian_det, inverse_jacobian);
540 c_vector<double, SPACE_DIM + 1> circumsphere = rMesh.GetElement(i)->CalculateCircumsphere(jacobian, inverse_jacobian);
541
542 c_vector<double, SPACE_DIM> circumcentre;
543 for (unsigned j = 0; j < SPACE_DIM; j++)
544 {
545 circumcentre(j) = circumsphere(j);
546 }
547
548 // Create a node in the Voronoi mesh at the location of this circumcentre
549 this->mNodes.push_back(new Node<SPACE_DIM>(i, circumcentre));
550 }
551}
553template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
554double VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)
555{
556 assert(SPACE_DIM == 2); // LCOV_EXCL_LINE - code will be removed at compile time
557
558 std::set<unsigned> node_indices_1;
559 for (unsigned i = 0; i < mElements[elementIndex1]->GetNumNodes(); i++)
560 {
561 node_indices_1.insert(mElements[elementIndex1]->GetNodeGlobalIndex(i));
562 }
563 std::set<unsigned> node_indices_2;
564 for (unsigned i = 0; i < mElements[elementIndex2]->GetNumNodes(); i++)
565 {
566 node_indices_2.insert(mElements[elementIndex2]->GetNodeGlobalIndex(i));
567 }
568
569 std::set<unsigned> shared_nodes;
570 std::set_intersection(node_indices_1.begin(), node_indices_1.end(),
571 node_indices_2.begin(), node_indices_2.end(),
572 std::inserter(shared_nodes, shared_nodes.begin()));
573
574 if (shared_nodes.size() == 1)
575 {
576 // It's possible that these two elements are actually infinite but are on the edge of the domain
577 EXCEPTION("Elements " << elementIndex1 << " and " << elementIndex2 << " share only one node.");
578 }
579 assert(shared_nodes.size() == 2);
580
581 unsigned index1 = *(shared_nodes.begin());
582 unsigned index2 = *(++(shared_nodes.begin()));
583
584 double edge_length = this->GetDistanceBetweenNodes(index1, index2);
585 return edge_length;
586}
587
588template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
591 assert(SPACE_DIM == 2); // LCOV_EXCL_LINE - code will be removed at compile time
592
593 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
594
595 double discriminant = sqrt((moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2));
596
597 // Note that as the matrix of second moments of area is symmetric, both its eigenvalues are real
598 double largest_eigenvalue = (moments(0) + moments(1) + discriminant) * 0.5;
599 double smallest_eigenvalue = (moments(0) + moments(1) - discriminant) * 0.5;
600
601 double elongation_shape_factor = sqrt(largest_eigenvalue / smallest_eigenvalue);
602 return elongation_shape_factor;
603}
604
605template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
607{
608 mpDelaunayMesh = nullptr;
609 this->mMeshChangesDuringSimulation = false;
610 Clear();
611}
613template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
618
619template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
621{
622 assert(index < this->mNodes.size());
623 return index;
624}
625
626template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
628{
629 assert(index < this->mElements.size());
630 return index;
631}
632
633template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
635{
637 // assert(index < this->mBoundaryElements.size() );
638 return index;
639}
640
641template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
643{
644 unsigned node_index = UNSIGNED_UNSET;
645
646 if (mVoronoiElementIndexMap.empty())
647 {
648 node_index = elementIndex;
649 }
650 else
651 {
652 for (std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.begin();
653 iter != mVoronoiElementIndexMap.end();
654 ++iter)
655 {
656 if (iter->second == elementIndex)
657 {
658 node_index = iter->first;
659 break;
660 }
661 }
662 }
663 assert(node_index != UNSIGNED_UNSET);
664 return node_index;
665}
666
667template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
669{
670 unsigned element_index = UNSIGNED_UNSET;
671
672 if (mVoronoiElementIndexMap.empty())
673 {
674 element_index = nodeIndex;
675 }
676 else
677 {
678 std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.find(nodeIndex);
679
680 if (iter == mVoronoiElementIndexMap.end())
681 {
682 EXCEPTION("This index does not correspond to a VertexElement");
683 }
684 else
685 {
686 element_index = iter->second;
687 }
688 }
689 assert(element_index != UNSIGNED_UNSET);
690 return element_index;
691}
692
693template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
695{
696 assert(SPACE_DIM == 2 || SPACE_DIM == 3); // LCOV_EXCL_LINE - code will be removed at compile time
697
698 // Get pointer to this element
699 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
700
701 // Loop over nodes in the current element and find which is contained in the most elements
702 unsigned rosette_rank = 0;
703 for (unsigned node_idx = 0; node_idx < p_element->GetNumNodes(); node_idx++)
704 {
705 unsigned num_elems_this_node = p_element->GetNode(node_idx)->rGetContainingElementIndices().size();
706
707 if (num_elems_this_node > rosette_rank)
708 {
709 rosette_rank = num_elems_this_node;
710 }
711 }
712
713 // Return the rosette rank
714 return rosette_rank;
715}
716
717template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
719{
720 // Delete elements
721 for (unsigned i = 0; i < mElements.size(); i++)
722 {
723 delete mElements[i];
724 }
725 mElements.clear();
726
727 // Delete faces
728 for (unsigned i = 0; i < mFaces.size(); i++)
729 {
730 delete mFaces[i];
731 }
732 mFaces.clear();
733
734
735 // Delete nodes
736 for (unsigned i = 0; i < this->mNodes.size(); i++)
737 {
738 delete this->mNodes[i];
739 }
740 this->mNodes.clear();
741}
742
743template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
745{
746 return this->mNodes.size();
747}
748
749template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
751{
752 return mElements.size();
753}
754
755template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
757{
758 return mElements.size();
759}
760
761template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
763{
764 return mFaces.size();
765}
766
767template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
769{
770 assert(index < mElements.size());
771 return mElements[index];
772}
773
774template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
775VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetFace(unsigned index) const
776{
777 assert(index < mFaces.size());
778 return mFaces[index];
779}
780
781template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
782c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCentroidOfElement(unsigned index)
783{
784 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
785 unsigned num_nodes = p_element->GetNumNodes();
786
787 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
788
789 switch (SPACE_DIM)
790 {
791 case 1:
792 {
793 centroid = 0.5 * (p_element->GetNodeLocation(0) + p_element->GetNodeLocation(1));
794 }
795 break;
796 case 2:
797 {
798 double centroid_x = 0;
799 double centroid_y = 0;
800
801 // Note that we cannot use GetVolumeOfElement() below as it returns the absolute, rather than signed, area
802 double element_signed_area = 0.0;
803
804 // Map the first vertex to the origin and employ GetVectorFromAtoB() to allow for periodicity
805 c_vector<double, SPACE_DIM> first_node_location = p_element->GetNodeLocation(0);
806 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
807
808 // Loop over vertices
809 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
810 {
811 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation((local_index + 1) % num_nodes);
812 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
813
814 double this_x = pos_1[0];
815 double this_y = pos_1[1];
816 double next_x = pos_2[0];
817 double next_y = pos_2[1];
818
819 double signed_area_term = this_x * next_y - this_y * next_x;
820
821 centroid_x += (this_x + next_x) * signed_area_term;
822 centroid_y += (this_y + next_y) * signed_area_term;
823 element_signed_area += 0.5 * signed_area_term;
824
825 pos_1 = pos_2;
826 }
827
828 assert(element_signed_area != 0.0);
829
830 // Finally, map back and employ GetVectorFromAtoB() to allow for periodicity
831 centroid = first_node_location;
832 centroid(0) += centroid_x / (6.0 * element_signed_area);
833 centroid(1) += centroid_y / (6.0 * element_signed_area);
834 }
835 break;
836 case 3:
837 {
839 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
840 {
841 centroid += p_element->GetNodeLocation(local_index);
842 }
843 centroid /= ((double)num_nodes);
844 }
845 break;
846 default:
848 }
849 return centroid;
850}
851
852template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
854{
855 // Create a set of neighbouring node indices
856 std::set<unsigned> neighbouring_node_indices;
857
858 // Find the indices of the elements owned by this node
859 std::set<unsigned> containing_elem_indices = this->GetNode(nodeIndex)->rGetContainingElementIndices();
860
861 // Iterate over these elements
862 for (std::set<unsigned>::iterator elem_iter = containing_elem_indices.begin();
863 elem_iter != containing_elem_indices.end();
864 ++elem_iter)
865 {
866 // Find the local index of this node in this element
867 unsigned local_index = GetElement(*elem_iter)->GetNodeLocalIndex(nodeIndex);
868
869 // Find the global indices of the preceding and successive nodes in this element
870 unsigned num_nodes = GetElement(*elem_iter)->GetNumNodes();
871 unsigned previous_local_index = (local_index + num_nodes - 1) % num_nodes;
872 unsigned next_local_index = (local_index + 1) % num_nodes;
873
874 // Add the global indices of these two nodes to the set of neighbouring node indices
875 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(previous_local_index));
876 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(next_local_index));
877 }
878
879 return neighbouring_node_indices;
880}
881
882template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
883std::set<unsigned> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNeighbouringNodeNotAlsoInElement(unsigned nodeIndex, unsigned elemIndex)
884{
885 // Get a pointer to this element
886 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elemIndex);
887
888 // Get the indices of the nodes contained in this element
889 std::set<unsigned> node_indices_in_this_element;
890 for (unsigned local_index = 0; local_index < p_element->GetNumNodes(); local_index++)
891 {
892 unsigned global_index = p_element->GetNodeGlobalIndex(local_index);
893 node_indices_in_this_element.insert(global_index);
894 }
895
896 // Check that the node is in fact contained in the element
897 if (node_indices_in_this_element.find(nodeIndex) == node_indices_in_this_element.end())
898 {
899 EXCEPTION("The given node is not contained in the given element.");
900 }
901
902 // Create a set of neighbouring node indices
903 std::set<unsigned> neighbouring_node_indices_not_in_this_element;
904
905 // Get the indices of this node's neighbours
906 std::set<unsigned> node_neighbours = GetNeighbouringNodeIndices(nodeIndex);
907
908 // Check if each neighbour is also in this element; if not, add it to the set
909 for (std::set<unsigned>::iterator iter = node_neighbours.begin();
910 iter != node_neighbours.end();
911 ++iter)
912 {
913 if (node_indices_in_this_element.find(*iter) == node_indices_in_this_element.end())
914 {
915 neighbouring_node_indices_not_in_this_element.insert(*iter);
916 }
917 }
918
919 return neighbouring_node_indices_not_in_this_element;
920}
921
922template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
924{
925 // Get a pointer to this element
926 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
927
928 // Create a set of neighbouring element indices
929 std::set<unsigned> neighbouring_element_indices;
930
931 // Loop over nodes owned by this element
932 for (unsigned local_index = 0; local_index < p_element->GetNumNodes(); local_index++)
933 {
934 // Get a pointer to this node
935 Node<SPACE_DIM>* p_node = p_element->GetNode(local_index);
936
937 // Find the indices of the elements owned by this node
938 std::set<unsigned> containing_elem_indices = p_node->rGetContainingElementIndices();
939
940 // Form the union of this set with the current set of neighbouring element indices
941 std::set<unsigned> all_elements;
942 std::set_union(neighbouring_element_indices.begin(), neighbouring_element_indices.end(),
943 containing_elem_indices.begin(), containing_elem_indices.end(),
944 std::inserter(all_elements, all_elements.begin()));
945
946 // Update the set of neighbouring element indices
947 neighbouring_element_indices = all_elements;
948 }
949
950 // Lastly remove this element's index from the set of neighbouring element indices
951 neighbouring_element_indices.erase(elementIndex);
952
953 return neighbouring_element_indices;
954}
955
956template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
961
963template <>
966{
967 EXCEPTION("VertexMesh<1,1>::ConstructFromMeshReader() is not implemented");
968}
969
971template <>
974{
975 EXCEPTION("VertexMesh<1,2>::ConstructFromMeshReader() is not implemented");
976}
977
979template <>
982{
983 EXCEPTION("VertexMesh<1,3>::ConstructFromMeshReader() is not implemented");
984}
985
987template <>
990{
991 EXCEPTION("VertexMesh<2,3>::ConstructFromMeshReader() is not implemented");
992}
993
995template <>
998{
999 assert(rMeshReader.HasNodePermutation() == false);
1000 // Store numbers of nodes and elements
1001 unsigned num_nodes = rMeshReader.GetNumNodes();
1002 unsigned num_elements = rMeshReader.GetNumElements();
1003
1004 // Reserve memory for nodes
1005 this->mNodes.reserve(num_nodes);
1006
1007 rMeshReader.Reset();
1008
1009 // Add nodes
1010 std::vector<double> node_data;
1011 for (unsigned i = 0; i < num_nodes; i++)
1012 {
1013 node_data = rMeshReader.GetNextNode();
1014 unsigned is_boundary_node = (bool)node_data[2];
1015 node_data.pop_back();
1016 this->mNodes.push_back(new Node<2>(i, node_data, is_boundary_node));
1017 }
1018
1019 rMeshReader.Reset();
1020
1021 // Reserve memory for elements
1022 mElements.reserve(rMeshReader.GetNumElements());
1023
1024 // Add elements
1025 for (unsigned elem_index = 0; elem_index < num_elements; elem_index++)
1026 {
1027 // Get the data for this element
1028 ElementData element_data = rMeshReader.GetNextElementData();
1029
1030 // Get the nodes owned by this element
1031 std::vector<Node<2>*> nodes;
1032 unsigned num_nodes_in_element = element_data.NodeIndices.size();
1033 for (unsigned j = 0; j < num_nodes_in_element; j++)
1034 {
1035 assert(element_data.NodeIndices[j] < this->mNodes.size());
1036 nodes.push_back(this->mNodes[element_data.NodeIndices[j]]);
1037 }
1038
1039 // Use nodes and index to construct this element
1040 VertexElement<2, 2>* p_element = new VertexElement<2, 2>(elem_index, nodes);
1041 mElements.push_back(p_element);
1042
1043 if (rMeshReader.GetNumElementAttributes() > 0)
1044 {
1045 assert(rMeshReader.GetNumElementAttributes() == 1);
1046 unsigned attribute_value = (unsigned)element_data.AttributeValue;
1047 p_element->SetAttribute(attribute_value);
1048 }
1049 }
1050 GenerateEdgesFromElements(mElements);
1051}
1052
1054template <>
1057{
1058 assert(rMeshReader.HasNodePermutation() == false);
1059
1060 // Store numbers of nodes and elements
1061 unsigned num_nodes = rMeshReader.GetNumNodes();
1062 unsigned num_elements = rMeshReader.GetNumElements();
1063
1064 // Reserve memory for nodes
1065 this->mNodes.reserve(num_nodes);
1066
1067 rMeshReader.Reset();
1068
1069 // Add nodes
1070 std::vector<double> node_data;
1071 for (unsigned i = 0; i < num_nodes; i++)
1072 {
1073 node_data = rMeshReader.GetNextNode();
1074 unsigned is_boundary_node = (bool)node_data[3];
1075 node_data.pop_back();
1076 this->mNodes.push_back(new Node<3>(i, node_data, is_boundary_node));
1077 }
1078
1079 rMeshReader.Reset();
1080
1081 // Reserve memory for nodes
1082 mElements.reserve(rMeshReader.GetNumElements());
1083
1084 // Use a std::set to keep track of which faces have been added to mFaces
1085 std::set<unsigned> faces_counted;
1086
1087 // Add elements
1088 for (unsigned elem_index = 0; elem_index < num_elements; elem_index++)
1089 {
1091 typedef VertexMeshReader<3, 3> VERTEX_MESH_READER;
1092 assert(dynamic_cast<VERTEX_MESH_READER*>(&rMeshReader) != nullptr);
1093
1094 // Get the data for this element
1095 VertexElementData element_data = static_cast<VERTEX_MESH_READER*>(&rMeshReader)->GetNextElementDataWithFaces();
1096
1097 // Get the nodes owned by this element
1098 std::vector<Node<3>*> nodes;
1099 unsigned num_nodes_in_element = element_data.NodeIndices.size();
1100 for (unsigned j = 0; j < num_nodes_in_element; j++)
1101 {
1102 assert(element_data.NodeIndices[j] < this->mNodes.size());
1103 nodes.push_back(this->mNodes[element_data.NodeIndices[j]]);
1104 }
1105
1106 // Get the faces owned by this element
1107 std::vector<VertexElement<2, 3>*> faces;
1108 unsigned num_faces_in_element = element_data.Faces.size();
1109 for (unsigned i = 0; i < num_faces_in_element; i++)
1110 {
1111 // Get the data for this face
1112 ElementData face_data = element_data.Faces[i];
1113
1114 // Get the face index
1115 unsigned face_index = (unsigned)face_data.AttributeValue;
1116
1117 // Get the nodes owned by this face
1118 std::vector<Node<3>*> nodes_in_face;
1119 unsigned num_nodes_in_face = face_data.NodeIndices.size();
1120 for (unsigned j = 0; j < num_nodes_in_face; j++)
1121 {
1122 assert(face_data.NodeIndices[j] < this->mNodes.size());
1123 nodes_in_face.push_back(this->mNodes[face_data.NodeIndices[j]]);
1124 }
1125
1126 // If this face index is not already encountered, create a new face and update faces_counted...
1127 if (faces_counted.find(face_index) == faces_counted.end())
1128 {
1129 // Use nodes and index to construct this face
1130 VertexElement<2, 3>* p_face = new VertexElement<2, 3>(face_index, nodes_in_face);
1131 mFaces.push_back(p_face);
1132 faces_counted.insert(face_index);
1133 faces.push_back(p_face);
1134 }
1135 else
1136 {
1137 // ... otherwise use the member of mFaces with this index
1138 bool face_added = false;
1139 for (unsigned k = 0; k < mFaces.size(); k++)
1140 {
1141 if (mFaces[k]->GetIndex() == face_index)
1142 {
1143 faces.push_back(mFaces[k]);
1144 face_added = true;
1145 break;
1146 }
1147 }
1148 UNUSED_OPT(face_added);
1149 assert(face_added == true);
1150 }
1151 }
1152
1154 std::vector<bool> orientations = std::vector<bool>(num_faces_in_element, true);
1155
1156 // Use faces and index to construct this element
1157 VertexElement<3, 3>* p_element = new VertexElement<3, 3>(elem_index, faces, orientations, nodes);
1158 mElements.push_back(p_element);
1159
1160 if (rMeshReader.GetNumElementAttributes() > 0)
1161 {
1162 assert(rMeshReader.GetNumElementAttributes() == 1);
1163 unsigned attribute_value = element_data.AttributeValue;
1164 p_element->SetAttribute(attribute_value);
1165 }
1166 }
1167}
1168
1169template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1171 const c_vector<double, SPACE_DIM>& rLocationA, const c_vector<double, SPACE_DIM>& rLocationB)
1172{
1173 c_vector<double, SPACE_DIM> vector;
1174 if (mpDelaunayMesh)
1175 {
1176 vector = mpDelaunayMesh->GetVectorFromAtoB(rLocationA, rLocationB);
1177 }
1178 else
1179 {
1180 vector = AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetVectorFromAtoB(rLocationA, rLocationB);
1181 }
1182 return vector;
1183}
1184
1185template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1187{
1188 assert(SPACE_DIM == 2 || SPACE_DIM == 3); // LCOV_EXCL_LINE - code will be removed at compile time
1189
1190 // Get pointer to this element
1191 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1192
1193 double element_volume = 0.0;
1194 if (SPACE_DIM == 2)
1195 {
1196 // We map the first vertex to the origin and employ GetVectorFromAtoB() to allow for periodicity
1197 c_vector<double, SPACE_DIM> first_node_location = p_element->GetNodeLocation(0);
1198 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
1199
1200 unsigned num_nodes = p_element->GetNumNodes();
1201 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1202 {
1203 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation((local_index + 1) % num_nodes);
1204 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
1205
1206 double this_x = pos_1[0];
1207 double this_y = pos_1[1];
1208 double next_x = pos_2[0];
1209 double next_y = pos_2[1];
1210
1211 element_volume += 0.5 * (this_x * next_y - next_x * this_y);
1212
1213 pos_1 = pos_2;
1214 }
1215 }
1216 else
1217 {
1218 // Loop over faces and add up pyramid volumes
1219 c_vector<double, SPACE_DIM> pyramid_apex = p_element->GetNodeLocation(0);
1220 for (unsigned face_index = 0; face_index < p_element->GetNumFaces(); face_index++)
1221 {
1222 // Get pointer to face
1223 VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* p_face = p_element->GetFace(face_index);
1224
1225 // Calculate the area of the face and get unit normal to this face
1226 c_vector<double, SPACE_DIM> unit_normal;
1227 double face_area = CalculateUnitNormalToFaceWithArea(p_face, unit_normal);
1228
1229 // Calculate the perpendicular distance from the plane of the face to the chosen apex
1230 c_vector<double, SPACE_DIM> base_to_apex = GetVectorFromAtoB(p_face->GetNodeLocation(0), pyramid_apex);
1231 double perpendicular_distance = fabs(inner_prod(base_to_apex, unit_normal));
1232
1233 // Use these to calculate the volume of the pyramid formed by the face and the point pyramid_apex
1234 element_volume += face_area * perpendicular_distance / 3;
1235 }
1236 }
1237
1238 // We take the absolute value just in case the nodes were really oriented clockwise
1239 return fabs(element_volume);
1240}
1241
1242template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1244{
1245 assert(SPACE_DIM == 2 || SPACE_DIM == 3); // LCOV_EXCL_LINE - code will be removed at compile time
1246
1247 // Get pointer to this element
1248 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1249
1250 double surface_area = 0.0;
1251 if (SPACE_DIM == 2)
1252 {
1253 unsigned num_nodes = p_element->GetNumNodes();
1254 unsigned this_node_index = p_element->GetNodeGlobalIndex(0);
1255 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1256 {
1257 unsigned next_node_index = p_element->GetNodeGlobalIndex((local_index + 1) % num_nodes);
1258
1259 surface_area += this->GetDistanceBetweenNodes(this_node_index, next_node_index);
1260 this_node_index = next_node_index;
1261 }
1262 }
1263 else
1264 {
1265 // Loop over faces and add up areas
1266 for (unsigned face_index = 0; face_index < p_element->GetNumFaces(); face_index++)
1267 {
1268 surface_area += CalculateAreaOfFace(p_element->GetFace(face_index));
1269 }
1270 }
1271 return surface_area;
1272}
1273
1275// 2D-specific methods //
1277template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1279 [[maybe_unused]] const c_vector<double, SPACE_DIM>& rTestPoint, [[maybe_unused]] unsigned elementIndex)
1280{
1281 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1282 {
1283 // Get the element
1284 VertexElement<2, 2>* p_element = GetElement(elementIndex);
1285 const unsigned num_nodes = p_element->GetNumNodes();
1286
1287 // Initialise boolean
1288 bool element_includes_point = false;
1289
1291 // // Initialise boolean
1292 // bool element_includes_point = true;
1293 //
1294 // unsigned winding_number = 0;
1295 //
1296 // c_vector<double, SPACE_DIM> first_node_location = p_element->GetNodeLocation(0);
1297 // c_vector<double, SPACE_DIM> test_point = this->GetVectorFromAtoB(first_node_location, rTestPoint);
1298 // c_vector<double, SPACE_DIM> this_node_location = zero_vector<double>(SPACE_DIM);
1299 //
1300 // // Loop over edges of the element
1301 // for (unsigned local_index=0; local_index<num_nodes; local_index++)
1302 // {
1303 // c_vector<double, SPACE_DIM> untransformed_vector = p_element->GetNodeLocation((local_index+1)%num_nodes);
1304 // c_vector<double, SPACE_DIM> next_node_location = this->GetVectorFromAtoB(first_node_location, untransformed_vector);
1305 //
1306 // // If this edge is crossing upward...
1307 // if (this_node_location[1] <= test_point[1])
1308 // {
1309 // if (next_node_location[1] > test_point[1])
1310 // {
1311 // double is_left = (next_node_location[0] - this_node_location[0])*(test_point[1] - this_node_location[1])
1312 // - (test_point[0] - this_node_location[0])*(next_node_location[1] - this_node_location[1]);
1313 //
1314 // // ...and the test point is to the left of the edge...
1315 // if (is_left > DBL_EPSILON)
1316 // {
1317 // // ...then there is a valid upward edge-ray intersection to the right of the test point
1318 // winding_number++;
1319 // }
1320 // }
1321 // }
1322 // else
1323 // {
1324 // // ...otherwise, if the edge is crossing downward
1325 // if (next_node_location[1] <= test_point[1])
1326 // {
1327 // double is_left = (next_node_location[0] - this_node_location[0])*(test_point[1] - this_node_location[1])
1328 // - (test_point[0] - this_node_location[0])*(next_node_location[1] - this_node_location[1]);
1329 //
1330 // // ...and the test point is to the right of the edge...
1331 // if (is_left < -DBL_EPSILON)
1332 // {
1333 // // ...then there is a valid downward edge-ray intersection to the right of the test point
1334 // winding_number--;
1335 // }
1336 // }
1337 // }
1338 // this_node_location = next_node_location;
1339 // }
1340 //
1341 // if (winding_number == 0)
1342 // {
1343 // element_includes_point = false;
1344 // }
1346
1347 // Remap the origin to the first vertex to allow alternative distance metrics to be used in subclasses
1348 const c_vector<double, 2> first_vertex = p_element->GetNodeLocation(0);
1349 c_vector<double, 2> test_point = GetVectorFromAtoB(first_vertex, rTestPoint);
1350
1351 // Loop over edges of the element
1352 c_vector<double, 2> vertexA = zero_vector<double>(2);
1353 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1354 {
1355 // Check if this edge crosses the ray running out horizontally (increasing x, fixed y) from the test point
1356 c_vector<double, 2> vector_a_to_point = GetVectorFromAtoB(vertexA, test_point);
1357
1358 // Pathological case - test point coincides with vertexA
1359 // (we will check vertexB next time we go through the for loop)
1360 if (norm_2(vector_a_to_point) < DBL_EPSILON)
1361 {
1362 return false;
1363 }
1364
1365 c_vector<double, 2> vertexB = GetVectorFromAtoB(first_vertex, p_element->GetNodeLocation((local_index + 1) % num_nodes));
1366 c_vector<double, 2> vector_b_to_point = GetVectorFromAtoB(vertexB, test_point);
1367 c_vector<double, 2> vector_a_to_b = GetVectorFromAtoB(vertexA, vertexB);
1368
1369 // Pathological case - ray coincides with horizontal edge
1370 if ((fabs(vector_a_to_b[1]) < DBL_EPSILON) && (fabs(vector_a_to_point[1]) < DBL_EPSILON) && (fabs(vector_b_to_point[1]) < DBL_EPSILON))
1371 {
1372 if ((vector_a_to_point[0] > 0) != (vector_b_to_point[0] > 0))
1373 {
1374 return false;
1375 }
1376 }
1377
1378 // Non-pathological case
1379 // A and B on different sides of the line y = test_point[1]
1380 if ((vertexA[1] > test_point[1]) != (vertexB[1] > test_point[1]))
1381 {
1382 // Intersection of y=test_point[1] and vector_a_to_b is on the right of test_point
1383 if (test_point[0] < vertexA[0] + vector_a_to_b[0] * vector_a_to_point[1] / vector_a_to_b[1])
1384 {
1385 element_includes_point = !element_includes_point;
1386 }
1387 }
1388
1389 vertexA = vertexB;
1390 }
1391 return element_includes_point;
1392 }
1393 else
1394 {
1396 }
1397}
1398
1399template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1401 [[maybe_unused]] const c_vector<double, SPACE_DIM>& rTestPoint, [[maybe_unused]] unsigned elementIndex)
1402{
1403 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1404 {
1405 // Get the element
1406 VertexElement<2, 2>* p_element = GetElement(elementIndex);
1407 const unsigned num_nodes = p_element->GetNumNodes();
1408
1409 double min_squared_normal_distance = DBL_MAX;
1410 unsigned min_distance_edge_index = UINT_MAX;
1411
1412 // Loop over edges of the element
1413 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1414 {
1415 // Get the end points of this edge
1416 const c_vector<double, 2> vertexA = p_element->GetNodeLocation(local_index);
1417 const c_vector<double, 2> vertexB = p_element->GetNodeLocation((local_index + 1) % num_nodes);
1418
1419 const c_vector<double, 2> vector_a_to_point = this->GetVectorFromAtoB(vertexA, rTestPoint);
1420 const c_vector<double, 2> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
1421 const double distance_a_to_b = norm_2(vector_a_to_b);
1422
1423 const c_vector<double, 2> edge_ab_unit_vector = vector_a_to_b / norm_2(vector_a_to_b);
1424 double distance_parallel_to_edge = inner_prod(vector_a_to_point, edge_ab_unit_vector);
1425
1426 double squared_distance_normal_to_edge = SmallPow(norm_2(vector_a_to_point), 2) - SmallPow(distance_parallel_to_edge, 2);
1427
1428 /*
1429 * If the point lies almost bang on the supporting line of the edge, then snap to the line.
1430 * This allows us to do floating point tie-breaks when line is exactly at a node.
1431 * We adopt a similar approach if the point is at the same position as a point in the
1432 * element.
1433 */
1434 if (squared_distance_normal_to_edge < DBL_EPSILON)
1435 {
1436 squared_distance_normal_to_edge = 0.0;
1437 }
1438
1439 if (fabs(distance_parallel_to_edge) < DBL_EPSILON)
1440 {
1441 distance_parallel_to_edge = 0.0;
1442 }
1443 else if (fabs(distance_parallel_to_edge - distance_a_to_b) < DBL_EPSILON)
1444 {
1445 distance_parallel_to_edge = distance_a_to_b;
1446 }
1447
1448 // Make sure node is within the confines of the edge and is the nearest edge to the node \this breaks for convex elements
1449 if (squared_distance_normal_to_edge < min_squared_normal_distance && distance_parallel_to_edge >= 0 && distance_parallel_to_edge <= distance_a_to_b)
1450 {
1451 min_squared_normal_distance = squared_distance_normal_to_edge;
1452 min_distance_edge_index = local_index;
1453 }
1454 }
1455
1456 assert(min_distance_edge_index < num_nodes);
1457 return min_distance_edge_index;
1458 }
1459 else
1460 {
1462 }
1463}
1464
1465template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1466c_vector<double, 3> VertexMesh<ELEMENT_DIM, SPACE_DIM>::CalculateMomentsOfElement([[maybe_unused]] unsigned index)
1467{
1468 if constexpr (SPACE_DIM == 2)
1469 {
1470 // Define helper variables
1471 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1472 unsigned num_nodes = p_element->GetNumNodes();
1473 c_vector<double, 3> moments = zero_vector<double>(3);
1474
1475 // Since we compute I_xx, I_yy and I_xy about the centroid, we must shift each vertex accordingly
1476 c_vector<double, SPACE_DIM> centroid = GetCentroidOfElement(index);
1477
1478 c_vector<double, SPACE_DIM> this_node_location = p_element->GetNodeLocation(0);
1479 c_vector<double, SPACE_DIM> pos_1 = this->GetVectorFromAtoB(centroid, this_node_location);
1480
1481 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1482 {
1483 unsigned next_index = (local_index + 1) % num_nodes;
1484 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation(next_index);
1485 c_vector<double, SPACE_DIM> pos_2 = this->GetVectorFromAtoB(centroid, next_node_location);
1486
1487 double signed_area_term = pos_1(0) * pos_2(1) - pos_2(0) * pos_1(1);
1488 // Ixx
1489 moments(0) += (pos_1(1) * pos_1(1) + pos_1(1) * pos_2(1) + pos_2(1) * pos_2(1)) * signed_area_term;
1490
1491 // Iyy
1492 moments(1) += (pos_1(0) * pos_1(0) + pos_1(0) * pos_2(0) + pos_2(0) * pos_2(0)) * signed_area_term;
1493
1494 // Ixy
1495 moments(2) += (pos_1(0) * pos_2(1) + 2 * pos_1(0) * pos_1(1) + 2 * pos_2(0) * pos_2(1) + pos_2(0) * pos_1(1)) * signed_area_term;
1496
1497 pos_1 = pos_2;
1498 }
1499
1500 moments(0) /= 12;
1501 moments(1) /= 12;
1502 moments(2) /= 24;
1503
1504 /*
1505 * If the nodes owned by the element were supplied in a clockwise rather
1506 * than anticlockwise manner, or if this arose as a result of enforcing
1507 * periodicity, then our computed quantities will be the wrong sign, so
1508 * we need to fix this.
1509 */
1510 if (moments(0) < 0.0)
1511 {
1512 moments(0) = -moments(0);
1513 moments(1) = -moments(1);
1514 moments(2) = -moments(2);
1515 }
1516 return moments;
1517 }
1518 else
1519 {
1521 }
1522}
1523
1524template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1525c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetShortAxisOfElement([[maybe_unused]] unsigned index)
1526{
1527 if constexpr (SPACE_DIM == 2)
1528 {
1529
1530 c_vector<double, SPACE_DIM> short_axis = zero_vector<double>(SPACE_DIM);
1531
1532 // Calculate the moments of the element about its centroid (recall that I_xx and I_yy must be non-negative)
1533 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
1534
1535 // Normalise the moments vector to remove problem of a very small discriminant (see #2874)
1536 moments /= norm_2(moments);
1537
1538 // If the principal moments are equal...
1539 double discriminant = (moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2);
1540 if (fabs(discriminant) < DBL_EPSILON)
1541 {
1542 // ...then every axis through the centroid is a principal axis, so return a random unit vector
1543 short_axis(0) = RandomNumberGenerator::Instance()->ranf();
1544 short_axis(1) = sqrt(1.0 - short_axis(0) * short_axis(0));
1545 }
1546 else
1547 {
1548 // If the product of inertia is zero, then the coordinate axes are the principal axes
1549 if (fabs(moments(2)) < DBL_EPSILON)
1550 {
1551 if (moments(0) < moments(1))
1552 {
1553 short_axis(0) = 0.0;
1554 short_axis(1) = 1.0;
1555 }
1556 else
1557 {
1558 short_axis(0) = 1.0;
1559 short_axis(1) = 0.0;
1560 }
1561 }
1562 else
1563 {
1564 // Otherwise we find the eigenvector of the inertia matrix corresponding to the largest eigenvalue
1565 double lambda = 0.5 * (moments(0) + moments(1) + sqrt(discriminant));
1566
1567 short_axis(0) = 1.0;
1568 short_axis(1) = (moments(0) - lambda) / moments(2);
1569
1570 // Normalise the short axis before returning it
1571 short_axis /= norm_2(short_axis);
1572 }
1573 }
1574
1575 return short_axis;
1576 }
1577 else
1578 {
1580 }
1581}
1582
1583template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1585 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement, [[maybe_unused]] unsigned localIndex)
1586{
1587 if constexpr (SPACE_DIM == 2)
1588 {
1589 unsigned num_nodes_in_element = pElement->GetNumNodes();
1590 unsigned next_local_index = (localIndex + 1) % num_nodes_in_element;
1591
1592 // We add an extra num_nodes_in_element in the line below as otherwise this term can be negative, which breaks the % operator
1593 unsigned previous_local_index = (num_nodes_in_element + localIndex - 1) % num_nodes_in_element;
1594
1595 c_vector<double, SPACE_DIM> previous_node_location = pElement->GetNodeLocation(previous_local_index);
1596 c_vector<double, SPACE_DIM> next_node_location = pElement->GetNodeLocation(next_local_index);
1597 c_vector<double, SPACE_DIM> difference_vector = this->GetVectorFromAtoB(previous_node_location, next_node_location);
1598
1599 c_vector<double, SPACE_DIM> area_gradient;
1600
1601 area_gradient[0] = 0.5 * difference_vector[1];
1602 area_gradient[1] = -0.5 * difference_vector[0];
1603
1604 return area_gradient;
1605 }
1606 else
1607 {
1609 }
1610}
1611
1612template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1614 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement, [[maybe_unused]] unsigned localIndex)
1615{
1616 if constexpr (SPACE_DIM == 2)
1617 {
1618 assert(SPACE_DIM == 2); // LCOV_EXCL_LINE - code will be removed at compile time
1619
1620 const unsigned num_nodes_in_element = pElement->GetNumNodes();
1621
1622 // We add an extra num_nodes_in_element-1 in the line below as otherwise this term can be negative, which breaks the % operator
1623 const unsigned previous_local_index = (num_nodes_in_element + localIndex - 1) % num_nodes_in_element;
1624
1625 const unsigned this_global_index = pElement->GetNodeGlobalIndex(localIndex);
1626 const unsigned previous_global_index = pElement->GetNodeGlobalIndex(previous_local_index);
1627
1628 const double previous_edge_length = this->GetDistanceBetweenNodes(this_global_index, previous_global_index);
1629 assert(previous_edge_length > DBL_EPSILON);
1630
1631 const c_vector<double, SPACE_DIM> previous_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(previous_local_index), pElement->GetNodeLocation(localIndex)) / previous_edge_length;
1632
1633 return previous_edge_gradient;
1634 }
1635 else
1636 {
1638 }
1639}
1640
1641template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1643 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement, [[maybe_unused]] unsigned localIndex)
1644{
1645 if constexpr (SPACE_DIM == 2)
1646 {
1647 const unsigned next_local_index = (localIndex + 1) % (pElement->GetNumNodes());
1648
1649 const unsigned this_global_index = pElement->GetNodeGlobalIndex(localIndex);
1650 const unsigned next_global_index = pElement->GetNodeGlobalIndex(next_local_index);
1651
1652 const double next_edge_length = this->GetDistanceBetweenNodes(this_global_index, next_global_index);
1653 assert(next_edge_length > DBL_EPSILON);
1654
1655 const c_vector<double, SPACE_DIM> next_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(next_local_index), pElement->GetNodeLocation(localIndex)) / next_edge_length;
1656
1657 return next_edge_gradient;
1658 }
1659 else
1660 {
1662 }
1663}
1664
1665template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1667 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement, [[maybe_unused]] unsigned localIndex)
1668{
1669 if constexpr (SPACE_DIM == 2)
1670 {
1671
1672 c_vector<double, SPACE_DIM> previous_edge_gradient = GetPreviousEdgeGradientOfElementAtNode(pElement, localIndex);
1673 c_vector<double, SPACE_DIM> next_edge_gradient = GetNextEdgeGradientOfElementAtNode(pElement, localIndex);
1674
1675 return previous_edge_gradient + next_edge_gradient;
1676 }
1677 else
1678 {
1680 }
1681}
1682
1684// 3D-specific methods //
1686
1687template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1689 [[maybe_unused]] VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* pFace,
1690 [[maybe_unused]] c_vector<double, SPACE_DIM>& rNormal)
1691{
1692 if constexpr (SPACE_DIM == 3)
1693 {
1694 assert(SPACE_DIM == 3); // LCOV_EXCL_LINE - code will be removed at compile time
1695
1696 // As we are in 3D, the face must have at least three vertices
1697 assert(pFace->GetNumNodes() >= 3u);
1698
1699 // Reset the answer
1700 rNormal = zero_vector<double>(SPACE_DIM);
1701
1702 c_vector<double, SPACE_DIM> v_minus_v0 = this->GetVectorFromAtoB(pFace->GetNode(0)->rGetLocation(), pFace->GetNode(1)->rGetLocation());
1703 for (unsigned local_index = 2; local_index < pFace->GetNumNodes(); local_index++)
1704 {
1705 c_vector<double, SPACE_DIM> vnext_minus_v0 = this->GetVectorFromAtoB(pFace->GetNode(0)->rGetLocation(), pFace->GetNode(local_index)->rGetLocation());
1706 rNormal += VectorProduct(v_minus_v0, vnext_minus_v0);
1707 v_minus_v0 = vnext_minus_v0;
1708 }
1709 double magnitude = norm_2(rNormal);
1710 if (magnitude != 0.0)
1711 {
1712 // Normalize the normal vector
1713 rNormal /= magnitude;
1714 // If all points are co-located, then magnitude==0.0 and there is potential for a floating point exception
1715 // here if we divide by zero, so we'll move on.
1716 }
1717 return magnitude / 2.0;
1718 }
1719 else
1720 {
1722 }
1723}
1724
1725template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1727 [[maybe_unused]] VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* pFace)
1728{
1729 if constexpr (SPACE_DIM == 3)
1730 {
1731 // Get the unit normal to the plane of this face
1732 c_vector<double, SPACE_DIM> unit_normal;
1733 return CalculateUnitNormalToFaceWithArea(pFace, unit_normal);
1734 }
1735 else
1736 {
1738 }
1739}
1740
1741
1742
1744#if defined(__xlC__)
1745template <>
1747{
1749}
1750#endif
1751
1752// Explicit instantiation
1753template class VertexMesh<1, 1>;
1754template class VertexMesh<1, 2>;
1755template class VertexMesh<1, 3>;
1756template class VertexMesh<2, 2>;
1757template class VertexMesh<2, 3>;
1758template class VertexMesh<3, 3>;
1759
1760// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define UNUSED_OPT(var)
#define NEVER_REACHED
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
void SetAttribute(double attribute)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetIndex() const
virtual void Reset()=0
virtual unsigned GetNumElements() const =0
virtual ElementData GetNextElementData()=0
virtual unsigned GetNumElementAttributes() const
virtual std::vector< double > GetNextNode()=0
virtual bool HasNodePermutation()
virtual unsigned GetNumNodes() const =0
bool mMeshChangesDuringSimulation
virtual c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocationA, const c_vector< double, SPACE_DIM > &rLocationB)
std::vector< Node< SPACE_DIM > * > mNodes
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
virtual unsigned GetNumElements() const
Definition Edge.hpp:52
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
Definition Node.hpp:59
std::set< unsigned > & rGetContainingElementIndices()
Definition Node.cpp:300
void AddElement(unsigned index)
Definition Node.cpp:268
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
bool IsBoundaryNode() const
Definition Node.cpp:164
unsigned GetIndex() const
Definition Node.cpp:158
static RandomNumberGenerator * Instance()
virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian) const
VertexElement< ELEMENT_DIM-1, SPACE_DIM > * GetFace(unsigned index) const
void AddFace(VertexElement< ELEMENT_DIM-1, SPACE_DIM > *pFace)
unsigned GetNumFaces() const
std::vector< VertexElement< ELEMENT_DIM - 1, SPACE_DIM > * > mFaces
virtual c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocationA, const c_vector< double, SPACE_DIM > &rLocationB)
double CalculateUnitNormalToFaceWithArea(VertexElement< ELEMENT_DIM - 1, SPACE_DIM > *pFace, c_vector< double, SPACE_DIM > &rNormal)
VertexElement< ELEMENT_DIM - 1, SPACE_DIM > * GetFace(unsigned index) const
virtual unsigned GetNumNodes() const
virtual double GetSurfaceAreaOfElement(unsigned index)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned nodeIndex)
std::set< unsigned > GetNeighbouringNodeNotAlsoInElement(unsigned nodeIndex, unsigned elemIndex)
void GenerateEdgesFromElements(std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > &rElements)
virtual double CalculateAreaOfFace(VertexElement< ELEMENT_DIM - 1, SPACE_DIM > *pFace)
const EdgeHelper< SPACE_DIM > & rGetEdgeHelper() const
c_vector< double, SPACE_DIM > GetNextEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
virtual unsigned GetNumFaces() const
VertexElement< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
unsigned GetLocalIndexForElementEdgeClosestToPoint(const c_vector< double, SPACE_DIM > &rTestPoint, unsigned elementIndex)
double GetElongationShapeFactorOfElement(unsigned elementIndex)
unsigned GetNumAllElements() const
unsigned SolveNodeMapping(unsigned index) const
TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpDelaunayMesh
bool ElementIncludesPoint(const c_vector< double, SPACE_DIM > &rTestPoint, unsigned elementIndex)
c_vector< double, SPACE_DIM > GetAreaGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
std::set< unsigned > GetNeighbouringElementIndices(unsigned elementIndex)
virtual ~VertexMesh()
c_vector< double, SPACE_DIM > GetPerimeterGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
double GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)
c_vector< double, SPACE_DIM > GetShortAxisOfElement(unsigned index)
unsigned GetRosetteRankOfElement(unsigned index)
Edge< SPACE_DIM > * GetEdge(unsigned index) const
unsigned SolveBoundaryElementMapping(unsigned index) const
unsigned SolveElementMapping(unsigned index) const
unsigned GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(unsigned nodeIndex)
virtual c_vector< double, SPACE_DIM > GetCentroidOfElement(unsigned index)
std::map< unsigned, unsigned > mVoronoiElementIndexMap
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
unsigned GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(unsigned elementIndex)
virtual void Clear()
unsigned GetNumEdges() const
c_vector< double, SPACE_DIM > GetPreviousEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > mElements
virtual c_vector< double, 3 > CalculateMomentsOfElement(unsigned index)
virtual VertexMesh< ELEMENT_DIM, SPACE_DIM > * GetMeshForVtk()
virtual double GetVolumeOfElement(unsigned index)
virtual unsigned GetNumElements() const
void GenerateVerticesFromElementCircumcentres(TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
std::vector< unsigned > NodeIndices
std::vector< ElementData > Faces
std::vector< unsigned > NodeIndices