53 : mWidth(rMesh.GetWidth(0)),
67 this->
mNodes.reserve(num_nodes);
72 for (
unsigned elem_index=0; elem_index<num_elements; elem_index++)
82 for (
unsigned i=0; i<
mNodes.size(); i++)
88 for (
unsigned i=0; i<num_nodes; i++)
92 for (
unsigned local_index=0; local_index<3; local_index++)
94 unsigned elem_index =
mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
95 unsigned num_nodes_in_elem =
mElements[elem_index]->GetNumNodes();
96 unsigned end_index = num_nodes_in_elem>0 ? num_nodes_in_elem-1 : 0;
105 std::vector<Node<2> *> nodes;
110 nodes.push_back(
new Node<2>(node_iter->GetIndex(), node_iter->rGetLocation(),node_iter->IsBoundaryNode()));
119 for (
unsigned j=0; j<3; j++)
127 std::set<unsigned> shared_elements;
128 std::set_intersection(node_a_element_indices.begin(),
129 node_a_element_indices.end(),
130 node_b_element_indices.begin(),
131 node_b_element_indices.end(),
132 std::inserter(shared_elements, shared_elements.begin()));
143 if (shared_elements.size() == 1)
146 c_vector<double,2> normal_vector;
148 normal_vector[0]= edge[1];
149 normal_vector[1]= -edge[0];
151 double dij = norm_2(normal_vector);
153 normal_vector /= dij;
166 c_vector<double,2> new_node_location = normal_vector + p_node_a->
rGetLocation() + ratio * edge;
169 bool node_clear =
true;
170 double node_clearance = 0.01;
172 for (
unsigned i=0; i<nodes.size(); i++)
174 double distance = norm_2(
mpDelaunayMesh->GetVectorFromAtoB(nodes[i]->rGetLocation(), new_node_location));
175 if (distance < node_clearance)
184 nodes.push_back(
new Node<2>(new_node_index, new_node_location));
196 bool bad_element =
false;
197 double edge_threshold = 1.5;
199 for (
unsigned j=0; j<3; j++)
212 for (
unsigned j=0; j<3; j++)
218 double edge_length = norm_2(edge);
221 if (edge_length<edge_threshold)
234 c_vector<double,2> normal_vector;
237 normal_vector[0]= edge[1];
238 normal_vector[1]= -edge[0];
240 double dij = norm_2(normal_vector);
242 normal_vector /= dij;
244 double bound_offset = 1.0;
245 c_vector<double,2> new_node_location = -bound_offset*normal_vector + p_node_a->
rGetLocation() + 0.5 * edge;
247 nodes.push_back(
new Node<2>(new_node_index, new_node_location));
255 for (
unsigned i=0; i<nodes.size(); i++)
266 this->
mNodes.reserve(num_nodes);
270 for (
unsigned elem_index = 0; elem_index < num_elements; elem_index++)
280 for (
unsigned i=0; i<
mNodes.size(); i++)
285 for (
unsigned i = 0; i < num_nodes; i++)
289 for (
unsigned local_index = 0; local_index < 3; local_index++)
291 unsigned elem_index = extended_mesh.
GetElement(i)->GetNodeGlobalIndex(local_index);
293 if (elem_index < num_elements)
295 unsigned num_nodes_in_elem =
mElements[elem_index]->GetNumNodes();
296 unsigned end_index = num_nodes_in_elem > 0 ? num_nodes_in_elem - 1 : 0;
305 for (
unsigned elem_index=0; elem_index<
mElements.size(); elem_index++)
312 std::vector<std::pair<double, unsigned> > index_angle_list;
313 for (
unsigned local_index=0; local_index<
mElements[elem_index]->GetNumNodes(); local_index++)
315 c_vector<double, 2> vectorA =
mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
316 c_vector<double, 2> vectorB =
mElements[elem_index]->GetNodeLocation(local_index);
317 c_vector<double, 2> centre_to_vertex =
mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
319 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
320 unsigned global_index =
mElements[elem_index]->GetNodeGlobalIndex(local_index);
322 std::pair<double, unsigned> pair(angle, global_index);
323 index_angle_list.push_back(pair);
327 sort(index_angle_list.begin(), index_angle_list.end());
331 for (
unsigned count = 0; count < index_angle_list.size(); count++)
333 unsigned local_index = count>1 ? count-1 : 0;
334 p_new_element->
AddNode(
mNodes[index_angle_list[count].second], local_index);
448 std::vector<Node<2>*> temp_nodes(3 * num_nodes);
449 std::vector<VertexElement<2, 2>*> elements;
452 for (
unsigned index=0; index<num_nodes; index++)
454 c_vector<double, 2> location;
455 location =
GetNode(index)->rGetLocation();
458 Node<2>* p_node =
new Node<2>(index,
false, location[0], location[1]);
459 temp_nodes[index] = p_node;
462 p_node =
new Node<2>(num_nodes + index,
false, location[0] +
mWidth, location[1]);
463 temp_nodes[num_nodes + index] = p_node;
466 p_node =
new Node<2>(2 * num_nodes + index,
false, location[0] -
mWidth, location[1]);
467 temp_nodes[2*num_nodes + index] = p_node;
475 unsigned elem_index = elem_iter->
GetIndex();
476 unsigned num_nodes_in_elem = elem_iter->GetNumNodes();
478 std::vector<Node<2>*> elem_nodes;
481 bool element_straddles_left_right_boundary =
false;
483 const c_vector<double, 2>& r_this_node_location = elem_iter->GetNode(0)->rGetLocation();
484 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
486 const c_vector<double, 2>& r_next_node_location = elem_iter->GetNode((local_index+1)%num_nodes_in_elem)->rGetLocation();
487 c_vector<double, 2> vector;
488 vector = r_next_node_location - r_this_node_location;
490 if (fabs(vector[0]) > 0.5 *
mWidth)
492 element_straddles_left_right_boundary =
true;
499 bool element_centre_on_right =
true;
503 double element_centre_x_location = this->
mpDelaunayMesh->GetNode(delaunay_index)->rGetLocation()[0];
504 if (element_centre_x_location < 0.5 *
mWidth)
506 element_centre_on_right =
false;
511 for (
unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
513 unsigned this_node_index = elem_iter->GetNodeGlobalIndex(local_index);
516 if (element_straddles_left_right_boundary)
519 bool node_is_right_of_centre = (elem_iter->GetNode(local_index)->rGetLocation()[0] - 0.5 *
mWidth > 0);
520 if (!node_is_right_of_centre && element_centre_on_right)
523 this_node_index += num_nodes;
525 else if (node_is_right_of_centre && !element_centre_on_right)
528 this_node_index += 2 * num_nodes;
532 elem_nodes.push_back(temp_nodes[this_node_index]);
536 elements.push_back(p_element);
540 std::vector<Node<2>*> nodes;
542 for (
unsigned index=0; index<temp_nodes.size(); index++)
544 unsigned num_elems_containing_this_node = temp_nodes[index]->rGetContainingElementIndices().size();
546 if (num_elems_containing_this_node == 0)
549 delete temp_nodes[index];
553 temp_nodes[index]->SetIndex(count);
554 nodes.push_back(temp_nodes[index]);