rectangle_with_moving_cylinder_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#ifndef OOMPH_RECTANGLE_WITH_MOVING_CYLINDER_MESH_TEMPLATE_HEADER
27#define OOMPH_RECTANGLE_WITH_MOVING_CYLINDER_MESH_TEMPLATE_HEADER
28
29#ifndef OOMPH_RECTANGLE_WITH_MOVING_CYLINDER_MESH_HEADER
30#error __FILE__ should only be included from rectangle_with_moving_cylinder_mesh.h.
31#endif // OOMPH_RECTANGLE_WITH_MOVING_CYLINDER_MESH_HEADER
32
33// PML mesh helpers
34#include "generic/pml_meshes.h"
35
36// Namespace extension
37namespace oomph
38{
39 //=============================================================================
40 /// Rectangular domain with circular whole
41 //=============================================================================
42 void RectangleWithHoleAndAnnularRegionDomain::
43 project_point_on_cylinder_to_annular_boundary(const unsigned& time,
44 const Vector<double>& xi,
45 Vector<double>& r)
46 {
47 // Lagrangian coordinate corresponds to theta=pi on the circle
48 Vector<double> xi_left(1, 4.0 * atan(1.0));
49
50 // Point on circle corresponding to theta=pi
51 Vector<double> point_left(2, 0.0);
52
53 // Get the position of the point on the circle
54 Cylinder_pt->position(time, xi_left, point_left);
55
56 // Lagrangian coordinate corresponds to theta=0 on the circle
57 Vector<double> xi_right(1, 0.0);
58
59 // Point on circle corresponding to theta=0
60 Vector<double> point_right(2, 0.0);
61
62 // Get the position of the point on the circle
63 Cylinder_pt->position(time, xi_right, point_right);
64
65 // Storage for the coordinates of the centre
66 Vector<double> point_centre(2, 0.0);
67
68 // Loop over the coordinates
69 for (unsigned i = 0; i < 2; i++)
70 {
71 // Calculate the i-th coordinate of the central point
72 point_centre[i] = 0.5 * (point_left[i] + point_right[i]);
73 }
74
75 // Calculate the x-coordinate of the projectd point
76 r[0] = point_centre[0] + Annular_region_radius * cos(xi[0]);
77
78 // Calculate the y-coordinate of the projectd point
79 r[1] = point_centre[1] + Annular_region_radius * sin(xi[0]);
80 } // End of project_point_on_cylinder_to_annular_boundary
81
82
83 /// Helper function that, given the Lagrangian coordinate, xi,
84 /// (associated with a point on the cylinder), returns the corresponding
85 /// point on the outer boundary of the annular region (where the inner
86 /// boundary is prescribed by the boundary of the cylinder)
87 void RectangleWithHoleAndAnnularRegionDomain::
88 project_point_on_cylinder_to_annular_boundary(const double& time,
89 const Vector<double>& xi,
90 Vector<double>& r)
91 {
92 // Lagrangian coordinate corresponds to theta=pi on the circle
93 Vector<double> xi_left(1, 4.0 * atan(1.0));
94
95 // Point on circle corresponding to theta=pi
96 Vector<double> point_left(2, 0.0);
97
98 // Get the position of the point on the circle
99 Cylinder_pt->position(time, xi_left, point_left);
100
101 // Lagrangian coordinate corresponds to theta=0 on the circle
102 Vector<double> xi_right(1, 0.0);
103
104 // Point on circle corresponding to theta=0
105 Vector<double> point_right(2, 0.0);
106
107 // Get the position of the point on the circle
108 Cylinder_pt->position(time, xi_right, point_right);
109
110 // Storage for the coordinates of the centre
111 Vector<double> point_centre(2, 0.0);
112
113 // Loop over the coordinates
114 for (unsigned i = 0; i < 2; i++)
115 {
116 // Calculate the i-th coordinate of the central point
117 point_centre[i] = 0.5 * (point_left[i] + point_right[i]);
118 }
119
120 // Calculate the x-coordinate of the projectd point
121 r[0] = point_centre[0] + Annular_region_radius * cos(xi[0]);
122
123 // Calculate the y-coordinate of the projectd point
124 r[1] = point_centre[1] + Annular_region_radius * sin(xi[0]);
125 } // End of project_point_on_cylinder_to_annular_boundary
126
127
128 /// Parametrisation of macro element boundaries: f(s) is the position
129 /// vector to macro-element m's boundary in the specified direction [N/S/E/W]
130 /// at the specified discrete time level (time=0: present; time>0: previous)
131 void RectangleWithHoleAndAnnularRegionDomain::macro_element_boundary(
132 const double& time,
133 const unsigned& m,
134 const unsigned& direction,
135 const Vector<double>& s,
136 Vector<double>& f)
137 {
138#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
139 // Create warning about time argument being moved to the front
140 std::string error_message_string =
141 "Order of function arguments has changed ";
142
143 // Finish the string off
144 error_message_string += "between versions 0.8 and 0.85";
145
146 // Output a warning
147 OomphLibWarning(
148 error_message_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
149#endif
150
151 // Lagrangian coordinate along surface of cylinder
152 Vector<double> xi(1, 0.0);
153
154 // Point on circle
155 Vector<double> point_on_circle(2, 0.0);
156
157 // Point on the outer boundary of the annular region
158 Vector<double> point_on_annular_ring(2, 0.0);
159
160 // Use the QuadTree enumeration for face directions
161 using namespace QuadTreeNames;
162
163 // Switch on the macro element
164 switch (m)
165 {
166 // Macro element 0 is immediately to the left of the cylinder, outside
167 // the inner annular region
168 case 0:
169
170 // Switch on the direction
171 switch (direction)
172 {
173 case N:
174 xi[0] = 3.0 * atan(1.0);
175 project_point_on_cylinder_to_annular_boundary(
176 time, xi, point_on_annular_ring);
177 linear_interpolate(Upper_mid_left, point_on_annular_ring, s[0], f);
178 break;
179
180 case S:
181 xi[0] = -3.0 * atan(1.0);
182 project_point_on_cylinder_to_annular_boundary(
183 time, xi, point_on_annular_ring);
184 linear_interpolate(Lower_mid_left, point_on_annular_ring, s[0], f);
185 break;
186
187 case W:
188 linear_interpolate(Lower_mid_left, Upper_mid_left, s[0], f);
189 break;
190
191 case E:
192 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
193 project_point_on_cylinder_to_annular_boundary(time, xi, f);
194 break;
195
196 default:
197
198 std::ostringstream error_stream;
199 error_stream << "Direction is incorrect: " << direction
200 << std::endl;
201 throw OomphLibError(error_stream.str(),
202 OOMPH_CURRENT_FUNCTION,
203 OOMPH_EXCEPTION_LOCATION);
204 }
205
206 break;
207
208 // Macro element 1 is immediately above the cylinder, outside
209 // the inner annular region
210 case 1:
211
212 switch (direction)
213 {
214 case N:
215 linear_interpolate(Upper_mid_left, Upper_mid_right, s[0], f);
216 break;
217
218 case S:
219 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
220 project_point_on_cylinder_to_annular_boundary(time, xi, f);
221 break;
222
223 case W:
224 xi[0] = 3.0 * atan(1.0);
225 project_point_on_cylinder_to_annular_boundary(
226 time, xi, point_on_annular_ring);
227 linear_interpolate(point_on_annular_ring, Upper_mid_left, s[0], f);
228 break;
229
230 case E:
231 xi[0] = 1.0 * atan(1.0);
232 project_point_on_cylinder_to_annular_boundary(
233 time, xi, point_on_annular_ring);
234 linear_interpolate(point_on_annular_ring, Upper_mid_right, s[0], f);
235 break;
236
237 default:
238
239 std::ostringstream error_stream;
240 error_stream << "Direction is incorrect: " << direction
241 << std::endl;
242 throw OomphLibError(error_stream.str(),
243 OOMPH_CURRENT_FUNCTION,
244 OOMPH_EXCEPTION_LOCATION);
245 }
246
247 break;
248
249 // Macro element 2 is immediately to the right of the cylinder, outside
250 // the inner annular region
251 case 2:
252
253 switch (direction)
254 {
255 case N:
256 xi[0] = 1.0 * atan(1.0);
257 project_point_on_cylinder_to_annular_boundary(
258 time, xi, point_on_annular_ring);
259 linear_interpolate(point_on_annular_ring, Upper_mid_right, s[0], f);
260 break;
261
262 case S:
263 xi[0] = -1.0 * atan(1.0);
264 project_point_on_cylinder_to_annular_boundary(
265 time, xi, point_on_annular_ring);
266 linear_interpolate(point_on_annular_ring, Lower_mid_right, s[0], f);
267 break;
268
269 case W:
270 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
271 project_point_on_cylinder_to_annular_boundary(time, xi, f);
272 break;
273
274 case E:
275 linear_interpolate(Lower_mid_right, Upper_mid_right, s[0], f);
276 break;
277
278 default:
279
280 std::ostringstream error_stream;
281 error_stream << "Direction is incorrect: " << direction
282 << std::endl;
283 throw OomphLibError(error_stream.str(),
284 OOMPH_CURRENT_FUNCTION,
285 OOMPH_EXCEPTION_LOCATION);
286 }
287
288 break;
289
290 // Macro element 3 is immediately below cylinder, outside
291 // the inner annular region
292 case 3:
293
294 switch (direction)
295 {
296 case N:
297 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
298 project_point_on_cylinder_to_annular_boundary(time, xi, f);
299 break;
300
301 case S:
302 linear_interpolate(Lower_mid_left, Lower_mid_right, s[0], f);
303 break;
304
305 case W:
306 xi[0] = -3.0 * atan(1.0);
307 project_point_on_cylinder_to_annular_boundary(
308 time, xi, point_on_annular_ring);
309 linear_interpolate(Lower_mid_left, point_on_annular_ring, s[0], f);
310 break;
311
312 case E:
313 xi[0] = -1.0 * atan(1.0);
314 project_point_on_cylinder_to_annular_boundary(
315 time, xi, point_on_annular_ring);
316 linear_interpolate(Lower_mid_right, point_on_annular_ring, s[0], f);
317 break;
318
319 default:
320
321 std::ostringstream error_stream;
322 error_stream << "Direction is incorrect: " << direction
323 << std::endl;
324 throw OomphLibError(error_stream.str(),
325 OOMPH_CURRENT_FUNCTION,
326 OOMPH_EXCEPTION_LOCATION);
327 }
328
329 break;
330
331 // Macro element 4, is immediately to the left of the cylinder
332 // lying in the inner annular region
333 case 4:
334
335 // Switch on the face of the m-th macro element
336 switch (direction)
337 {
338 case N:
339 xi[0] = 3.0 * atan(1.0);
340 Cylinder_pt->position(time, xi, point_on_circle);
341 project_point_on_cylinder_to_annular_boundary(
342 time, xi, point_on_annular_ring);
343 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
344 break;
345
346 case S:
347 xi[0] = -3.0 * atan(1.0);
348 Cylinder_pt->position(time, xi, point_on_circle);
349 project_point_on_cylinder_to_annular_boundary(
350 time, xi, point_on_annular_ring);
351 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
352 break;
353
354 case W:
355 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
356 project_point_on_cylinder_to_annular_boundary(time, xi, f);
357 break;
358
359 case E:
360 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
361 Cylinder_pt->position(time, xi, f);
362 break;
363
364 default:
365
366 std::ostringstream error_stream;
367 error_stream << "Direction is incorrect: " << direction
368 << std::endl;
369 throw OomphLibError(error_stream.str(),
370 OOMPH_CURRENT_FUNCTION,
371 OOMPH_EXCEPTION_LOCATION);
372 }
373
374 break;
375
376 // Macro element 5, is immediately above the cylinder in the
377 // inner annular region
378 case 5:
379
380 switch (direction)
381 {
382 case N:
383 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
384 project_point_on_cylinder_to_annular_boundary(time, xi, f);
385 break;
386
387 case S:
388 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
389 Cylinder_pt->position(time, xi, f);
390 break;
391
392 case W:
393 xi[0] = 3.0 * atan(1.0);
394 Cylinder_pt->position(time, xi, point_on_circle);
395 project_point_on_cylinder_to_annular_boundary(
396 time, xi, point_on_annular_ring);
397 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
398 break;
399
400 case E:
401 xi[0] = 1.0 * atan(1.0);
402 Cylinder_pt->position(time, xi, point_on_circle);
403 project_point_on_cylinder_to_annular_boundary(
404 time, xi, point_on_annular_ring);
405 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
406 break;
407
408 default:
409
410 std::ostringstream error_stream;
411 error_stream << "Direction is incorrect: " << direction
412 << std::endl;
413 throw OomphLibError(error_stream.str(),
414 OOMPH_CURRENT_FUNCTION,
415 OOMPH_EXCEPTION_LOCATION);
416 }
417
418 break;
419
420 // Macro element 6, is immediately to the right of the cylinder
421 // lying in the inner annular region
422 case 6:
423
424 switch (direction)
425 {
426 case N:
427 xi[0] = 1.0 * atan(1.0);
428 Cylinder_pt->position(time, xi, point_on_circle);
429 project_point_on_cylinder_to_annular_boundary(
430 time, xi, point_on_annular_ring);
431 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
432 break;
433
434 case S:
435 xi[0] = -1.0 * atan(1.0);
436 Cylinder_pt->position(time, xi, point_on_circle);
437 project_point_on_cylinder_to_annular_boundary(
438 time, xi, point_on_annular_ring);
439 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
440 break;
441
442 case W:
443 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
444 Cylinder_pt->position(time, xi, f);
445 break;
446
447 case E:
448 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
449 project_point_on_cylinder_to_annular_boundary(time, xi, f);
450 break;
451
452 default:
453
454 std::ostringstream error_stream;
455 error_stream << "Direction is incorrect: " << direction
456 << std::endl;
457 throw OomphLibError(error_stream.str(),
458 OOMPH_CURRENT_FUNCTION,
459 OOMPH_EXCEPTION_LOCATION);
460 }
461
462 break;
463
464 // Macro element 7, is immediately below the cylinder in the
465 // inner annular region
466 case 7:
467
468 switch (direction)
469 {
470 case N:
471 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
472 Cylinder_pt->position(time, xi, f);
473 break;
474
475 case S:
476 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
477 project_point_on_cylinder_to_annular_boundary(time, xi, f);
478 break;
479
480 case W:
481 xi[0] = -3.0 * atan(1.0);
482 Cylinder_pt->position(time, xi, point_on_circle);
483 project_point_on_cylinder_to_annular_boundary(
484 time, xi, point_on_annular_ring);
485 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
486 break;
487
488 case E:
489 xi[0] = -1.0 * atan(1.0);
490 Cylinder_pt->position(time, xi, point_on_circle);
491 project_point_on_cylinder_to_annular_boundary(
492 time, xi, point_on_annular_ring);
493 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
494 break;
495
496 default:
497
498 std::ostringstream error_stream;
499 error_stream << "Direction is incorrect: " << direction
500 << std::endl;
501 throw OomphLibError(error_stream.str(),
502 OOMPH_CURRENT_FUNCTION,
503 OOMPH_EXCEPTION_LOCATION);
504 }
505
506 break;
507
508 default:
509
510 std::ostringstream error_stream;
511 error_stream << "Wrong macro element number" << m << std::endl;
512 throw OomphLibError(
513 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
514 }
515 } // End of macro_element_boundary
516
517
518 /// Parametrisation of macro element boundaries: f(s) is the position
519 /// vector to macro-element m's boundary in the specified direction [N/S/E/W]
520 /// at the specified discrete time level (time=0: present; time>0: previous)
521 void RectangleWithHoleAndAnnularRegionDomain::macro_element_boundary(
522 const unsigned& time,
523 const unsigned& m,
524 const unsigned& direction,
525 const Vector<double>& s,
526 Vector<double>& f)
527 {
528#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
529 // Warn about time argument being moved to the front
530 OomphLibWarning(
531 "Order of function arguments has changed between versions 0.8 and 0.85",
532 "RectangleWithHoleAndAnnularRegionDomain::macro_element_boundary(...)",
533 OOMPH_EXCEPTION_LOCATION);
534#endif
535
536 // Lagrangian coordinate along surface of cylinder
537 Vector<double> xi(1, 0.0);
538
539 // Point on circle
540 Vector<double> point_on_circle(2, 0.0);
541
542 // Point on the outer boundary of the annular region
543 Vector<double> point_on_annular_ring(2, 0.0);
544
545 // Use the QuadTree enumeration for face directions
546 using namespace QuadTreeNames;
547
548 // Switch on the macro element
549 switch (m)
550 {
551 // Macro element 0 is immediately to the left of the cylinder, outside
552 // the inner annular region
553 case 0:
554
555 // Switch on the direction
556 switch (direction)
557 {
558 case N:
559 xi[0] = 3.0 * atan(1.0);
560 project_point_on_cylinder_to_annular_boundary(
561 time, xi, point_on_annular_ring);
562 linear_interpolate(Upper_mid_left, point_on_annular_ring, s[0], f);
563 break;
564
565 case S:
566 xi[0] = -3.0 * atan(1.0);
567 project_point_on_cylinder_to_annular_boundary(
568 time, xi, point_on_annular_ring);
569 linear_interpolate(Lower_mid_left, point_on_annular_ring, s[0], f);
570 break;
571
572 case W:
573 linear_interpolate(Lower_mid_left, Upper_mid_left, s[0], f);
574 break;
575
576 case E:
577 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
578 project_point_on_cylinder_to_annular_boundary(time, xi, f);
579 break;
580
581 default:
582
583 std::ostringstream error_stream;
584 error_stream << "Direction is incorrect: " << direction
585 << std::endl;
586 throw OomphLibError(error_stream.str(),
587 OOMPH_CURRENT_FUNCTION,
588 OOMPH_EXCEPTION_LOCATION);
589 }
590
591 break;
592
593 // Macro element 1 is immediately above the cylinder, outside
594 // the inner annular region
595 case 1:
596
597 switch (direction)
598 {
599 case N:
600 linear_interpolate(Upper_mid_left, Upper_mid_right, s[0], f);
601 break;
602
603 case S:
604 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
605 project_point_on_cylinder_to_annular_boundary(time, xi, f);
606 break;
607
608 case W:
609 xi[0] = 3.0 * atan(1.0);
610 project_point_on_cylinder_to_annular_boundary(
611 time, xi, point_on_annular_ring);
612 linear_interpolate(point_on_annular_ring, Upper_mid_left, s[0], f);
613 break;
614
615 case E:
616 xi[0] = 1.0 * atan(1.0);
617 project_point_on_cylinder_to_annular_boundary(
618 time, xi, point_on_annular_ring);
619 linear_interpolate(point_on_annular_ring, Upper_mid_right, s[0], f);
620 break;
621
622 default:
623
624 std::ostringstream error_stream;
625 error_stream << "Direction is incorrect: " << direction
626 << std::endl;
627 throw OomphLibError(error_stream.str(),
628 OOMPH_CURRENT_FUNCTION,
629 OOMPH_EXCEPTION_LOCATION);
630 }
631
632 break;
633
634 // Macro element 2 is immediately to the right of the cylinder, outside
635 // the inner annular region
636 case 2:
637
638 switch (direction)
639 {
640 case N:
641 xi[0] = 1.0 * atan(1.0);
642 project_point_on_cylinder_to_annular_boundary(
643 time, xi, point_on_annular_ring);
644 linear_interpolate(point_on_annular_ring, Upper_mid_right, s[0], f);
645 break;
646
647 case S:
648 xi[0] = -1.0 * atan(1.0);
649 project_point_on_cylinder_to_annular_boundary(
650 time, xi, point_on_annular_ring);
651 linear_interpolate(point_on_annular_ring, Lower_mid_right, s[0], f);
652 break;
653
654 case W:
655 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
656 project_point_on_cylinder_to_annular_boundary(time, xi, f);
657 break;
658
659 case E:
660 linear_interpolate(Lower_mid_right, Upper_mid_right, s[0], f);
661 break;
662
663 default:
664
665 std::ostringstream error_stream;
666 error_stream << "Direction is incorrect: " << direction
667 << std::endl;
668 throw OomphLibError(error_stream.str(),
669 OOMPH_CURRENT_FUNCTION,
670 OOMPH_EXCEPTION_LOCATION);
671 }
672
673 break;
674
675 // Macro element 3 is immediately below cylinder, outside
676 // the inner annular region
677 case 3:
678
679 switch (direction)
680 {
681 case N:
682 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
683 project_point_on_cylinder_to_annular_boundary(time, xi, f);
684 break;
685
686 case S:
687 linear_interpolate(Lower_mid_left, Lower_mid_right, s[0], f);
688 break;
689
690 case W:
691 xi[0] = -3.0 * atan(1.0);
692 project_point_on_cylinder_to_annular_boundary(
693 time, xi, point_on_annular_ring);
694 linear_interpolate(Lower_mid_left, point_on_annular_ring, s[0], f);
695 break;
696
697 case E:
698 xi[0] = -1.0 * atan(1.0);
699 project_point_on_cylinder_to_annular_boundary(
700 time, xi, point_on_annular_ring);
701 linear_interpolate(Lower_mid_right, point_on_annular_ring, s[0], f);
702 break;
703
704 default:
705
706 std::ostringstream error_stream;
707 error_stream << "Direction is incorrect: " << direction
708 << std::endl;
709 throw OomphLibError(error_stream.str(),
710 OOMPH_CURRENT_FUNCTION,
711 OOMPH_EXCEPTION_LOCATION);
712 }
713
714 break;
715
716 // Macro element 4, is immediately to the left of the cylinder
717 // lying in the inner annular region
718 case 4:
719
720 // Switch on the face of the m-th macro element
721 switch (direction)
722 {
723 case N:
724 xi[0] = 3.0 * atan(1.0);
725 Cylinder_pt->position(time, xi, point_on_circle);
726 project_point_on_cylinder_to_annular_boundary(
727 time, xi, point_on_annular_ring);
728 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
729 break;
730
731 case S:
732 xi[0] = -3.0 * atan(1.0);
733 Cylinder_pt->position(time, xi, point_on_circle);
734 project_point_on_cylinder_to_annular_boundary(
735 time, xi, point_on_annular_ring);
736 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
737 break;
738
739 case W:
740 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
741 project_point_on_cylinder_to_annular_boundary(time, xi, f);
742 break;
743
744 case E:
745 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
746 Cylinder_pt->position(time, xi, f);
747 break;
748
749 default:
750
751 std::ostringstream error_stream;
752 error_stream << "Direction is incorrect: " << direction
753 << std::endl;
754 throw OomphLibError(error_stream.str(),
755 OOMPH_CURRENT_FUNCTION,
756 OOMPH_EXCEPTION_LOCATION);
757 }
758
759 break;
760
761 // Macro element 5, is immediately above the cylinder in the
762 // inner annular region
763 case 5:
764
765 switch (direction)
766 {
767 case N:
768 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
769 project_point_on_cylinder_to_annular_boundary(time, xi, f);
770 break;
771
772 case S:
773 xi[0] = 3.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
774 Cylinder_pt->position(time, xi, f);
775 break;
776
777 case W:
778 xi[0] = 3.0 * atan(1.0);
779 Cylinder_pt->position(time, xi, point_on_circle);
780 project_point_on_cylinder_to_annular_boundary(
781 time, xi, point_on_annular_ring);
782 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
783 break;
784
785 case E:
786 xi[0] = 1.0 * atan(1.0);
787 Cylinder_pt->position(time, xi, point_on_circle);
788 project_point_on_cylinder_to_annular_boundary(
789 time, xi, point_on_annular_ring);
790 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
791 break;
792
793 default:
794
795 std::ostringstream error_stream;
796 error_stream << "Direction is incorrect: " << direction
797 << std::endl;
798 throw OomphLibError(error_stream.str(),
799 OOMPH_CURRENT_FUNCTION,
800 OOMPH_EXCEPTION_LOCATION);
801 }
802
803 break;
804
805 // Macro element 6, is immediately to the right of the cylinder
806 // lying in the inner annular region
807 case 6:
808
809 switch (direction)
810 {
811 case N:
812 xi[0] = 1.0 * atan(1.0);
813 Cylinder_pt->position(time, xi, point_on_circle);
814 project_point_on_cylinder_to_annular_boundary(
815 time, xi, point_on_annular_ring);
816 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
817 break;
818
819 case S:
820 xi[0] = -1.0 * atan(1.0);
821 Cylinder_pt->position(time, xi, point_on_circle);
822 project_point_on_cylinder_to_annular_boundary(
823 time, xi, point_on_annular_ring);
824 linear_interpolate(point_on_circle, point_on_annular_ring, s[0], f);
825 break;
826
827 case W:
828 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
829 Cylinder_pt->position(time, xi, f);
830 break;
831
832 case E:
833 xi[0] = -atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
834 project_point_on_cylinder_to_annular_boundary(time, xi, f);
835 break;
836
837 default:
838
839 std::ostringstream error_stream;
840 error_stream << "Direction is incorrect: " << direction
841 << std::endl;
842 throw OomphLibError(error_stream.str(),
843 OOMPH_CURRENT_FUNCTION,
844 OOMPH_EXCEPTION_LOCATION);
845 }
846
847 break;
848
849 // Macro element 7, is immediately below the cylinder in the
850 // inner annular region
851 case 7:
852
853 switch (direction)
854 {
855 case N:
856 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
857 Cylinder_pt->position(time, xi, f);
858 break;
859
860 case S:
861 xi[0] = -3.0 * atan(1.0) + 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
862 project_point_on_cylinder_to_annular_boundary(time, xi, f);
863 break;
864
865 case W:
866 xi[0] = -3.0 * atan(1.0);
867 Cylinder_pt->position(time, xi, point_on_circle);
868 project_point_on_cylinder_to_annular_boundary(
869 time, xi, point_on_annular_ring);
870 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
871 break;
872
873 case E:
874 xi[0] = -1.0 * atan(1.0);
875 Cylinder_pt->position(time, xi, point_on_circle);
876 project_point_on_cylinder_to_annular_boundary(
877 time, xi, point_on_annular_ring);
878 linear_interpolate(point_on_annular_ring, point_on_circle, s[0], f);
879 break;
880
881 default:
882
883 std::ostringstream error_stream;
884 error_stream << "Direction is incorrect: " << direction
885 << std::endl;
886 throw OomphLibError(error_stream.str(),
887 OOMPH_CURRENT_FUNCTION,
888 OOMPH_EXCEPTION_LOCATION);
889 }
890
891 break;
892
893 default:
894
895 std::ostringstream error_stream;
896 error_stream << "Wrong macro element number" << m << std::endl;
897 throw OomphLibError(
898 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
899 }
900 } // End of macro_element_boundary
901
902 ////////////////////////////////////////////////////////////////////////
903 ////////////////////////////////////////////////////////////////////////
904 ////////////////////////////////////////////////////////////////////////
905
906 //=============================================================================
907 /// Domain-based mesh for rectangular mesh with circular hole
908 //=============================================================================
909 /// Constructor: Pass pointer to geometric object that
910 /// represents the cylinder, the length and height of the domain.
911 /// The GeomObject must be parametrised such that
912 /// \f$\zeta \in [0,2\pi]\f$ sweeps around the circumference
913 /// in anticlockwise direction. Timestepper defaults to Steady
914 /// default timestepper.
915 template<class ELEMENT>
916 RectangleWithHoleAndAnnularRegionMesh<ELEMENT>::
917 RectangleWithHoleAndAnnularRegionMesh(GeomObject* cylinder_pt,
918 const double& annular_region_radius,
919 const double& length,
920 TimeStepper* time_stepper_pt)
921 {
922 // Create the domain
923 Domain_pt = new RectangleWithHoleAndAnnularRegionDomain(
924 cylinder_pt, annular_region_radius, length);
925
926 // Initialise the node counter
927 unsigned long node_count = 0;
928
929 // Vectors used to get data from domains
930 Vector<double> s(2), r(2);
931
932 // Setup temporary storage for the Nodes
933 Vector<Node*> tmp_node_pt;
934
935 // Number of macro elements in the domain
936 unsigned nmacro_element = Domain_pt->nmacro_element();
937
938 // Now blindly loop over the macro elements and associate a finite
939 // element with each
940 for (unsigned e = 0; e < nmacro_element; e++)
941 {
942 // Create the FiniteElement and add to the Element_pt Vector
943 Element_pt.push_back(new ELEMENT);
944
945 // Read out the number of linear points in the element
946 unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(e))->nnode_1d();
947
948 // Loop over nodes in the column
949 for (unsigned l1 = 0; l1 < np; l1++)
950 {
951 // Loop over the nodes in the row
952 for (unsigned l2 = 0; l2 < np; l2++)
953 {
954 // Allocate the memory for the node
955 tmp_node_pt.push_back(finite_element_pt(e)->construct_node(
956 l1 * np + l2, time_stepper_pt));
957
958 // Read out the position of the node from the macro element
959 s[0] = -1.0 + 2.0 * double(l2) / double(np - 1);
960 s[1] = -1.0 + 2.0 * double(l1) / double(np - 1);
961
962 // Use the macro element map from local coordinates to global
963 // coordinates
964 Domain_pt->macro_element_pt(e)->macro_map(s, r);
965
966 // Set the position of the node
967 tmp_node_pt[node_count]->x(0) = r[0];
968 tmp_node_pt[node_count]->x(1) = r[1];
969
970 // Increment the node number
971 node_count++;
972 }
973 } // for (unsigned l1=0;l1<np;l1++)
974 } // for (unsigned e=0;e<nmacro_element;e++)
975
976 //----------------------------------------------------------------
977 // Now the elements have been created, but there will be nodes in
978 // common, need to loop over the common edges and sort it, by
979 // reassigning pointers and the deleting excess nodes.
980 //----------------------------------------------------------------
981 // Read out the number of linear points in the element
982 unsigned np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
983
984 // Edge between Elements 0 and 1 (miss out the node which coincides
985 // with 4 elements; this will be dealt with separately later)
986 for (unsigned n = 1; n < np; n++)
987 {
988 // Set the nodes in element 1 to be the same as in element 0
989 finite_element_pt(1)->node_pt(n * np) =
990 finite_element_pt(0)->node_pt((np * np - 1) - n);
991
992 // Remove the nodes in element 1 from the temporary node list
993 delete tmp_node_pt[np * np + n * np];
994
995 // Make the corresponding pointer a null pointer
996 tmp_node_pt[np * np + n * np] = 0;
997 }
998
999 // Edge between Elements 0 and 3 (miss out the node which coincides
1000 // with 4 elements; this will be dealt with separately later)
1001 for (unsigned n = 0; n < np - 1; n++)
1002 {
1003 // Set the nodes in element 3 to be the same as in element 0
1004 finite_element_pt(3)->node_pt(n * np) = finite_element_pt(0)->node_pt(n);
1005
1006 // Remove the nodes in element 3 from the temporary node list
1007 delete tmp_node_pt[3 * np * np + n * np];
1008
1009 // Make the corresponding pointer a null pointer
1010 tmp_node_pt[3 * np * np + n * np] = 0;
1011 }
1012
1013 // Edge between Element 1 and 2 (miss out the node which coincides
1014 // with 4 elements; this will be dealt with separately later)
1015 for (unsigned n = 1; n < np; n++)
1016 {
1017 // Set the nodes in element 2 to be the same as in element 1
1018 finite_element_pt(2)->node_pt(np * (np - 1) + n) =
1019 finite_element_pt(1)->node_pt((np - 1) + n * np);
1020
1021 // Remove the nodes in element 2 from the temporary node list
1022 delete tmp_node_pt[2 * np * np + np * (np - 1) + n];
1023
1024 // Make the corresponding pointer a null pointer
1025 tmp_node_pt[2 * np * np + np * (np - 1) + n] = 0;
1026 }
1027
1028 // Edge between Element 3 and 2 (miss out the node which coincides
1029 // with 4 elements; this will be dealt with separately later)
1030 for (unsigned n = 1; n < np; n++)
1031 {
1032 // Set the nodes in element 2 to be the same as in element 3
1033 finite_element_pt(2)->node_pt(n) =
1034 finite_element_pt(3)->node_pt((np * np - 1) - n * np);
1035
1036 // Remove the nodes in element 2 from the temporary node list
1037 delete tmp_node_pt[2 * np * np + n];
1038
1039 // Make the corresponding pointer a null pointer
1040 tmp_node_pt[2 * np * np + n] = 0;
1041 }
1042
1043 // Edge between Elements 4 and 5 (miss out the node which coincides
1044 // with 4 elements; this will be dealt with separately later)
1045 for (unsigned n = 0; n < np - 1; n++)
1046 {
1047 // Set the nodes in element 1 to be the same as in element 0
1048 finite_element_pt(5)->node_pt(n * np) =
1049 finite_element_pt(4)->node_pt((np * np - 1) - n);
1050
1051 // Remove the nodes in element 1 from the temporary node list
1052 delete tmp_node_pt[5 * np * np + n * np];
1053
1054 // Make the corresponding pointer a null pointer
1055 tmp_node_pt[5 * np * np + n * np] = 0;
1056 }
1057
1058 // Edge between Elements 4 and 7 (miss out the node which coincides
1059 // with 4 elements; this will be dealt with separately later)
1060 for (unsigned n = 1; n < np; n++)
1061 {
1062 // Set the nodes in element 3 to be the same as in element 0
1063 finite_element_pt(7)->node_pt(n * np) = finite_element_pt(4)->node_pt(n);
1064
1065 // Remove the nodes in element 3 from the temporary node list
1066 delete tmp_node_pt[7 * np * np + n * np];
1067
1068 // Make the corresponding pointer a null pointer
1069 tmp_node_pt[7 * np * np + n * np] = 0;
1070 }
1071
1072 // Edge between Element 5 and 6 (miss out the node which coincides
1073 // with 4 elements; this will be dealt with separately later)
1074 for (unsigned n = 0; n < np - 1; n++)
1075 {
1076 // Set the nodes in element 2 to be the same as in element 1
1077 finite_element_pt(6)->node_pt(np * (np - 1) + n) =
1078 finite_element_pt(5)->node_pt((n + 1) * np - 1);
1079
1080 // Remove the nodes in element 2 from the temporary node list
1081 delete tmp_node_pt[6 * np * np + np * (np - 1) + n];
1082
1083 // Make the corresponding pointer a null pointer
1084 tmp_node_pt[6 * np * np + np * (np - 1) + n] = 0;
1085 }
1086
1087 // Edge between Element 7 and 6 (miss out the node which coincides
1088 // with 4 elements; this will be dealt with separately later)
1089 for (unsigned n = 0; n < np - 1; n++)
1090 {
1091 // Set the nodes in element 2 to be the same as in element 3
1092 finite_element_pt(6)->node_pt(n) =
1093 finite_element_pt(7)->node_pt((np * np - 1) - n * np);
1094
1095 // Remove the nodes in element 2 from the temporary node list
1096 delete tmp_node_pt[6 * np * np + n];
1097
1098 // Make the corresponding pointer a null pointer
1099 tmp_node_pt[6 * np * np + n] = 0;
1100 }
1101
1102 // Edge between Elements 0 and 4 (miss out the corner nodes as they have
1103 // coincide with 4 elements rather than just 2)
1104 for (unsigned n = 1; n < np - 1; n++)
1105 {
1106 // Set the nodes in element 1 to be the same as in element 0
1107 finite_element_pt(0)->node_pt((np - 1) + n * np) =
1108 finite_element_pt(4)->node_pt(n * np);
1109
1110 // Remove the nodes in element 1 from the temporary node list
1111 delete tmp_node_pt[(n + 1) * np - 1];
1112
1113 // Make the corresponding pointer a null pointer
1114 tmp_node_pt[(n + 1) * np - 1] = 0;
1115 }
1116
1117 // Edge between Elements 1 and 5 (miss out the corner nodes as they have
1118 // coincide with 4 elements rather than just 2)
1119 for (unsigned n = 1; n < np - 1; n++)
1120 {
1121 // Set the nodes in element 3 to be the same as in element 0
1122 finite_element_pt(1)->node_pt(n) =
1123 finite_element_pt(5)->node_pt(np * (np - 1) + n);
1124
1125 // Remove the nodes in element 1 from the temporary node list
1126 delete tmp_node_pt[np * np + n];
1127
1128 // Make the corresponding pointer a null pointer
1129 tmp_node_pt[np * np + n] = 0;
1130 }
1131
1132 // Edge between Element 2 and 6 (miss out the corner nodes as they have
1133 // coincide with 4 elements rather than just 2)
1134 for (unsigned n = 1; n < np - 1; n++)
1135 {
1136 // Set the nodes in element 2 to be the same as in element 1
1137 finite_element_pt(2)->node_pt(np * (np - 1) - n * np) =
1138 finite_element_pt(6)->node_pt((np * np - 1) - n * np);
1139
1140 // Remove the nodes in element 1 from the temporary node list
1141 delete tmp_node_pt[2 * np * np + np * (np - 1) - n * np];
1142
1143 // Make the corresponding pointer a null pointer
1144 tmp_node_pt[2 * np * np + np * (np - 1) - n * np] = 0;
1145 }
1146
1147 // Edge between Element 3 and 7 (miss out the corner nodes as they have
1148 // coincide with 4 elements rather than just 2)
1149 for (unsigned n = 1; n < np - 1; n++)
1150 {
1151 // Set the nodes in element 2 to be the same as in element 3
1152 finite_element_pt(3)->node_pt((np * np - 1) - n) =
1153 finite_element_pt(7)->node_pt((np - 1) - n);
1154
1155 // Remove the nodes in element 1 from the temporary node list
1156 delete tmp_node_pt[3 * np * np + (np * np - 1) - n];
1157
1158 // Make the corresponding pointer a null pointer
1159 tmp_node_pt[3 * np * np + (np * np - 1) - n] = 0;
1160 }
1161
1162 //--------------------------------------------------------------------
1163 // Now we'll deal with the corner nodes which lie in 4 elements rather
1164 // than just two elements. There are four of these to deal with.
1165 // Specifically:
1166 // Node 0: Lies in elements 0, 3, 4 and 7 (copy from element 0)
1167 // Node 1: Lies in elements 0, 1, 4 and 5 (copy from element 0)
1168 // Node 2: Lies in elements 1, 2, 5 and 6 (copy from element 1)
1169 // Node 3: Lies in elements 2, 3, 6 and 7 (copy from element 2)
1170 //--------------------------------------------------------------------
1171 //---------------
1172 // Corner node 0:
1173 //---------------
1174 // Set the corner node in element 3 to be the same as in element 0
1175 finite_element_pt(3)->node_pt(np * (np - 1)) =
1176 finite_element_pt(0)->node_pt(np - 1);
1177
1178 // Set the corner node in element 4 to be the same as in element 0
1179 finite_element_pt(4)->node_pt(0) = finite_element_pt(0)->node_pt(np - 1);
1180
1181 // Set the corner node in element 7 to be the same as in element 0
1182 finite_element_pt(7)->node_pt(0) = finite_element_pt(0)->node_pt(np - 1);
1183
1184 // Remove the overwritten node in element 3 from the temporary node list
1185 delete tmp_node_pt[3 * np * np + np * (np - 1)];
1186
1187 // Make the corresponding pointer a null pointer
1188 tmp_node_pt[3 * np * np + np * (np - 1)] = 0;
1189
1190 // Remove the overwritten node in element 4 from the temporary node list
1191 delete tmp_node_pt[4 * np * np];
1192
1193 // Make the corresponding pointer a null pointer
1194 tmp_node_pt[4 * np * np] = 0;
1195
1196 // Remove the overwritten node in element 7 from the temporary node list
1197 delete tmp_node_pt[7 * np * np];
1198
1199 // Make the corresponding pointer a null pointer
1200 tmp_node_pt[7 * np * np] = 0;
1201
1202 //---------------
1203 // Corner node 1:
1204 //---------------
1205 // Set the corner node in element 3 to be the same as in element 0
1206 finite_element_pt(1)->node_pt(0) =
1207 finite_element_pt(0)->node_pt(np * np - 1);
1208
1209 // Set the corner node in element 4 to be the same as in element 0
1210 finite_element_pt(4)->node_pt(np * (np - 1)) =
1211 finite_element_pt(0)->node_pt(np * np - 1);
1212
1213 // Set the corner node in element 7 to be the same as in element 0
1214 finite_element_pt(5)->node_pt(np * (np - 1)) =
1215 finite_element_pt(0)->node_pt(np * np - 1);
1216
1217 // Remove the overwritten node in element 1 from the temporary node list
1218 delete tmp_node_pt[np * np];
1219
1220 // Make the corresponding pointer a null pointer
1221 tmp_node_pt[np * np] = 0;
1222
1223 // Remove the overwritten node in element 4 from the temporary node list
1224 delete tmp_node_pt[4 * np * np + np * (np - 1)];
1225
1226 // Make the corresponding pointer a null pointer
1227 tmp_node_pt[4 * np * np + np * (np - 1)] = 0;
1228
1229 // Remove the overwritten node in element 5 from the temporary node list
1230 delete tmp_node_pt[5 * np * np + np * (np - 1)];
1231
1232 // Make the corresponding pointer a null pointer
1233 tmp_node_pt[5 * np * np + np * (np - 1)] = 0;
1234
1235 //---------------
1236 // Corner node 2:
1237 //---------------
1238 // Set the corner node in element 2 to be the same as in element 1
1239 finite_element_pt(2)->node_pt(np * (np - 1)) =
1240 finite_element_pt(1)->node_pt(np - 1);
1241
1242 // Set the corner node in element 5 to be the same as in element 1
1243 finite_element_pt(5)->node_pt(np * np - 1) =
1244 finite_element_pt(1)->node_pt(np - 1);
1245
1246 // Set the corner node in element 6 to be the same as in element 1
1247 finite_element_pt(6)->node_pt(np * np - 1) =
1248 finite_element_pt(1)->node_pt(np - 1);
1249
1250 // Remove the overwritten node in element 2 from the temporary node list
1251 delete tmp_node_pt[2 * np * np + np * (np - 1)];
1252
1253 // Make the corresponding pointer a null pointer
1254 tmp_node_pt[2 * np * np + np * (np - 1)] = 0;
1255
1256 // Remove the overwritten node in element 5 from the temporary node list
1257 delete tmp_node_pt[5 * np * np + np * np - 1];
1258
1259 // Make the corresponding pointer a null pointer
1260 tmp_node_pt[5 * np * np + np * np - 1] = 0;
1261
1262 // Remove the overwritten node in element 6 from the temporary node list
1263 delete tmp_node_pt[6 * np * np + np * np - 1];
1264
1265 // Make the corresponding pointer a null pointer
1266 tmp_node_pt[6 * np * np + np * np - 1] = 0;
1267
1268 //---------------
1269 // Corner node 3:
1270 //---------------
1271 // Set the corner node in element 3 to be the same as in element 2
1272 finite_element_pt(3)->node_pt(np * np - 1) =
1273 finite_element_pt(2)->node_pt(0);
1274
1275 // Set the corner node in element 4 to be the same as in element 2
1276 finite_element_pt(6)->node_pt(np - 1) = finite_element_pt(2)->node_pt(0);
1277
1278 // Set the corner node in element 7 to be the same as in element 2
1279 finite_element_pt(7)->node_pt(np - 1) = finite_element_pt(2)->node_pt(0);
1280
1281 // Remove the overwritten node in element 2 from the temporary node list
1282 delete tmp_node_pt[3 * np * np + np * np - 1];
1283
1284 // Make the corresponding pointer a null pointer
1285 tmp_node_pt[3 * np * np + np * np - 1] = 0;
1286
1287 // Remove the overwritten node in element 5 from the temporary node list
1288 delete tmp_node_pt[6 * np * np + np - 1];
1289
1290 // Make the corresponding pointer a null pointer
1291 tmp_node_pt[6 * np * np + np - 1] = 0;
1292
1293 // Remove the overwritten node in element 6 from the temporary node list
1294 delete tmp_node_pt[7 * np * np + np - 1];
1295
1296 // Make the corresponding pointer a null pointer
1297 tmp_node_pt[7 * np * np + np - 1] = 0;
1298
1299 //----------------------------------------------------------------------
1300 // Now all the nodes have been set up and the nodes which coincided with
1301 // each other have been deleted and the correct pointers have been set.
1302 // All that remains is to copy the nodes in the temporary storage into
1303 // the global storage vector Node_pt.
1304 //----------------------------------------------------------------------
1305 // Now set the actual true nodes
1306 for (unsigned long n = 0; n < node_count; n++)
1307 {
1308 // If it's not a null pointer
1309 if (tmp_node_pt[n] != 0)
1310 {
1311 // Add it to the global Node list
1312 Node_pt.push_back(tmp_node_pt[n]);
1313 }
1314 } // for (unsigned long n=0;n<node_count;n++)
1315
1316 //-------------------------------------------------------------
1317 // Now sort out which boundary each node lies on. The outer box
1318 // has boundary numbering 0 for the South face, 1 for the East
1319 // face, 2 for the North face and 3 for the West face. Finally,
1320 // the cylinder boundary corresponds to boundary 4.
1321 //-------------------------------------------------------------
1322 // Finally set the nodes on the boundaries
1323 set_nboundary(5);
1324
1325 // Loop over the nodes along an edge
1326 for (unsigned n = 0; n < np; n++)
1327 {
1328 // West face
1329 Node* nod_pt = finite_element_pt(0)->node_pt(n * np);
1330 convert_to_boundary_node(nod_pt);
1331 add_boundary_node(3, nod_pt);
1332
1333 // North face
1334 nod_pt = finite_element_pt(1)->node_pt(np * (np - 1) + n);
1335 convert_to_boundary_node(nod_pt);
1336 add_boundary_node(2, nod_pt);
1337
1338 // East face
1339 nod_pt = finite_element_pt(2)->node_pt(n * np + np - 1);
1340 convert_to_boundary_node(nod_pt);
1341 add_boundary_node(1, nod_pt);
1342
1343 // South face
1344 nod_pt = finite_element_pt(3)->node_pt(n);
1345 convert_to_boundary_node(nod_pt);
1346 add_boundary_node(0, nod_pt);
1347 }
1348
1349 // Loop over all but one node on an edge
1350 for (unsigned n = 0; n < np - 1; n++)
1351 {
1352 // First part of hole
1353 Node* nod_pt = finite_element_pt(4)->node_pt((np - 1) + n * np);
1354 convert_to_boundary_node(nod_pt);
1355 add_boundary_node(4, nod_pt);
1356
1357 // Next part of hole
1358 nod_pt = finite_element_pt(5)->node_pt(n);
1359 convert_to_boundary_node(nod_pt);
1360 add_boundary_node(4, nod_pt);
1361
1362 // Next part of hole
1363 nod_pt = finite_element_pt(6)->node_pt(np * (np - 1) - n * np);
1364 convert_to_boundary_node(nod_pt);
1365 add_boundary_node(4, nod_pt);
1366
1367 // Final part of hole boundary
1368 nod_pt = finite_element_pt(7)->node_pt((np * np - 1) - n);
1369 convert_to_boundary_node(nod_pt);
1370 add_boundary_node(4, nod_pt);
1371 }
1372
1373#ifdef PARANOID
1374 // Make sure there are no duplicate nodes (i.e. two or more nodes that
1375 // occupy the same Eulerian position)
1376 if (check_for_repeated_nodes())
1377 {
1378 // Throw a warning
1379 OomphLibWarning("The mesh contains repeated nodes!\n",
1380 OOMPH_CURRENT_FUNCTION,
1381 OOMPH_EXCEPTION_LOCATION);
1382 }
1383#endif
1384 } // End of RectangleWithHoleAndAnnularRegionMesh
1385
1386 ////////////////////////////////////////////////////////////////////////
1387 ////////////////////////////////////////////////////////////////////////
1388 ////////////////////////////////////////////////////////////////////////
1389
1390 //=============================================================================
1391 /// Constructor. Pass pointer to geometric object that
1392 /// represents the cylinder; hierher the length and height of the domain.
1393 /// The GeomObject must be parametrised such that
1394 /// \f$\zeta \in [0,2\pi]\f$ sweeps around the circumference
1395 /// in anticlockwise direction. Timestepper defaults to Steady
1396 /// default timestepper.
1397 //=============================================================================
1398 template<class ELEMENT>
1399 RefineableQuadMeshWithMovingCylinder<ELEMENT>::
1400 RefineableQuadMeshWithMovingCylinder(GeomObject* cylinder_pt,
1401 const double& annular_region_radius,
1402 const double& length_of_central_box,
1403 const double& x_left,
1404 const double& x_right,
1405 const double& height,
1406 TimeStepper* time_stepper_pt)
1407 {
1408 // Do you want a coarse problem?
1409 Coarse_problem = false; // true; // hierher
1410
1411 // Vector for pointers to all the temporary meshes
1412 Vector<Mesh*> all_temp_mesh_pt;
1413
1414 // Build central mesh
1415 Central_mesh_pt =
1416 new RefineableRectangleWithHoleAndAnnularRegionMesh<ELEMENT>(
1417 cylinder_pt,
1418 annular_region_radius,
1419 length_of_central_box,
1420 time_stepper_pt);
1421
1422 // Add this mesh to the list of temproary meshes
1423 all_temp_mesh_pt.push_back(Central_mesh_pt);
1424
1425 // Bulk mesh left boundary id
1426 unsigned int left_boundary_id = 3;
1427
1428 // Bulk mesh top boundary id
1429 unsigned int top_boundary_id = 2;
1430
1431 // Bulk mesh right boundary id
1432 unsigned int right_boundary_id = 1;
1433
1434 // Bulk mesh bottom boundary id
1435 unsigned int bottom_boundary_id = 0;
1436
1437 // Build the surrounding meshes
1438
1439 // Right mesh
1440 double l_right = x_right - 0.5 * length_of_central_box;
1441 unsigned n_right =
1442 std::max(unsigned(1), unsigned(l_right / length_of_central_box));
1443 if (Coarse_problem)
1444 {
1445 n_right = 1;
1446 }
1447
1448 Mesh* surrounding_right_mesh_pt =
1449 TwoDimensionalPMLHelper::create_right_pml_mesh<PMLLayerElement<ELEMENT>>(
1450 Central_mesh_pt, right_boundary_id, n_right, l_right, time_stepper_pt);
1451
1452 all_temp_mesh_pt.push_back(surrounding_right_mesh_pt);
1453
1454 // Remove nodes from top and bottom boundaries (0 and 2)
1455 for (unsigned ibound_rm = 0; ibound_rm <= 2; ibound_rm += 2)
1456 {
1457 unsigned num_nod = surrounding_right_mesh_pt->nboundary_node(ibound_rm);
1458 for (unsigned inod = 0; inod < num_nod; inod++)
1459 {
1460 Node* nod_pt =
1461 surrounding_right_mesh_pt->boundary_node_pt(ibound_rm, inod);
1462 if (nod_pt->is_on_boundary(ibound_rm))
1463 {
1464 nod_pt->remove_from_boundary(ibound_rm);
1465 }
1466 }
1467 }
1468
1469 // Top mesh
1470 double l_top = 0.5 * height - 0.5 * length_of_central_box;
1471 unsigned n_top =
1472 std::max(unsigned(1), unsigned(l_top / length_of_central_box));
1473 if (Coarse_problem)
1474 {
1475 n_top = 1;
1476 }
1477 Mesh* surrounding_top_mesh_pt =
1478 TwoDimensionalPMLHelper::create_top_pml_mesh<PMLLayerElement<ELEMENT>>(
1479 Central_mesh_pt, top_boundary_id, n_top, l_top, time_stepper_pt);
1480 all_temp_mesh_pt.push_back(surrounding_top_mesh_pt);
1481
1482 // Remove nodes from right and left boundaries (1 and 3)
1483 for (unsigned ibound_rm = 1; ibound_rm <= 3; ibound_rm += 2)
1484 {
1485 unsigned num_nod = surrounding_top_mesh_pt->nboundary_node(ibound_rm);
1486 for (unsigned inod = 0; inod < num_nod; inod++)
1487 {
1488 Node* nod_pt =
1489 surrounding_top_mesh_pt->boundary_node_pt(ibound_rm, inod);
1490 if (nod_pt->is_on_boundary(ibound_rm))
1491 {
1492 nod_pt->remove_from_boundary(ibound_rm);
1493 }
1494 }
1495 }
1496
1497 // Left mesh
1498 double l_left = -0.5 * length_of_central_box - x_left;
1499 unsigned n_left =
1500 std::max(unsigned(1), unsigned(l_left / length_of_central_box));
1501 if (Coarse_problem)
1502 {
1503 n_left = 1;
1504 }
1505 Mesh* surrounding_left_mesh_pt =
1506 TwoDimensionalPMLHelper::create_left_pml_mesh<PMLLayerElement<ELEMENT>>(
1507 Central_mesh_pt, left_boundary_id, n_left, l_left, time_stepper_pt);
1508 all_temp_mesh_pt.push_back(surrounding_left_mesh_pt);
1509
1510 // Remove nodes from top and bottom boundaries (0 and 2)
1511 for (unsigned ibound_rm = 0; ibound_rm <= 2; ibound_rm += 2)
1512 {
1513 unsigned num_nod = surrounding_left_mesh_pt->nboundary_node(ibound_rm);
1514 for (unsigned inod = 0; inod < num_nod; inod++)
1515 {
1516 Node* nod_pt =
1517 surrounding_left_mesh_pt->boundary_node_pt(ibound_rm, inod);
1518 if (nod_pt->is_on_boundary(ibound_rm))
1519 {
1520 nod_pt->remove_from_boundary(ibound_rm);
1521 }
1522 }
1523 }
1524
1525 // Bottom mesh
1526 double l_bottom = 0.5 * height - 0.5 * length_of_central_box;
1527 unsigned n_bottom =
1528 std::max(unsigned(1), unsigned(l_bottom / length_of_central_box));
1529 if (Coarse_problem)
1530 {
1531 n_bottom = 1;
1532 }
1533 Mesh* surrounding_bottom_mesh_pt =
1534 TwoDimensionalPMLHelper::create_bottom_pml_mesh<PMLLayerElement<ELEMENT>>(
1535 Central_mesh_pt,
1536 bottom_boundary_id,
1537 n_bottom,
1538 l_bottom,
1539 time_stepper_pt);
1540 all_temp_mesh_pt.push_back(surrounding_bottom_mesh_pt);
1541
1542 // Remove nodes from right and left boundaries (1 and 3)
1543 for (unsigned ibound_rm = 1; ibound_rm <= 3; ibound_rm += 2)
1544 {
1545 unsigned num_nod = surrounding_bottom_mesh_pt->nboundary_node(ibound_rm);
1546 for (unsigned inod = 0; inod < num_nod; inod++)
1547 {
1548 Node* nod_pt =
1549 surrounding_bottom_mesh_pt->boundary_node_pt(ibound_rm, inod);
1550 if (nod_pt->is_on_boundary(ibound_rm))
1551 {
1552 nod_pt->remove_from_boundary(ibound_rm);
1553 }
1554 }
1555 }
1556
1557 // Build corner Surrounding meshes
1558 Mesh* surrounding_top_right_mesh_pt =
1559 TwoDimensionalPMLHelper::create_top_right_pml_mesh<
1560 PMLLayerElement<ELEMENT>>(surrounding_right_mesh_pt,
1561 surrounding_top_mesh_pt,
1562 Central_mesh_pt,
1563 right_boundary_id,
1564 time_stepper_pt);
1565 all_temp_mesh_pt.push_back(surrounding_top_right_mesh_pt);
1566
1567 Mesh* surrounding_bottom_right_mesh_pt =
1568 TwoDimensionalPMLHelper::create_bottom_right_pml_mesh<
1569 PMLLayerElement<ELEMENT>>(surrounding_right_mesh_pt,
1570 surrounding_bottom_mesh_pt,
1571 Central_mesh_pt,
1572 right_boundary_id,
1573 time_stepper_pt);
1574 all_temp_mesh_pt.push_back(surrounding_bottom_right_mesh_pt);
1575
1576 Mesh* surrounding_top_left_mesh_pt =
1577 TwoDimensionalPMLHelper::create_top_left_pml_mesh<
1578 PMLLayerElement<ELEMENT>>(surrounding_left_mesh_pt,
1579 surrounding_top_mesh_pt,
1580 Central_mesh_pt,
1581 left_boundary_id,
1582 time_stepper_pt);
1583 all_temp_mesh_pt.push_back(surrounding_top_left_mesh_pt);
1584
1585 Mesh* surrounding_bottom_left_mesh_pt =
1586 TwoDimensionalPMLHelper::create_bottom_left_pml_mesh<
1587 PMLLayerElement<ELEMENT>>(surrounding_left_mesh_pt,
1588 surrounding_bottom_mesh_pt,
1589 Central_mesh_pt,
1590 left_boundary_id,
1591 time_stepper_pt);
1592 all_temp_mesh_pt.push_back(surrounding_bottom_left_mesh_pt);
1593
1594 // Copy all elements across
1595 unsigned n_mesh = all_temp_mesh_pt.size();
1596 unsigned nel = 0;
1597 unsigned nnod = 0;
1598 for (unsigned i = 0; i < n_mesh; i++)
1599 {
1600 nel += all_temp_mesh_pt[i]->nelement();
1601 nnod += all_temp_mesh_pt[i]->nnode();
1602 }
1603 this->Element_pt.resize(nel);
1604 this->Node_pt.resize(nnod);
1605
1606 unsigned count_el = 0;
1607 unsigned count_nod = 0;
1608 for (unsigned i = 0; i < n_mesh; i++)
1609 {
1610 unsigned nel = all_temp_mesh_pt[i]->nelement();
1611 for (unsigned e = 0; e < nel; e++)
1612 {
1613 this->Element_pt[count_el] = all_temp_mesh_pt[i]->element_pt(e);
1614 count_el++;
1615 }
1616 unsigned nnod = all_temp_mesh_pt[i]->nnode();
1617 for (unsigned j = 0; j < nnod; j++)
1618 {
1619 this->Node_pt[count_nod] = all_temp_mesh_pt[i]->node_pt(j);
1620 count_nod++;
1621 }
1622 }
1623
1624 // Remove nodes on outer boundaries (0 to 3) of central mesh
1625 // from their respective boundaries (otherwise they get added
1626 // to the mesh boundaries during adaptation)
1627 for (unsigned b = 0; b < 4; b++)
1628 {
1629 unsigned nnod = Central_mesh_pt->nboundary_node(b);
1630 for (unsigned j = 0; j < nnod; j++)
1631 {
1632 Node* nod_pt = Central_mesh_pt->boundary_node_pt(b, j);
1633 nod_pt->remove_from_boundary(b);
1634 }
1635 }
1636
1637 // Setup boundaries:
1638 //==================
1639 this->set_nboundary(5);
1640
1641 Vector<Mesh*> pml_mesh_pt;
1642
1643 // Top meshes:
1644 //------------
1645 {
1646 // Upper boundary in top mesh (same number in this mesh)
1647 unsigned ibound = 2;
1648
1649 // Loop over relevant meshes
1650 Vector<Mesh*> all_mesh_pt;
1651 all_mesh_pt.push_back(surrounding_top_left_mesh_pt);
1652 all_mesh_pt.push_back(surrounding_top_mesh_pt);
1653 all_mesh_pt.push_back(surrounding_top_right_mesh_pt);
1654 unsigned n_mesh = all_mesh_pt.size();
1655 for (unsigned m = 0; m < n_mesh; m++)
1656 {
1657 pml_mesh_pt.push_back(all_mesh_pt[m]);
1658 unsigned num_nod = all_mesh_pt[m]->nboundary_node(ibound);
1659 for (unsigned inod = 0; inod < num_nod; inod++)
1660 {
1661 Node* nod_pt = all_mesh_pt[m]->boundary_node_pt(ibound, inod);
1662 this->add_boundary_node(ibound, nod_pt);
1663 }
1664 }
1665 }
1666
1667 // Left meshes:
1668 //------------
1669 {
1670 // Left boundary (same number in this mesh)
1671 unsigned ibound = 3;
1672
1673 pml_mesh_pt.push_back(surrounding_left_mesh_pt);
1674
1675 // Loop over relevant meshes
1676 Vector<Mesh*> all_mesh_pt;
1677 all_mesh_pt.push_back(surrounding_top_left_mesh_pt);
1678 all_mesh_pt.push_back(surrounding_left_mesh_pt);
1679 all_mesh_pt.push_back(surrounding_bottom_left_mesh_pt);
1680 unsigned n_mesh = all_mesh_pt.size();
1681 for (unsigned m = 0; m < n_mesh; m++)
1682 {
1683 unsigned num_nod = all_mesh_pt[m]->nboundary_node(ibound);
1684 for (unsigned inod = 0; inod < num_nod; inod++)
1685 {
1686 Node* nod_pt = all_mesh_pt[m]->boundary_node_pt(ibound, inod);
1687 this->add_boundary_node(ibound, nod_pt);
1688 }
1689 }
1690 }
1691
1692 // Bottom meshes:
1693 //---------------
1694 {
1695 // Bottom boundary (same number in this mesh)
1696 unsigned ibound = 0;
1697
1698 // Loop over relevant meshes
1699 Vector<Mesh*> all_mesh_pt;
1700 all_mesh_pt.push_back(surrounding_bottom_left_mesh_pt);
1701 all_mesh_pt.push_back(surrounding_bottom_mesh_pt);
1702 all_mesh_pt.push_back(surrounding_bottom_right_mesh_pt);
1703 unsigned n_mesh = all_mesh_pt.size();
1704 for (unsigned m = 0; m < n_mesh; m++)
1705 {
1706 pml_mesh_pt.push_back(all_mesh_pt[m]);
1707 unsigned num_nod = all_mesh_pt[m]->nboundary_node(ibound);
1708 for (unsigned inod = 0; inod < num_nod; inod++)
1709 {
1710 Node* nod_pt = all_mesh_pt[m]->boundary_node_pt(ibound, inod);
1711 this->add_boundary_node(ibound, nod_pt);
1712 }
1713 }
1714 }
1715
1716 // Right meshes:
1717 //--------------
1718 {
1719 // Right boundary (same number in this mesh)
1720 unsigned ibound = 1;
1721
1722 pml_mesh_pt.push_back(surrounding_right_mesh_pt);
1723
1724 // Loop over relevant meshes
1725 Vector<Mesh*> all_mesh_pt;
1726 all_mesh_pt.push_back(surrounding_bottom_right_mesh_pt);
1727 all_mesh_pt.push_back(surrounding_right_mesh_pt);
1728 all_mesh_pt.push_back(surrounding_top_right_mesh_pt);
1729 unsigned n_mesh = all_mesh_pt.size();
1730 for (unsigned m = 0; m < n_mesh; m++)
1731 {
1732 unsigned num_nod = all_mesh_pt[m]->nboundary_node(ibound);
1733 for (unsigned inod = 0; inod < num_nod; inod++)
1734 {
1735 Node* nod_pt = all_mesh_pt[m]->boundary_node_pt(ibound, inod);
1736 this->add_boundary_node(ibound, nod_pt);
1737 }
1738 }
1739 }
1740
1741 // No slip on boundary of inner mesh
1742 //----------------------------------
1743 {
1744 // Inner boundary
1745 unsigned ibound = 4;
1746
1747 // Pin both velocities
1748 unsigned num_nod = Central_mesh_pt->nboundary_node(ibound);
1749 for (unsigned inod = 0; inod < num_nod; inod++)
1750 {
1751 Node* nod_pt = Central_mesh_pt->boundary_node_pt(ibound, inod);
1752 this->add_boundary_node(ibound, nod_pt);
1753 }
1754 }
1755
1756 // Setup quadtree forest for mesh refinement
1757 this->setup_quadtree_forest();
1758
1759 // Setup boundary element lookup schemes
1760 this->setup_boundary_element_info();
1761
1762 // Cleanup. NOTE: Can't delete Central_mesh_pt as it's responsible for
1763 // deleting Domain_pt but
1764 unsigned n_pml_mesh = pml_mesh_pt.size();
1765 for (unsigned j = 0; j < n_pml_mesh; j++)
1766 {
1767 pml_mesh_pt[j]->flush_element_and_node_storage();
1768 delete pml_mesh_pt[j];
1769 pml_mesh_pt[j] = 0;
1770 }
1771 Central_mesh_pt->flush_element_and_node_storage();
1772 } // End of RefineableQuadMeshWithMovingCylinder
1773} // namespace oomph
1774#endif
e
e
Definition cfortran.h:571
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::FiniteElement::size
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
oomph::FiniteElement::Node_pt
Node ** Node_pt
Storage for pointers to the nodes in the element.
Definition elements.h:1323
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
oomph::FiniteElement::nnode_1d
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
oomph::GeomObject
A geometric object is an object that provides a parametrised description of its shape via the functio...
Definition geom_objects.h:101
oomph::GeomObject::position
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition geom_objects.h:192
oomph::Mesh
A general mesh class.
Definition mesh.h:67
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
oomph::Node::resize
void resize(const unsigned &n_value)
Resize the number of equations.
Definition nodes.cc:2167
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition oomph_definitions.h:274
oomph::RectangleWithHoleAndAnnularRegionDomain
Rectangular domain with circular whole DRAIG: This looks like a redefinition of the RectangleWithHole...
Definition rectangle_with_moving_cylinder_mesh.h:54
oomph::RectangleWithHoleAndAnnularRegionDomain::Upper_mid_left
Vector< double > Upper_mid_left
Where the "radial" line from circle meets upper boundary on left.
Definition rectangle_with_moving_cylinder_mesh.h:192
oomph::RectangleWithHoleAndAnnularRegionDomain::linear_interpolate
void linear_interpolate(const Vector< double > &left, const Vector< double > &right, const double &s, Vector< double > &f)
Helper function to interpolate linearly between the "right" and "left" points; .
Definition rectangle_with_moving_cylinder_mesh.h:138
oomph::RectangleWithHoleAndAnnularRegionDomain::macro_element_boundary
void macro_element_boundary(const double &time, const unsigned &m, const unsigned &direction, const Vector< double > &s, Vector< double > &f)
Parametrisation of macro element boundaries: f(s) is the position vector to macro-element m's boundar...
Definition rectangle_with_moving_cylinder_mesh.template.cc:131
oomph::RectangleWithHoleAndAnnularRegionDomain::Cylinder_pt
GeomObject * Cylinder_pt
Pointer to geometric object that represents the central cylinder.
Definition rectangle_with_moving_cylinder_mesh.h:198
oomph::RectangleWithHoleAndAnnularRegionDomain::project_point_on_cylinder_to_annular_boundary
void project_point_on_cylinder_to_annular_boundary(const unsigned &time, const Vector< double > &xi, Vector< double > &r)
Helper function that, given the Lagrangian coordinate, xi, (associated with a point on the cylinder),...
Definition rectangle_with_moving_cylinder_mesh.template.cc:43
oomph::RectangleWithHoleAndAnnularRegionDomain::Lower_mid_left
Vector< double > Lower_mid_left
Where the "radial" line from circle meets lower boundary on left.
Definition rectangle_with_moving_cylinder_mesh.h:180
oomph::RectangleWithHoleAndAnnularRegionDomain::Lower_mid_right
Vector< double > Lower_mid_right
Where the "radial" line from circle meets lower boundary on right.
Definition rectangle_with_moving_cylinder_mesh.h:183
oomph::RectangleWithHoleAndAnnularRegionDomain::Annular_region_radius
double Annular_region_radius
The radius of the outer boundary of the annular region whose inner boundary is described by Cylinder_...
Definition rectangle_with_moving_cylinder_mesh.h:204
oomph::RectangleWithHoleAndAnnularRegionDomain::Upper_mid_right
Vector< double > Upper_mid_right
Where the "radial" line from circle meets upper boundary on right.
Definition rectangle_with_moving_cylinder_mesh.h:195
oomph::RectangleWithHoleAndAnnularRegionMesh::RectangleWithHoleAndAnnularRegionMesh
RectangleWithHoleAndAnnularRegionMesh(GeomObject *cylinder_pt, const double &annular_region_radius, const double &length, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that represents the cylinder, the length and height of ...
Definition rectangle_with_moving_cylinder_mesh.template.cc:917
oomph::RefineableQuadMeshWithMovingCylinder::RefineableQuadMeshWithMovingCylinder
RefineableQuadMeshWithMovingCylinder(GeomObject *cylinder_pt, const double &annular_region_radius, const double &length_of_central_box, const double &x_left, const double &x_right, const double &height, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. Pass pointer to geometric object that represents the cylinder; hierher the length and he...
Definition rectangle_with_moving_cylinder_mesh.template.cc:1400
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TAdvectionDiffusionReactionElement::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Definition Tadvection_diffusion_reaction_elements.h:70
oomph::TimeStepper
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition timesteppers.h:231
oomph::TwoDimensionalPMLHelper::create_bottom_left_pml_mesh
Mesh * create_bottom_left_pml_mesh(Mesh *pml_left_mesh_pt, Mesh *pml_bottom_mesh_pt, Mesh *bulk_mesh_pt, const unsigned &left_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML bottom left corner mesh, aligned with the existing PML meshes
Definition pml_meshes.h:2366
oomph::TwoDimensionalPMLHelper::create_top_left_pml_mesh
Mesh * create_top_left_pml_mesh(Mesh *pml_left_mesh_pt, Mesh *pml_top_mesh_pt, Mesh *bulk_mesh_pt, const unsigned &left_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML top left corner mesh, aligned with the existing PML meshes
Definition pml_meshes.h:2233
oomph::TwoDimensionalPMLHelper::create_top_right_pml_mesh
Mesh * create_top_right_pml_mesh(Mesh *pml_right_mesh_pt, Mesh *pml_top_mesh_pt, Mesh *bulk_mesh_pt, const unsigned &right_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML top right corner mesh, aligned with the existing PML meshes
Definition pml_meshes.h:1967
oomph::TwoDimensionalPMLHelper::create_bottom_right_pml_mesh
Mesh * create_bottom_right_pml_mesh(Mesh *pml_right_mesh_pt, Mesh *pml_bottom_mesh_pt, Mesh *bulk_mesh_pt, const unsigned &right_boundary_id, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML bottom right corner mesh, aligned with the existing PML meshes
Definition pml_meshes.h:2102
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30
pml_meshes.h