full_circle_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-2025 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
27#ifndef OOMPH_FULL_CIRCLE_MESH_TEMPLATE_HEADER
28#define OOMPH_FULL_CIRCLE_MESH_TEMPLATE_HEADER
29
30#ifndef OOMPH_FULL_CIRCLE_MESH_HEADER
31#error __FILE__ should only be included from full_circle_mesh.h.
32#endif // OOMPH_FULL_CIRCLE_MESH_HEADER
33
34namespace oomph
35{
36 //====================================================================
37 /// Constructor for deformable quarter tube mesh class.
38 /// The domain is specified by the GeomObject that
39 /// identifies the entire volume.
40 //====================================================================
41 template<class ELEMENT>
42 FullCircleMesh<ELEMENT>::FullCircleMesh(GeomObject* area_pt,
43 const Vector<double>& theta_positions,
44 const Vector<double>& radius_box,
45 TimeStepper* time_stepper_pt)
46 : Area_pt(area_pt)
47 {
48// Check that the vectors are the correct sizes.
49#ifdef PARANOID
50 if (radius_box.size() != 4 || theta_positions.size() != 4)
51 {
52 std::string err =
53 "This mesh is hard coded to only work for the case when there are 5 "
54 "elements: the central square and 4 surrounding elements, but you gave "
55 "vectors inconsistent with this.";
56 throw OomphLibError(
57 err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
58 }
59#endif
60
61 // Mesh can only be built with 2D Qelements.
62 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
63
64 // Build macro element-based domain
65 Domain_pt = new FullCircleDomain(area_pt, theta_positions, radius_box);
66
67 // Set the number of boundaries
68 set_nboundary(1);
69
70 // We have only bothered to parametrise the only boundary (boundary 0)
71 set_boundary_coordinate_exists(0);
72
73 // Allocate the store for the elements
74 const unsigned nelem = 5;
75 Element_pt.resize(nelem);
76
77 // Create dummy element so we can determine the number of nodes
78 ELEMENT* dummy_el_pt = new ELEMENT;
79
80 // Read out the number of linear points in the element
81 unsigned n_p = dummy_el_pt->nnode_1d();
82
83 // Kill the element
84 delete dummy_el_pt;
85
86 // Can now allocate the store for the nodes
87 unsigned nnodes_total = (n_p * n_p + 4 * (n_p - 1) * (n_p - 1));
88 Node_pt.resize(nnodes_total);
89
90 Vector<double> s(2);
91 Vector<double> r(2);
92
93 // Storage for the intrinsic boundary coordinate
94 Vector<double> zeta(1);
95
96 // Loop over elements and create all nodes
97 for (unsigned ielem = 0; ielem < nelem; ielem++)
98 {
99 // Create element
100 Element_pt[ielem] = new ELEMENT;
101
102 // Loop over rows in y/s_1-direction
103 for (unsigned i1 = 0; i1 < n_p; i1++)
104 {
105 // Loop over rows in x/s_0-direction
106 for (unsigned i0 = 0; i0 < n_p; i0++)
107 {
108 // Local node number
109 unsigned jnod_local = i0 + i1 * n_p;
110
111 // Create the node
112 Node* node_pt = finite_element_pt(ielem)->construct_node(
113 jnod_local, time_stepper_pt);
114
115 // Set the position of the node from macro element mapping
116 s[0] = -1.0 + 2.0 * double(i0) / double(n_p - 1);
117 s[1] = -1.0 + 2.0 * double(i1) / double(n_p - 1);
118 Domain_pt->macro_element_pt(ielem)->macro_map(s, r);
119
120 node_pt->x(0) = r[0];
121 node_pt->x(1) = r[1];
122 }
123 }
124 }
125
126 // Initialise number of global nodes
127 unsigned node_count = 0;
128
129 // Tolerance for node killing:
130 double node_kill_tol = 1.0e-12;
131
132 // Check for error in node killing
133 bool stopit = false;
134
135 // Define pine
136 const double pi = MathematicalConstants::Pi;
137
138 // Loop over elements
139 for (unsigned ielem = 0; ielem < nelem; ielem++)
140 {
141 // Which macro element?
142 switch (ielem)
143 {
144 // Macro element 0: Central box, create all the nodes
145 //----------------------------------------------------
146 case 0:
147
148 // Loop over rows in y/s_1-direction
149 for (unsigned i1 = 0; i1 < n_p; i1++)
150 {
151 // Loop over rows in x/s_0-direction
152 for (unsigned i0 = 0; i0 < n_p; i0++)
153 {
154 // Local node number
155 unsigned jnod_local = i0 + i1 * n_p;
156
157 // Add the node to the mesh
158 Node_pt[node_count] =
159 finite_element_pt(ielem)->node_pt(jnod_local);
160
161 // Increment node counter
162 node_count++;
163 }
164 }
165
166 break;
167
168 // Macro element 1: Bottom
169 //---------------------------------
170 case 1:
171
172 // Loop over rows in y/s_1-direction
173 for (unsigned i1 = 0; i1 < n_p; i1++)
174 {
175 // Loop over rows in x/s_0-direction
176 for (unsigned i0 = 0; i0 < n_p; i0++)
177 {
178 // Local node number
179 unsigned jnod_local = i0 + i1 * n_p;
180
181 // Has the node been killed?
182 bool killed = false;
183
184 // Duplicate node: kill and set pointer to central element
185 if (i1 == (n_p - 1))
186 {
187 // Neighbour element
188 unsigned ielem_neigh = ielem - 1;
189
190 // Node in neighbour element
191 unsigned i0_neigh = i0;
192 unsigned i1_neigh = 0;
193 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
194
195 // Check:
196 for (unsigned i = 0; i < 2; i++)
197 {
198 double error = std::fabs(
199 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
200 finite_element_pt(ielem_neigh)
201 ->node_pt(jnod_local_neigh)
202 ->x(i));
203 if (error > node_kill_tol)
204 {
205 oomph_info << "Error in node killing for i " << i << " "
206 << error << std::endl;
207 stopit = true;
208 }
209 }
210
211 // Kill node
212 delete finite_element_pt(ielem)->node_pt(jnod_local);
213 killed = true;
214
215 // Set pointer to neighbour:
216 finite_element_pt(ielem)->node_pt(jnod_local) =
217 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
218 }
219
220 // No duplicate node: Copy across to mesh
221 if (!killed)
222 {
223 Node_pt[node_count] =
224 finite_element_pt(ielem)->node_pt(jnod_local);
225
226 // Set boundaries:
227
228 // Bottom: outer boundary
229 if (i1 == 0)
230 {
231 this->convert_to_boundary_node(Node_pt[node_count]);
232 add_boundary_node(0, Node_pt[node_count]);
233
234 // Get azimuthal boundary coordinate
235 zeta[0] = theta_positions[0] +
236 double(i1) / double(n_p - 1) * 0.5 *
237 (theta_positions[1] - theta_positions[0]);
238
239 Node_pt[node_count]->set_coordinates_on_boundary(0, zeta);
240 }
241
242 // Increment node counter
243 node_count++;
244 }
245 }
246 } // End of loop over nodes
247
248 break;
249
250 // Macro element 2: Right element
251 //--------------------------------
252 case 2:
253
254 // Loop over rows in y/s_1-direction
255 for (unsigned i1 = 0; i1 < n_p; i1++)
256 {
257 // Loop over rows in x/s_0-direction
258 for (unsigned i0 = 0; i0 < n_p; i0++)
259 {
260 // Local node number
261 unsigned jnod_local = i0 + i1 * n_p;
262
263 // Has the node been killed?
264 bool killed = false;
265
266 // Duplicate node: kill and set pointer to previous element
267 if (i1 == 0)
268 {
269 // Neighbour element
270 unsigned ielem_neigh = ielem - 1;
271
272 // Node in neighbour element
273 unsigned i0_neigh = n_p - 1;
274 unsigned i1_neigh = n_p - 1 - i0;
275
276 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
277
278 // Check:
279 for (unsigned i = 0; i < 2; i++)
280 {
281 double error = std::fabs(
282 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
283 finite_element_pt(ielem_neigh)
284 ->node_pt(jnod_local_neigh)
285 ->x(i));
286 if (error > node_kill_tol)
287 {
288 oomph_info << "Error in node killing for i " << i << " "
289 << error << std::endl;
290 stopit = true;
291 }
292 }
293
294 // Kill node
295 delete finite_element_pt(ielem)->node_pt(jnod_local);
296 killed = true;
297
298 // Set pointer to neighbour:
299 finite_element_pt(ielem)->node_pt(jnod_local) =
300 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
301 }
302
303 // Duplicate node: kill and set pointer to central element
304 if ((i0 == 0) && (i1 != 0))
305 {
306 // Neighbour element
307 unsigned ielem_neigh = ielem - 2;
308
309 // Node in neighbour element
310 unsigned i0_neigh = n_p - 1;
311 unsigned i1_neigh = i1;
312 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
313
314 // Check:
315 for (unsigned i = 0; i < 2; i++)
316 {
317 double error = std::fabs(
318 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
319 finite_element_pt(ielem_neigh)
320 ->node_pt(jnod_local_neigh)
321 ->x(i));
322 if (error > node_kill_tol)
323 {
324 oomph_info << "Error in node killing for i " << i << " "
325 << error << std::endl;
326 stopit = true;
327 }
328 }
329
330 // Kill node
331 delete finite_element_pt(ielem)->node_pt(jnod_local);
332 killed = true;
333
334 // Set pointer to neighbour:
335 finite_element_pt(ielem)->node_pt(jnod_local) =
336 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
337 }
338
339 // No duplicate node: Copy across to mesh
340 if (!killed)
341 {
342 Node_pt[node_count] =
343 finite_element_pt(ielem)->node_pt(jnod_local);
344
345 // Set boundaries:
346
347 // FullCircle boundary
348 if (i0 == n_p - 1)
349 {
350 this->convert_to_boundary_node(Node_pt[node_count]);
351 add_boundary_node(0, Node_pt[node_count]);
352
353 // Get azimuthal boundary coordinate
354 zeta[0] = theta_positions[1] +
355 double(i1) / double(n_p - 1) * 0.5 *
356 (theta_positions[2] - theta_positions[1]);
357
358 Node_pt[node_count]->set_coordinates_on_boundary(0, zeta);
359 }
360
361 // Increment node counter
362 node_count++;
363 }
364 }
365 }
366
367 break;
368
369 // Macro element 3: Top element
370 //--------------------------------
371 case 3:
372
373 // Loop over rows in y/s_1-direction
374 for (unsigned i1 = 0; i1 < n_p; i1++)
375 {
376 // Loop over rows in x/s_0-direction
377 for (unsigned i0 = 0; i0 < n_p; i0++)
378 {
379 // Local node number
380 unsigned jnod_local = i0 + i1 * n_p;
381
382 // Has the node been killed?
383 bool killed = false;
384
385 // Duplicate node: kill and set pointer to previous element
386 if (i0 == n_p - 1)
387 {
388 // Neighbour element
389 unsigned ielem_neigh = ielem - 1;
390
391 // Node in neighbour element
392 unsigned i0_neigh = i1;
393 unsigned i1_neigh = n_p - 1;
394 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
395
396 // Check:
397 for (unsigned i = 0; i < 2; i++)
398 {
399 double error = std::fabs(
400 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
401 finite_element_pt(ielem_neigh)
402 ->node_pt(jnod_local_neigh)
403 ->x(i));
404 if (error > node_kill_tol)
405 {
406 oomph_info << "Error in node killing for i " << i << " "
407 << error << std::endl;
408 stopit = true;
409 }
410 }
411
412 // Kill node
413 delete finite_element_pt(ielem)->node_pt(jnod_local);
414 killed = true;
415
416 // Set pointer to neighbour:
417 finite_element_pt(ielem)->node_pt(jnod_local) =
418 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
419 }
420
421 // Duplicate node: kill and set pointer to central element
422 if ((i1 == 0) && (i0 != n_p - 1))
423 {
424 // Neighbour element
425 unsigned ielem_neigh = ielem - 3;
426
427 // Node in neighbour element
428 unsigned i0_neigh = i0;
429 unsigned i1_neigh = n_p - 1;
430 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
431
432 // Check:
433 for (unsigned i = 0; i < 2; i++)
434 {
435 double error = std::fabs(
436 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
437 finite_element_pt(ielem_neigh)
438 ->node_pt(jnod_local_neigh)
439 ->x(i));
440 if (error > node_kill_tol)
441 {
442 oomph_info << "Error in node killing for i " << i << " "
443 << error << std::endl;
444 stopit = true;
445 }
446 }
447
448 // Kill node
449 delete finite_element_pt(ielem)->node_pt(jnod_local);
450 killed = true;
451
452 // Set pointer to neighbour:
453 finite_element_pt(ielem)->node_pt(jnod_local) =
454 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
455 }
456
457 // No duplicate node: Copy across to mesh
458 if (!killed)
459 {
460 Node_pt[node_count] =
461 finite_element_pt(ielem)->node_pt(jnod_local);
462
463 // Set boundaries:
464
465 // FullCircle boundary
466 if (i1 == n_p - 1)
467 {
468 this->convert_to_boundary_node(Node_pt[node_count]);
469 add_boundary_node(0, Node_pt[node_count]);
470
471 // Get azimuthal boundary coordinate
472 zeta[0] = theta_positions[3] +
473 double(i0) / double(n_p - 1) * 0.5 *
474 (theta_positions[2] - theta_positions[3]);
475
476 Node_pt[node_count]->set_coordinates_on_boundary(0, zeta);
477 }
478
479 // Increment node counter
480 node_count++;
481 }
482 }
483 }
484 break;
485
486 // Macro element 4: Left element
487 //--------------------------------
488 case 4:
489
490 // Loop over rows in y/s_1-direction
491 for (unsigned i1 = 0; i1 < n_p; i1++)
492 {
493 // Loop over rows in x/s_0-direction
494 for (unsigned i0 = 0; i0 < n_p; i0++)
495 {
496 // Local node number
497 unsigned jnod_local = i0 + i1 * n_p;
498
499 // Has the node been killed?
500 bool killed = false;
501
502 // Duplicate node: kill and set pointer to previous element
503 if (i1 == n_p - 1)
504 {
505 // Neighbour element
506 unsigned ielem_neigh = ielem - 1;
507
508 // Node in neighbour element
509 unsigned i0_neigh = 0;
510 unsigned i1_neigh = n_p - 1 - i0;
511 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
512
513 // Check:
514 for (unsigned i = 0; i < 2; i++)
515 {
516 double error = std::fabs(
517 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
518 finite_element_pt(ielem_neigh)
519 ->node_pt(jnod_local_neigh)
520 ->x(i));
521 if (error > node_kill_tol)
522 {
523 oomph_info << "Error in node killing for i " << i << " "
524 << error << std::endl;
525 stopit = true;
526 }
527 }
528
529 // Kill node
530 delete finite_element_pt(ielem)->node_pt(jnod_local);
531 killed = true;
532
533 // Set pointer to neighbour:
534 finite_element_pt(ielem)->node_pt(jnod_local) =
535 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
536 }
537
538 // Duplicate node: kill and set pointer to central element
539 if ((i0 == n_p - 1) && (i1 != n_p - 1))
540 {
541 // Neighbour element
542 unsigned ielem_neigh = ielem - 4;
543
544 // Node in neighbour element
545 unsigned i0_neigh = 0;
546 unsigned i1_neigh = i1;
547 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
548
549 // Check:
550 for (unsigned i = 0; i < 2; i++)
551 {
552 double error = std::fabs(
553 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
554 finite_element_pt(ielem_neigh)
555 ->node_pt(jnod_local_neigh)
556 ->x(i));
557 if (error > node_kill_tol)
558 {
559 oomph_info << "Error in node killing for i " << i << " "
560 << error << std::endl;
561 stopit = true;
562 }
563 }
564
565 // Kill node
566 delete finite_element_pt(ielem)->node_pt(jnod_local);
567 killed = true;
568
569 // Set pointer to neighbour:
570 finite_element_pt(ielem)->node_pt(jnod_local) =
571 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
572 }
573
574 // Duplicate node: kill and set pointer to other ring element
575 if ((i1 == 0) && (i0 != n_p - 1))
576 {
577 // Neighbour element
578 unsigned ielem_neigh = ielem - 3;
579
580 // Node in neighbour element
581 unsigned i0_neigh = 0;
582 unsigned i1_neigh = i0;
583 unsigned jnod_local_neigh = i0_neigh + i1_neigh * n_p;
584
585 // Check:
586 for (unsigned i = 0; i < 2; i++)
587 {
588 double error = std::fabs(
589 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
590 finite_element_pt(ielem_neigh)
591 ->node_pt(jnod_local_neigh)
592 ->x(i));
593 if (error > node_kill_tol)
594 {
595 oomph_info << "Error in node killing for i " << i << " "
596 << error << std::endl;
597 stopit = true;
598 }
599 }
600
601 // Kill node
602 delete finite_element_pt(ielem)->node_pt(jnod_local);
603 killed = true;
604
605 // Set pointer to neighbour:
606 finite_element_pt(ielem)->node_pt(jnod_local) =
607 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
608 }
609
610 // No duplicate node: Copy across to mesh
611 if (!killed)
612 {
613 Node_pt[node_count] =
614 finite_element_pt(ielem)->node_pt(jnod_local);
615
616 // Set boundaries:
617
618 // FullCircle boundary
619 if (i0 == 0)
620 {
621 this->convert_to_boundary_node(Node_pt[node_count]);
622 add_boundary_node(0, Node_pt[node_count]);
623
624 // Get azimuthal boundary coordinate
625 zeta[0] =
626 theta_positions[0] + 2.0 * pi +
627 double(i1) / double(n_p - 1) * 0.5 *
628 (theta_positions[3] - theta_positions[0] + 2.0 * pi);
629
630 Node_pt[node_count]->set_coordinates_on_boundary(0, zeta);
631 }
632
633 // Increment node counter
634 node_count++;
635 }
636 }
637 }
638 break;
639 }
640 }
641
642 // Terminate if there's been an error
643 if (stopit)
644 {
645 std::ostringstream error_stream;
646 error_stream << "Error in killing nodes\n"
647 << "The most probable cause is that the domain is not\n"
648 << "compatible with the mesh.\n"
649 << "For the FullCircleMesh, the domain must be\n"
650 << "topologically consistent with a quarter tube with a\n"
651 << "non-curved centreline.\n";
652 throw OomphLibError(
653 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
654 }
655
656 // Setup boundary element lookup schemes
657 setup_boundary_element_info();
658 }
659
660} // namespace oomph
661#endif
oomph::FullCircleDomain
Topologically circular domain, e.g. a tube cross section. The entire domain must be defined by a Geom...
Definition full_circle_domain.h:69
oomph::FullCircleMesh::Domain_pt
FullCircleDomain * Domain_pt
Pointer to domain.
Definition full_circle_mesh.h:88
oomph::FullCircleMesh::area_pt
GeomObject *& area_pt()
Access function to GeomObject representing wall.
Definition full_circle_mesh.h:69
oomph::FullCircleMesh::FullCircleMesh
FullCircleMesh(GeomObject *wall_pt, const Vector< double > &theta_positions, const Vector< double > &radius_box, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the area; values of theta at which divid...
Definition full_circle_mesh.template.cc:42
oomph::MacroElementNodeUpdateCollapsibleChannelMesh
Collapsible channel mesh with MacroElement-based node update. The collapsible segment is represented ...
Definition collapsible_channel_mesh.h:236