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-2023 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
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//====================================================================
27
28
29 // PML mesh helpers
30 #include "../generic/pml_meshes.h"
31
32 // Namespace extension
33 namespace oomph
34 {
35  //=============================================================================
36  /// Rectangular domain with circular whole
37  //=============================================================================
40  const Vector<double>& xi,
41  Vector<double>& r)
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,
86  Vector<double>& r)
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,
132  Vector<double>& f)
133  {
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
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,
522  Vector<double>& f)
523  {
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
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>
915  const double& length,
916  TimeStepper* time_stepper_pt)
917  {
918  // Create the domain
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);
1328
1329  // North face
1330  nod_pt = finite_element_pt(1)->node_pt(np * (np - 1) + n);
1331  convert_to_boundary_node(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);
1338
1339  // South face
1340  nod_pt = finite_element_pt(3)->node_pt(n);
1341  convert_to_boundary_node(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);
1352
1353  // Next part of hole
1354  nod_pt = finite_element_pt(5)->node_pt(n);
1355  convert_to_boundary_node(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);
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);
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>
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,
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);
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);
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);
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);
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);
1753  }
1754  }
1755
1756  // Setup quadtree forest for mesh refinement
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
unsigned nmacro_element()
Number of macro elements in domain.
Definition: domain.h:123
/////////////////////////////////////////////////////////////////////
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
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.
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.
RectangleWithHoleAndAnnularRegionDomain(GeomObject *cylinder_pt, const double &annular_region_radius, const double &length)
Constructor. Pass pointer to geometric object that represents the cylinder, the length of the (square...
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_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
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_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_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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...