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>
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(
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
66
67 // Set the number of boundaries
69
70 // We have only bothered to parametrise the only boundary (boundary 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
92
93 // Storage for the intrinsic boundary coordinate
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
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);
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
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(
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
213 killed = true;
214
215 // Set pointer to neighbour:
218 }
219
220 // No duplicate node: Copy across to mesh
221 if (!killed)
222 {
225
226 // Set boundaries:
227
228 // Bottom: outer boundary
229 if (i1 == 0)
230 {
231 this->convert_to_boundary_node(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 *
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(
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
296 killed = true;
297
298 // Set pointer to neighbour:
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(
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
332 killed = true;
333
334 // Set pointer to neighbour:
337 }
338
339 // No duplicate node: Copy across to mesh
340 if (!killed)
341 {
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]);
352
353 // Get azimuthal boundary coordinate
354 zeta[0] = theta_positions[1] +
355 double(i1) / double(n_p - 1) * 0.5 *
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(
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
414 killed = true;
415
416 // Set pointer to neighbour:
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(
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
450 killed = true;
451
452 // Set pointer to neighbour:
455 }
456
457 // No duplicate node: Copy across to mesh
458 if (!killed)
459 {
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]);
470
471 // Get azimuthal boundary coordinate
472 zeta[0] = theta_positions[3] +
473 double(i0) / double(n_p - 1) * 0.5 *
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(
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
531 killed = true;
532
533 // Set pointer to neighbour:
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(
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
567 killed = true;
568
569 // Set pointer to neighbour:
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(
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
603 killed = true;
604
605 // Set pointer to neighbour:
608 }
609
610 // No duplicate node: Copy across to mesh
611 if (!killed)
612 {
615
616 // Set boundaries:
617
618 // FullCircle boundary
619 if (i0 == 0)
620 {
621 this->convert_to_boundary_node(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(
654 }
655
656 // Setup boundary element lookup schemes
658 }
659
660} // namespace oomph
661#endif
static char t char * s
Definition cfortran.h:568
cstr elem_len * i
Definition cfortran.h:603
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition domain.h:116
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition elements.h:2513
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition elements.h:2222
Topologically circular domain, e.g. a tube cross section. The entire domain must be defined by a Geom...
FullCircleDomain * Domain_pt
Pointer to domain.
GeomObject *& area_pt()
Access function to GeomObject representing wall.
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...
A geometric object is an object that provides a parametrised description of its shape via the functio...
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition mesh.h:183
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition mesh.h:477
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition mesh.h:509
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition mesh.h:440
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this,...
Definition mesh.cc:2590
void set_boundary_coordinate_exists(const unsigned &i)
Set boundary coordinate on the i-th boundary to be existing.
Definition mesh.h:584
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition mesh.h:186
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
void setup_boundary_element_info()
Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to ...
Definition quad_mesh.h:73
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
const double Pi
50 digits from maple
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...