36#include "VertexMesh.hpp"
37#include "MutableMesh.hpp"
38#include "RandomNumberGenerator.hpp"
41template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
44 : mpDelaunayMesh(nullptr)
51 for (
unsigned node_index = 0; node_index < nodes.size(); node_index++)
54 this->
mNodes.push_back(p_temp_node);
56 for (
unsigned elem_index = 0; elem_index < vertexElements.size(); elem_index++)
59 mElements.push_back(p_temp_vertex_element);
66 std::set<unsigned> faces_counted;
69 for (
unsigned elem_index = 0; elem_index <
mElements.size(); elem_index++)
72 for (
unsigned face_index = 0; face_index <
mElements[elem_index]->GetNumFaces(); face_index++)
75 unsigned global_index = p_face->
GetIndex();
78 if (faces_counted.find(global_index) == faces_counted.end())
81 faces_counted.insert(global_index);
88 for (
unsigned index = 0; index <
mElements.size(); index++)
92 unsigned element_index = p_element->
GetIndex();
93 unsigned num_nodes_in_element = p_element->
GetNumNodes();
95 for (
unsigned node_index = 0; node_index < num_nodes_in_element; node_index++)
97 p_element->
GetNode(node_index)->AddElement(element_index);
106template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
110 : mpDelaunayMesh(nullptr)
116 for (
unsigned node_index = 0; node_index < nodes.size(); node_index++)
119 this->
mNodes.push_back(p_temp_node);
122 for (
unsigned face_index = 0; face_index < faces.size(); face_index++)
124 VertexElement<ELEMENT_DIM - 1, SPACE_DIM>* p_temp_face = faces[face_index];
125 mFaces.push_back(p_temp_face);
128 for (
unsigned elem_index = 0; elem_index < vertexElements.size(); elem_index++)
131 mElements.push_back(p_temp_vertex_element);
135 for (
unsigned index = 0; index <
mElements.size(); index++)
138 for (
unsigned node_index = 0; node_index < p_temp_vertex_element->
GetNumNodes(); node_index++)
154 : mpDelaunayMesh(&rMesh)
167 this->
mNodes.reserve(num_nodes);
171 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
181 for (
unsigned i = 0; i < num_nodes; i++)
185 for (
unsigned local_index = 0; local_index < 3; local_index++)
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;
198 std::vector<Node<2> *> nodes;
203 nodes.push_back(
new Node<2>(node_iter->GetIndex(), node_iter->rGetLocation(),node_iter->IsBoundaryNode()));
212 for (
unsigned j=0; j<3; j++)
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()));
234 if (shared_elements.size() == 1)
237 double edge_length = norm_2(edge);
238 c_vector<double,2> normal_vector;
240 normal_vector[0]= edge[1];
241 normal_vector[1]= -edge[0];
243 double dij = norm_2(normal_vector);
245 normal_vector /= dij;
247 double extra_node_scaling = 1.0;
249 int num_sections = ceil(edge_length*extra_node_scaling);
250 for (
int section=0; section<=num_sections; section++)
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));
263 unsigned num_nodes = extended_mesh.GetNumAllElements();
266 this->
mNodes.reserve(num_nodes);
270 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
280 for (
unsigned i = 0; i < num_nodes; i++)
284 for (
unsigned local_index = 0; local_index < 3; local_index++)
286 unsigned elem_index = extended_mesh.GetElement(i)->GetNodeGlobalIndex(local_index);
288 if (elem_index < num_elements)
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;
300 for (
unsigned elem_index = 0; elem_index <
mElements.size(); elem_index++)
307 std::vector<std::pair<double, unsigned> > index_angle_list;
308 for (
unsigned local_index = 0; local_index <
mElements[elem_index]->GetNumNodes(); local_index++)
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);
314 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
315 unsigned global_index =
mElements[elem_index]->GetNodeGlobalIndex(local_index);
317 std::pair<double, unsigned> pair(angle, global_index);
318 index_angle_list.push_back(pair);
322 sort(index_angle_list.begin(), index_angle_list.end());
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);
350 : mpDelaunayMesh(&rMesh)
358 this->
mNodes.reserve(num_nodes);
363 std::map<unsigned, VertexElement<3, 3>*> index_element_map;
364 unsigned face_index = 0;
365 unsigned element_index = 0;
372 Node<3>* p_node_a = edge_iterator.GetNodeA();
373 Node<3>* p_node_b = edge_iterator.GetNodeB();
379 std::set<unsigned> edge_element_indices;
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()));
387 c_vector<double, 3> edge_vector;
390 c_vector<double, 3> mid_edge;
393 unsigned element0_index = *(edge_element_indices.begin());
395 c_vector<double, 3> basis_vector1;
396 basis_vector1 =
mNodes[element0_index]->rGetLocation() - mid_edge;
398 c_vector<double, 3> basis_vector2;
406 std::vector<std::pair<double, unsigned> > index_angle_list;
409 for (std::set<unsigned>::iterator index_iter = edge_element_indices.begin();
410 index_iter != edge_element_indices.end();
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);
419 double angle = atan2(local_vertex_dot_basis_vector2, local_vertex_dot_basis_vector1);
421 std::pair<double, unsigned> pair(angle, *index_iter);
422 index_angle_list.push_back(pair);
426 sort(index_angle_list.begin(), index_angle_list.end());
431 for (
unsigned count = 0; count < index_angle_list.size(); count++)
433 unsigned local_index = count > 1 ? count - 1 : 0;
434 p_face->
AddNode(
mNodes[index_angle_list[count].second], local_index);
443 unsigned node_a_index = p_node_a->
GetIndex();
444 if (index_element_map[node_a_index])
447 index_element_map[node_a_index]->AddFace(p_face);
456 index_element_map[node_a_index] = p_element;
461 unsigned node_b_index = p_node_b->
GetIndex();
462 if (index_element_map[node_b_index])
465 index_element_map[node_b_index]->AddFace(p_face);
474 index_element_map[node_b_index] = p_element;
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();
486 mElements.push_back(element_iter->second);
498template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
503 for (
auto elem : rElements)
505 elem->SetEdgeHelper(&mEdgeHelper);
510template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
513 return mEdgeHelper.GetNumEdges();
516template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
519 return mEdgeHelper.GetEdge(index);
522template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
528template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
531 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
532 c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
540 c_vector<double, SPACE_DIM + 1> circumsphere = rMesh.
GetElement(i)->CalculateCircumsphere(jacobian, inverse_jacobian);
542 c_vector<double, SPACE_DIM> circumcentre;
543 for (
unsigned j = 0; j < SPACE_DIM; j++)
545 circumcentre(j) = circumsphere(j);
553template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
556 assert(SPACE_DIM == 2);
558 std::set<unsigned> node_indices_1;
559 for (
unsigned i = 0; i < mElements[elementIndex1]->GetNumNodes(); i++)
561 node_indices_1.insert(mElements[elementIndex1]->GetNodeGlobalIndex(i));
563 std::set<unsigned> node_indices_2;
564 for (
unsigned i = 0; i < mElements[elementIndex2]->GetNumNodes(); i++)
566 node_indices_2.insert(mElements[elementIndex2]->GetNodeGlobalIndex(i));
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()));
574 if (shared_nodes.size() == 1)
577 EXCEPTION(
"Elements " << elementIndex1 <<
" and " << elementIndex2 <<
" share only one node.");
579 assert(shared_nodes.size() == 2);
581 unsigned index1 = *(shared_nodes.begin());
582 unsigned index2 = *(++(shared_nodes.begin()));
584 double edge_length = this->GetDistanceBetweenNodes(index1, index2);
588template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
591 assert(SPACE_DIM == 2);
593 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
595 double discriminant = sqrt((moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2));
598 double largest_eigenvalue = (moments(0) + moments(1) + discriminant) * 0.5;
599 double smallest_eigenvalue = (moments(0) + moments(1) - discriminant) * 0.5;
601 double elongation_shape_factor = sqrt(largest_eigenvalue / smallest_eigenvalue);
602 return elongation_shape_factor;
605template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
608 mpDelaunayMesh =
nullptr;
609 this->mMeshChangesDuringSimulation =
false;
613template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
619template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
622 assert(index < this->mNodes.size());
626template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
629 assert(index < this->mElements.size());
633template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
641template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
646 if (mVoronoiElementIndexMap.empty())
648 node_index = elementIndex;
652 for (std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.begin();
653 iter != mVoronoiElementIndexMap.end();
656 if (iter->second == elementIndex)
658 node_index = iter->first;
667template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
672 if (mVoronoiElementIndexMap.empty())
674 element_index = nodeIndex;
678 std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.find(nodeIndex);
680 if (iter == mVoronoiElementIndexMap.end())
682 EXCEPTION(
"This index does not correspond to a VertexElement");
686 element_index = iter->second;
690 return element_index;
693template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
696 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
702 unsigned rosette_rank = 0;
703 for (
unsigned node_idx = 0; node_idx < p_element->
GetNumNodes(); node_idx++)
705 unsigned num_elems_this_node = p_element->
GetNode(node_idx)->rGetContainingElementIndices().size();
707 if (num_elems_this_node > rosette_rank)
709 rosette_rank = num_elems_this_node;
717template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
721 for (
unsigned i = 0; i < mElements.size(); i++)
728 for (
unsigned i = 0; i < mFaces.size(); i++)
736 for (
unsigned i = 0; i < this->mNodes.size(); i++)
738 delete this->mNodes[i];
740 this->mNodes.clear();
743template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
746 return this->mNodes.size();
749template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
752 return mElements.size();
755template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
758 return mElements.size();
761template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
764 return mFaces.size();
767template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
770 assert(index < mElements.size());
771 return mElements[index];
774template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
777 assert(index < mFaces.size());
778 return mFaces[index];
781template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
787 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
798 double centroid_x = 0;
799 double centroid_y = 0;
802 double element_signed_area = 0.0;
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);
809 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
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);
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];
819 double signed_area_term = this_x * next_y - this_y * next_x;
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;
828 assert(element_signed_area != 0.0);
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);
839 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
843 centroid /= ((
double)num_nodes);
852template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
856 std::set<unsigned> neighbouring_node_indices;
859 std::set<unsigned> containing_elem_indices = this->GetNode(nodeIndex)->rGetContainingElementIndices();
862 for (std::set<unsigned>::iterator elem_iter = containing_elem_indices.begin();
863 elem_iter != containing_elem_indices.end();
867 unsigned local_index = GetElement(*elem_iter)->GetNodeLocalIndex(nodeIndex);
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;
875 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(previous_local_index));
876 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(next_local_index));
879 return neighbouring_node_indices;
882template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
889 std::set<unsigned> node_indices_in_this_element;
890 for (
unsigned local_index = 0; local_index < p_element->
GetNumNodes(); local_index++)
893 node_indices_in_this_element.insert(global_index);
897 if (node_indices_in_this_element.find(nodeIndex) == node_indices_in_this_element.end())
899 EXCEPTION(
"The given node is not contained in the given element.");
903 std::set<unsigned> neighbouring_node_indices_not_in_this_element;
906 std::set<unsigned> node_neighbours = GetNeighbouringNodeIndices(nodeIndex);
909 for (std::set<unsigned>::iterator iter = node_neighbours.begin();
910 iter != node_neighbours.end();
913 if (node_indices_in_this_element.find(*iter) == node_indices_in_this_element.end())
915 neighbouring_node_indices_not_in_this_element.insert(*iter);
919 return neighbouring_node_indices_not_in_this_element;
922template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
929 std::set<unsigned> neighbouring_element_indices;
932 for (
unsigned local_index = 0; local_index < p_element->
GetNumNodes(); local_index++)
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()));
947 neighbouring_element_indices = all_elements;
951 neighbouring_element_indices.erase(elementIndex);
953 return neighbouring_element_indices;
956template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
967 EXCEPTION(
"VertexMesh<1,1>::ConstructFromMeshReader() is not implemented");
975 EXCEPTION(
"VertexMesh<1,2>::ConstructFromMeshReader() is not implemented");
983 EXCEPTION(
"VertexMesh<1,3>::ConstructFromMeshReader() is not implemented");
991 EXCEPTION(
"VertexMesh<2,3>::ConstructFromMeshReader() is not implemented");
1005 this->mNodes.reserve(num_nodes);
1007 rMeshReader.
Reset();
1010 std::vector<double> node_data;
1011 for (
unsigned i = 0; i < num_nodes; i++)
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));
1019 rMeshReader.
Reset();
1025 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
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++)
1035 assert(element_data.
NodeIndices[j] < this->mNodes.size());
1036 nodes.push_back(this->mNodes[element_data.
NodeIndices[j]]);
1041 mElements.push_back(p_element);
1050 GenerateEdgesFromElements(mElements);
1065 this->mNodes.reserve(num_nodes);
1067 rMeshReader.
Reset();
1070 std::vector<double> node_data;
1071 for (
unsigned i = 0; i < num_nodes; i++)
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));
1079 rMeshReader.
Reset();
1085 std::set<unsigned> faces_counted;
1088 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
1092 assert(
dynamic_cast<VERTEX_MESH_READER*
>(&rMeshReader) !=
nullptr);
1095 VertexElementData element_data =
static_cast<VERTEX_MESH_READER*
>(&rMeshReader)->GetNextElementDataWithFaces();
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++)
1102 assert(element_data.
NodeIndices[j] < this->mNodes.size());
1103 nodes.push_back(this->mNodes[element_data.
NodeIndices[j]]);
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++)
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++)
1122 assert(face_data.NodeIndices[j] < this->mNodes.size());
1123 nodes_in_face.push_back(this->mNodes[face_data.NodeIndices[j]]);
1127 if (faces_counted.find(face_index) == faces_counted.end())
1131 mFaces.push_back(p_face);
1132 faces_counted.insert(face_index);
1133 faces.push_back(p_face);
1138 bool face_added =
false;
1139 for (
unsigned k = 0; k < mFaces.size(); k++)
1141 if (mFaces[k]->GetIndex() == face_index)
1143 faces.push_back(mFaces[k]);
1149 assert(face_added ==
true);
1154 std::vector<bool> orientations = std::vector<bool>(num_faces_in_element,
true);
1158 mElements.push_back(p_element);
1169template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1171 const c_vector<double, SPACE_DIM>& rLocationA,
const c_vector<double, SPACE_DIM>& rLocationB)
1173 c_vector<double, SPACE_DIM> vector;
1176 vector = mpDelaunayMesh->GetVectorFromAtoB(rLocationA, rLocationB);
1185template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1188 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
1193 double element_volume = 0.0;
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);
1201 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
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);
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];
1211 element_volume += 0.5 * (this_x * next_y - next_x * this_y);
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++)
1226 c_vector<double, SPACE_DIM> unit_normal;
1227 double face_area = CalculateUnitNormalToFaceWithArea(p_face, unit_normal);
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));
1234 element_volume += face_area * perpendicular_distance / 3;
1239 return fabs(element_volume);
1242template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1245 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
1250 double surface_area = 0.0;
1255 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
1257 unsigned next_node_index = p_element->
GetNodeGlobalIndex((local_index + 1) % num_nodes);
1259 surface_area += this->GetDistanceBetweenNodes(this_node_index, next_node_index);
1260 this_node_index = next_node_index;
1266 for (
unsigned face_index = 0; face_index < p_element->
GetNumFaces(); face_index++)
1268 surface_area += CalculateAreaOfFace(p_element->
GetFace(face_index));
1271 return surface_area;
1277template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1279 [[maybe_unused]]
const c_vector<double, SPACE_DIM>& rTestPoint, [[maybe_unused]]
unsigned elementIndex)
1281 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1285 const unsigned num_nodes = p_element->
GetNumNodes();
1288 bool element_includes_point =
false;
1348 const c_vector<double, 2> first_vertex = p_element->
GetNodeLocation(0);
1349 c_vector<double, 2> test_point = GetVectorFromAtoB(first_vertex, rTestPoint);
1352 c_vector<double, 2> vertexA = zero_vector<double>(2);
1353 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
1356 c_vector<double, 2> vector_a_to_point = GetVectorFromAtoB(vertexA, test_point);
1360 if (norm_2(vector_a_to_point) < DBL_EPSILON)
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);
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))
1372 if ((vector_a_to_point[0] > 0) != (vector_b_to_point[0] > 0))
1380 if ((vertexA[1] > test_point[1]) != (vertexB[1] > test_point[1]))
1383 if (test_point[0] < vertexA[0] + vector_a_to_b[0] * vector_a_to_point[1] / vector_a_to_b[1])
1385 element_includes_point = !element_includes_point;
1391 return element_includes_point;
1399template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1401 [[maybe_unused]]
const c_vector<double, SPACE_DIM>& rTestPoint, [[maybe_unused]]
unsigned elementIndex)
1403 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1407 const unsigned num_nodes = p_element->
GetNumNodes();
1409 double min_squared_normal_distance = DBL_MAX;
1410 unsigned min_distance_edge_index = UINT_MAX;
1413 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
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);
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);
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);
1426 double squared_distance_normal_to_edge = SmallPow(norm_2(vector_a_to_point), 2) - SmallPow(distance_parallel_to_edge, 2);
1434 if (squared_distance_normal_to_edge < DBL_EPSILON)
1436 squared_distance_normal_to_edge = 0.0;
1439 if (fabs(distance_parallel_to_edge) < DBL_EPSILON)
1441 distance_parallel_to_edge = 0.0;
1443 else if (fabs(distance_parallel_to_edge - distance_a_to_b) < DBL_EPSILON)
1445 distance_parallel_to_edge = distance_a_to_b;
1449 if (squared_distance_normal_to_edge < min_squared_normal_distance && distance_parallel_to_edge >= 0 && distance_parallel_to_edge <= distance_a_to_b)
1451 min_squared_normal_distance = squared_distance_normal_to_edge;
1452 min_distance_edge_index = local_index;
1456 assert(min_distance_edge_index < num_nodes);
1457 return min_distance_edge_index;
1465template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1468 if constexpr (SPACE_DIM == 2)
1473 c_vector<double, 3> moments = zero_vector<double>(3);
1476 c_vector<double, SPACE_DIM> centroid = GetCentroidOfElement(index);
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);
1481 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
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);
1487 double signed_area_term = pos_1(0) * pos_2(1) - pos_2(0) * pos_1(1);
1489 moments(0) += (pos_1(1) * pos_1(1) + pos_1(1) * pos_2(1) + pos_2(1) * pos_2(1)) * signed_area_term;
1492 moments(1) += (pos_1(0) * pos_1(0) + pos_1(0) * pos_2(0) + pos_2(0) * pos_2(0)) * signed_area_term;
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;
1510 if (moments(0) < 0.0)
1512 moments(0) = -moments(0);
1513 moments(1) = -moments(1);
1514 moments(2) = -moments(2);
1524template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1527 if constexpr (SPACE_DIM == 2)
1530 c_vector<double, SPACE_DIM> short_axis = zero_vector<double>(SPACE_DIM);
1533 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
1536 moments /= norm_2(moments);
1539 double discriminant = (moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2);
1540 if (fabs(discriminant) < DBL_EPSILON)
1544 short_axis(1) = sqrt(1.0 - short_axis(0) * short_axis(0));
1549 if (fabs(moments(2)) < DBL_EPSILON)
1551 if (moments(0) < moments(1))
1553 short_axis(0) = 0.0;
1554 short_axis(1) = 1.0;
1558 short_axis(0) = 1.0;
1559 short_axis(1) = 0.0;
1565 double lambda = 0.5 * (moments(0) + moments(1) + sqrt(discriminant));
1567 short_axis(0) = 1.0;
1568 short_axis(1) = (moments(0) - lambda) / moments(2);
1571 short_axis /= norm_2(short_axis);
1583template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1587 if constexpr (SPACE_DIM == 2)
1589 unsigned num_nodes_in_element = pElement->GetNumNodes();
1590 unsigned next_local_index = (localIndex + 1) % num_nodes_in_element;
1593 unsigned previous_local_index = (num_nodes_in_element + localIndex - 1) % num_nodes_in_element;
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);
1599 c_vector<double, SPACE_DIM> area_gradient;
1601 area_gradient[0] = 0.5 * difference_vector[1];
1602 area_gradient[1] = -0.5 * difference_vector[0];
1604 return area_gradient;
1612template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1616 if constexpr (SPACE_DIM == 2)
1618 assert(SPACE_DIM == 2);
1620 const unsigned num_nodes_in_element = pElement->GetNumNodes();
1623 const unsigned previous_local_index = (num_nodes_in_element + localIndex - 1) % num_nodes_in_element;
1625 const unsigned this_global_index = pElement->GetNodeGlobalIndex(localIndex);
1626 const unsigned previous_global_index = pElement->GetNodeGlobalIndex(previous_local_index);
1628 const double previous_edge_length = this->GetDistanceBetweenNodes(this_global_index, previous_global_index);
1629 assert(previous_edge_length > DBL_EPSILON);
1631 const c_vector<double, SPACE_DIM> previous_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(previous_local_index), pElement->GetNodeLocation(localIndex)) / previous_edge_length;
1633 return previous_edge_gradient;
1641template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1645 if constexpr (SPACE_DIM == 2)
1647 const unsigned next_local_index = (localIndex + 1) % (pElement->GetNumNodes());
1649 const unsigned this_global_index = pElement->GetNodeGlobalIndex(localIndex);
1650 const unsigned next_global_index = pElement->GetNodeGlobalIndex(next_local_index);
1652 const double next_edge_length = this->GetDistanceBetweenNodes(this_global_index, next_global_index);
1653 assert(next_edge_length > DBL_EPSILON);
1655 const c_vector<double, SPACE_DIM> next_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(next_local_index), pElement->GetNodeLocation(localIndex)) / next_edge_length;
1657 return next_edge_gradient;
1665template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1669 if constexpr (SPACE_DIM == 2)
1672 c_vector<double, SPACE_DIM> previous_edge_gradient = GetPreviousEdgeGradientOfElementAtNode(pElement, localIndex);
1673 c_vector<double, SPACE_DIM> next_edge_gradient = GetNextEdgeGradientOfElementAtNode(pElement, localIndex);
1675 return previous_edge_gradient + next_edge_gradient;
1687template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1690 [[maybe_unused]] c_vector<double, SPACE_DIM>& rNormal)
1692 if constexpr (SPACE_DIM == 3)
1694 assert(SPACE_DIM == 3);
1697 assert(pFace->GetNumNodes() >= 3u);
1700 rNormal = zero_vector<double>(SPACE_DIM);
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++)
1705 c_vector<double, SPACE_DIM> vnext_minus_v0 = this->GetVectorFromAtoB(pFace->GetNode(0)->rGetLocation(), pFace->GetNode(local_index)->rGetLocation());
1707 v_minus_v0 = vnext_minus_v0;
1709 double magnitude = norm_2(rNormal);
1710 if (magnitude != 0.0)
1713 rNormal /= magnitude;
1717 return magnitude / 2.0;
1725template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1729 if constexpr (SPACE_DIM == 3)
1732 c_vector<double, SPACE_DIM> unit_normal;
1733 return CalculateUnitNormalToFaceWithArea(pFace, unit_normal);
#define EXCEPTION(message)
const unsigned UNSIGNED_UNSET
#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 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
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
std::set< unsigned > & rGetContainingElementIndices()
void AddElement(unsigned index)
const c_vector< double, SPACE_DIM > & rGetLocation() const
bool IsBoundaryNode() const
unsigned GetIndex() const
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)
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)
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