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