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