36 #include "VertexMesh.hpp"
37 #include "RandomNumberGenerator.hpp"
40 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
43 : mpDelaunayMesh(NULL)
50 for (
unsigned node_index=0; node_index<nodes.size(); node_index++)
53 this->
mNodes.push_back(p_temp_node);
55 for (
unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
58 mElements.push_back(p_temp_vertex_element);
65 std::set<unsigned> faces_counted;
68 for (
unsigned elem_index=0; elem_index<
mElements.size(); elem_index++)
71 for (
unsigned face_index=0; face_index<
mElements[elem_index]->GetNumFaces(); face_index++)
74 unsigned global_index = p_face->GetIndex();
77 if (faces_counted.find(global_index) == faces_counted.end())
80 faces_counted.insert(global_index);
87 for (
unsigned index=0; index<
mElements.size(); index++)
91 unsigned element_index = p_element->
GetIndex();
92 unsigned num_nodes_in_element = p_element->
GetNumNodes();
94 for (
unsigned node_index=0; node_index<num_nodes_in_element; node_index++)
96 p_element->
GetNode(node_index)->AddElement(element_index);
104 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
108 : mpDelaunayMesh(NULL)
114 for (
unsigned node_index=0; node_index<nodes.size(); node_index++)
117 this->
mNodes.push_back(p_temp_node);
120 for (
unsigned face_index=0; face_index<faces.size(); face_index++)
122 VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_temp_face = faces[face_index];
123 mFaces.push_back(p_temp_face);
126 for (
unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
129 mElements.push_back(p_temp_vertex_element);
133 for (
unsigned index=0; index<
mElements.size(); index++)
136 for (
unsigned node_index=0; node_index<p_temp_vertex_element->
GetNumNodes(); node_index++)
155 : mpDelaunayMesh(&rMesh)
166 this->
mNodes.reserve(num_nodes);
170 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
180 for (
unsigned i=0; i<num_nodes; i++)
184 for (
unsigned local_index=0; local_index<3; local_index++)
186 unsigned elem_index =
mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
187 unsigned num_nodes_in_elem =
mElements[elem_index]->GetNumNodes();
188 unsigned end_index = num_nodes_in_elem>0 ? num_nodes_in_elem-1 : 0;
195 for (
unsigned elem_index=0; elem_index<
mElements.size(); elem_index++)
202 std::vector<std::pair<double, unsigned> > index_angle_list;
203 for (
unsigned local_index=0; local_index<
mElements[elem_index]->GetNumNodes(); local_index++)
205 c_vector<double, 2> vectorA =
mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
206 c_vector<double, 2> vectorB =
mElements[elem_index]->GetNodeLocation(local_index);
207 c_vector<double, 2> centre_to_vertex =
mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
209 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
210 unsigned global_index =
mElements[elem_index]->GetNodeGlobalIndex(local_index);
212 std::pair<double, unsigned> pair(angle, global_index);
213 index_angle_list.push_back(pair);
217 sort(index_angle_list.begin(), index_angle_list.end());
221 for (
unsigned count = 0; count < index_angle_list.size(); count++)
223 unsigned local_index = count>1 ? count-1 : 0;
224 p_new_element->
AddNode(
mNodes[index_angle_list[count].second], local_index);
245 : mpDelaunayMesh(&rMesh)
253 this->
mNodes.reserve(num_nodes);
258 std::map<unsigned, VertexElement<3,3>*> index_element_map;
259 unsigned face_index = 0;
260 unsigned element_index = 0;
267 Node<3>* p_node_a = edge_iterator.GetNodeA();
268 Node<3>* p_node_b = edge_iterator.GetNodeB();
274 std::set<unsigned> edge_element_indices;
276 std::set_intersection(node_a_element_indices.begin(),
277 node_a_element_indices.end(),
278 node_b_element_indices.begin(),
279 node_b_element_indices.end(),
280 std::inserter(edge_element_indices, edge_element_indices.begin()));
283 c_vector<double,3> mid_edge = edge_vector*0.5 + p_node_a->
rGetLocation();
285 unsigned element0_index = *(edge_element_indices.begin());
287 c_vector<double,3> basis_vector1 =
mNodes[element0_index]->rGetLocation() - mid_edge;
289 c_vector<double,3> basis_vector2 =
VectorProduct(edge_vector, basis_vector1);
296 std::vector<std::pair<double, unsigned> > index_angle_list;
299 for (std::set<unsigned>::iterator index_iter = edge_element_indices.begin();
300 index_iter != edge_element_indices.end();
304 c_vector<double, 3> vertex_vector =
mNodes[*index_iter]->rGetLocation() - mid_edge;
306 double local_vertex_dot_basis_vector1 = inner_prod(vertex_vector, basis_vector1);
307 double local_vertex_dot_basis_vector2 = inner_prod(vertex_vector, basis_vector2);
309 double angle = atan2(local_vertex_dot_basis_vector2, local_vertex_dot_basis_vector1);
311 std::pair<double, unsigned> pair(angle, *index_iter);
312 index_angle_list.push_back(pair);
316 sort(index_angle_list.begin(), index_angle_list.end());
321 for (
unsigned count = 0; count < index_angle_list.size(); count++)
323 unsigned local_index = count>1 ? count-1 : 0;
324 p_face->AddNode(
mNodes[index_angle_list[count].second], local_index);
333 unsigned node_a_index = p_node_a->
GetIndex();
334 if (index_element_map[node_a_index])
337 index_element_map[node_a_index]->AddFace(p_face);
346 index_element_map[node_a_index] = p_element;
351 unsigned node_b_index = p_node_b->
GetIndex();
352 if (index_element_map[node_b_index])
355 index_element_map[node_b_index]->AddFace(p_face);
364 index_element_map[node_b_index] = p_element;
371 unsigned elem_count = 0;
372 for (std::map<
unsigned,
VertexElement<3,3>*>::iterator element_iter = index_element_map.begin();
373 element_iter != index_element_map.end();
376 mElements.push_back(element_iter->second);
385 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
388 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
389 c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
397 c_vector<double, SPACE_DIM+1> circumsphere = rMesh.
GetElement(i)->CalculateCircumsphere(jacobian, inverse_jacobian);
399 c_vector<double, SPACE_DIM> circumcentre;
400 for (
unsigned j=0; j<SPACE_DIM; j++)
402 circumcentre(j) = circumsphere(j);
411 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
414 assert(SPACE_DIM == 2);
416 std::set<unsigned> node_indices_1;
417 for (
unsigned i=0; i<mElements[elementIndex1]->GetNumNodes(); i++)
419 node_indices_1.insert(mElements[elementIndex1]->GetNodeGlobalIndex(i));
421 std::set<unsigned> node_indices_2;
422 for (
unsigned i=0; i<mElements[elementIndex2]->GetNumNodes(); i++)
424 node_indices_2.insert(mElements[elementIndex2]->GetNodeGlobalIndex(i));
427 std::set<unsigned> shared_nodes;
428 std::set_intersection(node_indices_1.begin(), node_indices_1.end(),
429 node_indices_2.begin(), node_indices_2.end(),
430 std::inserter(shared_nodes, shared_nodes.begin()));
432 if (shared_nodes.size() == 1)
435 EXCEPTION(
"Elements "<< elementIndex1 <<
" and "<< elementIndex2<<
" share only one node.");
437 assert(shared_nodes.size() == 2);
439 unsigned index1 = *(shared_nodes.begin());
440 unsigned index2 = *(++(shared_nodes.begin()));
442 double edge_length = this->GetDistanceBetweenNodes(index1, index2);
446 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
449 assert(SPACE_DIM == 2);
451 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
453 double discriminant = sqrt((moments(0) - moments(1))*(moments(0) - moments(1)) + 4.0*moments(2)*moments(2));
456 double largest_eigenvalue = (moments(0) + moments(1) + discriminant)*0.5;
457 double smallest_eigenvalue = (moments(0) + moments(1) - discriminant)*0.5;
459 double elongation_shape_factor = sqrt(largest_eigenvalue/smallest_eigenvalue);
460 return elongation_shape_factor;
464 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
467 mpDelaunayMesh = NULL;
468 this->mMeshChangesDuringSimulation =
false;
473 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
480 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
483 assert(index < this->mNodes.size());
488 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
491 assert(index < this->mElements.size());
496 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
505 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
510 if (mVoronoiElementIndexMap.empty())
512 node_index = elementIndex;
516 for (std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.begin();
517 iter != mVoronoiElementIndexMap.end();
520 if (iter->second == elementIndex)
522 node_index = iter->first;
532 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
537 if (mVoronoiElementIndexMap.empty())
539 element_index = nodeIndex;
543 std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.find(nodeIndex);
545 if (iter == mVoronoiElementIndexMap.end())
547 EXCEPTION(
"This index does not correspond to a VertexElement");
551 element_index = iter->second;
555 return element_index;
559 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
562 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
568 unsigned rosette_rank = 0;
569 for (
unsigned node_idx = 0 ; node_idx < p_element->
GetNumNodes() ; node_idx++)
571 unsigned num_elems_this_node = p_element->
GetNode(node_idx)->rGetContainingElementIndices().size();
573 if (num_elems_this_node > rosette_rank)
575 rosette_rank = num_elems_this_node;
584 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
588 for (
unsigned i=0; i<mElements.size(); i++)
595 for (
unsigned i=0; i<mFaces.size(); i++)
602 for (
unsigned i=0; i<this->mNodes.size(); i++)
604 delete this->mNodes[i];
606 this->mNodes.clear();
610 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
613 return this->mNodes.size();
617 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
620 return mElements.size();
624 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
627 return mElements.size();
631 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
634 return mFaces.size();
638 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
641 assert(index < mElements.size());
642 return mElements[index];
646 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
649 assert(index < mFaces.size());
650 return mFaces[index];
654 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
660 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
671 double centroid_x = 0;
672 double centroid_y = 0;
675 double element_signed_area = 0.0;
678 c_vector<double, SPACE_DIM> first_node_location = p_element->
GetNodeLocation(0);
679 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
682 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
684 c_vector<double, SPACE_DIM> next_node_location = p_element->
GetNodeLocation((local_index+1)%num_nodes);
685 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
687 double this_x = pos_1[0];
688 double this_y = pos_1[1];
689 double next_x = pos_2[0];
690 double next_y = pos_2[1];
692 double signed_area_term = this_x*next_y - this_y*next_x;
694 centroid_x += (this_x + next_x)*signed_area_term;
695 centroid_y += (this_y + next_y)*signed_area_term;
696 element_signed_area += 0.5*signed_area_term;
701 assert(element_signed_area != 0.0);
704 centroid = first_node_location;
705 centroid(0) += centroid_x / (6.0*element_signed_area);
706 centroid(1) += centroid_y / (6.0*element_signed_area);
712 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
716 centroid /= ((
double) num_nodes);
726 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
730 std::set<unsigned> neighbouring_node_indices;
733 std::set<unsigned> containing_elem_indices = this->GetNode(nodeIndex)->rGetContainingElementIndices();
736 for (std::set<unsigned>::iterator elem_iter = containing_elem_indices.begin();
737 elem_iter != containing_elem_indices.end();
741 unsigned local_index = GetElement(*elem_iter)->GetNodeLocalIndex(nodeIndex);
744 unsigned num_nodes = GetElement(*elem_iter)->GetNumNodes();
745 unsigned previous_local_index = (local_index + num_nodes - 1)%num_nodes;
746 unsigned next_local_index = (local_index + 1)%num_nodes;
749 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(previous_local_index));
750 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(next_local_index));
753 return neighbouring_node_indices;
757 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
764 std::set<unsigned> node_indices_in_this_element;
765 for (
unsigned local_index=0; local_index<p_element->
GetNumNodes(); local_index++)
768 node_indices_in_this_element.insert(global_index);
772 if (node_indices_in_this_element.find(nodeIndex) == node_indices_in_this_element.end())
774 EXCEPTION(
"The given node is not contained in the given element.");
778 std::set<unsigned> neighbouring_node_indices_not_in_this_element;
781 std::set<unsigned> node_neighbours = GetNeighbouringNodeIndices(nodeIndex);
784 for (std::set<unsigned>::iterator iter = node_neighbours.begin();
785 iter != node_neighbours.end();
788 if (node_indices_in_this_element.find(*iter) == node_indices_in_this_element.end())
790 neighbouring_node_indices_not_in_this_element.insert(*iter);
794 return neighbouring_node_indices_not_in_this_element;
798 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
805 std::set<unsigned> neighbouring_element_indices;
808 for (
unsigned local_index=0; local_index<p_element->
GetNumNodes(); local_index++)
817 std::set<unsigned> all_elements;
818 std::set_union(neighbouring_element_indices.begin(), neighbouring_element_indices.end(),
819 containing_elem_indices.begin(), containing_elem_indices.end(),
820 std::inserter(all_elements, all_elements.begin()));
823 neighbouring_element_indices = all_elements;
827 neighbouring_element_indices.erase(elementIndex);
829 return neighbouring_element_indices;
872 this->mNodes.reserve(num_nodes);
877 std::vector<double> node_data;
878 for (
unsigned i=0; i<num_nodes; i++)
881 unsigned is_boundary_node = (
bool) node_data[2];
882 node_data.pop_back();
883 this->mNodes.push_back(
new Node<2>(i, node_data, is_boundary_node));
892 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
898 std::vector<Node<2>*> nodes;
899 unsigned num_nodes_in_element = element_data.
NodeIndices.size();
900 for (
unsigned j=0; j<num_nodes_in_element; j++)
902 assert(element_data.
NodeIndices[j] < this->mNodes.size());
903 nodes.push_back(this->mNodes[element_data.
NodeIndices[j]]);
908 mElements.push_back(p_element);
931 this->mNodes.reserve(num_nodes);
936 std::vector<double> node_data;
937 for (
unsigned i=0; i<num_nodes; i++)
940 unsigned is_boundary_node = (
bool) node_data[3];
941 node_data.pop_back();
942 this->mNodes.push_back(
new Node<3>(i, node_data, is_boundary_node));
951 std::set<unsigned> faces_counted;
954 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
958 assert(dynamic_cast<VERTEX_MESH_READER*>(&rMeshReader) != NULL);
961 VertexElementData element_data =
static_cast<VERTEX_MESH_READER*
>(&rMeshReader)->GetNextElementDataWithFaces();
964 std::vector<Node<3>*> nodes;
965 unsigned num_nodes_in_element = element_data.
NodeIndices.size();
966 for (
unsigned j=0; j<num_nodes_in_element; j++)
968 assert(element_data.
NodeIndices[j] < this->mNodes.size());
969 nodes.push_back(this->mNodes[element_data.
NodeIndices[j]]);
973 std::vector<VertexElement<2,3>*> faces;
974 unsigned num_faces_in_element = element_data.
Faces.size();
975 for (
unsigned i=0; i<num_faces_in_element; i++)
984 std::vector<
Node<3>*> nodes_in_face;
985 unsigned num_nodes_in_face = face_data.
NodeIndices.size();
986 for (
unsigned j=0; j<num_nodes_in_face; j++)
988 assert(face_data.
NodeIndices[j] < this->mNodes.size());
989 nodes_in_face.push_back(this->mNodes[face_data.
NodeIndices[j]]);
993 if (faces_counted.find(face_index) == faces_counted.end())
997 mFaces.push_back(p_face);
998 faces_counted.insert(face_index);
999 faces.push_back(p_face);
1004 bool face_added =
false;
1005 for (
unsigned k=0; k<mFaces.size(); k++)
1007 if (mFaces[k]->GetIndex() == face_index)
1009 faces.push_back(mFaces[k]);
1015 assert(face_added ==
true);
1020 std::vector<bool> orientations = std::vector<bool>(num_faces_in_element,
true);
1024 mElements.push_back(p_element);
1030 p_element->SetAttribute(attribute_value);
1036 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1038 const c_vector<double, SPACE_DIM>& rLocationA,
const c_vector<double, SPACE_DIM>& rLocationB)
1040 c_vector<double, SPACE_DIM> vector;
1043 vector = mpDelaunayMesh->GetVectorFromAtoB(rLocationA, rLocationB);
1053 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1056 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
1061 double element_volume = 0.0;
1065 c_vector<double, SPACE_DIM> first_node_location = p_element->
GetNodeLocation(0);
1066 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
1069 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
1071 c_vector<double, SPACE_DIM> next_node_location = p_element->
GetNodeLocation((local_index+1)%num_nodes);
1072 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
1074 double this_x = pos_1[0];
1075 double this_y = pos_1[1];
1076 double next_x = pos_2[0];
1077 double next_y = pos_2[1];
1079 element_volume += 0.5*(this_x*next_y - next_x*this_y);
1087 c_vector<double, SPACE_DIM> pyramid_apex = p_element->
GetNodeLocation(0);
1088 for (
unsigned face_index=0; face_index<p_element->
GetNumFaces(); face_index++)
1095 c_vector<double, SPACE_DIM> unit_normal;
1096 double face_area = CalculateUnitNormalToFaceWithArea(p_face, unit_normal);
1099 c_vector<double, SPACE_DIM> base_to_apex = GetVectorFromAtoB(p_face->GetNodeLocation(0), pyramid_apex);
1100 double perpendicular_distance = fabs(inner_prod(base_to_apex, unit_normal));
1104 element_volume += face_area * perpendicular_distance / 3;
1109 return fabs(element_volume);
1113 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1116 assert(SPACE_DIM == 2 || SPACE_DIM == 3);
1121 double surface_area = 0.0;
1126 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
1130 surface_area += this->GetDistanceBetweenNodes(this_node_index, next_node_index);
1131 this_node_index = next_node_index;
1137 for (
unsigned face_index=0; face_index<p_element->
GetNumFaces(); face_index++)
1139 surface_area += CalculateAreaOfFace(p_element->
GetFace(face_index));
1142 return surface_area;
1150 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1153 assert(SPACE_DIM == 2);
1154 assert(ELEMENT_DIM == SPACE_DIM);
1161 bool element_includes_point =
false;
1221 c_vector<double, SPACE_DIM> first_vertex = p_element->
GetNodeLocation(0);
1222 c_vector<double, SPACE_DIM> test_point = GetVectorFromAtoB(first_vertex, rTestPoint);
1225 c_vector<double, SPACE_DIM> vertexA = zero_vector<double>(SPACE_DIM);
1226 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
1229 c_vector<double, SPACE_DIM> vector_a_to_point = GetVectorFromAtoB(vertexA, test_point);
1233 if (norm_2(vector_a_to_point) < DBL_EPSILON)
1238 c_vector<double, SPACE_DIM> vertexB = GetVectorFromAtoB(first_vertex, p_element->
GetNodeLocation((local_index+1)%num_nodes));
1239 c_vector<double, SPACE_DIM> vector_b_to_point = GetVectorFromAtoB(vertexB, test_point);
1240 c_vector<double, SPACE_DIM> vector_a_to_b = GetVectorFromAtoB(vertexA, vertexB);
1243 if ( (fabs(vector_a_to_b[1]) < DBL_EPSILON) &&
1244 (fabs(vector_a_to_point[1]) < DBL_EPSILON) &&
1245 (fabs(vector_b_to_point[1]) < DBL_EPSILON) )
1247 if ( (vector_a_to_point[0]>0) != (vector_b_to_point[0]>0) )
1255 if ( (vertexA[1] > test_point[1]) != (vertexB[1] > test_point[1]) )
1258 if (test_point[0] < vertexA[0] + vector_a_to_b[0]*vector_a_to_point[1]/vector_a_to_b[1])
1260 element_includes_point = !element_includes_point;
1266 return element_includes_point;
1270 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1274 assert(SPACE_DIM == 2);
1275 assert(ELEMENT_DIM == SPACE_DIM);
1281 double min_squared_normal_distance = DBL_MAX;
1282 unsigned min_distance_edge_index = UINT_MAX;
1285 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
1288 c_vector<double, SPACE_DIM> vertexA = p_element->
GetNodeLocation(local_index);
1289 c_vector<double, SPACE_DIM> vertexB = p_element->
GetNodeLocation((local_index+1)%num_nodes);
1291 c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, rTestPoint);
1292 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
1293 double distance_a_to_b = norm_2(vector_a_to_b);
1295 c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
1296 double distance_parallel_to_edge = inner_prod(vector_a_to_point, edge_ab_unit_vector);
1298 double squared_distance_normal_to_edge =
SmallPow(norm_2(vector_a_to_point), 2) -
SmallPow(distance_parallel_to_edge, 2);
1306 if (squared_distance_normal_to_edge < DBL_EPSILON)
1308 squared_distance_normal_to_edge = 0.0;
1311 if (fabs(distance_parallel_to_edge) < DBL_EPSILON)
1313 distance_parallel_to_edge = 0.0;
1315 else if (fabs(distance_parallel_to_edge-distance_a_to_b) < DBL_EPSILON)
1317 distance_parallel_to_edge = distance_a_to_b;
1321 if (squared_distance_normal_to_edge < min_squared_normal_distance &&
1322 distance_parallel_to_edge >=0 &&
1323 distance_parallel_to_edge <= distance_a_to_b)
1325 min_squared_normal_distance = squared_distance_normal_to_edge;
1326 min_distance_edge_index = local_index;
1330 assert(min_distance_edge_index < num_nodes);
1331 return min_distance_edge_index;
1335 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1338 assert(SPACE_DIM == 2);
1343 c_vector<double, 3> moments = zero_vector<double>(3);
1346 c_vector<double, SPACE_DIM> centroid = GetCentroidOfElement(index);
1348 c_vector<double, SPACE_DIM> this_node_location = p_element->
GetNodeLocation(0);
1349 c_vector<double, SPACE_DIM> pos_1 = this->GetVectorFromAtoB(centroid, this_node_location);
1351 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
1353 unsigned next_index = (local_index+1)%num_nodes;
1354 c_vector<double, SPACE_DIM> next_node_location = p_element->
GetNodeLocation(next_index);
1355 c_vector<double, SPACE_DIM> pos_2 = this->GetVectorFromAtoB(centroid, next_node_location);
1357 double signed_area_term = pos_1(0)*pos_2(1) - pos_2(0)*pos_1(1);
1359 moments(0) += (pos_1(1)*pos_1(1) + pos_1(1)*pos_2(1) + pos_2(1)*pos_2(1) ) * signed_area_term;
1362 moments(1) += (pos_1(0)*pos_1(0) + pos_1(0)*pos_2(0) + pos_2(0)*pos_2(0)) * signed_area_term;
1365 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;
1380 if (moments(0) < 0.0)
1382 moments(0) = -moments(0);
1383 moments(1) = -moments(1);
1384 moments(2) = -moments(2);
1390 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1393 assert(SPACE_DIM == 2);
1395 c_vector<double, SPACE_DIM> short_axis = zero_vector<double>(SPACE_DIM);
1398 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
1401 double discriminant = (moments(0) - moments(1))*(moments(0) - moments(1)) + 4.0*moments(2)*moments(2);
1402 if (fabs(discriminant) < 1e-10)
1406 short_axis(1) = sqrt(1.0 - short_axis(0)*short_axis(0));
1411 if (moments(2) == 0.0)
1413 if (moments(0) < moments(1))
1415 short_axis(0) = 0.0;
1416 short_axis(1) = 1.0;
1420 short_axis(0) = 1.0;
1421 short_axis(1) = 0.0;
1427 double lambda = 0.5*(moments(0) + moments(1) + sqrt(discriminant));
1429 short_axis(0) = 1.0;
1430 short_axis(1) = (moments(0) - lambda)/moments(2);
1432 double magnitude = norm_2(short_axis);
1433 short_axis = short_axis / magnitude;
1441 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1444 assert(SPACE_DIM==2);
1446 unsigned num_nodes_in_element = pElement->
GetNumNodes();
1447 unsigned next_local_index = (localIndex+1)%num_nodes_in_element;
1450 unsigned previous_local_index = (num_nodes_in_element+localIndex-1)%num_nodes_in_element;
1452 c_vector<double, SPACE_DIM> previous_node_location = pElement->
GetNodeLocation(previous_local_index);
1453 c_vector<double, SPACE_DIM> next_node_location = pElement->
GetNodeLocation(next_local_index);
1454 c_vector<double, SPACE_DIM> difference_vector = this->GetVectorFromAtoB(previous_node_location, next_node_location);
1456 c_vector<double, SPACE_DIM> area_gradient;
1458 area_gradient[0] = 0.5*difference_vector[1];
1459 area_gradient[1] = -0.5*difference_vector[0];
1461 return area_gradient;
1465 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1468 assert(SPACE_DIM==2);
1470 unsigned num_nodes_in_element = pElement->
GetNumNodes();
1473 unsigned previous_local_index = (num_nodes_in_element+localIndex-1)%num_nodes_in_element;
1478 double previous_edge_length = this->GetDistanceBetweenNodes(this_global_index, previous_global_index);
1479 assert(previous_edge_length > DBL_EPSILON);
1481 c_vector<double, SPACE_DIM> previous_edge_gradient = this->GetVectorFromAtoB(pElement->
GetNodeLocation(previous_local_index), pElement->
GetNodeLocation(localIndex))/previous_edge_length;
1483 return previous_edge_gradient;
1487 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1490 assert(SPACE_DIM==2);
1492 unsigned next_local_index = (localIndex+1)%(pElement->
GetNumNodes());
1497 double next_edge_length = this->GetDistanceBetweenNodes(this_global_index, next_global_index);
1498 assert(next_edge_length > DBL_EPSILON);
1500 c_vector<double, SPACE_DIM> next_edge_gradient = this->GetVectorFromAtoB(pElement->
GetNodeLocation(next_local_index), pElement->
GetNodeLocation(localIndex))/next_edge_length;
1502 return next_edge_gradient;
1506 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1509 assert(SPACE_DIM==2);
1511 c_vector<double, SPACE_DIM> previous_edge_gradient = GetPreviousEdgeGradientOfElementAtNode(pElement, localIndex);
1512 c_vector<double, SPACE_DIM> next_edge_gradient = GetNextEdgeGradientOfElementAtNode(pElement, localIndex);
1514 return previous_edge_gradient + next_edge_gradient;
1523 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1526 assert(SPACE_DIM == 3);
1531 rNormal = zero_vector<double>(SPACE_DIM);
1533 c_vector<double, SPACE_DIM> v_minus_v0 = this->GetVectorFromAtoB(pFace->
GetNode(0)->rGetLocation(), pFace->
GetNode(1)->rGetLocation());
1534 for (
unsigned local_index=2; local_index<pFace->
GetNumNodes(); local_index++)
1537 c_vector<double, SPACE_DIM> vnext_minus_v0 = this->GetVectorFromAtoB(pFace->
GetNode(0)->rGetLocation(), pFace->
GetNode(local_index)->rGetLocation());
1539 v_minus_v0 = vnext_minus_v0;
1541 double magnitude = norm_2(rNormal);
1542 if ( magnitude != 0.0 )
1545 rNormal /= magnitude;
1549 return magnitude/2.0;
1553 template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
1556 assert(SPACE_DIM == 3);
1559 c_vector<double, SPACE_DIM> unit_normal;
1560 return CalculateUnitNormalToFaceWithArea(pFace, unit_normal);
1563 #if defined(__xlC__)
virtual c_vector< double, 3 > CalculateMomentsOfElement(unsigned index)
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
void SetAttribute(double attribute)
VertexElement< ELEMENT_DIM-1, SPACE_DIM > * GetFace(unsigned index) const
unsigned GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(unsigned elementIndex)
virtual unsigned GetNumFaces() const
virtual c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocationA, const c_vector< double, SPACE_DIM > &rLocationB)
virtual ElementData GetNextElementData()=0
std::set< unsigned > GetNeighbouringNodeIndices(unsigned nodeIndex)
virtual unsigned GetNumElements() const
double SmallPow(double x, unsigned exponent)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetNumAllElements() const
bool IsBoundaryNode() const
virtual unsigned GetNumElements() const
virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian) const
#define EXCEPTION(message)
virtual unsigned GetNumNodes() const
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
void GenerateVerticesFromElementCircumcentres(TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
std::set< unsigned > & rGetContainingElementIndices()
c_vector< double, SPACE_DIM > GetNextEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
std::set< unsigned > GetNeighbouringNodeNotAlsoInElement(unsigned nodeIndex, unsigned elemIndex)
VertexElement< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
c_vector< double, SPACE_DIM > GetPreviousEdgeGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
unsigned GetNumFaces() const
std::vector< ElementData > Faces
std::vector< unsigned > NodeIndices
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
virtual bool HasNodePermutation()
VertexElement< ELEMENT_DIM-1, SPACE_DIM > * GetFace(unsigned index) const
bool mMeshChangesDuringSimulation
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
std::vector< Node< SPACE_DIM > * > mNodes
double CalculateUnitNormalToFaceWithArea(VertexElement< ELEMENT_DIM-1, SPACE_DIM > *pFace, c_vector< double, SPACE_DIM > &rNormal)
virtual unsigned GetNumElements() const =0
virtual double CalculateAreaOfFace(VertexElement< ELEMENT_DIM-1, SPACE_DIM > *pFace)
unsigned GetNumNodes() const
const unsigned UNSIGNED_UNSET
double GetElongationShapeFactorOfElement(unsigned elementIndex)
void AddFace(VertexElement< ELEMENT_DIM-1, SPACE_DIM > *pFace)
static RandomNumberGenerator * Instance()
virtual double GetVolumeOfElement(unsigned index)
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
virtual unsigned GetNumElementAttributes() const
std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > mElements
unsigned SolveNodeMapping(unsigned index) const
c_vector< double, SPACE_DIM > GetPerimeterGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
const c_vector< double, SPACE_DIM > & rGetLocation() const
void AddElement(unsigned index)
unsigned GetIndex() const
unsigned SolveBoundaryElementMapping(unsigned index) const
virtual c_vector< double, SPACE_DIM > GetCentroidOfElement(unsigned index)
virtual c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocationA, const c_vector< double, SPACE_DIM > &rLocationB)
std::set< unsigned > GetNeighbouringElementIndices(unsigned elementIndex)
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
virtual std::vector< double > GetNextNode()=0
c_vector< double, SPACE_DIM > GetAreaGradientOfElementAtNode(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned localIndex)
unsigned GetIndex() const
std::vector< VertexElement< ELEMENT_DIM-1, SPACE_DIM > * > mFaces
double GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)
unsigned GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(unsigned nodeIndex)
virtual double GetSurfaceAreaOfElement(unsigned index)
std::vector< unsigned > NodeIndices
unsigned SolveElementMapping(unsigned index) const
virtual unsigned GetNumNodes() const =0
bool ElementIncludesPoint(const c_vector< double, SPACE_DIM > &rTestPoint, unsigned elementIndex)
unsigned GetLocalIndexForElementEdgeClosestToPoint(const c_vector< double, SPACE_DIM > &rTestPoint, unsigned elementIndex)
c_vector< double, SPACE_DIM > GetShortAxisOfElement(unsigned index)
TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpDelaunayMesh
std::map< unsigned, unsigned > mVoronoiElementIndexMap
unsigned GetRosetteRankOfElement(unsigned index)