pml_meshes.h
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#ifndef OOMPH_PML_MESH_HEADER
27#define OOMPH_PML_MESH_HEADER
28
29// Config header generated by autoconfig
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34#include "mesh.h"
35#include "../meshes/rectangular_quadmesh.template.h"
36#include "../meshes/rectangular_quadmesh.template.cc"
37
38
39namespace oomph
40{
41 //=======================================================================
42 /// General definition of policy class defining the elements to
43 /// be used in the actual PML layers. Has to be instantiated for
44 /// each specific "bulk" PML element type.
45 //=======================================================================
46 template<class ELEMENT>
48 {
49 };
50
51 /// /////////////////////////////////////////////////////////////////////
52 /// /////////////////////////////////////////////////////////////////////
53 /// /////////////////////////////////////////////////////////////////////
54
55 //==============================================================
56 /// Base class for elements with pml capabilities
57 //==============================================================
58 template<unsigned DIM>
60 {
61 public:
62 /// Constructor
64 : Pml_is_enabled(false),
65 Pml_direction_active(DIM, false),
66 Pml_inner_boundary(DIM, 0.0),
67 Pml_outer_boundary(DIM, 0.0)
68 {
69 }
70
71 /// Virtual destructor
72 virtual ~PMLElementBase() {}
73
74 /// Disable pml. Ensures the PML-ification in all directions
75 /// has been deactivated
77 {
78 // Disable the PML-ification
79 Pml_is_enabled = false;
80
81 // Loop over the entries in Pml_direction_active and deactivate the
82 // PML in this direction
83 for (unsigned direction = 0; direction < DIM; direction++)
84 {
85 // Disable the PML in this direction
86 Pml_direction_active[direction] = false;
87 }
88 } // End of disable_pml
89
90 /// Enable pml. Specify the coordinate direction along which
91 /// pml boundary is constant, as well as the coordinate
92 /// along the dimension for the interface between the physical and
93 /// artificial domains and the coordinate for the outer boundary. All of
94 /// these are used to adjust the perfectly matched layer mechanism. Needs to
95 /// be called separately for each pml-ified direction (if needed -- e.g. in
96 /// corner elements)
97 void enable_pml(const int& direction,
98 const double& interface_border_value,
99 const double& outer_domain_border_value)
100 {
101 Pml_is_enabled = true;
102 Pml_direction_active[direction] = true;
103 Pml_inner_boundary[direction] = interface_border_value;
104 Pml_outer_boundary[direction] = outer_domain_border_value;
105 }
106
107 /// Pure virtual function in which we have to specify the
108 /// values to be pinned (and set to zero) on the outer edge of
109 /// the pml layer. This is usually all of the nodal values
110 /// (values 0 and 1 (real and imag part) for Helmholtz;
111 /// values 0,1,2 and 3 (real and imag part of x- and y-displacement
112 /// for 2D time-harmonic linear elasticity; etc.). Vector
113 /// must be resized internally!
115 Vector<unsigned>& values_to_pin) = 0;
116
117 protected:
118 /// Boolean indicating if element is used in pml mode
120
121 /// Coordinate direction along which pml boundary is constant;
122 /// alternatively: coordinate direction in which coordinate stretching
123 /// is performed.
124 std::vector<bool> Pml_direction_active;
125
126 /// Coordinate of inner pml boundary
127 /// (Storage is provided for any coordinate
128 /// direction; only the entries for "active" directions is used.)
130
131 /// Coordinate of outer pml boundary
132 /// (Storage is provided for any coordinate
133 /// direction; only the entries for "active" directions is used.)
135 };
136
137 /// ///////////////////////////////////////////////////////////////////
138 /// ///////////////////////////////////////////////////////////////////
139 /// ///////////////////////////////////////////////////////////////////
140
141 //===================================================================
142 /// All helper routines for 2D bulk boundary mesh usage in order to
143 /// generate PML meshes aligned to the main mesh
144 //===================================================================
145 namespace TwoDimensionalPMLHelper
146 {
147 /// helper function for sorting the right boundary nodes
148 extern bool sorter_right_boundary(Node* nod_i_pt, Node* nod_j_pt);
149
150 /// helper function for sorting the top boundary nodes
151 extern bool sorter_top_boundary(Node* nod_i_pt, Node* nod_j_pt);
152
153 /// helper function for sorting the left boundary nodes
154 extern bool sorter_left_boundary(Node* nod_i_pt, Node* nod_j_pt);
155
156 /// helper function for sorting the bottom boundary nodes
157 extern bool sorter_bottom_boundary(Node* nod_i_pt, Node* nod_j_pt);
158
159 } // namespace TwoDimensionalPMLHelper
160
161 /// ///////////////////////////////////////////////////////////////
162 /// ///////////////////////////////////////////////////////////////
163 /// ///////////////////////////////////////////////////////////////
164
165 //====================================================================
166 /// PML mesh base class. Contains a pure virtual locate_zeta function
167 /// to be uploaded in PMLQuadMesh and PMLBrickMesh (once the code for
168 /// it has been written)
169 //====================================================================
171 {
172 public:
173 /// Pure virtual function to provide an optimised version of
174 /// the locate_zeta function for PML meshes. As the PML meshes are
175 /// rectangular (in 2D, or rectangular prisms in 3D) the location of
176 /// a point in the mesh is simpler than in a more complex mesh
177 virtual void pml_locate_zeta(const Vector<double>& x,
178 FiniteElement*& el_pt) = 0;
179 };
180
181 //====================================================================
182 /// PML mesh class. Policy class for 2D PML meshes
183 //====================================================================
184 template<class ELEMENT>
185 class PMLQuadMeshBase : public RectangularQuadMesh<ELEMENT>,
186 public PMLMeshBase
187 {
188 public:
189 /// Constructor: Create the underlying rectangular quad mesh
190 PMLQuadMeshBase(const unsigned& n_pml_x,
191 const unsigned& n_pml_y,
192 const double& x_pml_start,
193 const double& x_pml_end,
194 const double& y_pml_start,
195 const double& y_pml_end,
196 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
197 : RectangularQuadMesh<ELEMENT>(n_pml_x,
198 n_pml_y,
199 x_pml_start,
200 x_pml_end,
201 y_pml_start,
202 y_pml_end,
203 time_stepper_pt)
204 {
205 }
206
207 /// Overloaded function to allow the user to locate an element
208 /// in mesh given the (Eulerian) position of a point. Note, the result
209 /// may be nonunique if the given point lies on the boundary of two
210 /// elements in the mesh. Additionally, we assume that the PML mesh is
211 /// axis aligned when deciding if the given point lies inside the mesh
213 FiniteElement*& coarse_mesh_el_pt)
214 {
215 //------------------------------------------
216 // Make sure the point lies inside the mesh:
217 //------------------------------------------
218 // Get the lower x limit of the mesh
219 const double x_min = this->x_min();
220
221 // Get the upper x limit of the mesh
222 const double x_max = this->x_max();
223
224 // Get the lower y limit of the mesh
225 const double y_min = this->y_min();
226
227 // Get the upper y limit of the mesh
228 const double y_max = this->y_max();
229
230 // Due to roundoff error we will choose a small tolerance which will be
231 // used to determine whether or not the point lies in the mesh
232 double tol = 1.0e-12;
233
234 // Assuming the mesh is axis-aligned we merely need to check if the
235 // x-coordinate of the input lies in the range [x_min,x_max] and the
236 // y-coordinate of the input lies in the range [y_min,y_max]
237 if ((x[0] < x_min - tol) || (x[0] > x_max + tol) ||
238 (x[1] < y_min - tol) || (x[1] > y_max + tol))
239 {
240 // Open an output string stream
241 std::ostringstream error_message_stream;
242
243 // Create an error message
244 error_message_stream << "Point does not lie in the mesh." << std::endl;
245
246 // Throw the error message
247 throw OomphLibError(error_message_stream.str(),
248 OOMPH_CURRENT_FUNCTION,
249 OOMPH_EXCEPTION_LOCATION);
250 }
251
252 //-----------------------------------------
253 // Collect some information about the mesh:
254 //-----------------------------------------
255 // Find out how many elements there are in the x-direction
256 const unsigned nx = this->nx();
257
258 // Find out how many elements there are in the y-direction
259 const unsigned ny = this->ny();
260
261 // Find out how many nodes there are in one direction in each element
262 unsigned nnode_1d = this->finite_element_pt(0)->nnode_1d();
263
264 // Find out how many nodes there are in each element
265 unsigned nnode = this->finite_element_pt(0)->nnode();
266
267 // Create a pointer to store the element pointer as a FiniteElement
268 FiniteElement* el_pt = 0;
269
270 // Vector to store the coordinate of each element corner node which
271 // lies on the bottom boundary of the mesh and varies, i.e. if we're
272 // on the bottom boundary of the PML mesh then the y-coordinate will
273 // be fixed therefore we need only store the x-coordinates of the
274 // corner nodes of each element attached to this boundary
275 Vector<double> bottom_boundary_x_coordinates(nx + 1, 0.0);
276
277 // Vector to store the coordinate of each element corner node which lies
278 // on the right boundary of the mesh and varies
279 Vector<double> right_boundary_y_coordinates(ny + 1, 0.0);
280
281 // Recall, we already know the start and end coordinates of the mesh
282 // in the x-direction so we can assign these entries straight away.
283 // Note, these values are exactly the same as those had we instead
284 // considered the upper boundary so it is only necessary to store
285 // the information of the one of these boundaries. A similar argument
286 // holds for the left and right boundaries.
287
288 // The first entry of bottom_boundary_x_coordinates is:
289 bottom_boundary_x_coordinates[0] = x_min;
290
291 // The last entry is:
292 bottom_boundary_x_coordinates[nx] = x_max;
293
294 // The first entry of right_boundary_y_coordinates is:
295 right_boundary_y_coordinates[0] = y_min;
296
297 // The last entry is:
298 right_boundary_y_coordinates[ny] = y_max;
299
300 // To collect the remaining entries we need to loop over all of the
301 // elements on the bottom boundary and collect the x-coordinate of
302 // the bottom right node in the element. To avoid assigning the
303 // last entry twice we ignore the last element.
304
305 // Store the lower boundary ID
306 unsigned lower_boundary_id = 0;
307
308 // Store the right boundary ID
309 unsigned right_boundary_id = 1;
310
311 // Loop over the elements on the bottom boundary
312 for (unsigned i = 0; i < nx; i++)
313 {
314 // Assign the element pointer
315 el_pt = this->boundary_element_pt(lower_boundary_id, i);
316
317 // Store the x-coordinate of the bottom right node in this element
318 bottom_boundary_x_coordinates[i + 1] =
319 el_pt->node_pt(nnode_1d - 1)->x(0);
320 }
321
322 // To collect the remaining entries we need to loop over all of the
323 // elements on the right boundary and collect the y-coordinate of
324 // the top right node in the element. To avoid assigning the
325 // last entry twice we ignore the last element.
326
327 // Loop over the elements on the bottom boundary
328 for (unsigned i = 0; i < ny; i++)
329 {
330 // Assign the element pointer
331 el_pt = this->boundary_element_pt(right_boundary_id, i);
332
333 // Store the y-coordinate of the top right node in this element
334 right_boundary_y_coordinates[i + 1] = el_pt->node_pt(nnode - 1)->x(1);
335 }
336
337 //---------------------------
338 // Find the matching element:
339 //---------------------------
340 // Boolean variable to indicate that the ID of the element in the
341 // x-direction has been found. Note, the value of this ID must lie
342 // in the range [0,nx]
343 bool element_x_id_has_been_found = false;
344
345 // Boolean variable to indicate that the ID of the element in the
346 // y-direction has been found. Note, the value of this ID must lie
347 // in the range [0,ny]
348 bool element_y_id_has_been_found = false;
349
350 // Initialise the ID of the element in the x-direction
351 unsigned element_x_id = 0;
352
353 // Initialise the ID of the element in the y-direction
354 unsigned element_y_id = 0;
355
356 // Loop over all of the entries in bottom_boundary_x_coordinates
357 for (unsigned i = 0; i < nx; i++)
358 {
359 // Check if the x-coordinate of the input lies in this interval
360 if ((x[0] >= bottom_boundary_x_coordinates[i]) &&
361 (x[0] <= bottom_boundary_x_coordinates[i + 1]))
362 {
363 // The point lies in the i-th interval so the element ID is i
364 element_x_id = i;
365
366 // Indicate that the element ID has been found
367 element_x_id_has_been_found = true;
368 }
369 } // for (unsigned i=0;i<nx;i++)
370
371 // If the element ID hasn't been found
372 if (!element_x_id_has_been_found)
373 {
374 // Open an output string stream
375 std::ostringstream error_message_stream;
376
377 // Create an error message
378 error_message_stream << "The ID of the element in the x-direction "
379 << "has not been found." << std::endl;
380
381 // Throw the error message
382 throw OomphLibError(error_message_stream.str(),
383 OOMPH_CURRENT_FUNCTION,
384 OOMPH_EXCEPTION_LOCATION);
385 }
386
387 // Loop over all of the entries in right_boundary_y_coordinates
388 for (unsigned i = 0; i < ny; i++)
389 {
390 // Check if the y-coordinate of the input lies in this interval
391 if ((x[1] >= right_boundary_y_coordinates[i]) &&
392 (x[1] <= right_boundary_y_coordinates[i + 1]))
393 {
394 // The point lies in the i-th interval so the element ID is i
395 element_y_id = i;
396
397 // Indicate that the element ID has been found
398 element_y_id_has_been_found = true;
399 }
400 } // for (unsigned i=0;i<ny;i++)
401
402 // If the element ID hasn't been found
403 if (!element_y_id_has_been_found)
404 {
405 // Open an output string stream
406 std::ostringstream error_message_stream;
407
408 // Create an error message
409 error_message_stream << "The ID of the element in the y-direction "
410 << "has not been found." << std::endl;
411
412 // Throw the error message
413 throw OomphLibError(error_message_stream.str(),
414 OOMPH_CURRENT_FUNCTION,
415 OOMPH_EXCEPTION_LOCATION);
416 }
417
418 // Calculate the element number in the mesh
419 unsigned el_id = element_y_id * nx + element_x_id;
420
421 // Set the pointer to this element
422 coarse_mesh_el_pt = dynamic_cast<FiniteElement*>(this->element_pt(el_id));
423 } // End of pml_locate_zeta
424 };
425
426 //====================================================================
427 /// PML mesh, derived from RectangularQuadMesh.
428 //====================================================================
429 template<class ELEMENT>
430 class PMLQuadMesh : public PMLQuadMeshBase<ELEMENT>
431 {
432 public:
433 /// Constructor: Pass pointer to "bulk" mesh,
434 /// the boundary ID of axis aligned boundary to which the
435 /// mesh is supposed to be attached and the boundary ID in the
436 /// rectangular quad mesh that contains the pml elements.
437 /// 0: constant y, bulk mesh below;
438 /// 1: constant x, bulk mesh to the right;
439 /// 2: constant y, bulk mesh above;
440 /// 3: constant x, bulk mesh to left.
441 PMLQuadMesh(Mesh* bulk_mesh_pt,
442 const unsigned& boundary_id,
443 const unsigned& quad_boundary_id,
444 const unsigned& n_pml_x,
445 const unsigned& n_pml_y,
446 const double& x_pml_start,
447 const double& x_pml_end,
448 const double& y_pml_start,
449 const double& y_pml_end,
450 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
451 : PMLQuadMeshBase<ELEMENT>(n_pml_x,
452 n_pml_y,
453 x_pml_start,
454 x_pml_end,
455 y_pml_start,
456 y_pml_end,
457 time_stepper_pt)
458 {
459 unsigned n_boundary_node = bulk_mesh_pt->nboundary_node(boundary_id);
460
461 // Create a vector of ordered boundary nodes
462 Vector<Node*> ordered_boundary_node_pt(n_boundary_node);
463
464 // Fill the vector with the nodes on the respective boundary
465 for (unsigned n = 0; n < n_boundary_node; n++)
466 {
467 ordered_boundary_node_pt[n] =
468 bulk_mesh_pt->boundary_node_pt(boundary_id, n);
469 }
470
471 /// Sort them depending on the boundary being used
472
473 // Right boundary
474 if (quad_boundary_id == 3)
475 {
476 std::sort(ordered_boundary_node_pt.begin(),
477 ordered_boundary_node_pt.end(),
479 }
480
481 /// Top boundary
482 if (quad_boundary_id == 0)
483 {
484 std::sort(ordered_boundary_node_pt.begin(),
485 ordered_boundary_node_pt.end(),
487 }
488
489 /// Left boundary
490 if (quad_boundary_id == 1)
491 {
492 std::sort(ordered_boundary_node_pt.begin(),
493 ordered_boundary_node_pt.end(),
495 }
496
497 /// Bottom boundary
498 if (quad_boundary_id == 2)
499 {
500 std::sort(ordered_boundary_node_pt.begin(),
501 ordered_boundary_node_pt.end(),
503 }
504
505 unsigned nnode_1d = this->finite_element_pt(0)->nnode_1d();
506
507 /// Simple interior node numbering helpers
508 /// to be precomputed before the main loop
509
510 /// Top left node in element
511 unsigned interior_node_nr_helper_1 = nnode_1d * (nnode_1d - 1);
512 /// Lower right node in element
513 unsigned interior_node_nr_helper_2 = nnode_1d - 1;
514 /// Used to find nodes in top row
515 unsigned interior_node_nr_helper_3 = nnode_1d * (nnode_1d - 1) - 1;
516
517 /// Set all nodes from the PML mesh that must disappear
518 /// after merging as obsolete
519 unsigned nnod = this->nboundary_node(quad_boundary_id);
520 for (unsigned j = 0; j < nnod; j++)
521 {
522 this->boundary_node_pt(quad_boundary_id, j)->set_obsolete();
523 }
524
525 // Kill the obsolete nodes
526 this->prune_dead_nodes();
527
528 // Find the number of elements inside the PML mesh
529 unsigned n_pml_element = this->nelement();
530
531 /// Simple interior element numbering helpers
532
533 /// Last element in mesh (top right)
534 unsigned interior_element_nr_helper_1 = n_pml_element - 1;
535
536
537 // Connect the elements in the pml mesh to the ones
538 // in the triangular mesh at element level
539 unsigned count = 0;
540
541 // Each boundary requires a specific treatment
542 // Right boundary
543 if (quad_boundary_id == 3)
544 {
545 for (unsigned e = 0; e < n_pml_element; e++)
546 {
547 // If element is on the right boundary
548 if ((e % n_pml_x) == 0)
549 {
550 // Upcast from GeneralisedElement to bulk element
551 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
552
553 // Loop over all nodes in element
554 unsigned n_node = el_pt->nnode();
555 for (unsigned inod = 0; inod < n_node; inod++)
556 {
557 if (inod % nnode_1d == 0)
558 {
559 // Get the pointer from the triangular mesh
560 el_pt->node_pt(inod) = ordered_boundary_node_pt[count];
561 count++;
562
563 // Node between two elements
564 if (inod == interior_node_nr_helper_1)
565 {
566 count--;
567 }
568 }
569 }
570 }
571 }
572 }
573
574 // Top boundary
575 if (quad_boundary_id == 0)
576 {
577 for (unsigned e = 0; e < n_pml_element; e++)
578 {
579 // If element is on the right boundary
580 if ((int)(e / n_pml_x) == 0)
581 {
582 // Upcast from GeneralisedElement to bulk element
583 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
584
585 // Loop over all nodes in element
586 unsigned n_node = el_pt->nnode();
587 for (unsigned inod = 0; inod < n_node; inod++)
588 {
589 if ((int)(inod / nnode_1d) == 0)
590 {
591 // Get the pointer from the triangular mesh
592 el_pt->node_pt(inod) = ordered_boundary_node_pt[count];
593 count++;
594
595 // Node between two elements
596 if (inod == interior_node_nr_helper_2)
597 {
598 count--;
599 }
600 }
601 }
602 }
603 }
604 }
605
606 // Left boundary
607 if (quad_boundary_id == 1)
608 {
609 for (unsigned e = interior_element_nr_helper_1; e < n_pml_element; e--)
610 {
611 // If element is on the right boundary
612 if ((e % n_pml_x) == (n_pml_x - 1))
613 {
614 // Upcast from GeneralisedElement to bulk element
615 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
616
617 // Loop over all nodes in element
618 unsigned n_node = el_pt->nnode();
619 unsigned starter = n_node - 1;
620 for (unsigned inod = starter; inod <= starter; inod--)
621 {
622 if (inod % nnode_1d == interior_node_nr_helper_2)
623 {
624 // Get the pointer from the triangular mesh
625 el_pt->node_pt(inod) = ordered_boundary_node_pt[count];
626 count++;
627
628 // Node between two elements
629 if (inod == interior_node_nr_helper_2)
630 {
631 count--;
632 }
633 }
634 }
635 }
636 }
637 }
638
639 // Bottom boundary
640 if (quad_boundary_id == 2)
641 {
642 for (unsigned e = interior_element_nr_helper_1; e < n_pml_element; e--)
643 {
644 // If element is on the top boundary
645 if (e >= (n_pml_x * (n_pml_y - 1)))
646 {
647 // Upcast from GeneralisedElement to bulk element
648 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
649
650 // Loop over all nodes in element
651 unsigned n_node = el_pt->nnode();
652 unsigned starter = n_node - 1;
653 for (unsigned inod = starter; inod <= starter; inod--)
654 {
655 if (inod > interior_node_nr_helper_3)
656 {
657 // Get the pointer from the triangular mesh
658 el_pt->node_pt(inod) = ordered_boundary_node_pt[count];
659 count++;
660
661 // Node between two elements
662 if (inod == interior_node_nr_helper_1)
663 {
664 count--;
665 }
666 }
667 }
668 }
669 }
670 }
671
672 /// The alignment is done individually for each boundary
673 /// and depends on the ordering of the nodes, in this case of the
674 /// RectangularQuadMesh<2,3> for each boundary. There are operations
675 /// with mod and div 3 necessary in this case, as well as specific
676 /// mechanisms to loop over the boundary in a certain way in order
677 /// to obtain the convenient coordinates.
678
679 // Loop over all elements and make sure the coordinates are aligned
680 for (unsigned e = 0; e < n_pml_element; e++)
681 {
682 // Upcast from GeneralisedElement to bulk element
683 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
684 unsigned n_node = el_pt->nnode();
685
686 // Loop over all nodes in element
687 double temp_coordinate = 0.0;
688 for (unsigned inod = 0; inod < n_node; inod++)
689 {
690 // Check if we are looking at the left boundary of the quad mesh
691 if (quad_boundary_id == 3)
692 {
693 // If it is one of the ones on the left boundary
694 if (inod % nnode_1d == 0)
695 {
696 // Get the y-coordinate of the leftmost node in that element
697 temp_coordinate = el_pt->node_pt(inod)->x(1);
698 }
699
700 // Each node's y-coordinate is reset to be the one of the leftmost
701 // node in that element on its x-coordinate
702 el_pt->node_pt(inod)->x(1) = temp_coordinate;
703 }
704 // End of left quad boundary check
705
706 // Check if we are looking at the top boundary of the quad mesh
707 if (quad_boundary_id == 0)
708 {
709 // If it is one of the ones on the bottom boundary
710 if (inod > interior_node_nr_helper_2)
711 {
712 // Get the y-coordinate of the leftmost node in that element
713 el_pt->node_pt(inod)->x(0) =
714 el_pt->node_pt(inod - nnode_1d)->x(0);
715 }
716 }
717 // End of top quad boundary check
718 }
719 }
720
721 for (unsigned e = interior_element_nr_helper_1; e < n_pml_element; e--)
722 {
723 // Upcast from GeneralisedElement to bulk element
724 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
725 unsigned n_node = el_pt->nnode();
726
727 // Loop over all nodes in element
728 double temp_coordinate = 0.0;
729 unsigned starter = n_node - 1;
730 for (unsigned inod = starter; inod <= starter; inod--)
731 {
732 // Check if we are looking at the right boundary of the quad mesh
733 if (quad_boundary_id == 1)
734 {
735 // If it is one of the ones on the left boundary
736 if (inod % nnode_1d == interior_node_nr_helper_2)
737 {
738 // Get the y-coordinate of the leftmost node in that element
739 temp_coordinate = el_pt->node_pt(inod)->x(1);
740 }
741
742 // Each node's y-coordinate is reset to be the one of the leftmost
743 // node in that element on its x-coordinate
744 el_pt->node_pt(inod)->x(1) = temp_coordinate;
745 }
746 // End of right quad boundary check
747
748 // Check if we are looking at the top boundary of the quad mesh
749 if (quad_boundary_id == 2)
750 {
751 // If it is one of the ones on the bottom boundary
752 if (inod < interior_node_nr_helper_1)
753 {
754 // Get the y-coordinate of the leftmost node in that element
755 el_pt->node_pt(inod)->x(0) =
756 el_pt->node_pt(inod + nnode_1d)->x(0);
757 }
758 }
759 // End of top quad boundary check
760 }
761 }
762 // End of alignment
763 }
764 };
765
766
767 /// ///////////////////////////////////////////////////////////////////
768 /// ///////////////////////////////////////////////////////////////////
769 /// ///////////////////////////////////////////////////////////////////
770
771
772 //====================================================================
773 /// PML mesh, derived from RectangularQuadMesh.
774 //====================================================================
775 template<class ELEMENT>
776 class PMLCornerQuadMesh : public PMLQuadMeshBase<ELEMENT>
777 {
778 public:
779 /// Constructor: Pass pointer to "bulk" mesh
780 /// and the two existing PML meshes in order to construct the corner
781 /// PML mesh in between them based on their element number
782 /// and coordinates.
783 PMLCornerQuadMesh(Mesh* PMLQuad_mesh_x_pt,
784 Mesh* PMLQuad_mesh_y_pt,
785 Mesh* bulk_mesh_pt,
786 Node* special_corner_node_pt,
787 const unsigned& parent_boundary_x_id,
788 const unsigned& parent_boundary_y_id,
789 const unsigned& current_boundary_x_id,
790 const unsigned& current_boundary_y_id,
791 const unsigned& n_pml_x,
792 const unsigned& n_pml_y,
793 const double& x_pml_start,
794 const double& x_pml_end,
795 const double& y_pml_start,
796 const double& y_pml_end,
797 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
798 : PMLQuadMeshBase<ELEMENT>(n_pml_x,
799 n_pml_y,
800 x_pml_start,
801 x_pml_end,
802 y_pml_start,
803 y_pml_end,
804 time_stepper_pt)
805 {
806 unsigned nnode_1d = this->finite_element_pt(0)->nnode_1d();
807
808 /// Simple interior node numbering helpers
809 /// to be precomputed before the main loop
810
811 /// Top left node in element
812 unsigned interior_node_nr_helper_1 = nnode_1d * (nnode_1d - 1);
813 /// Lower right node in element
814 unsigned interior_node_nr_helper_2 = nnode_1d - 1;
815 /// Top right node in element
816 unsigned interior_node_nr_helper_3 = nnode_1d * nnode_1d - 1;
817
818 // Set up top right corner element
819 //--------------------------------
820 if ((parent_boundary_x_id == 2) && (parent_boundary_y_id == 1))
821 {
822 // Get the number of nodes to be connected on the horizontal boundary
823 unsigned n_boundary_x_node =
824 PMLQuad_mesh_x_pt->nboundary_node(parent_boundary_x_id);
825
826 // Create a vector of ordered boundary nodes
827 Vector<Node*> ordered_boundary_x_node_pt(n_boundary_x_node);
828
829 // Fill the vector with the nodes on the respective boundary
830 for (unsigned n = 0; n < n_boundary_x_node; n++)
831 {
832 ordered_boundary_x_node_pt[n] =
833 PMLQuad_mesh_x_pt->boundary_node_pt(parent_boundary_x_id, n);
834 }
835
836 // Sort them from lowest to highest (in x coordinate)
837 if (parent_boundary_x_id == 2)
838 {
839 std::sort(ordered_boundary_x_node_pt.begin(),
840 ordered_boundary_x_node_pt.end(),
842 }
843
844 // Get the number of nodes to be connected on the vertical boundary
845 unsigned n_boundary_y_node =
846 PMLQuad_mesh_y_pt->nboundary_node(parent_boundary_y_id);
847
848 // Create a vector of ordered boundary nodes
849 Vector<Node*> ordered_boundary_y_node_pt(n_boundary_y_node);
850
851 // Fill the vector with the nodes on the respective boundary
852 for (unsigned n = 0; n < n_boundary_y_node; n++)
853 {
854 ordered_boundary_y_node_pt[n] =
855 PMLQuad_mesh_y_pt->boundary_node_pt(parent_boundary_y_id, n);
856 }
857
858 // Sort them
859 if (parent_boundary_y_id == 1)
860 {
861 std::sort(ordered_boundary_y_node_pt.begin(),
862 ordered_boundary_y_node_pt.end(),
864 }
865
866 unsigned x_nnod = this->nboundary_node(current_boundary_x_id);
867 for (unsigned j = 0; j < x_nnod; j++)
868 {
869 this->boundary_node_pt(current_boundary_x_id, j)->set_obsolete();
870 }
871
872 unsigned y_nnod = this->nboundary_node(current_boundary_y_id);
873 for (unsigned j = 0; j < y_nnod; j++)
874 {
875 this->boundary_node_pt(current_boundary_y_id, j)->set_obsolete();
876 }
877
878 // Kill the obsolete nodes
879 this->prune_dead_nodes();
880
881 unsigned n_pml_element = this->nelement();
882
883 // Connect the elements in the pml mesh to the ones
884 // in the triangular mesh at element level
885 unsigned count = 0;
886
887 if (parent_boundary_y_id == 1)
888 {
889 for (unsigned e = 0; e < n_pml_element; e++)
890 {
891 // If element is on the right boundary
892 if ((e % n_pml_x) == 0)
893 {
894 // Upcast from GeneralisedElement to bulk element
895 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
896
897 // Loop over all nodes in element
898 unsigned n_node = el_pt->nnode();
899 for (unsigned inod = 0; inod < n_node; inod++)
900 {
901 // If it is one of the ones on the left boundary
902 if (e == 0)
903 {
904 if (inod == 0) el_pt->node_pt(inod) = special_corner_node_pt;
905 if ((inod % nnode_1d == 0) && (inod > 0))
906 {
907 // Get the pointer from the triangular mesh
908 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
909 count++;
910
911 // Node between two elements
912 if (inod == interior_node_nr_helper_1)
913 {
914 count--;
915 }
916 }
917 }
918 else
919 {
920 if ((inod % nnode_1d) == 0)
921 {
922 // Get the pointer from the triangular mesh
923 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
924 count++;
925
926 // Node between two elements
927 if (inod == interior_node_nr_helper_1)
928 {
929 count--;
930 }
931 }
932 }
933 }
934 }
935 }
936 }
937
938 count = 0;
939
940 if (parent_boundary_x_id == 2)
941 {
942 for (unsigned e = 0; e < n_pml_element; e++)
943 {
944 // If element is on the right boundary
945 if ((int)(e / n_pml_x) == 0)
946 {
947 // Upcast from GeneralisedElement to bulk element
948 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
949
950 // Loop over all nodes in element
951 unsigned n_node = el_pt->nnode();
952 for (unsigned inod = 0; inod < n_node; inod++)
953 {
954 if (e == 0)
955 {
956 if (((int)(inod / nnode_1d) == 0) && (inod > 0))
957 {
958 // Get the pointer from the triangular mesh
959 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
960 count++;
961
962 // Node between two elements
963 if (inod == interior_node_nr_helper_2)
964 {
965 count--;
966 }
967 }
968 }
969 else
970 {
971 if ((int)(inod / nnode_1d) == 0)
972 {
973 // Get the pointer from the triangular mesh
974 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
975 count++;
976
977 // Node between two elements
978 if (inod == interior_node_nr_helper_2)
979 {
980 count--;
981 }
982 }
983 }
984 }
985 }
986 }
987 }
988 }
989
990 // Set up bottom right corner element
991 //--------------------------------
992 if ((parent_boundary_x_id == 0) && (parent_boundary_y_id == 1))
993 {
994 // Get the number of nodes to be connected on the horizontal boundary
995 unsigned n_boundary_x_node =
996 PMLQuad_mesh_x_pt->nboundary_node(parent_boundary_x_id);
997
998 // Create a vector of ordered boundary nodes
999 Vector<Node*> ordered_boundary_x_node_pt(n_boundary_x_node);
1000
1001 // Fill the vector with the nodes on the respective boundary
1002 for (unsigned n = 0; n < n_boundary_x_node; n++)
1003 {
1004 ordered_boundary_x_node_pt[n] =
1005 PMLQuad_mesh_x_pt->boundary_node_pt(parent_boundary_x_id, n);
1006 }
1007
1008 // Sort them from lowest to highest (in x coordinate)
1009 if (parent_boundary_x_id == 0)
1010 {
1011 std::sort(ordered_boundary_x_node_pt.begin(),
1012 ordered_boundary_x_node_pt.end(),
1014 }
1015
1016 // Get the number of nodes to be connected on the vertical boundary
1017 unsigned n_boundary_y_node =
1018 PMLQuad_mesh_y_pt->nboundary_node(parent_boundary_y_id);
1019
1020 // Create a vector of ordered boundary nodes
1021 Vector<Node*> ordered_boundary_y_node_pt(n_boundary_y_node);
1022
1023 // Fill the vector with the nodes on the respective boundary
1024 for (unsigned n = 0; n < n_boundary_y_node; n++)
1025 {
1026 ordered_boundary_y_node_pt[n] =
1027 PMLQuad_mesh_y_pt->boundary_node_pt(parent_boundary_y_id, n);
1028 }
1029
1030 // Sort them
1031 if (parent_boundary_y_id == 1)
1032 {
1033 std::sort(ordered_boundary_y_node_pt.begin(),
1034 ordered_boundary_y_node_pt.end(),
1036 }
1037
1038 unsigned x_nnod = this->nboundary_node(current_boundary_x_id);
1039 for (unsigned j = 0; j < x_nnod; j++)
1040 {
1041 this->boundary_node_pt(current_boundary_x_id, j)->set_obsolete();
1042 }
1043
1044 unsigned y_nnod = this->nboundary_node(current_boundary_y_id);
1045 for (unsigned j = 0; j < y_nnod; j++)
1046 {
1047 this->boundary_node_pt(current_boundary_y_id, j)->set_obsolete();
1048 }
1049
1050 // Kill the obsolete nodes
1051 this->prune_dead_nodes();
1052
1053 // Get the number of elements in the PML mesh
1054 unsigned n_pml_element = this->nelement();
1055
1056 // Connect the elements in the pml mesh to the ones
1057 // in the triangular mesh at element level
1058 unsigned count = 0;
1059
1060 if (parent_boundary_y_id == 1)
1061 {
1062 for (unsigned e = 0; e < n_pml_element; e++)
1063 {
1064 // If element is on the right boundary
1065 if ((e % n_pml_x) == 0)
1066 {
1067 // Upcast from GeneralisedElement to bulk element
1068 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1069
1070 // Loop over all nodes in element
1071 unsigned n_node = el_pt->nnode();
1072 for (unsigned inod = 0; inod < n_node; inod++)
1073 {
1074 if (e == ((n_pml_x) * (n_pml_y - 1)))
1075 {
1076 if (inod == interior_node_nr_helper_1)
1077 {
1078 el_pt->node_pt(inod) = special_corner_node_pt;
1079 }
1080 if ((inod % nnode_1d == 0) &&
1081 (inod < interior_node_nr_helper_1))
1082 {
1083 // Get the pointer from the triangular mesh
1084 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1085 count++;
1086
1087 // Node between two elements
1088 if (inod == interior_node_nr_helper_1)
1089 {
1090 count--;
1091 }
1092 }
1093 }
1094 else
1095 {
1096 if ((inod % nnode_1d) == 0)
1097 {
1098 // Get the pointer from the triangular mesh
1099 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1100 count++;
1101
1102 // Node between two elements
1103 if (inod == interior_node_nr_helper_1)
1104 {
1105 count--;
1106 }
1107 }
1108 }
1109 }
1110 }
1111 }
1112 }
1113
1114 count = 0;
1115
1116 if (parent_boundary_x_id == 0)
1117 {
1118 for (unsigned e = 0; e < n_pml_element; e++)
1119 {
1120 // If element is on the right boundary
1121 if (e >= ((n_pml_x - 0) * (n_pml_y - 1)))
1122 {
1123 // Upcast from GeneralisedElement to bulk element
1124 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1125
1126 // Loop over all nodes in element
1127 unsigned n_node = el_pt->nnode();
1128 for (unsigned inod = 0; inod < n_node; inod++)
1129 {
1130 if (e == ((n_pml_x) * (n_pml_y - 1)))
1131 {
1132 if (((unsigned)(inod / nnode_1d) ==
1133 interior_node_nr_helper_2) &&
1134 (inod > interior_node_nr_helper_1))
1135 {
1136 // Get the pointer from the triangular mesh
1137 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1138 count++;
1139
1140 // Node between two elements
1141 if (inod == interior_node_nr_helper_3)
1142 {
1143 count--;
1144 }
1145 }
1146 }
1147 else
1148 {
1149 if ((unsigned)(inod / nnode_1d) == interior_node_nr_helper_2)
1150 {
1151 // Get the pointer from the triangular mesh
1152 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1153 count++;
1154
1155 // Node between two elements
1156 if (inod == interior_node_nr_helper_3)
1157 {
1158 count--;
1159 }
1160 }
1161 }
1162 }
1163 }
1164 }
1165 }
1166 }
1167
1168 // Set up top left corner element
1169 //--------------------------------
1170 if ((parent_boundary_x_id == 2) && (parent_boundary_y_id == 3))
1171 {
1172 // Get the number of nodes to be connected on the horizontal boundary
1173 unsigned n_boundary_x_node =
1174 PMLQuad_mesh_x_pt->nboundary_node(parent_boundary_x_id);
1175
1176 // Create a vector of ordered boundary nodes
1177 Vector<Node*> ordered_boundary_x_node_pt(n_boundary_x_node);
1178
1179 // Fill the vector with the nodes on the respective boundary
1180 for (unsigned n = 0; n < n_boundary_x_node; n++)
1181 {
1182 ordered_boundary_x_node_pt[n] =
1183 PMLQuad_mesh_x_pt->boundary_node_pt(parent_boundary_x_id, n);
1184 }
1185
1186 // Sort them from lowest to highest (in x coordinate)
1187 if (parent_boundary_x_id == 2)
1188 {
1189 std::sort(ordered_boundary_x_node_pt.begin(),
1190 ordered_boundary_x_node_pt.end(),
1192 }
1193
1194 // Get the number of nodes to be connected on the vertical boundary
1195 unsigned n_boundary_y_node =
1196 PMLQuad_mesh_y_pt->nboundary_node(parent_boundary_y_id);
1197
1198 // Create a vector of ordered boundary nodes
1199 Vector<Node*> ordered_boundary_y_node_pt(n_boundary_y_node);
1200
1201 // Fill the vector with the nodes on the respective boundary
1202 for (unsigned n = 0; n < n_boundary_y_node; n++)
1203 {
1204 ordered_boundary_y_node_pt[n] =
1205 PMLQuad_mesh_y_pt->boundary_node_pt(parent_boundary_y_id, n);
1206 }
1207
1208 // Sort them from lowest to highest (in x coordinate)
1209 if (parent_boundary_y_id == 1)
1210 {
1211 std::sort(ordered_boundary_y_node_pt.begin(),
1212 ordered_boundary_y_node_pt.end(),
1214 }
1215
1216 unsigned x_nnod = this->nboundary_node(current_boundary_x_id);
1217 for (unsigned j = 0; j < x_nnod; j++)
1218 {
1219 this->boundary_node_pt(current_boundary_x_id, j)->set_obsolete();
1220 }
1221
1222 unsigned y_nnod = this->nboundary_node(current_boundary_y_id);
1223 for (unsigned j = 0; j < y_nnod; j++)
1224 {
1225 this->boundary_node_pt(current_boundary_y_id, j)->set_obsolete();
1226 }
1227
1228 // Kill the obsolete nodes
1229 this->prune_dead_nodes();
1230
1231 // Get the number of elements in the PML mesh
1232 unsigned n_pml_element = this->nelement();
1233
1234 // Connect the elements in the pml mesh to the ones
1235 // in the triangular mesh at element level
1236 unsigned count = 0;
1237
1238 if (parent_boundary_y_id == 3)
1239 {
1240 for (unsigned e = 0; e < n_pml_element; e++)
1241 {
1242 // If element is on the right boundary
1243 if ((e % n_pml_x) == (n_pml_x - 1))
1244 {
1245 // Upcast from GeneralisedElement to bulk element
1246 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1247
1248 // Loop over all nodes in element
1249 unsigned n_node = el_pt->nnode();
1250 for (unsigned inod = 0; inod < n_node; inod++)
1251 {
1252 if (e == (n_pml_x - 1))
1253 {
1254 if (inod == interior_node_nr_helper_2)
1255 el_pt->node_pt(inod) = special_corner_node_pt;
1256 if ((inod % nnode_1d == interior_node_nr_helper_2) &&
1257 (inod > (nnode_1d - 1)))
1258 {
1259 // Get the pointer from the triangular mesh
1260 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1261 count++;
1262
1263 // Node between two elements
1264 if (inod == interior_node_nr_helper_3)
1265 {
1266 count--;
1267 }
1268 }
1269 }
1270 else
1271 {
1272 if ((inod % nnode_1d) == interior_node_nr_helper_2)
1273 {
1274 // Get the pointer from the triangular mesh
1275 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1276 count++;
1277
1278 // Node between two elements
1279 if (inod == interior_node_nr_helper_3)
1280 {
1281 count--;
1282 }
1283 }
1284 }
1285 }
1286 }
1287 }
1288 }
1289
1290 count = 0;
1291
1292 if (parent_boundary_x_id == 2)
1293 {
1294 for (unsigned e = 0; e < n_pml_element; e++)
1295 {
1296 // If element is on the right boundary
1297 if ((int)(e / n_pml_x) == 0)
1298 {
1299 // Upcast from GeneralisedElement to bulk element
1300 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1301
1302 // Loop over all nodes in element
1303 unsigned n_node = el_pt->nnode();
1304 for (unsigned inod = 0; inod < n_node; inod++)
1305 {
1306 // If it is one of the ones on the left boundary
1307 if (e == (n_pml_x - 1))
1308 {
1309 if (((int)(inod / nnode_1d) == 0) &&
1310 (inod < interior_node_nr_helper_2))
1311 {
1312 // Get the pointer from the triangular mesh
1313 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1314 count++;
1315
1316 // Node between two elements
1317 if (inod == (nnode_1d - 1))
1318 {
1319 count--;
1320 }
1321 }
1322 }
1323 else
1324 {
1325 if ((int)(inod / nnode_1d) == 0)
1326 {
1327 // Get the pointer from the triangular mesh
1328 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1329 count++;
1330
1331 // Node between two elements
1332 if (inod == interior_node_nr_helper_2)
1333 {
1334 count--;
1335 }
1336 }
1337 }
1338 }
1339 }
1340 }
1341 }
1342 }
1343
1344 // Set up bottom left corner element
1345 //--------------------------------
1346 if ((parent_boundary_x_id == 0) && (parent_boundary_y_id == 3))
1347 {
1348 // Get the number of nodes to be connected on the horizontal boundary
1349 unsigned n_boundary_x_node =
1350 PMLQuad_mesh_x_pt->nboundary_node(parent_boundary_x_id);
1351
1352 // Create a vector of ordered boundary nodes
1353 Vector<Node*> ordered_boundary_x_node_pt(n_boundary_x_node);
1354
1355 // Fill the vector with the nodes on the respective boundary
1356 for (unsigned n = 0; n < n_boundary_x_node; n++)
1357 {
1358 ordered_boundary_x_node_pt[n] =
1359 PMLQuad_mesh_x_pt->boundary_node_pt(parent_boundary_x_id, n);
1360 }
1361
1362 // Sort them
1363 if (parent_boundary_x_id == 0)
1364 {
1365 std::sort(ordered_boundary_x_node_pt.begin(),
1366 ordered_boundary_x_node_pt.end(),
1368 }
1369
1370 // Get the number of nodes to be connected on the vertical boundary
1371 unsigned n_boundary_y_node =
1372 PMLQuad_mesh_y_pt->nboundary_node(parent_boundary_y_id);
1373
1374 // Create a vector of ordered boundary nodes
1375 Vector<Node*> ordered_boundary_y_node_pt(n_boundary_y_node);
1376
1377 // Fill the vector with the nodes on the respective boundary
1378 for (unsigned n = 0; n < n_boundary_y_node; n++)
1379 {
1380 ordered_boundary_y_node_pt[n] =
1381 PMLQuad_mesh_y_pt->boundary_node_pt(parent_boundary_y_id, n);
1382 }
1383
1384 // Sort them
1385 if (parent_boundary_y_id == 3)
1386 {
1387 std::sort(ordered_boundary_y_node_pt.begin(),
1388 ordered_boundary_y_node_pt.end(),
1390 }
1391
1392 unsigned x_nnod = this->nboundary_node(current_boundary_x_id);
1393 for (unsigned j = 0; j < x_nnod; j++)
1394 {
1395 this->boundary_node_pt(current_boundary_x_id, j)->set_obsolete();
1396 }
1397
1398 unsigned y_nnod = this->nboundary_node(current_boundary_y_id);
1399 for (unsigned j = 0; j < y_nnod; j++)
1400 {
1401 this->boundary_node_pt(current_boundary_y_id, j)->set_obsolete();
1402 }
1403
1404 // Kill the obsolete nodes
1405 this->prune_dead_nodes();
1406
1407 unsigned n_pml_element = this->nelement();
1408
1409 // Connect the elements in the pml mesh to the ones
1410 // in the triangular mesh at element level
1411 unsigned count = 0;
1412
1413 if (parent_boundary_y_id == 3)
1414 {
1415 for (unsigned e = 0; e < n_pml_element; e++)
1416 {
1417 // If element is on the right boundary
1418 if ((e % n_pml_x) == (n_pml_x - 1))
1419 {
1420 // Upcast from GeneralisedElement to bulk element
1421 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1422
1423 // Loop over all nodes in element
1424 unsigned n_node = el_pt->nnode();
1425 for (unsigned inod = 0; inod < n_node; inod++)
1426 {
1427 if (e == (n_pml_element - 1))
1428 {
1429 if (inod == interior_node_nr_helper_3)
1430 {
1431 el_pt->node_pt(inod) = special_corner_node_pt;
1432 }
1433 if ((inod % nnode_1d == interior_node_nr_helper_2) &&
1434 (inod < interior_node_nr_helper_3))
1435 {
1436 // Get the pointer from the triangular mesh
1437 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1438 count++;
1439
1440 // Node between two elements
1441 if (inod == interior_node_nr_helper_3)
1442 {
1443 count--;
1444 }
1445 }
1446 }
1447 else
1448 {
1449 if ((inod % nnode_1d) == interior_node_nr_helper_2)
1450 {
1451 // Get the pointer from the triangular mesh
1452 el_pt->node_pt(inod) = ordered_boundary_y_node_pt[count];
1453 count++;
1454
1455 // Node between two elements
1456 if (inod == interior_node_nr_helper_3)
1457 {
1458 count--;
1459 }
1460 }
1461 }
1462 }
1463 }
1464 }
1465 }
1466
1467 count = 0;
1468
1469 if (parent_boundary_x_id == 0)
1470 {
1471 for (unsigned e = 0; e < n_pml_element; e++)
1472 {
1473 // If element is on the right boundary
1474 if (e >= ((n_pml_x) * (n_pml_y - 1)))
1475 {
1476 // Upcast from GeneralisedElement to bulk element
1477 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->element_pt(e));
1478
1479 // Loop over all nodes in element
1480 unsigned n_node = el_pt->nnode();
1481 for (unsigned inod = 0; inod < n_node; inod++)
1482 {
1483 if (e == (n_pml_element - 1))
1484 {
1485 if (((unsigned)(inod / nnode_1d) ==
1486 interior_node_nr_helper_2) &&
1487 (inod < interior_node_nr_helper_3))
1488 {
1489 // Get the pointer from the triangular mesh
1490 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1491 count++;
1492
1493 // Node between two elements
1494 if (inod == interior_node_nr_helper_3)
1495 {
1496 count--;
1497 }
1498 }
1499 }
1500 else
1501 {
1502 if ((unsigned)(inod / nnode_1d) == interior_node_nr_helper_2)
1503 {
1504 // Get the pointer from the triangular mesh
1505 el_pt->node_pt(inod) = ordered_boundary_x_node_pt[count];
1506 count++;
1507
1508 // Node between two elements
1509 if (inod == interior_node_nr_helper_3)
1510 {
1511 count--;
1512 }
1513 }
1514 }
1515 }
1516 }
1517 }
1518 }
1519 }
1520 }
1521 };
1522
1523 /// /////////////////////////////////////////////////////////////
1524 /// /////////////////////////////////////////////////////////////
1525 /// /////////////////////////////////////////////////////////////
1526
1527
1528 //===================================================================
1529 /// All helper routines for 2D bulk boundary mesh usage in order to
1530 /// generate PML meshes aligned to the main mesh
1531 //===================================================================
1532 namespace TwoDimensionalPMLHelper
1533 {
1534 //============================================================================
1535 /// "Constructor" for PML mesh,aligned with the right physical domain
1536 /// boundary
1537 //============================================================================
1538 template<class ASSOCIATED_PML_QUAD_ELEMENT>
1540 Mesh* bulk_mesh_pt,
1541 const unsigned& right_boundary_id,
1542 const unsigned& n_x_right_pml,
1543 const double& width_x_right_pml,
1544 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1545 {
1546 // Look at the right boundary of the triangular mesh
1547 unsigned n_right_boundary_node =
1548 bulk_mesh_pt->nboundary_node(right_boundary_id);
1549
1550 // Create a vector of ordered boundary nodes
1551 Vector<Node*> ordered_right_boundary_node_pt(n_right_boundary_node);
1552
1553 // Fill the vector with the nodes on the respective boundary
1554 for (unsigned n = 0; n < n_right_boundary_node; n++)
1555 {
1556 ordered_right_boundary_node_pt[n] =
1557 bulk_mesh_pt->boundary_node_pt(right_boundary_id, n);
1558 }
1559
1560 // Sort them from lowest to highest (in y coordinate)
1561 std::sort(ordered_right_boundary_node_pt.begin(),
1562 ordered_right_boundary_node_pt.end(),
1564
1565 // The number of elements in y is taken from the triangular mesh
1566 unsigned n_y_right_pml =
1567 bulk_mesh_pt->nboundary_element(right_boundary_id);
1568
1569 // Specific PML sizes needed, taken directly from physical domain
1570 double l_pml_right_x_start = ordered_right_boundary_node_pt[0]->x(0);
1571 /// PML layer with added to the bulk mesh coordinate
1572 double l_pml_right_x_end =
1573 width_x_right_pml + ordered_right_boundary_node_pt[0]->x(0);
1574 double l_pml_right_y_start = ordered_right_boundary_node_pt[0]->x(1);
1575 double l_pml_right_y_end =
1576 ordered_right_boundary_node_pt[n_right_boundary_node - 1]->x(1);
1577
1578 // Rectangular boundary id to be merged with triangular mesh
1579 unsigned right_quadPML_boundary_id = 3;
1580
1581 // Create the mesh to be designated to the PML
1582 Mesh* pml_right_mesh_pt = 0;
1583
1584 // Build the right one
1585 pml_right_mesh_pt =
1587 right_boundary_id,
1588 right_quadPML_boundary_id,
1589 n_x_right_pml,
1590 n_y_right_pml,
1591 l_pml_right_x_start,
1592 l_pml_right_x_end,
1593 l_pml_right_y_start,
1594 l_pml_right_y_end,
1595 time_stepper_pt);
1596
1597 // Enable PML damping on the entire mesh
1598 unsigned n_element_pml_right = pml_right_mesh_pt->nelement();
1599 for (unsigned e = 0; e < n_element_pml_right; e++)
1600 {
1601 // Upcast
1602 PMLElementBase<2>* el_pt =
1603 dynamic_cast<PMLElementBase<2>*>(pml_right_mesh_pt->element_pt(e));
1604 el_pt->enable_pml(0, l_pml_right_x_start, l_pml_right_x_end);
1605 }
1606
1607 // Get the values to be pinned from the first element (which we
1608 // assume exists!)
1609 PMLElementBase<2>* el_pt =
1610 dynamic_cast<PMLElementBase<2>*>(pml_right_mesh_pt->element_pt(0));
1611 Vector<unsigned> values_to_pin;
1612 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
1613 unsigned npin = values_to_pin.size();
1614
1615 // Exterior boundary needs to be set to Dirichlet 0
1616 // in both real and imaginary parts
1617 unsigned n_bound_pml_right = pml_right_mesh_pt->nboundary();
1618 for (unsigned b = 0; b < n_bound_pml_right; b++)
1619 {
1620 unsigned n_node = pml_right_mesh_pt->nboundary_node(b);
1621 for (unsigned n = 0; n < n_node; n++)
1622 {
1623 Node* nod_pt = pml_right_mesh_pt->boundary_node_pt(b, n);
1624 if (b == 1)
1625 {
1626 for (unsigned j = 0; j < npin; j++)
1627 {
1628 unsigned j_val = values_to_pin[j];
1629 nod_pt->pin(j_val);
1630 nod_pt->set_value(j_val, 0.0);
1631 }
1632 }
1633 }
1634 }
1635
1636 /// Return the finalized mesh, with PML enabled
1637 /// and boundary conditions added
1638 return pml_right_mesh_pt;
1639 }
1640
1641 //===========================================================================
1642 /// "Constructor" for PML mesh, aligned with the top physical domain
1643 /// boundary
1644 //===========================================================================
1645 template<class ASSOCIATED_PML_QUAD_ELEMENT>
1647 Mesh* bulk_mesh_pt,
1648 const unsigned& top_boundary_id,
1649 const unsigned& n_y_top_pml,
1650 const double& width_y_top_pml,
1651 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1652 {
1653 // Look at the top boundary of the triangular mesh
1654 unsigned n_top_boundary_node =
1655 bulk_mesh_pt->nboundary_node(top_boundary_id);
1656
1657 // Create a vector of ordered boundary nodes
1658 Vector<Node*> ordered_top_boundary_node_pt(n_top_boundary_node);
1659
1660 // Fill the vector with the nodes on the respective boundary
1661 for (unsigned n = 0; n < n_top_boundary_node; n++)
1662 {
1663 ordered_top_boundary_node_pt[n] =
1664 bulk_mesh_pt->boundary_node_pt(top_boundary_id, n);
1665 }
1666
1667 // Sort them from lowest to highest (in x coordinate)
1668 std::sort(ordered_top_boundary_node_pt.begin(),
1669 ordered_top_boundary_node_pt.end(),
1671
1672 // The number of elements in x is taken from the triangular mesh
1673 unsigned n_x_top_pml = bulk_mesh_pt->nboundary_element(top_boundary_id);
1674
1675 // Specific PML sizes needed, taken directly from physical domain
1676 double l_pml_top_x_start = ordered_top_boundary_node_pt[0]->x(0);
1677 double l_pml_top_x_end =
1678 ordered_top_boundary_node_pt[n_top_boundary_node - 1]->x(0);
1679 double l_pml_top_y_start = ordered_top_boundary_node_pt[0]->x(1);
1680 /// PML layer width added to the bulk mesh coordinate
1681 double l_pml_top_y_end =
1682 width_y_top_pml + ordered_top_boundary_node_pt[0]->x(1);
1683
1684 unsigned top_quadPML_boundary_id = 0;
1685
1686 // Create the mesh to be designated to the PML
1687 Mesh* pml_top_mesh_pt = 0;
1688
1689 // Build the top PML mesh
1690 pml_top_mesh_pt =
1692 top_boundary_id,
1693 top_quadPML_boundary_id,
1694 n_x_top_pml,
1695 n_y_top_pml,
1696 l_pml_top_x_start,
1697 l_pml_top_x_end,
1698 l_pml_top_y_start,
1699 l_pml_top_y_end,
1700 time_stepper_pt);
1701
1702 // Enable PML damping on the entire mesh
1703 unsigned n_element_pml_top = pml_top_mesh_pt->nelement();
1704 for (unsigned e = 0; e < n_element_pml_top; e++)
1705 {
1706 // Upcast
1707 PMLElementBase<2>* el_pt =
1708 dynamic_cast<PMLElementBase<2>*>(pml_top_mesh_pt->element_pt(e));
1709 el_pt->enable_pml(1, l_pml_top_y_start, l_pml_top_y_end);
1710 }
1711
1712 // Get the values to be pinned from the first element (which we
1713 // assume exists!)
1714 PMLElementBase<2>* el_pt =
1715 dynamic_cast<PMLElementBase<2>*>(pml_top_mesh_pt->element_pt(0));
1716 Vector<unsigned> values_to_pin;
1717 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
1718 unsigned npin = values_to_pin.size();
1719
1720 // Exterior boundary needs to be set to Dirichlet 0
1721 // for both real and imaginary parts of all fields
1722 // in the problem
1723 unsigned n_bound_pml_top = pml_top_mesh_pt->nboundary();
1724 for (unsigned b = 0; b < n_bound_pml_top; b++)
1725 {
1726 unsigned n_node = pml_top_mesh_pt->nboundary_node(b);
1727 for (unsigned n = 0; n < n_node; n++)
1728 {
1729 Node* nod_pt = pml_top_mesh_pt->boundary_node_pt(b, n);
1730 if (b == 2)
1731 {
1732 for (unsigned j = 0; j < npin; j++)
1733 {
1734 unsigned j_val = values_to_pin[j];
1735 nod_pt->pin(j_val);
1736 nod_pt->set_value(j_val, 0.0);
1737 }
1738 }
1739 }
1740 }
1741
1742 /// Return the finalized mesh, with PML enabled
1743 /// and boundary conditions added
1744 return pml_top_mesh_pt;
1745 }
1746
1747 //============================================================================
1748 /// "Constructor" for PML mesh, aligned with the left physical domain
1749 /// boundary
1750 //============================================================================
1751 template<class ASSOCIATED_PML_QUAD_ELEMENT>
1753 Mesh* bulk_mesh_pt,
1754 const unsigned& left_boundary_id,
1755 const unsigned& n_x_left_pml,
1756 const double& width_x_left_pml,
1757 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1758 {
1759 // Look at the left boundary of the triangular mesh
1760 unsigned n_left_boundary_node =
1761 bulk_mesh_pt->nboundary_node(left_boundary_id);
1762
1763 // Create a vector of ordered boundary nodes
1764 Vector<Node*> ordered_left_boundary_node_pt(n_left_boundary_node);
1765
1766 // Fill the vector with the nodes on the respective boundary
1767 for (unsigned n = 0; n < n_left_boundary_node; n++)
1768 {
1769 ordered_left_boundary_node_pt[n] =
1770 bulk_mesh_pt->boundary_node_pt(left_boundary_id, n);
1771 }
1772
1773 // Sort them from lowest to highest (in y coordinate)
1774 std::sort(ordered_left_boundary_node_pt.begin(),
1775 ordered_left_boundary_node_pt.end(),
1777
1778 // The number of elements in y is taken from the triangular mesh
1779 unsigned n_y_left_pml = bulk_mesh_pt->nboundary_element(left_boundary_id);
1780
1781 // Specific PML sizes needed, taken directly from physical domain
1782 /// PML layer width subtracted from left bulk mesh coordinate
1783 double l_pml_left_x_start =
1784 -width_x_left_pml +
1785 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(0);
1786 double l_pml_left_x_end =
1787 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(0);
1788 double l_pml_left_y_start =
1789 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(1);
1790 double l_pml_left_y_end = ordered_left_boundary_node_pt[0]->x(1);
1791
1792 unsigned left_quadPML_boundary_id = 1;
1793
1794 // Create the mesh to be designated to the PML
1795 Mesh* pml_left_mesh_pt = 0;
1796
1797 // Build the left PML mesh
1798 pml_left_mesh_pt =
1800 left_boundary_id,
1801 left_quadPML_boundary_id,
1802 n_x_left_pml,
1803 n_y_left_pml,
1804 l_pml_left_x_start,
1805 l_pml_left_x_end,
1806 l_pml_left_y_start,
1807 l_pml_left_y_end,
1808 time_stepper_pt);
1809
1810 // Enable PML damping on the entire mesh
1811 unsigned n_element_pml_left = pml_left_mesh_pt->nelement();
1812 for (unsigned e = 0; e < n_element_pml_left; e++)
1813 {
1814 // Upcast
1815 PMLElementBase<2>* el_pt =
1816 dynamic_cast<PMLElementBase<2>*>(pml_left_mesh_pt->element_pt(e));
1817 el_pt->enable_pml(0, l_pml_left_x_end, l_pml_left_x_start);
1818 }
1819
1820 // Get the values to be pinned from the first element (which we
1821 // assume exists!)
1822 PMLElementBase<2>* el_pt =
1823 dynamic_cast<PMLElementBase<2>*>(pml_left_mesh_pt->element_pt(0));
1824 Vector<unsigned> values_to_pin;
1825 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
1826 unsigned npin = values_to_pin.size();
1827
1828 // Exterior boundary needs to be set to Dirichlet 0
1829 // for both real and imaginary parts of all fields
1830 // in the problem
1831 unsigned n_bound_pml_left = pml_left_mesh_pt->nboundary();
1832 for (unsigned b = 0; b < n_bound_pml_left; b++)
1833 {
1834 unsigned n_node = pml_left_mesh_pt->nboundary_node(b);
1835 for (unsigned n = 0; n < n_node; n++)
1836 {
1837 Node* nod_pt = pml_left_mesh_pt->boundary_node_pt(b, n);
1838 if (b == 3)
1839 {
1840 for (unsigned j = 0; j < npin; j++)
1841 {
1842 unsigned j_val = values_to_pin[j];
1843 nod_pt->pin(j_val);
1844 nod_pt->set_value(j_val, 0.0);
1845 }
1846 }
1847 }
1848 }
1849
1850 /// Return the finalized mesh, with PML enabled
1851 /// and boundary conditions added
1852 return pml_left_mesh_pt;
1853 }
1854
1855 //============================================================================
1856 /// "Constructor" for PML mesh,aligned with the bottom physical domain
1857 /// boundary
1858 //============================================================================
1859 template<class ASSOCIATED_PML_QUAD_ELEMENT>
1861 Mesh* bulk_mesh_pt,
1862 const unsigned& bottom_boundary_id,
1863 const unsigned& n_y_bottom_pml,
1864 const double& width_y_bottom_pml,
1865 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1866 {
1867 // Look at the bottom boundary of the triangular mesh
1868 unsigned n_bottom_boundary_node =
1869 bulk_mesh_pt->nboundary_node(bottom_boundary_id);
1870
1871 // Create a vector of ordered boundary nodes
1872 Vector<Node*> ordered_bottom_boundary_node_pt(n_bottom_boundary_node);
1873
1874 // Fill the vector with the nodes on the respective boundary
1875 for (unsigned n = 0; n < n_bottom_boundary_node; n++)
1876 {
1877 ordered_bottom_boundary_node_pt[n] =
1878 bulk_mesh_pt->boundary_node_pt(bottom_boundary_id, n);
1879 }
1880
1881 // Sort them from highest to lowest (in x coordinate)
1882 std::sort(ordered_bottom_boundary_node_pt.begin(),
1883 ordered_bottom_boundary_node_pt.end(),
1885
1886 // The number of elements in y is taken from the triangular mesh
1887 unsigned n_x_bottom_pml =
1888 bulk_mesh_pt->nboundary_element(bottom_boundary_id);
1889
1890 // Specific PML sizes needed, taken directly from physical domain
1891 double l_pml_bottom_x_start =
1892 ordered_bottom_boundary_node_pt[n_bottom_boundary_node - 1]->x(0);
1893 double l_pml_bottom_x_end = ordered_bottom_boundary_node_pt[0]->x(0);
1894 /// PML layer width subtracted from the bulk mesh lower
1895 /// boundary coordinate
1896 double l_pml_bottom_y_start =
1897 -width_y_bottom_pml + ordered_bottom_boundary_node_pt[0]->x(1);
1898 double l_pml_bottom_y_end = ordered_bottom_boundary_node_pt[0]->x(1);
1899
1900 unsigned bottom_quadPML_boundary_id = 2;
1901
1902 // Create the mesh to be designated to the PML
1903 Mesh* pml_bottom_mesh_pt = 0;
1904
1905 // Build the bottom PML mesh
1906 pml_bottom_mesh_pt =
1908 bottom_boundary_id,
1909 bottom_quadPML_boundary_id,
1910 n_x_bottom_pml,
1911 n_y_bottom_pml,
1912 l_pml_bottom_x_start,
1913 l_pml_bottom_x_end,
1914 l_pml_bottom_y_start,
1915 l_pml_bottom_y_end,
1916 time_stepper_pt);
1917
1918 // Enable PML damping on the entire mesh
1919 unsigned n_element_pml_bottom = pml_bottom_mesh_pt->nelement();
1920 for (unsigned e = 0; e < n_element_pml_bottom; e++)
1921 {
1922 // Upcast
1923 PMLElementBase<2>* el_pt =
1924 dynamic_cast<PMLElementBase<2>*>(pml_bottom_mesh_pt->element_pt(e));
1925 el_pt->enable_pml(1, l_pml_bottom_y_end, l_pml_bottom_y_start);
1926 }
1927
1928 // Get the values to be pinned from the first element (which we
1929 // assume exists!)
1930 PMLElementBase<2>* el_pt =
1931 dynamic_cast<PMLElementBase<2>*>(pml_bottom_mesh_pt->element_pt(0));
1932 Vector<unsigned> values_to_pin;
1933 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
1934 unsigned npin = values_to_pin.size();
1935
1936 // Exterior boundary needs to be set to Dirichlet 0
1937 // for both real and imaginary parts of all fields
1938 // in the problem
1939 unsigned n_bound_pml_bottom = pml_bottom_mesh_pt->nboundary();
1940 for (unsigned b = 0; b < n_bound_pml_bottom; b++)
1941 {
1942 unsigned n_node = pml_bottom_mesh_pt->nboundary_node(b);
1943 for (unsigned n = 0; n < n_node; n++)
1944 {
1945 Node* nod_pt = pml_bottom_mesh_pt->boundary_node_pt(b, n);
1946 if (b == 0)
1947 {
1948 for (unsigned j = 0; j < npin; j++)
1949 {
1950 unsigned j_val = values_to_pin[j];
1951 nod_pt->pin(j_val);
1952 nod_pt->set_value(j_val, 0.0);
1953 }
1954 }
1955 }
1956 }
1957
1958 /// Return the finalized mesh, with PML enabled
1959 /// and boundary conditions added
1960 return pml_bottom_mesh_pt;
1961 }
1962
1963 //==========================================================================
1964 /// "Constructor" for PML top right corner mesh,
1965 /// aligned with the existing PML meshes
1966 //==========================================================================
1967 template<class ASSOCIATED_PML_QUAD_ELEMENT>
1969 Mesh* pml_right_mesh_pt,
1970 Mesh* pml_top_mesh_pt,
1971 Mesh* bulk_mesh_pt,
1972 const unsigned& right_boundary_id,
1973 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
1974 {
1975 /// Relevant boundary id's to be used in construction
1976 /// Parent id refers to already existing PML meshes
1977 unsigned parent_boundary_x_id = 2;
1978 unsigned parent_boundary_y_id = 1;
1979 // Current id refers to the mesh that is to be constructed
1980 unsigned current_boundary_x_id = 0;
1981 unsigned current_boundary_y_id = 3;
1982
1983 // Look at the right boundary of the triangular mesh
1984 unsigned n_right_boundary_node =
1985 bulk_mesh_pt->nboundary_node(right_boundary_id);
1986
1987 // Create a vector of ordered boundary nodes
1988 Vector<Node*> ordered_right_boundary_node_pt(n_right_boundary_node);
1989
1990 // Fill the vector with the nodes on the respective boundary
1991 for (unsigned n = 0; n < n_right_boundary_node; n++)
1992 {
1993 ordered_right_boundary_node_pt[n] =
1994 bulk_mesh_pt->boundary_node_pt(right_boundary_id, n);
1995 }
1996
1997 // Sort them from lowest to highest (in y coordinate)
1998 std::sort(ordered_right_boundary_node_pt.begin(),
1999 ordered_right_boundary_node_pt.end(),
2001
2002 /// Number of elements and boundary nodes to be acted upon during
2003 /// construction are extracted from the 'parent' PML meshes
2004 unsigned n_x_right_pml =
2005 pml_right_mesh_pt->nboundary_element(parent_boundary_x_id);
2006 unsigned n_y_top_pml =
2007 pml_top_mesh_pt->nboundary_element(parent_boundary_y_id);
2008 unsigned n_x_boundary_nodes =
2009 pml_right_mesh_pt->nboundary_node(parent_boundary_x_id);
2010 unsigned n_y_boundary_nodes =
2011 pml_top_mesh_pt->nboundary_node(parent_boundary_y_id);
2012
2013 /// Specific PML sizes needed, taken directly from physical domain
2014 /// and existing PML meshes
2015 double l_pml_right_x_start =
2016 ordered_right_boundary_node_pt[n_right_boundary_node - 1]->x(0);
2017 double l_pml_right_x_end =
2018 pml_right_mesh_pt
2019 ->boundary_node_pt(parent_boundary_x_id, n_x_boundary_nodes - 1)
2020 ->x(0);
2021 double l_pml_top_y_start =
2022 ordered_right_boundary_node_pt[n_right_boundary_node - 1]->x(1);
2023 double l_pml_top_y_end =
2024 pml_top_mesh_pt
2025 ->boundary_node_pt(parent_boundary_y_id, n_y_boundary_nodes - 1)
2026 ->x(1);
2027
2028 // Create the mesh to be designated to the PML
2029 Mesh* pml_top_right_mesh_pt = 0;
2030
2031 // Build the top right corner PML mesh
2032 pml_top_right_mesh_pt =
2034 pml_right_mesh_pt,
2035 pml_top_mesh_pt,
2036 bulk_mesh_pt,
2037 ordered_right_boundary_node_pt[n_right_boundary_node - 1],
2038 parent_boundary_x_id,
2039 parent_boundary_y_id,
2040 current_boundary_x_id,
2041 current_boundary_y_id,
2042 n_x_right_pml,
2043 n_y_top_pml,
2044 l_pml_right_x_start,
2045 l_pml_right_x_end,
2046 l_pml_top_y_start,
2047 l_pml_top_y_end,
2048 time_stepper_pt);
2049
2050 // Enable PML damping on the entire mesh
2051 /// The enabling must be perfromed in both x- and y-directions
2052 /// as this is a corner PML mesh
2053 unsigned n_element_pml_top_right = pml_top_right_mesh_pt->nelement();
2054 for (unsigned e = 0; e < n_element_pml_top_right; e++)
2055 {
2056 // Upcast
2057 PMLElementBase<2>* el_pt = dynamic_cast<PMLElementBase<2>*>(
2058 pml_top_right_mesh_pt->element_pt(e));
2059 el_pt->enable_pml(0, l_pml_right_x_start, l_pml_right_x_end);
2060 el_pt->enable_pml(1, l_pml_top_y_start, l_pml_top_y_end);
2061 }
2062
2063 // Get the values to be pinned from the first element (which we
2064 // assume exists!)
2065 PMLElementBase<2>* el_pt =
2066 dynamic_cast<PMLElementBase<2>*>(pml_top_right_mesh_pt->element_pt(0));
2067 Vector<unsigned> values_to_pin;
2068 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
2069 unsigned npin = values_to_pin.size();
2070
2071 // Exterior boundary needs to be set to Dirichlet 0
2072 // for both real and imaginary parts of all fields
2073 // in the problem
2074 unsigned n_bound_pml_top_right = pml_top_right_mesh_pt->nboundary();
2075 for (unsigned b = 0; b < n_bound_pml_top_right; b++)
2076 {
2077 unsigned n_node = pml_top_right_mesh_pt->nboundary_node(b);
2078 for (unsigned n = 0; n < n_node; n++)
2079 {
2080 Node* nod_pt = pml_top_right_mesh_pt->boundary_node_pt(b, n);
2081 if ((b == 1) || (b == 2))
2082 {
2083 for (unsigned j = 0; j < npin; j++)
2084 {
2085 unsigned j_val = values_to_pin[j];
2086 nod_pt->pin(j_val);
2087 nod_pt->set_value(j_val, 0.0);
2088 }
2089 }
2090 }
2091 }
2092
2093 /// Return the finalized mesh, with PML enabled
2094 /// and boundary conditions added
2095 return pml_top_right_mesh_pt;
2096 }
2097
2098 //==========================================================================
2099 /// "Constructor" for PML bottom right corner mesh,
2100 /// aligned with the existing PML meshes
2101 //==========================================================================
2102 template<class ASSOCIATED_PML_QUAD_ELEMENT>
2104 Mesh* pml_right_mesh_pt,
2105 Mesh* pml_bottom_mesh_pt,
2106 Mesh* bulk_mesh_pt,
2107 const unsigned& right_boundary_id,
2108 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2109 {
2110 /// Relevant boundary id's to be used in construction
2111 /// Parent id refers to already existing PML meshes
2112 unsigned parent_boundary_x_id = 0;
2113 unsigned parent_boundary_y_id = 1;
2114 // Current id refers to the mesh that is to be constructed
2115 unsigned current_boundary_x_id = 2;
2116 unsigned current_boundary_y_id = 3;
2117
2118 // Look at the right boundary of the triangular mesh
2119 unsigned n_right_boundary_node =
2120 bulk_mesh_pt->nboundary_node(right_boundary_id);
2121
2122 // Create a vector of ordered boundary nodes
2123 Vector<Node*> ordered_right_boundary_node_pt(n_right_boundary_node);
2124
2125 // Fill the vector with the nodes on the respective boundary
2126 for (unsigned n = 0; n < n_right_boundary_node; n++)
2127 {
2128 ordered_right_boundary_node_pt[n] =
2129 bulk_mesh_pt->boundary_node_pt(right_boundary_id, n);
2130 }
2131
2132 // Sort them from lowest to highest (in y coordinate)
2133 std::sort(ordered_right_boundary_node_pt.begin(),
2134 ordered_right_boundary_node_pt.end(),
2136
2137 /// Number of elements and boundary nodes to be acted upon during
2138 /// construction are extracted from the 'parent' PML meshes
2139 unsigned n_x_right_pml =
2140 pml_right_mesh_pt->nboundary_element(parent_boundary_x_id);
2141 unsigned n_y_bottom_pml =
2142 pml_bottom_mesh_pt->nboundary_element(parent_boundary_y_id);
2143 unsigned n_x_boundary_nodes =
2144 pml_right_mesh_pt->nboundary_node(parent_boundary_x_id);
2145
2146 /// Specific PML sizes needed, taken directly from physical domain
2147 /// and existing PML meshes
2148 double l_pml_right_x_start = ordered_right_boundary_node_pt[0]->x(0);
2149 double l_pml_right_x_end =
2150 pml_right_mesh_pt
2151 ->boundary_node_pt(parent_boundary_x_id, n_x_boundary_nodes - 1)
2152 ->x(0);
2153 double l_pml_bottom_y_start =
2154 pml_bottom_mesh_pt->boundary_node_pt(parent_boundary_y_id, 0)->x(1);
2155 double l_pml_bottom_y_end = ordered_right_boundary_node_pt[0]->x(1);
2156
2157 // Create the mesh to be designated to the PML
2158 Mesh* pml_bottom_right_mesh_pt = 0;
2159
2160 // Build the bottom right corner PML mesh
2161 pml_bottom_right_mesh_pt =
2163 pml_right_mesh_pt,
2164 pml_bottom_mesh_pt,
2165 bulk_mesh_pt,
2166 ordered_right_boundary_node_pt[0],
2167 parent_boundary_x_id,
2168 parent_boundary_y_id,
2169 current_boundary_x_id,
2170 current_boundary_y_id,
2171 n_x_right_pml,
2172 n_y_bottom_pml,
2173 l_pml_right_x_start,
2174 l_pml_right_x_end,
2175 l_pml_bottom_y_start,
2176 l_pml_bottom_y_end,
2177 time_stepper_pt);
2178
2179 // Enable PML damping on the entire mesh
2180 /// The enabling must be perfromed in both x- and y-directions
2181 /// as this is a corner PML mesh
2182 unsigned n_element_pml_bottom_right =
2183 pml_bottom_right_mesh_pt->nelement();
2184
2185 for (unsigned e = 0; e < n_element_pml_bottom_right; e++)
2186 {
2187 // Upcast
2188 PMLElementBase<2>* el_pt = dynamic_cast<PMLElementBase<2>*>(
2189 pml_bottom_right_mesh_pt->element_pt(e));
2190 el_pt->enable_pml(0, l_pml_right_x_start, l_pml_right_x_end);
2191 el_pt->enable_pml(1, l_pml_bottom_y_end, l_pml_bottom_y_start);
2192 }
2193
2194 // Get the values to be pinned from the first element (which we
2195 // assume exists!)
2196 PMLElementBase<2>* el_pt = dynamic_cast<PMLElementBase<2>*>(
2197 pml_bottom_right_mesh_pt->element_pt(0));
2198 Vector<unsigned> values_to_pin;
2199 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
2200 unsigned npin = values_to_pin.size();
2201
2202 // Exterior boundary needs to be set to Dirichlet 0
2203 // for both real and imaginary parts of all fields
2204 // in the problem
2205 unsigned n_bound_pml_bottom_right = pml_bottom_right_mesh_pt->nboundary();
2206 for (unsigned b = 0; b < n_bound_pml_bottom_right; b++)
2207 {
2208 unsigned n_node = pml_bottom_right_mesh_pt->nboundary_node(b);
2209 for (unsigned n = 0; n < n_node; n++)
2210 {
2211 Node* nod_pt = pml_bottom_right_mesh_pt->boundary_node_pt(b, n);
2212 if ((b == 0) || (b == 1))
2213 {
2214 for (unsigned j = 0; j < npin; j++)
2215 {
2216 unsigned j_val = values_to_pin[j];
2217 nod_pt->pin(j_val);
2218 nod_pt->set_value(j_val, 0.0);
2219 }
2220 }
2221 }
2222 }
2223
2224 /// Return the finalized mesh, with PML enabled
2225 /// and boundary conditions added
2226 return pml_bottom_right_mesh_pt;
2227 }
2228
2229 //==========================================================================
2230 /// "Constructor" for PML top left corner mesh,
2231 /// aligned with the existing PML meshes
2232 //==========================================================================
2233 template<class ASSOCIATED_PML_QUAD_ELEMENT>
2235 Mesh* pml_left_mesh_pt,
2236 Mesh* pml_top_mesh_pt,
2237 Mesh* bulk_mesh_pt,
2238 const unsigned& left_boundary_id,
2239 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2240 {
2241 /// Relevant boundary id's to be used in construction
2242 /// Parent id refers to already existing PML meshes
2243 unsigned parent_boundary_x_id = 2;
2244 unsigned parent_boundary_y_id = 3;
2245 // Current id refers to the mesh that is to be constructed
2246 unsigned current_boundary_x_id = 0;
2247 unsigned current_boundary_y_id = 1;
2248
2249 // Look at the left boundary of the triangular mesh
2250 unsigned n_left_boundary_node =
2251 bulk_mesh_pt->nboundary_node(left_boundary_id);
2252
2253 // Create a vector of ordered boundary nodes
2254 Vector<Node*> ordered_left_boundary_node_pt(n_left_boundary_node);
2255
2256 // Fill the vector with the nodes on the respective boundary
2257 for (unsigned n = 0; n < n_left_boundary_node; n++)
2258 {
2259 ordered_left_boundary_node_pt[n] =
2260 bulk_mesh_pt->boundary_node_pt(left_boundary_id, n);
2261 }
2262
2263 /// Sort them from lowest to highest (in y coordinate)
2264 /// sorter_right_boundary is still functional, as the sorting
2265 /// is performed by the same criterion
2266 std::sort(ordered_left_boundary_node_pt.begin(),
2267 ordered_left_boundary_node_pt.end(),
2269
2270 /// Number of elements and boundary nodes to be acted upon during
2271 /// construction are extracted from the 'parent' PML meshes
2272 unsigned n_x_left_pml =
2273 pml_left_mesh_pt->nboundary_element(parent_boundary_x_id);
2274 unsigned n_y_top_pml =
2275 pml_top_mesh_pt->nboundary_element(parent_boundary_y_id);
2276 unsigned n_y_boundary_nodes =
2277 pml_top_mesh_pt->nboundary_node(parent_boundary_y_id);
2278
2279 /// Specific PML sizes needed, taken directly from physical domain
2280 /// and existing PML meshes
2281 double l_pml_left_x_start =
2282 pml_left_mesh_pt->boundary_node_pt(parent_boundary_x_id, 0)->x(0);
2283 double l_pml_left_x_end =
2284 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(0);
2285 double l_pml_top_y_start =
2286 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(1);
2287 double l_pml_top_y_end =
2288 pml_top_mesh_pt
2289 ->boundary_node_pt(parent_boundary_y_id, n_y_boundary_nodes - 1)
2290 ->x(1);
2291
2292 // Create the mesh to be designated to the PML
2293 Mesh* pml_top_left_mesh_pt = 0;
2294
2295 // Build the top left corner PML mesh
2296 pml_top_left_mesh_pt = new PMLCornerQuadMesh<ASSOCIATED_PML_QUAD_ELEMENT>(
2297 pml_left_mesh_pt,
2298 pml_top_mesh_pt,
2299 bulk_mesh_pt,
2300 ordered_left_boundary_node_pt[n_left_boundary_node - 1],
2301 parent_boundary_x_id,
2302 parent_boundary_y_id,
2303 current_boundary_x_id,
2304 current_boundary_y_id,
2305 n_x_left_pml,
2306 n_y_top_pml,
2307 l_pml_left_x_start,
2308 l_pml_left_x_end,
2309 l_pml_top_y_start,
2310 l_pml_top_y_end,
2311 time_stepper_pt);
2312
2313 // Enable PML damping on the entire mesh
2314 /// The enabling must be perfromed in both x- and y-directions
2315 /// as this is a corner PML mesh
2316 unsigned n_element_pml_top_left = pml_top_left_mesh_pt->nelement();
2317
2318 for (unsigned e = 0; e < n_element_pml_top_left; e++)
2319 {
2320 // Upcast
2321 PMLElementBase<2>* el_pt =
2322 dynamic_cast<PMLElementBase<2>*>(pml_top_left_mesh_pt->element_pt(e));
2323 el_pt->enable_pml(0, l_pml_left_x_end, l_pml_left_x_start);
2324 el_pt->enable_pml(1, l_pml_top_y_start, l_pml_top_y_end);
2325 }
2326
2327 // Get the values to be pinned from the first element (which we
2328 // assume exists!)
2329 PMLElementBase<2>* el_pt =
2330 dynamic_cast<PMLElementBase<2>*>(pml_top_left_mesh_pt->element_pt(0));
2331 Vector<unsigned> values_to_pin;
2332 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
2333 unsigned npin = values_to_pin.size();
2334
2335 // Exterior boundary needs to be set to Dirichlet 0
2336 // for both real and imaginary parts of all fields
2337 // in the problem
2338 unsigned n_bound_pml_top_left = pml_top_left_mesh_pt->nboundary();
2339 for (unsigned b = 0; b < n_bound_pml_top_left; b++)
2340 {
2341 unsigned n_node = pml_top_left_mesh_pt->nboundary_node(b);
2342 for (unsigned n = 0; n < n_node; n++)
2343 {
2344 Node* nod_pt = pml_top_left_mesh_pt->boundary_node_pt(b, n);
2345 if ((b == 2) || (b == 3))
2346 {
2347 for (unsigned j = 0; j < npin; j++)
2348 {
2349 unsigned j_val = values_to_pin[j];
2350 nod_pt->pin(j_val);
2351 nod_pt->set_value(j_val, 0.0);
2352 }
2353 }
2354 }
2355 }
2356
2357 /// Return the finalized mesh, with PML enabled
2358 /// and boundary conditions added
2359 return pml_top_left_mesh_pt;
2360 }
2361
2362 //==========================================================================
2363 /// "Constructor" for PML bottom left corner mesh,
2364 /// aligned with the existing PML meshes
2365 //==========================================================================
2366 template<class ASSOCIATED_PML_QUAD_ELEMENT>
2368 Mesh* pml_left_mesh_pt,
2369 Mesh* pml_bottom_mesh_pt,
2370 Mesh* bulk_mesh_pt,
2371 const unsigned& left_boundary_id,
2372 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2373 {
2374 /// Relevant boundary id's to be used in construction
2375 /// Parent id refers to already existing PML meshes
2376 unsigned parent_boundary_x_id = 0;
2377 unsigned parent_boundary_y_id = 3;
2378 // Current id refers to the mesh that is to be constructed
2379 unsigned current_boundary_x_id = 2;
2380 unsigned current_boundary_y_id = 1;
2381
2382 // Look at the left boundary of the triangular mesh
2383 unsigned n_left_boundary_node =
2384 bulk_mesh_pt->nboundary_node(left_boundary_id);
2385
2386 // Create a vector of ordered boundary nodes
2387 Vector<Node*> ordered_left_boundary_node_pt(n_left_boundary_node);
2388
2389 // Fill the vector with the nodes on the respective boundary
2390 for (unsigned n = 0; n < n_left_boundary_node; n++)
2391 {
2392 ordered_left_boundary_node_pt[n] =
2393 bulk_mesh_pt->boundary_node_pt(left_boundary_id, n);
2394 }
2395
2396 /// Sort them from lowest to highest (in y coordinate)
2397 /// sorter_right_boundary is still functional, as the sorting
2398 /// is performed by the same criterion
2399 std::sort(ordered_left_boundary_node_pt.begin(),
2400 ordered_left_boundary_node_pt.end(),
2402
2403 /// Number of elements and boundary nodes to be acted upon during
2404 /// construction are extracted from the 'parent' PML meshes
2405 unsigned n_x_left_pml =
2406 pml_left_mesh_pt->nboundary_element(parent_boundary_x_id);
2407 unsigned n_y_bottom_pml =
2408 pml_bottom_mesh_pt->nboundary_element(parent_boundary_y_id);
2409
2410 /// Specific PML sizes needed, taken directly from physical domain
2411 /// and existing PML meshes
2412 double l_pml_left_x_start =
2413 pml_left_mesh_pt->boundary_node_pt(parent_boundary_x_id, 0)->x(0);
2414 double l_pml_left_x_end =
2415 ordered_left_boundary_node_pt[n_left_boundary_node - 1]->x(0);
2416 double l_pml_bottom_y_start =
2417 pml_bottom_mesh_pt->boundary_node_pt(parent_boundary_y_id, 0)->x(1);
2418 double l_pml_bottom_y_end = ordered_left_boundary_node_pt[0]->x(1);
2419
2420 // Create the mesh to be designated to the PML
2421 Mesh* pml_bottom_left_mesh_pt = 0;
2422
2423 // Build the bottom left corner PML mesh
2424 pml_bottom_left_mesh_pt =
2426 pml_left_mesh_pt,
2427 pml_bottom_mesh_pt,
2428 bulk_mesh_pt,
2429 ordered_left_boundary_node_pt[0],
2430 parent_boundary_x_id,
2431 parent_boundary_y_id,
2432 current_boundary_x_id,
2433 current_boundary_y_id,
2434 n_x_left_pml,
2435 n_y_bottom_pml,
2436 l_pml_left_x_start,
2437 l_pml_left_x_end,
2438 l_pml_bottom_y_start,
2439 l_pml_bottom_y_end,
2440 time_stepper_pt);
2441
2442 // Enable PML damping on the entire mesh
2443 /// The enabling must be perfromed in both x- and y-directions
2444 /// as this is a corner PML mesh
2445 unsigned n_element_pml_bottom_left = pml_bottom_left_mesh_pt->nelement();
2446 for (unsigned e = 0; e < n_element_pml_bottom_left; e++)
2447 {
2448 // Upcast
2449 PMLElementBase<2>* el_pt = dynamic_cast<PMLElementBase<2>*>(
2450 pml_bottom_left_mesh_pt->element_pt(e));
2451 el_pt->enable_pml(0, l_pml_left_x_end, l_pml_left_x_start);
2452 el_pt->enable_pml(1, l_pml_bottom_y_end, l_pml_bottom_y_start);
2453 }
2454
2455 // Get the values to be pinned from the first element (which we
2456 // assume exists!)
2457 PMLElementBase<2>* el_pt = dynamic_cast<PMLElementBase<2>*>(
2458 pml_bottom_left_mesh_pt->element_pt(0));
2459 Vector<unsigned> values_to_pin;
2460 el_pt->values_to_be_pinned_on_outer_pml_boundary(values_to_pin);
2461 unsigned npin = values_to_pin.size();
2462
2463 // Exterior boundary needs to be set to Dirichlet 0
2464 // for both real and imaginary parts of all fields
2465 // in the problem
2466 unsigned n_bound_pml_bottom_left = pml_bottom_left_mesh_pt->nboundary();
2467 for (unsigned b = 0; b < n_bound_pml_bottom_left; b++)
2468 {
2469 unsigned n_node = pml_bottom_left_mesh_pt->nboundary_node(b);
2470 for (unsigned n = 0; n < n_node; n++)
2471 {
2472 Node* nod_pt = pml_bottom_left_mesh_pt->boundary_node_pt(b, n);
2473 if ((b == 0) || (b == 3))
2474 {
2475 for (unsigned j = 0; j < npin; j++)
2476 {
2477 unsigned j_val = values_to_pin[j];
2478 nod_pt->pin(j_val);
2479 nod_pt->set_value(j_val, 0.0);
2480 }
2481 }
2482 }
2483 }
2484
2485 /// Return the finalized mesh, with PML enabled
2486 /// and boundary conditions added
2487 return pml_bottom_left_mesh_pt;
2488 }
2489
2490 } // namespace TwoDimensionalPMLHelper
2491
2492 /// ///////////////////////////////////////////////////////////////////
2493 /// ///////////////////////////////////////////////////////////////////
2494 /// ///////////////////////////////////////////////////////////////////
2495
2496} // namespace oomph
2497
2498#endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
A general Finite Element class.
Definition: elements.h:1313
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
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
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
Vector< Node * > prune_dead_nodes()
Prune nodes. Nodes that have been marked as obsolete are removed from the mesh (and its boundary-node...
Definition: mesh.cc:966
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
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
void set_obsolete()
Mark node as obsolete.
Definition: nodes.h:1436
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: pml_meshes.h:777
PMLCornerQuadMesh(Mesh *PMLQuad_mesh_x_pt, Mesh *PMLQuad_mesh_y_pt, Mesh *bulk_mesh_pt, Node *special_corner_node_pt, const unsigned &parent_boundary_x_id, const unsigned &parent_boundary_y_id, const unsigned &current_boundary_x_id, const unsigned &current_boundary_y_id, const unsigned &n_pml_x, const unsigned &n_pml_y, const double &x_pml_start, const double &x_pml_end, const double &y_pml_start, const double &y_pml_end, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to "bulk" mesh and the two existing PML meshes in order to construct the co...
Definition: pml_meshes.h:783
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: pml_meshes.h:60
bool Pml_is_enabled
Boolean indicating if element is used in pml mode.
Definition: pml_meshes.h:119
virtual ~PMLElementBase()
Virtual destructor.
Definition: pml_meshes.h:72
void enable_pml(const int &direction, const double &interface_border_value, const double &outer_domain_border_value)
Enable pml. Specify the coordinate direction along which pml boundary is constant,...
Definition: pml_meshes.h:97
void disable_pml()
Disable pml. Ensures the PML-ification in all directions has been deactivated.
Definition: pml_meshes.h:76
std::vector< bool > Pml_direction_active
Coordinate direction along which pml boundary is constant; alternatively: coordinate direction in whi...
Definition: pml_meshes.h:124
Vector< double > Pml_outer_boundary
Coordinate of outer pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:134
Vector< double > Pml_inner_boundary
Coordinate of inner pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:129
PMLElementBase()
Constructor.
Definition: pml_meshes.h:63
virtual void values_to_be_pinned_on_outer_pml_boundary(Vector< unsigned > &values_to_pin)=0
Pure virtual function in which we have to specify the values to be pinned (and set to zero) on the ou...
General definition of policy class defining the elements to be used in the actual PML layers....
Definition: pml_meshes.h:48
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
Definition: pml_meshes.h:171
virtual void pml_locate_zeta(const Vector< double > &x, FiniteElement *&el_pt)=0
Pure virtual function to provide an optimised version of the locate_zeta function for PML meshes....
PML mesh class. Policy class for 2D PML meshes.
Definition: pml_meshes.h:187
void pml_locate_zeta(const Vector< double > &x, FiniteElement *&coarse_mesh_el_pt)
Overloaded function to allow the user to locate an element in mesh given the (Eulerian) position of a...
Definition: pml_meshes.h:212
PMLQuadMeshBase(const unsigned &n_pml_x, const unsigned &n_pml_y, const double &x_pml_start, const double &x_pml_end, const double &y_pml_start, const double &y_pml_end, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Create the underlying rectangular quad mesh.
Definition: pml_meshes.h:190
PML mesh, derived from RectangularQuadMesh.
Definition: pml_meshes.h:431
PMLQuadMesh(Mesh *bulk_mesh_pt, const unsigned &boundary_id, const unsigned &quad_boundary_id, const unsigned &n_pml_x, const unsigned &n_pml_y, const double &x_pml_start, const double &x_pml_end, const double &y_pml_start, const double &y_pml_end, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to "bulk" mesh, the boundary ID of axis aligned boundary to which the mesh ...
Definition: pml_meshes.h:441
RectangularQuadMesh is a two-dimensional mesh of Quad elements with Nx elements in the "x" (horizonal...
const unsigned & ny() const
Return number of elements in y direction.
const double y_min() const
Return the minimum value of y coordinate.
const double x_max() const
Return the maximum value of x coordinate.
const double y_max() const
Return the maximum value of y coordinate.
const unsigned & nx() const
Return number of elements in x direction.
const double x_min() const
Return the minimum value of x coordinate.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
bool sorter_bottom_boundary(Node *nod_i_pt, Node *nod_j_pt)
helper function for sorting the bottom boundary nodes
Definition: pml_meshes.cc:57
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_bottom_pml_mesh(Mesh *bulk_mesh_pt, const unsigned &bottom_boundary_id, const unsigned &n_y_bottom_pml, const double &width_y_bottom_pml, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML mesh,aligned with the bottom physical domain boundary
Definition: pml_meshes.h:1860
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
bool sorter_top_boundary(Node *nod_i_pt, Node *nod_j_pt)
helper function for sorting the top boundary nodes
Definition: pml_meshes.cc:45
Mesh * create_right_pml_mesh(Mesh *bulk_mesh_pt, const unsigned &right_boundary_id, const unsigned &n_x_right_pml, const double &width_x_right_pml, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML mesh,aligned with the right physical domain boundary
Definition: pml_meshes.h:1539
Mesh * create_left_pml_mesh(Mesh *bulk_mesh_pt, const unsigned &left_boundary_id, const unsigned &n_x_left_pml, const double &width_x_left_pml, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML mesh, aligned with the left physical domain boundary
Definition: pml_meshes.h:1752
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
bool sorter_right_boundary(Node *nod_i_pt, Node *nod_j_pt)
helper function for sorting the right boundary nodes
Definition: pml_meshes.cc:39
Mesh * create_top_pml_mesh(Mesh *bulk_mesh_pt, const unsigned &top_boundary_id, const unsigned &n_y_top_pml, const double &width_y_top_pml, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
"Constructor" for PML mesh, aligned with the top physical domain boundary
Definition: pml_meshes.h:1646
bool sorter_left_boundary(Node *nod_i_pt, Node *nod_j_pt)
helper function for sorting the left boundary nodes
Definition: pml_meshes.cc:51
//////////////////////////////////////////////////////////////////// ////////////////////////////////...