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-2022 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
33using namespace std;
34using namespace oomph;
35
36
37/// /////////////////////////////////////////////////////////////////
38/// /////////////////////////////////////////////////////////////////
39/// /////////////////////////////////////////////////////////////////
40
41
42namespace 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
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
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
Definition: mesh.h:2692
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
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:1365
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
void build_from_scaffold(TriangleScaffoldMesh *tmp_mesh_pt, TimeStepper *time_stepper_pt, const bool &use_attributes)
Build the quad mesh from the given scaffold mesh.
void adapt(const Vector< double > &elem_error)
Overload the adapt function (to ensure nodes are snapped to the boundary)
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
void adapt(const Vector< double > &elemental_error)
Adapt mesh: Refine elements whose error is lager than err_max and (try to) unrefine those whose error...
Triangle Mesh that is based on input files generated by the triangle mesh generator Triangle.
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...