quad_from_triangle_mesh.template.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 0.90. August 3, 2009.
7 //LIC//
8 //LIC// Copyright (C) 2006-2009 Matthias Heil and Andrew Hazel
9 //LIC//
10 //LIC// This library is free software; you can redistribute it and/or
11 //LIC// modify it under the terms of the GNU Lesser General Public
12 //LIC// License as published by the Free Software Foundation; either
13 //LIC// version 2.1 of the License, or (at your option) any later version.
14 //LIC//
15 //LIC// This library is distributed in the hope that it will be useful,
16 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
17 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 //LIC// Lesser General Public License for more details.
19 //LIC//
20 //LIC// You should have received a copy of the GNU Lesser General Public
21 //LIC// License along with this library; if not, write to the Free Software
22 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
23 //LIC// 02110-1301 USA.
24 //LIC//
25 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
26 //LIC//
27 //LIC//====================================================================
28 //Driver for 2D moving block
29 #ifndef OOMPH_QUAD_FROM_TRIANGLE_MESH_TEMPLATE_CC
30 #define OOMPH_QUAD_FROM_TRIANGLE_MESH_TEMPLATE_CC
31 
32 // The mesh
34 
35 using namespace std;
36 using namespace oomph;
37 
38 
39 ////////////////////////////////////////////////////////////////////
40 ////////////////////////////////////////////////////////////////////
41 ////////////////////////////////////////////////////////////////////
42 
43 
44 
45 namespace oomph
46 {
47 
48 //======================================================================
49 /// \short Build the full mesh with the help of the scaffold mesh coming
50 /// from the triangle mesh generator, Triangle. To build this quad
51 /// element based mesh we make use of the fact that a triangle element
52 /// can be split as shown in the diagram below:
53 ///
54 /// N2
55 /// | | N0 : 1st scaffold element node
56 /// | | N1 : 2nd scaffold element node
57 /// | | N2 : 3rd scaffold element node
58 /// | |
59 /// C | Q_2 | B Edge 0 : N0 --> N1
60 /// | | | | Edge 1 : N1 --> N2
61 /// | | | | Edge 2 : N2 --> N0
62 /// | | | |
63 /// | | | | A : Midpoint of edge 0
64 /// | Q_0 | Q_1 | B : Midpoint of edge 1
65 /// | | | C : Midpoint of edge 2
66 /// | | |
67 /// N0 __________|__________ N1
68 /// A
69 ///
70 /// The intersection of all three quad elements is the centroid. Using
71 /// this splitting, the subsequent mesh will consist of quadrilaterals
72 /// whose shape which depend on the structure of the underlying mesh.
73 //======================================================================
74  template<class ELEMENT>
76  TriangleScaffoldMesh* tmp_mesh_pt,TimeStepper* time_stepper_pt,
77  const bool &use_attributes)
78  {
79  // Create space for elements
80  unsigned nelem=tmp_mesh_pt->nelement();
81 
82  // We will have 3 quad elements per scaffold element
83  Element_pt.resize(3*nelem);
84 
85  // Set number of boundaries
86  unsigned nbound=tmp_mesh_pt->nboundary();
87 
88  // Resize the boundary information (the number of boundaries doesn't
89  // change)
90  set_nboundary(nbound);
91 
92  // Stores each element attached to a boundary and the index of the
93  // face of the given element attached to the boundary
94  Boundary_element_pt.resize(nbound);
95  Face_index_at_boundary.resize(nbound);
96 
97  // Create a quad element for nodal data
98  ELEMENT* temp_el_pt=new ELEMENT;
99 
100  // Get the number of nodes in a quad element
101  unsigned nnode_el=temp_el_pt->nnode();
102 
103  // Find the number of nodes along one edge of a quad element
104  unsigned nnode_1d=temp_el_pt->nnode_1d();
105 
106  // Calculate the number of nodes that will lie along an edge of a
107  // triangle element in the scaffold mesh
108  unsigned nnode_edge=2*nnode_1d-1;
109 
110  // Delete the element pointer
111  delete temp_el_pt;
112 
113  // Make it a null pointer
114  temp_el_pt=0;
115 
116  // Create dummy linear quad for geometry
117  QElement<2,2> dummy_element;
118 
119  // The dimension of the element
120  unsigned n_dim=2;
121 
122  // The position type
123  unsigned n_position_type=1;
124 
125  // Don't assign memory for any values
126  unsigned initial_n_value=0;
127 
128  // Loop over the nodes of the element and make them
129  for (unsigned j=0;j<4;j++)
130  {
131  dummy_element.node_pt(j)=new Node(n_dim,n_position_type,initial_n_value);
132  }
133 
134  // Local node number of each quad element corner
135  unsigned corner_0=0;
136  unsigned corner_1=nnode_1d-1;
137  unsigned corner_2=nnode_el-nnode_1d;
138  unsigned corner_3=nnode_el-1;
139 
140  // Create a map to return a vector of pointers to nnode_1d nodes where
141  // the input is an edge. If the edge hasn't been set up then this will
142  // return a null pointer. Note: all node pointers on an edge will be
143  // stored in clockwise ordering. Therefore, to copy the data of an
144  // edge into the adjoining element we must proceed through the vector
145  // backwards (as progressing through an edge of an element in a clockwise
146  // manner is equivalent to proceeding through the edge of the neighbouring
147  // element in an anti-clockwise manner)
148  std::map<Edge,Vector<Node*> > edge_nodes_map;
149 
150  // Set up a map to check if the scaffold mesh node has been set up in the
151  // quad mesh. If the node has been set up this map will return a pointer
152  // to it otherwise it will return a null pointer
153  std::map<Node*,Node*> scaffold_to_quad_mesh_node;
154 
155  // Loop over elements in scaffold mesh
156  unsigned new_el_count=0;
157 
158  // Create storage for the coordinates of the centroid
159  Vector<double> centroid(2);
160 
161  // Create storage for the coordinates of the vertices of the triangle
162  Vector<Vector<double> > triangle_vertex(3);
163 
164  // Loop over all of the elements in the scaffold mesh
165  for (unsigned e=0;e<nelem;e++)
166  {
167  // Initialise centroid values for the e-th triangle element
168  centroid[0]=0.0;
169  centroid[1]=0.0;
170 
171  // Loop over the scaffold element nodes
172  for (unsigned j=0;j<3;j++)
173  {
174  // Resize the j-th triangle_vertex entry to contain the x and
175  // y-coordinate
176  triangle_vertex[j].resize(2);
177 
178  // Get the coordinates
179  double x=tmp_mesh_pt->finite_element_pt(e)->node_pt(j)->x(0);
180  double y=tmp_mesh_pt->finite_element_pt(e)->node_pt(j)->x(1);
181 
182  // Increment the centroid coordinates
183  centroid[0]+=x;
184  centroid[1]+=y;
185 
186  // Assign the triangle_vertex coordinates
187  triangle_vertex[j][0]=x;
188  triangle_vertex[j][1]=y;
189  }
190 
191  // Divide the centroid entries by 3 to get the centroid coordinates
192  centroid[0]/=3.0;
193  centroid[1]/=3.0;
194 
195  // Create element pointers and assign them to a vector
196  //----------------------------------------------------
197  // Make new quad elements of the type specified by the template parameter
198  ELEMENT* el0_pt=new ELEMENT;
199  ELEMENT* el1_pt=new ELEMENT;
200  ELEMENT* el2_pt=new ELEMENT;
201 
202  // Create a vector of ELEMENTs to store el0_pt, el1_pt and el2_pt
203  Vector<ELEMENT*> el_vector_pt(3,0);
204 
205  // Assign the entries to el_vector_pt
206  el_vector_pt[0]=el0_pt;
207  el_vector_pt[1]=el1_pt;
208  el_vector_pt[2]=el2_pt;
209 
210 
211  // Create the first node in each quad element and store in Node_pt.
212  // These correspond to the nodes of the simplex triangle stored in
213  // Tmp_mesh_pt. If they have already been set up then we do nothing:
214  //------------------------------------------------------------------
215  // Loop over the scaffold element nodes and check to see if they have
216  // been set up
217  for (unsigned j=0;j<3;j++)
218  {
219  // Pointer to node in the scaffold mesh
220  Node* scaffold_node_pt=tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
221 
222  // Check if the node has been set up yet
223  Node* qmesh_node_pt=scaffold_to_quad_mesh_node[scaffold_node_pt];
224 
225  // Haven't done this one yet
226  if (qmesh_node_pt==0)
227  {
228  // Get pointer to set of mesh boundaries that this
229  // scaffold node occupies; NULL if the node is not on any boundary
230  std::set<unsigned>* boundaries_pt;
231  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
232 
233  // Check to see if it's on any boundaries
234  if (boundaries_pt!=0)
235  {
236  // Create new boundary node. The scaffold element nodes are the corners
237  // of a simplex triangle and thus always correspond to the first node
238  // in each quad element
239  qmesh_node_pt=el_vector_pt[j]->construct_boundary_node(corner_0,
240  time_stepper_pt);
241 
242  // Add to boundaries
243  for (std::set<unsigned>::iterator it=boundaries_pt->begin();
244  it!=boundaries_pt->end();++it)
245  {
246  add_boundary_node(*it,qmesh_node_pt);
247  }
248  }
249  // Build normal node
250  else
251  {
252  // Create new normal node
253  qmesh_node_pt=el_vector_pt[j]->construct_node(corner_0,time_stepper_pt);
254  }
255 
256  // Add the mapping from the scaffold mesh node to the quad mesh node
257  scaffold_to_quad_mesh_node[scaffold_node_pt]=qmesh_node_pt;
258 
259  // Copy new node, created using the NEW element's construct_node
260  // function into global storage, using the same global
261  // node number that we've just associated with the
262  // corresponding node in the scaffold mesh
263  Node_pt.push_back(qmesh_node_pt);
264  }
265  // If this node has already been done we need to copy the data across
266  else
267  {
268  el_vector_pt[j]->node_pt(corner_0)=qmesh_node_pt;
269  }
270 
271 
272  // Set global coordinate
273  el_vector_pt[j]->node_pt(corner_0)->x(0)=triangle_vertex[j][0];
274  el_vector_pt[j]->node_pt(corner_0)->x(1)=triangle_vertex[j][1];
275  }
276 
277 
278  // Create the edges of the scaffold element and check to see if
279  // they've been set up yet or not. If they haven't:
280  //--------------------------------------------------------------
281  // Make the three edges of the triangle
282  Edge edge0(tmp_mesh_pt->finite_element_pt(e)->node_pt(0),
283  tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
284  Edge edge1(tmp_mesh_pt->finite_element_pt(e)->node_pt(1),
285  tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
286  Edge edge2(tmp_mesh_pt->finite_element_pt(e)->node_pt(2),
287  tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
288 
289  // Check if the edges have been set up (each will have size nnode_1d).
290  // If they have not been set up yet, this will
291  Vector<Node*> edge0_vector_pt=edge_nodes_map[edge0];
292  Vector<Node*> edge1_vector_pt=edge_nodes_map[edge1];
293  Vector<Node*> edge2_vector_pt=edge_nodes_map[edge2];
294 
295  // Bools to indicate whether or not the edges have been set up
296  bool edge0_setup=(edge0_vector_pt.size()!=0);
297  bool edge1_setup=(edge1_vector_pt.size()!=0);
298  bool edge2_setup=(edge2_vector_pt.size()!=0);
299 
300  // If edge 0 hasn't been set up (node 0 to node 1)
301  if (!edge0_setup)
302  {
303  // Resize the vector to have length nnode_1d
304  edge0_vector_pt.resize(nnode_edge,0);
305 
306  // First node along edge 0 is the first node of element 0
307  edge0_vector_pt[0]=el_vector_pt[0]->node_pt(0);
308 
309  // Last node along edge 0 is the first node of element 1
310  edge0_vector_pt[nnode_edge-1]=el_vector_pt[1]->node_pt(0);
311  }
312 
313  // If edge 1 hasn't been set up (node 1 to node 2)
314  if (!edge1_setup)
315  {
316  // Resize the vector to have length nnode_1d
317  edge1_vector_pt.resize(nnode_edge,0);
318 
319  // First node along edge 1 is the first node of element 1
320  edge1_vector_pt[0]=el_vector_pt[1]->node_pt(0);
321 
322  // Last node along edge 1 is the first node of element 2
323  edge1_vector_pt[nnode_edge-1]=el_vector_pt[2]->node_pt(0);
324  }
325 
326  // If edge 2 hasn't been set up (node 2 to node 0)
327  if (!edge2_setup)
328  {
329  // Resize the vector to have length nnode_1d
330  edge2_vector_pt.resize(nnode_edge,0);
331 
332  // First node along edge 2 is the first node of element 2
333  edge2_vector_pt[0]=el_vector_pt[2]->node_pt(0);
334 
335  // Last node along edge 2 is the first node of element 0
336  edge2_vector_pt[nnode_edge-1]=el_vector_pt[0]->node_pt(0);
337  }
338 
339 
340 #ifdef PARANOID
341  // If any of the edges have been set up, make sure that that the endpoints
342  // in the returned vectors have the same address as those on the vertices
343 
344  // Come back and finish this off.
345  // To check:
346  // - If two edges which have been set up have the same node in the middle
347  // - If an edge has already been set up then the map will return the
348  // same node as in the vector
349 #endif
350 
351  // Boundary IDs for bottom and left edge of quad
352  // from scaffold mesh (if these remain zero the edges
353  // are not on a boundary)
354  unsigned q0_lower_boundary_id=0;
355  unsigned q0_left_boundary_id=0;
356  unsigned q1_lower_boundary_id=0;
357  unsigned q1_left_boundary_id=0;
358  unsigned q2_lower_boundary_id=0;
359  unsigned q2_left_boundary_id=0;
360 
361  // Lower/left boundary IDs for quad 0; the lower edge in quad 0 is on
362  // edge 0 of the scaffold triangle and the left edge in quad is on edge
363  // 2 in scaffold triangle
364  q0_lower_boundary_id=tmp_mesh_pt->edge_boundary(e,0);
365  q0_left_boundary_id=tmp_mesh_pt->edge_boundary(e,2);
366 
367  // Lower/left boundary IDs for quad 1; the lower edge in quad 1 is on
368  // edge 1 of the scaffold triangle and the left edge in quad is on edge
369  // 0 of the scaffold triangle
370  q1_lower_boundary_id=tmp_mesh_pt->edge_boundary(e,1);
371  q1_left_boundary_id=tmp_mesh_pt->edge_boundary(e,0);
372 
373  // Lower/left boundary IDs for quad 2; the lower edge in quad 2 is on
374  // edge 2 of the scaffold triangle and the left edge in quad is on edge
375  // 1 of the scaffold triangle
376  q2_lower_boundary_id=tmp_mesh_pt->edge_boundary(e,2);
377  q2_left_boundary_id=tmp_mesh_pt->edge_boundary(e,1);
378 
379  // Store the boundary IDs as a vector; allows us to loop over them easily
380  Vector<unsigned> boundary_id_vector(6,0);
381  boundary_id_vector[0]=q0_lower_boundary_id;
382  boundary_id_vector[1]=q0_left_boundary_id;
383  boundary_id_vector[2]=q1_lower_boundary_id;
384  boundary_id_vector[3]=q1_left_boundary_id;
385  boundary_id_vector[4]=q2_lower_boundary_id;
386  boundary_id_vector[5]=q2_left_boundary_id;
387 
388  // Loop over the quad elements and store the boundary elements in the vector
389  // Boundary_element_pt
390  for (unsigned j=0;j<3;j++)
391  {
392  // Loop over the lower and the left boundary ID in the j'th element
393  for (unsigned k=0;k<2;k++)
394  {
395  // The quad element lies on a boundary of the mesh
396  if (boundary_id_vector[2*j+k]>0)
397  {
398  // Since the j'th quad element lies on a boundary of the mesh we add a
399  // pointer to the element to the appropriate entry of Boundary_element_pt
400  Boundary_element_pt[boundary_id_vector[2*j+k]-1].
401  push_back(el_vector_pt[j]);
402 
403  // If k=0 then the lower boundary of the quad element lies on
404  // the boundary of the mesh and if k=1 then the left boundary
405  // of the quad element lies on the mesh boundary. For quad elements
406  // the indices are as follows:
407  // North face: 2
408  // East face: 1
409  // South face: -2
410  // West face: -1
411  if (k==0)
412  {
413  Face_index_at_boundary[boundary_id_vector[2*j+k]-1].push_back(-2);
414  }
415  else
416  {
417  Face_index_at_boundary[boundary_id_vector[2*j+k]-1].push_back(-1);
418  } // if (k==0)
419  } // if (boundary_id_vector[2*j+k]>0)
420  } // for (unsigned k=0;k<2;k++)
421  } // for (unsigned j=0;j<3;j++)
422 
423 
424  // The upper right node is always the centroid. Note: The 'corner_3' node
425  // lies within each of the three quad elements so we simply share the pointer
426  // to it with each element:
427  //---------------------------------------------------------------------------
428  // Construct the centroid node
429  Node* nod_pt=el0_pt->construct_node(corner_3,time_stepper_pt);
430 
431  // Add the pointer to the vector of nodal pointers
432  Node_pt.push_back(nod_pt);
433 
434  // Quad 0
435  el0_pt->node_pt(corner_3)->x(0)=centroid[0];
436  el0_pt->node_pt(corner_3)->x(1)=centroid[1];
437 
438  // Quad 1
439  el1_pt->node_pt(corner_3)=el0_pt->node_pt(corner_3);
440 
441  // Quad 2
442  el2_pt->node_pt(corner_3)=el0_pt->node_pt(corner_3);
443 
444 
445  // Set the nodal positions of the dummy element to emulate the FIRST
446  // quad element (this allows for simple linear interpolation later):
447  //------------------------------------------------------------------
448  // Bottom-left corner
449  dummy_element.node_pt(0)->x(0)=triangle_vertex[0][0];
450  dummy_element.node_pt(0)->x(1)=triangle_vertex[0][1];
451 
452  // Bottom-right corner
453  dummy_element.node_pt(1)->x(0)=
454  0.5*(triangle_vertex[0][0]+triangle_vertex[1][0]);
455  dummy_element.node_pt(1)->x(1)=
456  0.5*(triangle_vertex[0][1]+triangle_vertex[1][1]);
457 
458  // Top-left corner
459  dummy_element.node_pt(2)->x(0)=
460  0.5*(triangle_vertex[0][0]+triangle_vertex[2][0]);
461  dummy_element.node_pt(2)->x(1)=
462  0.5*(triangle_vertex[0][1]+triangle_vertex[2][1]);
463 
464  // Top-right corner
465  dummy_element.node_pt(3)->x(0)=centroid[0];
466  dummy_element.node_pt(3)->x(1)=centroid[1];
467 
468 
469  // Set up all of the nodes in the first quad element (Q0):
470  //--------------------------------------------------------
471  // Local and global coordinate vectors for the nodes
472  Vector<double> s(2);
473  Vector<double> x(2);
474 
475  // Loop over all of nodes in Q0 noting that the lower left corner node
476  // (node 0) and the upper right corner node (centroid) have already
477  // been set up
478  for (unsigned j=1;j<corner_3;j++)
479  {
480  // Indicates whether or not the node has been set up yet
481  bool done=false;
482 
483  // On the lower edge
484  if (j<nnode_1d)
485  {
486  // If the lower edge has already been set up then we already know the
487  // nodes along this edge
488  if (edge0_setup)
489  {
490  // The j'th node along this edge is the (nnode_1d-j)'th node in the
491  // vector (remembering that the ordering is backwards since it has
492  // already been set up)
493  el0_pt->node_pt(j)=edge0_vector_pt[(nnode_edge-1)-j];
494 
495  // Since the node has already been set up we do not need to sort
496  // out its global coordinate data so skip to the next node
497  continue;
498  }
499  // If the edge hasn't been set up yet
500  else
501  {
502  // If the node lies on a boundary too then we need to construct a
503  // boundary node
504  if (q0_lower_boundary_id>0)
505  {
506  // Construct a boundary node
507  Node* nod_pt=el0_pt->construct_boundary_node(j,time_stepper_pt);
508 
509  // Add it to the list of boundary nodes
510  add_boundary_node(q0_lower_boundary_id-1,nod_pt);
511 
512  // Add the node into the vector of nodes on edge 0
513  edge0_vector_pt[j]=nod_pt;
514 
515  // Add it to the list of nodes in the mesh
516  Node_pt.push_back(nod_pt);
517 
518  // Indicate the j'th node has been set up
519  done=true;
520  }
521  }
522 
523  // Node is not on a mesh boundary but on the lower edge
524  if (!done)
525  {
526  // Construct a normal node
527  Node* nod_pt=el0_pt->construct_node(j,time_stepper_pt);
528 
529  // Add the node into the vector of nodes on edge 0
530  edge0_vector_pt[j]=nod_pt;
531 
532  // Add it to the list of nodes in the mesh
533  Node_pt.push_back(nod_pt);
534 
535  // Indicate the j'th node has been set up
536  done=true;
537  }
538  }
539  // On the left edge
540  else if (j%nnode_1d==0)
541  {
542  // If the left edge has already been set up then we already know the
543  // nodes along this edge
544  if (edge2_setup)
545  {
546  // The j'th node is the (j/nnode_1d)'th node along this edge and thus
547  // the (j/nnode_1d)'th entry in the edge vector
548  el0_pt->node_pt(j)=edge2_vector_pt[j/nnode_1d];
549 
550  // Since the node has already been set up we do not need to sort
551  // out its global coordinate data
552  continue;
553  }
554  // If the edge hasn't been set up yet
555  else
556  {
557  if (q0_left_boundary_id>0)
558  {
559  // Construct a boundary node
560  Node* nod_pt=el0_pt->construct_boundary_node(j,time_stepper_pt);
561 
562  // Add it to the list of boundary nodes
563  add_boundary_node(q0_left_boundary_id-1,nod_pt);
564 
565  // Add the node into the vector of nodes on edge 2 in clockwise order
566  edge2_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
567 
568  // Add it to the list of nodes in the mesh
569  Node_pt.push_back(nod_pt);
570 
571  // Indicate that the j'th node has been set up
572  done=true;
573  }
574  }
575 
576  // Node is not on a mesh boundary but on the left edge
577  if (!done)
578  {
579  // Construct a normal node
580  Node* nod_pt=el0_pt->construct_node(j,time_stepper_pt);
581 
582  // Add the node into the vector of nodes on edge 2 in clockwise order
583  edge2_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
584 
585  // Add it to the list of nodes in the mesh
586  Node_pt.push_back(nod_pt);
587 
588  // Indicate the j'th node has been set up
589  done=true;
590  }
591  }
592 
593  // Node is not on a mesh boundary or on the edge of the scaffold element
594  if (!done)
595  {
596  // Construct a normal node
597  Node* nod_pt=el0_pt->construct_node(j,time_stepper_pt);
598 
599  // Add it to the list of nodes in the mesh
600  Node_pt.push_back(nod_pt);
601  }
602 
603  // Get local coordinate
604  el0_pt->local_coordinate_of_node(j,s);
605 
606  // Interpolate position linearly between vertex nodes
607  dummy_element.interpolated_x(s,x);
608  el0_pt->node_pt(j)->x(0)=x[0];
609  el0_pt->node_pt(j)->x(1)=x[1];
610  }
611 
612 
613  // Set the nodal positions of the dummy element to now emulate the
614  // SECOND quad element:
615  //------------------------------------------------------------------
616  // Note: we do not need to change the top-right corner since it always
617  // coincides with the centroid of the triangle element
618 
619  // Bottom-left corner
620  dummy_element.node_pt(0)->x(0)=triangle_vertex[1][0];
621  dummy_element.node_pt(0)->x(1)=triangle_vertex[1][1];
622 
623  // Bottom-right corner
624  dummy_element.node_pt(1)->x(0)=
625  0.5*(triangle_vertex[1][0]+triangle_vertex[2][0]);
626  dummy_element.node_pt(1)->x(1)=
627  0.5*(triangle_vertex[1][1]+triangle_vertex[2][1]);
628 
629  // Top-left corner
630  dummy_element.node_pt(2)->x(0)=
631  0.5*(triangle_vertex[0][0]+triangle_vertex[1][0]);
632  dummy_element.node_pt(2)->x(1)=
633  0.5*(triangle_vertex[0][1]+triangle_vertex[1][1]);
634 
635 
636  // Set up all of the nodes in the second quad element (Q1):
637  //--------------------------------------------------------
638  // Here we need to notice that we have already set up the final nnode_1d
639  // nodes (the upper edge of Q1 coincides with the right edge of Q0)
640 
641  // Loop over nodes 1 to (corner_2-1) in Q1 noting that the lower left
642  // corner node (node 0) and the upper edge of Q1 contains nodes
643  // corner_2 to corner_3
644  for (unsigned j=1;j<corner_2;j++)
645  {
646  // Indicates whether or not the node has been set up yet
647  bool done=false;
648 
649  // On the lower edge
650  if (j<nnode_1d)
651  {
652  // If the lower edge has already been set up then we already know the
653  // nodes along this edge
654  if (edge1_setup)
655  {
656  // The j'th node along this edge is the (nnode_1d-j)'th node in the
657  // vector (remembering that the ordering is backwards if it has
658  // already been set up)
659  el1_pt->node_pt(j)=edge1_vector_pt[(nnode_edge-1)-j];
660 
661  // Since the node has already been set up we do not need to sort
662  // out its global coordinate data
663  continue;
664  }
665  // If the edge hasn't been set up yet
666  else
667  {
668  // If the node lies on a boundary too then we need to construct a
669  // boundary node
670  if (q1_lower_boundary_id>0)
671  {
672  // Construct a boundary node
673  Node* nod_pt=el1_pt->construct_boundary_node(j,time_stepper_pt);
674 
675  // Add it to the list of boundary nodes
676  add_boundary_node(q1_lower_boundary_id-1,nod_pt);
677 
678  // Add the node into the vector of nodes on edge 1
679  edge1_vector_pt[j]=nod_pt;
680 
681  // Add it to the list of nodes in the mesh
682  Node_pt.push_back(nod_pt);
683 
684  // Indicate the j'th node has been set up
685  done=true;
686  }
687  }
688 
689  // Node is not on a mesh boundary but on the lower edge
690  if (!done)
691  {
692  // Construct a normal node
693  Node* nod_pt=el1_pt->construct_node(j,time_stepper_pt);
694 
695  // Add the node into the vector of nodes on edge 1
696  edge1_vector_pt[j]=nod_pt;
697 
698  // Add it to the list of nodes in the mesh
699  Node_pt.push_back(nod_pt);
700 
701  // Indicate the j'th node has been set up
702  done=true;
703  }
704  }
705  // On the left edge
706  else if (j%nnode_1d==0)
707  {
708  // If the left edge has already been set up then we already know the
709  // nodes along this edge
710  if (edge0_setup)
711  {
712  // The j'th node along this edge is the (j/nnode_1d)'th node in the
713  // vector
714  el1_pt->node_pt(j)=edge0_vector_pt[j/nnode_1d];
715 
716  // Since the node has already been set up we do not need to sort
717  // out its global coordinate data
718  continue;
719  }
720  // If the edge hasn't been set up yet
721  else
722  {
723  if (q1_left_boundary_id>0)
724  {
725  // Construct a boundary node
726  Node* nod_pt=el1_pt->construct_boundary_node(j,time_stepper_pt);
727 
728  // Add it to the list of boundary nodes
729  add_boundary_node(q1_left_boundary_id-1,nod_pt);
730 
731  // Add the node into the vector of nodes on edge 0 in clockwise order
732  edge0_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
733 
734  // Add it to the list of nodes in the mesh
735  Node_pt.push_back(nod_pt);
736 
737  // Indicate that the j'th node has been set up
738  done=true;
739  }
740  }
741 
742  // Node is not on a mesh boundary but on the left edge
743  if (!done)
744  {
745  // Construct a normal node
746  Node* nod_pt=el1_pt->construct_node(j,time_stepper_pt);
747 
748  // Add the node into the vector of nodes on edge 0 in clockwise order
749  edge0_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
750 
751  // Add it to the list of nodes in the mesh
752  Node_pt.push_back(nod_pt);
753 
754  // Indicate the j'th node has been set up
755  done=true;
756  }
757  }
758 
759  // Node is not on a mesh boundary or the scaffold element boundary
760  if (!done)
761  {
762  // Construct a normal node
763  Node* nod_pt=el1_pt->construct_node(j,time_stepper_pt);
764 
765  // Add it to the list of nodes in the mesh
766  Node_pt.push_back(nod_pt);
767  }
768 
769  // Get local coordinate
770  el1_pt->local_coordinate_of_node(j,s);
771 
772  // Interpolate position linearly between vertex nodes
773  dummy_element.interpolated_x(s,x);
774  el1_pt->node_pt(j)->x(0)=x[0];
775  el1_pt->node_pt(j)->x(1)=x[1];
776  }
777 
778 
779  // We now need to loop over nodes corner_2 to (corner_3-1) and copy the
780  // given information from Q0. We do not need to set up the (corner_3)'th
781  // node since it coincides with the centroid which has already been set up
782  for (unsigned j=corner_2;j<corner_3;j++)
783  {
784  // The nodes along the upper edge of Q1 go from corner_2 to corner_3-1
785  // while the nodes along the right edge of Q0 go from corner_1 to
786  // (corner_3-nnode_1d) in increments of nnode_1d
787  el1_pt->node_pt(j)=el0_pt->node_pt(corner_1+(j-corner_2)*nnode_1d);
788  }
789 
790 
791  // Set the nodal positions of the dummy element to now emulate the
792  // THIRD quad element:
793  //------------------------------------------------------------------
794  // Note: we do not need to change the top-right corner since it always
795  // coincides with the centroid of the triangle element
796 
797  // Bottom-left corner
798  dummy_element.node_pt(0)->x(0)=triangle_vertex[2][0];
799  dummy_element.node_pt(0)->x(1)=triangle_vertex[2][1];
800 
801  // Bottom-right corner
802  dummy_element.node_pt(1)->x(0)=
803  0.5*(triangle_vertex[0][0]+triangle_vertex[2][0]);
804  dummy_element.node_pt(1)->x(1)=
805  0.5*(triangle_vertex[0][1]+triangle_vertex[2][1]);
806 
807  // Top-left corner
808  dummy_element.node_pt(2)->x(0)=
809  0.5*(triangle_vertex[1][0]+triangle_vertex[2][0]);
810  dummy_element.node_pt(2)->x(1)=
811  0.5*(triangle_vertex[1][1]+triangle_vertex[2][1]);
812 
813 
814  // Set up all of the nodes in the third quad element (Q2):
815  //--------------------------------------------------------
816  // Here we need to notice that we have already set up the final nnode_1d
817  // nodes (the upper edge of Q2 coincides with the right edge of Q1).
818  // We have also already set up the nodes on the right edge of Q2 (the
819  // right edge of Q2 coincides with the upper edge of Q0)
820 
821  // Loop over nodes 1 to (corner_2-1)
822  for (unsigned j=1;j<corner_2;j++)
823  {
824  // Indicates whether or not the node has been set up yet
825  bool done=false;
826 
827  // On the lower edge
828  if (j<nnode_1d-1)
829  {
830  // If the lower edge has already been set up then we already know the
831  // nodes along this edge
832  if (edge2_setup)
833  {
834  // The j'th node along this edge is the (nnode_1d-j)'th node in the
835  // vector (remembering that the ordering is backwards if it has
836  // already been set up)
837  el2_pt->node_pt(j)=edge2_vector_pt[(nnode_edge-1)-j];
838 
839  // Since the node has already been set up we do not need to sort
840  // out its global coordinate data
841  continue;
842  }
843  // If the edge hasn't been set up yet
844  else
845  {
846  // If the node lies on a boundary too then we need to construct a
847  // boundary node
848  if (q2_lower_boundary_id>0)
849  {
850  // Construct a boundary node
851  Node* nod_pt=el2_pt->construct_boundary_node(j,time_stepper_pt);
852 
853  // Add it to the list of boundary nodes
854  add_boundary_node(q2_lower_boundary_id-1,nod_pt);
855 
856  // Add the node into the vector of nodes on edge 2
857  edge2_vector_pt[j]=nod_pt;
858 
859  // Add it to the list of nodes in the mesh
860  Node_pt.push_back(nod_pt);
861 
862  // Indicate the j'th node has been set up
863  done=true;
864  }
865  }
866 
867  // Node is not on a mesh boundary but on the lower edge
868  if (!done)
869  {
870  // Construct a normal node
871  Node* nod_pt=el2_pt->construct_node(j,time_stepper_pt);
872 
873  // Add the node into the vector of nodes on edge 2
874  edge2_vector_pt[j]=nod_pt;
875 
876  // Add it to the list of nodes in the mesh
877  Node_pt.push_back(nod_pt);
878 
879  // Indicate the j'th node has been set up
880  done=true;
881  }
882  }
883  // On the right edge
884  else if ((j+1)%nnode_1d==0)
885  {
886  // Copy the data from the top edge of element 0 to element 2
887  el2_pt->node_pt(j)=el0_pt->node_pt((corner_2-1)+(j+1)/nnode_1d);
888 
889  // We don't need to set up the global coordinate data so just
890  // skip to the next node in the element
891  continue;
892  }
893  // On the left edge
894  else if (j%nnode_1d==0)
895  {
896  // If the left edge has already been set up then we already know the
897  // nodes along this edge
898  if (edge1_setup)
899  {
900  // The j'th node along this edge is the (j/nnode_1d)'th node in the
901  // vector
902  el2_pt->node_pt(j)=edge1_vector_pt[j/nnode_1d];
903 
904  // Since the node has already been set up we do not need to sort
905  // out its global coordinate data
906  continue;
907  }
908  // If the edge hasn't been set up yet
909  else
910  {
911  if (q2_left_boundary_id>0)
912  {
913  // Construct a boundary node
914  Node* nod_pt=el2_pt->construct_boundary_node(j,time_stepper_pt);
915 
916  // Add it to the list of boundary nodes
917  add_boundary_node(q2_left_boundary_id-1,nod_pt);
918 
919  // Add the node into the vector of nodes on edge 1 in clockwise order
920  edge1_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
921 
922  // Add it to the list of nodes in the mesh
923  Node_pt.push_back(nod_pt);
924 
925  // Indicate that the j'th node has been set up
926  done=true;
927  }
928  }
929 
930  // Node is not on a mesh boundary but on the left edge
931  if (!done)
932  {
933  // Construct a normal node
934  Node* nod_pt=el2_pt->construct_node(j,time_stepper_pt);
935 
936  // Add the node into the vector of nodes on edge 1 in clockwise order
937  edge1_vector_pt[(nnode_edge-1)-(j/nnode_1d)]=nod_pt;
938 
939  // Add it to the list of nodes in the mesh
940  Node_pt.push_back(nod_pt);
941 
942  // Indicate the j'th node has been set up
943  done=true;
944  }
945  }
946 
947  // Node is not on a mesh boundary
948  if (!done)
949  {
950  // Construct a normal node
951  Node* nod_pt=el2_pt->construct_node(j,time_stepper_pt);
952 
953  // Add it to the list of nodes in the mesh
954  Node_pt.push_back(nod_pt);
955  }
956 
957  // Get local coordinate
958  el2_pt->local_coordinate_of_node(j,s);
959 
960  // Interpolate position linearly between vertex nodes
961  dummy_element.interpolated_x(s,x);
962  el2_pt->node_pt(j)->x(0)=x[0];
963  el2_pt->node_pt(j)->x(1)=x[1];
964  }
965 
966  // We now need to loop over nodes corner_2 to (corner_3-1) and copy the
967  // given information from Q1. We do not need to set up the (corner_3)'th
968  // node since it coincides with the centroid which has already been set up
969  for (unsigned j=corner_2;j<corner_3;j++)
970  {
971  // The nodes along the upper edge of Q2 go from corner_2 to corner_3-1
972  // while the nodes along the right edge of Q1 go from corner_1 to
973  // (corner_3-nnode_1d) in increments of nnode_1d
974  el2_pt->node_pt(j)=el1_pt->node_pt(corner_1+(j-corner_2)*nnode_1d);
975  }
976 
977  // Add the element pointers to Element_pt and then increment the counter
978  Element_pt[new_el_count]=el0_pt;
979  Element_pt[new_el_count+1]=el1_pt;
980  Element_pt[new_el_count+2]=el2_pt;
981  new_el_count+=3;
982 
983  // Assign the edges to the edge map
984  edge_nodes_map[edge0]=edge0_vector_pt;
985  edge_nodes_map[edge1]=edge1_vector_pt;
986  edge_nodes_map[edge2]=edge2_vector_pt;
987  }
988 
989  // Indicate that the look up scheme has been set up
990  Lookup_for_elements_next_boundary_is_setup=true;
991 
992  // Clean the dummy element nodes
993  for (unsigned j=0;j<4;j++)
994  {
995  // Delete the j-th dummy element node
996  delete dummy_element.node_pt(j);
997 
998  // Make it a null pointer
999  dummy_element.node_pt(j)=0;
1000  }
1001 
1002  }
1003 
1004 
1005  ////////////////////////////////////////////////////////////////////
1006  ////////////////////////////////////////////////////////////////////
1007  ////////////////////////////////////////////////////////////////////
1008 
1009 
1010 
1011 //======================================================================
1012 /// Adapt problem based on specified elemental error estimates
1013 //======================================================================
1014  template <class ELEMENT>
1016  adapt(const Vector<double>& elem_error)
1017  {
1018  // Call the adapt function from the TreeBasedRefineableMeshBase base class
1019  TreeBasedRefineableMeshBase::adapt(elem_error);
1020 
1021 #ifdef OOMPH_HAS_TRIANGLE_LIB
1022  // Move the nodes on the new boundary onto the old curvilinear
1023  // boundary. If the boundary is straight this will do precisely
1024  // nothing but will be somewhat inefficient
1025  this->snap_nodes_onto_geometric_objects();
1026 #endif
1027  } // End of adapt
1028 } // End of namespace oomph
1029 
1030 #endif // OOMPH_QUAD_FROM_TRIANGLE_MESH_TEMPLATE_CC
Triangle Mesh that is based on input files generated by the triangle mesh generator Triangle...
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
unsigned edge_boundary(const unsigned &e, const unsigned &i) const
Return the boundary id of the i-th edge in the e-th element: This is zero-based as in triangle...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
e
Definition: cfortran.h:575
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1287
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
Unstructured refineable QuadFromTriangleMesh.
static char t char * s
Definition: cfortran.h:572
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Edge class.
Definition: mesh.h:2401
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219