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