466 EXCEPTION(
"RefineElement could not be started (point is not in element)");
475 for (
unsigned i = 0; i < ELEMENT_DIM; i++)
479 unsigned new_elt_index;
480 if (mDeletedElementIndices.empty())
482 new_elt_index = this->mElements.size();
486 new_elt_index = mDeletedElementIndices.back();
487 mDeletedElementIndices.pop_back();
494 p_new_element->
UpdateNode(ELEMENT_DIM-1-i, this->mNodes[new_node_index]);
497 if ((
unsigned) new_elt_index == this->mElements.size())
499 this->mElements.push_back(p_new_element);
503 delete this->mElements[new_elt_index];
504 this->mElements[new_elt_index] = p_new_element;
509 pElement->
UpdateNode(ELEMENT_DIM, this->mNodes[new_node_index]);
511 return new_node_index;
517 if (!this->mNodes[index]->IsBoundaryNode() )
519 EXCEPTION(
" You may only delete a boundary node ");
522 this->mNodes[index]->MarkAsDeleted();
523 mDeletedNodeIndices.push_back(index);
526 typename std::vector<Node<SPACE_DIM>*>::iterator b_node_iter
527 = std::find(this->mBoundaryNodes.begin(), this->mBoundaryNodes.end(), this->mNodes[index]);
528 this->mBoundaryNodes.erase(b_node_iter);
531 std::set<unsigned> boundary_element_indices = this->mNodes[index]->rGetContainingBoundaryElementIndices();
532 std::set<unsigned>::const_iterator boundary_element_indices_iterator = boundary_element_indices.begin();
533 while (boundary_element_indices_iterator != boundary_element_indices.end())
535 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element = this->GetBoundaryElement(*boundary_element_indices_iterator);
536 p_boundary_element->MarkAsDeleted();
537 mDeletedBoundaryElementIndices.push_back(*boundary_element_indices_iterator);
538 boundary_element_indices_iterator++;
542 std::set<unsigned> element_indices = this->mNodes[index]->rGetContainingElementIndices();
543 std::set<unsigned>::const_iterator element_indices_iterator = element_indices.begin();
544 while (element_indices_iterator != element_indices.end())
547 for (
unsigned i=0; i<p_element->
GetNumNodes(); i++)
554 this->mBoundaryNodes.push_back(p_node);
558 mDeletedElementIndices.push_back(p_element->
GetIndex());
559 element_indices_iterator++;
566 assert(!mAddedNodes);
567 map.
Resize(this->GetNumAllNodes());
569 std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> live_elements;
571 for (
unsigned i=0; i<this->mElements.size(); i++)
573 assert(i==this->mElements[i]->GetIndex());
574 if (this->mElements[i]->IsDeleted())
576 delete this->mElements[i];
580 live_elements.push_back(this->mElements[i]);
582 unsigned this_element_index = this->mElements[i]->GetIndex();
583 if (SPACE_DIM == ELEMENT_DIM)
585 this->mElementJacobians[live_elements.size()-1] = this->mElementJacobians[this_element_index];
586 this->mElementInverseJacobians[live_elements.size()-1] = this->mElementInverseJacobians[this_element_index];
590 this->mElementWeightedDirections[live_elements.size()-1] = this->mElementWeightedDirections[this_element_index];
592 this->mElementJacobianDeterminants[live_elements.size()-1] = this->mElementJacobianDeterminants[this_element_index];
596 assert(mDeletedElementIndices.size() == this->mElements.size()-live_elements.size());
597 mDeletedElementIndices.clear();
598 this->mElements = live_elements;
599 unsigned num_elements = this->mElements.size();
601 if (SPACE_DIM == ELEMENT_DIM)
603 this->mElementJacobians.resize(num_elements);
604 this->mElementInverseJacobians.resize(num_elements);
608 this->mElementWeightedDirections.resize(num_elements);
610 this->mElementJacobianDeterminants.resize(num_elements);
612 std::vector<Node<SPACE_DIM> *> live_nodes;
613 for (
unsigned i=0; i<this->mNodes.size(); i++)
615 if (this->mNodes[i]->IsDeleted())
617 delete this->mNodes[i];
622 live_nodes.push_back(this->mNodes[i]);
625 map.
SetNewIndex(i, (
unsigned)(live_nodes.size()-1));
629 assert(mDeletedNodeIndices.size() == this->mNodes.size()-live_nodes.size());
630 this->mNodes = live_nodes;
631 mDeletedNodeIndices.clear();
633 std::vector<
BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> live_boundary_elements;
634 for (
unsigned i=0; i<this->mBoundaryElements.size(); i++)
636 if (this->mBoundaryElements[i]->IsDeleted())
638 delete this->mBoundaryElements[i];
642 live_boundary_elements.push_back(this->mBoundaryElements[i]);
644 this->mBoundaryElementWeightedDirections[live_boundary_elements.size()-1] = this->mBoundaryElementWeightedDirections[this->mBoundaryElements[i]->GetIndex()];
645 this->mBoundaryElementJacobianDeterminants[live_boundary_elements.size()-1] = this->mBoundaryElementJacobianDeterminants[this->mBoundaryElements[i]->GetIndex()];
649 assert(mDeletedBoundaryElementIndices.size() == this->mBoundaryElements.size()-live_boundary_elements.size());
650 this->mBoundaryElements = live_boundary_elements;
651 mDeletedBoundaryElementIndices.clear();
653 unsigned num_boundary_elements = this->mBoundaryElements.size();
655 this->mBoundaryElementWeightedDirections.resize(num_boundary_elements);
656 this->mBoundaryElementJacobianDeterminants.resize(num_boundary_elements);
658 for (
unsigned i=0; i<this->mNodes.size(); i++)
660 this->mNodes[i]->SetIndex(i);
663 for (
unsigned i=0; i<this->mElements.size(); i++)
665 this->mElements[i]->ResetIndex(i);
668 for (
unsigned i=0; i<this->mBoundaryElements.size(); i++)
670 this->mBoundaryElements[i]->ResetIndex(i);
678 assert( ELEMENT_DIM == SPACE_DIM );
682 if (GetNumNodes() <= SPACE_DIM)
684 EXCEPTION(
"The number of nodes must exceed the spatial dimension.");
688 map.
Resize(this->GetNumAllNodes());
689 if (mAddedNodes || !mDeletedNodeIndices.empty())
692 if (this->mpDistributedVectorFactory)
694 delete this->mpDistributedVectorFactory;
701 std::vector<c_vector<double, SPACE_DIM> > old_node_locations;
702 unsigned new_index = 0;
703 for (
unsigned i=0; i<this->GetNumAllNodes(); i++)
705 if (this->mNodes[i]->IsDeleted())
712 old_node_locations.push_back(this->mNodes[i]->rGetLocation());
721 for (
unsigned node_index=0; node_index<old_node_locations.size(); node_index++)
724 bool is_boundary_node = (node_index==0 || node_index==old_node_locations.size()-1);
727 this->mNodes.push_back(p_node);
729 if (is_boundary_node)
731 this->mBoundaryNodes.push_back(p_node);
736 std::map<double, unsigned> location_index_map;
737 for (
unsigned i=0; i<this->mNodes.size(); i++)
739 location_index_map[this->mNodes[i]->rGetLocation()[0]] = this->mNodes[i]->GetIndex();
743 std::vector<unsigned> node_indices_ordered_spatially;
744 for (std::map<double, unsigned>::iterator iter = location_index_map.begin();
745 iter != location_index_map.end();
748 node_indices_ordered_spatially.push_back(iter->second);
752 this->mElements.reserve(old_node_locations.size()-1);
753 for (
unsigned element_index=0; element_index<old_node_locations.size()-1; element_index++)
755 std::vector<Node<SPACE_DIM>*> nodes;
756 for (
unsigned j=0; j<2; j++)
758 unsigned global_node_index = node_indices_ordered_spatially[element_index + j];
759 assert(global_node_index < this->mNodes.size());
760 nodes.push_back(this->mNodes[global_node_index]);
766 std::vector<Node<SPACE_DIM>*> nodes;
767 nodes.push_back(this->mNodes[0]);
771 nodes.push_back(this->mNodes[old_node_locations.size()-1]);
774 this->RefreshJacobianCachedData();
776 else if (SPACE_DIM==2)
778 struct triangulateio mesher_input, mesher_output;
779 this->InitialiseTriangulateIo(mesher_input);
780 this->InitialiseTriangulateIo(mesher_output);
782 this->ExportToMesher(map, mesher_input);
785 triangulate((
char*)
"Qze", &mesher_input, &mesher_output,
nullptr);
787 this->ImportFromMesher(mesher_output, mesher_output.numberoftriangles, mesher_output.trianglelist, mesher_output.numberofedges, mesher_output.edgelist, mesher_output.edgemarkerlist);
790 this->FreeTriangulateIo(mesher_input);
791 this->FreeTriangulateIo(mesher_output);
796 class tetgen::tetgenio mesher_input, mesher_output;
798 this->ExportToMesher(map, mesher_input);
801 tetgen::tetrahedralize((
char*)
"Qz", &mesher_input, &mesher_output);
803 this->ImportFromMesher(mesher_output, mesher_output.numberoftetrahedra, mesher_output.tetrahedronlist, mesher_output.numberoftrifaces, mesher_output.trifacelist,
nullptr);
817 assert(ELEMENT_DIM == 2);
818 assert(SPACE_DIM == 3);
820 std::vector<c_vector<unsigned, 5> > history;
822 bool long_edge_exists =
true;
824 while (long_edge_exists)
826 std::set<std::pair<unsigned, unsigned> > long_edges;
830 elem_iter != this->GetElementIteratorEnd();
833 unsigned num_nodes = ELEMENT_DIM+1;
836 for (
unsigned local_index=0; local_index<num_nodes; local_index++)
840 unsigned local_index_plus_one = (local_index+1)%num_nodes;
844 double distance_between_nodes = this->GetDistanceBetweenNodes(p_node_a->
GetIndex(), p_node_b->
GetIndex());
846 if (distance_between_nodes > cutoffLength)
850 std::pair<unsigned, unsigned> long_edge(p_node_a->
GetIndex(),p_node_b->
GetIndex());
851 long_edges.insert(long_edge);
855 std::pair<unsigned, unsigned> long_edge(p_node_b->
GetIndex(),p_node_a->
GetIndex());
856 long_edges.insert(long_edge);
862 if (long_edges.size() > 0)
864 while (long_edges.size() > 0)
866 double longest_edge = 0.0;
867 std::set<std::pair<unsigned, unsigned> >::iterator longest_edge_iter;
870 for (std::set<std::pair<unsigned, unsigned> >::iterator edge_iter = long_edges.begin();
871 edge_iter != long_edges.end();
874 unsigned node_a_global_index = edge_iter->first;
875 unsigned node_b_global_index = edge_iter->second;
877 double distance_between_nodes = this->GetDistanceBetweenNodes(node_a_global_index, node_b_global_index);
879 if (distance_between_nodes > longest_edge)
881 longest_edge = distance_between_nodes;
882 longest_edge_iter = edge_iter;
885 assert(longest_edge >0);
887 c_vector<unsigned, 3> new_node_index = SplitEdge(this->GetNode(longest_edge_iter->first), this->GetNode(longest_edge_iter->second));
889 c_vector<unsigned, 5> node_set;
890 node_set(0) = new_node_index[0];
891 node_set(1) = longest_edge_iter->first;
892 node_set(2) = longest_edge_iter->second;
893 node_set(3) = new_node_index[1];
894 node_set(4) = new_node_index[2];
895 history.push_back(node_set);
898 long_edges.erase(*longest_edge_iter);
903 long_edge_exists =
false;
913 c_vector<unsigned, 3> new_node_index_vector;
918 std::set<unsigned> intersection_elements;
919 std::set_intersection(elements_of_node_a.begin(), elements_of_node_a.end(),
920 elements_of_node_b.begin(), elements_of_node_b.end(),
921 std::inserter(intersection_elements, intersection_elements.begin()));
926 bool is_boundary_node = intersection_elements.size() == 1;
930 unsigned new_node_index = this->AddNode(p_new_node);
932 new_node_index_vector[0] = new_node_index;
934 unsigned counter = 1;
936 for (std::set<unsigned>::const_iterator it = intersection_elements.begin(); it != intersection_elements.end(); ++it)
938 unsigned elementIndex = *it;
946 AddElement(p_new_element);
949 p_new_element->
ReplaceNode(pNodeA, this->mNodes[new_node_index]);
952 p_original_element->
ReplaceNode(pNodeB, this->mNodes[new_node_index]);
976 new_node_index_vector[counter] = other_node_index;
988 return new_node_index_vector;
994 assert(ELEMENT_DIM == SPACE_DIM);
996 std::set<unsigned> neighbouring_elements_indices;
997 std::set< Element<ELEMENT_DIM,SPACE_DIM> *> neighbouring_elements;
998 std::set<unsigned> neighbouring_nodes_indices;
1001 for (
unsigned i=0; i<num_nodes; i++)
1005 for (std::set<unsigned>::const_iterator it = neighbouring_elements_indices.begin();
1006 it != neighbouring_elements_indices.end();
1009 neighbouring_elements.insert(this->GetElement(*it));
1012 neighbouring_elements.erase(pElement);
1015 typedef typename std::set<Element<ELEMENT_DIM,SPACE_DIM> *>::const_iterator
ElementIterator;
1018 it != neighbouring_elements.end();
1021 for (
unsigned i=0; i<num_nodes; i++)
1023 neighbouring_nodes_indices.insert((*it)->GetNodeGlobalIndex(i));
1028 for (
unsigned i = 0; i < num_nodes; i++)
1034 c_vector<double, SPACE_DIM+1> this_circum_centre = zero_vector<double>(SPACE_DIM+1);
1039 c_vector<double, ELEMENT_DIM> circum_centre = zero_vector<double>(ELEMENT_DIM);
1040 for (
unsigned i=0; i<ELEMENT_DIM; i++)
1042 circum_centre[i] = this_circum_centre[i];
1045 for (std::set<unsigned>::const_iterator it = neighbouring_nodes_indices.begin();
1046 it != neighbouring_nodes_indices.end();
1049 c_vector<double, ELEMENT_DIM> node_location = this->GetNode(*it)->rGetLocation();
1052 node_location -= circum_centre;
1055 double squared_distance = inner_prod(node_location, node_location);
1059 if (squared_distance < this_circum_centre[ELEMENT_DIM])
1062 double radius = sqrt(this_circum_centre[ELEMENT_DIM]);
1063 double distance = radius - sqrt(squared_distance);
1066 if (distance/radius > maxPenetration)