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