55 : mWidth(rMesh.GetWidth(0)),
56 mHeight(rMesh.GetWidth(1)),
70 this->
mNodes.reserve(num_nodes);
75 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
85 for (
unsigned i=0; i<num_nodes; i++)
91 for (
unsigned i=0; i<num_nodes; i++)
95 for (
unsigned local_index=0; local_index<3; local_index++)
97 unsigned elem_index =
mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
98 unsigned num_nodes_in_elem =
mElements[elem_index]->GetNumNodes();
99 unsigned end_index = num_nodes_in_elem>0 ? num_nodes_in_elem-1 : 0;
108 std::vector<Node<2> *> nodes;
113 nodes.push_back(
new Node<2>(node_iter->GetIndex(), node_iter->rGetLocation(),node_iter->IsBoundaryNode()));
122 bool bad_element =
false;
123 double edge_threshold = 1.5;
125 for (
unsigned j=0; j<3; j++)
138 for (
unsigned j=0; j<3; j++)
144 double edge_length = norm_2(edge);
147 if (edge_length<edge_threshold)
150 c_vector<double,2> normal_vector;
153 normal_vector[0]= edge[1];
154 normal_vector[1]= -edge[0];
156 double dij = norm_2(normal_vector);
158 normal_vector /= dij;
160 double bound_offset = 0.5;
161 c_vector<double,2> new_node_location = -bound_offset*normal_vector + p_node_a->
rGetLocation() + 0.5*edge;
163 nodes.push_back(
new Node<2>(new_node_index, new_node_location));
167 unsigned num_sections = 1;
168 for (
unsigned section=0; section <= num_sections; section++)
170 double ratio = (
double)section/(
double)num_sections;
171 c_vector<double,2> new_node_location = -normal_vector + ratio*p_node_a->
rGetLocation() + (1-ratio)*p_node_b->
rGetLocation();
174 bool node_clear =
true;
175 double node_clearance = 0.3;
177 for (
unsigned i=0; i<nodes.size(); i++)
179 double distance = norm_2(
mpDelaunayMesh->GetVectorFromAtoB(nodes[i]->rGetLocation(), new_node_location));
180 if (distance < node_clearance)
189 nodes.push_back(
new Node<2>(new_node_index, new_node_location));
199 for (
unsigned i=0; i<nodes.size(); i++)
210 this->
mNodes.reserve(num_nodes);
214 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
224 for (
unsigned i=0; i<num_nodes; i++)
231 for (
unsigned i = 0; i < num_nodes; i++)
235 for (
unsigned local_index = 0; local_index < 3; local_index++)
237 unsigned elem_index = extended_mesh.
GetElement(i)->GetNodeGlobalIndex(local_index);
239 if (elem_index < num_elements)
241 unsigned num_nodes_in_elem =
mElements[elem_index]->GetNumNodes();
242 unsigned end_index = num_nodes_in_elem > 0 ? num_nodes_in_elem - 1 : 0;
251 for (
unsigned elem_index=0; elem_index<
mElements.size(); elem_index++)
258 std::vector<std::pair<double, unsigned> > index_angle_list;
259 for (
unsigned local_index=0; local_index<
mElements[elem_index]->GetNumNodes(); local_index++)
261 c_vector<double, 2> vectorA =
mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
262 c_vector<double, 2> vectorB =
mElements[elem_index]->GetNodeLocation(local_index);
263 c_vector<double, 2> centre_to_vertex =
mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
265 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
266 unsigned global_index =
mElements[elem_index]->GetNodeGlobalIndex(local_index);
268 std::pair<double, unsigned> pair(angle, global_index);
269 index_angle_list.push_back(pair);
273 sort(index_angle_list.begin(), index_angle_list.end());
277 for (
unsigned count = 0; count < index_angle_list.size(); count++)
279 unsigned local_index = count>1 ? count-1 : 0;
280 p_new_element->
AddNode(
mNodes[index_angle_list[count].second], local_index);
427 std::vector<Node<2>*> temp_nodes;
428 std::vector<VertexElement<2, 2>*> elements;
432 temp_nodes.resize(4 * num_nodes);
435 for (
unsigned index = 0; index < num_nodes; index++)
437 c_vector<double, 2> location;
438 location =
GetNode(index)->rGetLocation();
441 Node<2>* p_node =
new Node<2>(index,
false, location[0], location[1]);
442 temp_nodes[index] = p_node;
445 p_node =
new Node<2>(num_nodes + index,
false, location[0] +
mWidth, location[1]);
446 temp_nodes[num_nodes + index] = p_node;
449 p_node =
new Node<2>(2*num_nodes + index,
false, location[0], location[1] +
mHeight);
450 temp_nodes[2*num_nodes + index] = p_node;
454 temp_nodes[3*num_nodes + index] = p_node;
462 unsigned elem_index = elem_iter->
GetIndex();
463 unsigned num_nodes_in_elem = elem_iter->GetNumNodes();
465 std::vector<Node<2>*> elem_nodes;
468 bool element_straddles_left_right_boundary =
false;
469 bool element_straddles_top_bottom_boundary =
false;
471 const c_vector<double, 2>& r_this_node_location = elem_iter->GetNode(0)->rGetLocation();
472 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
474 const c_vector<double, 2>& r_next_node_location = elem_iter->GetNode((local_index+1)%num_nodes_in_elem)->rGetLocation();
475 c_vector<double, 2> vector;
476 vector = r_next_node_location - r_this_node_location;
478 if (fabs(vector[0]) > 0.5*
mWidth)
480 element_straddles_left_right_boundary =
true;
482 if (fabs(vector[1]) > 0.5*
mHeight)
484 element_straddles_top_bottom_boundary =
true;
489 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
491 unsigned this_node_index = elem_iter->GetNodeGlobalIndex(local_index);
494 if (element_straddles_left_right_boundary)
497 bool node_is_right_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[0] - 0.5*
mWidth > 0);
498 if (!node_is_right_of_centre)
501 this_node_index += num_nodes;
506 if (element_straddles_top_bottom_boundary)
509 bool node_is_above_centre = (elem_iter->GetNode(local_index)->rGetLocation()[1] - 0.5*
mHeight > 0);
510 if (!node_is_above_centre)
513 this_node_index += 2*num_nodes;
517 elem_nodes.push_back(temp_nodes[this_node_index]);
521 elements.push_back(p_element);
526 temp_nodes.resize(9*num_nodes);
529 for (
unsigned index=0; index<num_nodes; index++)
531 c_vector<double, 2> location;
532 location =
GetNode(index)->rGetLocation();
535 Node<2>* p_node =
new Node<2>(index,
false, location[0], location[1]);
536 temp_nodes[index] = p_node;
539 p_node =
new Node<2>(num_nodes + index,
false, location[0] +
mWidth, location[1]);
540 temp_nodes[num_nodes + index] = p_node;
544 temp_nodes[2*num_nodes + index] = p_node;
547 p_node =
new Node<2>(3*num_nodes + index,
false, location[0], location[1] +
mHeight);
548 temp_nodes[3*num_nodes + index] = p_node;
552 temp_nodes[4*num_nodes + index] = p_node;
555 p_node =
new Node<2>(5*num_nodes + index,
false, location[0] -
mWidth, location[1]);
556 temp_nodes[5*num_nodes + index] = p_node;
560 temp_nodes[6*num_nodes + index] = p_node;
563 p_node =
new Node<2>(7*num_nodes + index,
false, location[0], location[1] -
mHeight);
564 temp_nodes[7*num_nodes + index] = p_node;
568 temp_nodes[8*num_nodes + index] = p_node;
576 unsigned elem_index = elem_iter->
GetIndex();
577 unsigned num_nodes_in_elem = elem_iter->GetNumNodes();
579 std::vector<Node<2>*> elem_nodes;
582 bool element_straddles_left_right_boundary =
false;
583 bool element_straddles_top_bottom_boundary =
false;
585 const c_vector<double, 2>& r_this_node_location = elem_iter->GetNode(0)->rGetLocation();
586 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
588 const c_vector<double, 2>& r_next_node_location = elem_iter->GetNode((local_index+1)%num_nodes_in_elem)->rGetLocation();
589 c_vector<double, 2> vector;
590 vector = r_next_node_location - r_this_node_location;
592 if (fabs(vector[0]) > 0.5*
mWidth)
594 element_straddles_left_right_boundary =
true;
596 if (fabs(vector[1]) > 0.5*
mHeight)
598 element_straddles_top_bottom_boundary =
true;
604 bool element_centre_on_right =
true;
605 bool element_centre_on_top =
true;
608 double element_centre_x_location = this->
mpDelaunayMesh->GetNode(delaunay_index)->rGetLocation()[0];
609 double element_centre_y_location = this->
mpDelaunayMesh->GetNode(delaunay_index)->rGetLocation()[1];
611 if (element_centre_x_location < 0.5*
mWidth)
613 element_centre_on_right =
false;
615 if (element_centre_y_location < 0.5*
mHeight)
617 element_centre_on_top =
false;
621 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
623 unsigned this_node_index = elem_iter->GetNodeGlobalIndex(local_index);
626 if (element_straddles_left_right_boundary && !element_straddles_top_bottom_boundary)
629 bool node_is_right_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[0] - 0.5*
mWidth > 0);
630 if (!node_is_right_of_centre && element_centre_on_right)
633 this_node_index += num_nodes;
635 else if (node_is_right_of_centre && !element_centre_on_right)
638 this_node_index += 5*num_nodes;
641 else if (!element_straddles_left_right_boundary && element_straddles_top_bottom_boundary)
644 bool node_is_top_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[1] - 0.5*
mHeight > 0);
645 if (!node_is_top_of_centre && element_centre_on_top)
648 this_node_index += 3*num_nodes;
650 else if (node_is_top_of_centre && !element_centre_on_top)
653 this_node_index += 7*num_nodes;
656 else if (element_straddles_left_right_boundary && element_straddles_top_bottom_boundary)
658 bool node_is_right_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[0] - 0.5*
mWidth > 0);
659 bool node_is_top_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[1] - 0.5*
mHeight > 0);
661 if (!node_is_top_of_centre && element_centre_on_top)
663 if (!node_is_right_of_centre && element_centre_on_right)
665 this_node_index += 2*num_nodes;
667 else if (node_is_right_of_centre && !element_centre_on_right)
669 this_node_index += 4*num_nodes;
673 this_node_index += 3*num_nodes;
676 else if (node_is_top_of_centre && !element_centre_on_top)
678 if (!node_is_right_of_centre && element_centre_on_right)
680 this_node_index += 8*num_nodes;
682 else if (node_is_right_of_centre && !element_centre_on_right)
684 this_node_index += 6*num_nodes;
688 this_node_index += 7*num_nodes;
693 if (!node_is_right_of_centre && element_centre_on_right)
695 this_node_index += num_nodes;
697 else if (node_is_right_of_centre && !element_centre_on_right)
699 this_node_index += 5*num_nodes;
716 elem_nodes.push_back(temp_nodes[this_node_index]);
720 elements.push_back(p_element);
725 std::vector<Node<2>*> nodes;
727 for (
unsigned index=0; index<temp_nodes.size(); index++)
729 unsigned num_elems_containing_this_node = temp_nodes[index]->rGetContainingElementIndices().size();
731 if (num_elems_containing_this_node == 0)
734 delete temp_nodes[index];
738 temp_nodes[index]->SetIndex(count);
739 nodes.push_back(temp_nodes[index]);