unstructured_two_d_helmholtz.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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 // Driver for a specific 2D Helmholtz problem with
31 // perfectly matched layer treatment for the exterior boundaries
32 
33 #include<fenv.h>
34 
35 #include "math.h"
36 #include <complex>
37 
38 // Generic routines
39 #include "generic.h"
40 
41 // The pml Helmholtz equations
42 #include "pml_helmholtz.h"
43 
44 // The meshes needed in the PML constructions
45 #include "meshes/triangle_mesh.h"
46 #include "meshes/rectangular_quadmesh.h"
47 
48 using namespace oomph;
49 using namespace std;
50 
51 /////////////////////////////////////////////////////////////////////
52 /////////////////////////////////////////////////////////////////////
53 /////////////////////////////////////////////////////////////////////
54 
55 //===== start_of_namespace=============================================
56 /// Namespace for the Helmholtz problem parameters
57 //=====================================================================
58 namespace GlobalParameters
59 {
60 
61  /// Wavenumber (also known as k), k=omega/c
62  double Wavenumber = sqrt(50.0);
63 
64  /// Square of the wavenumber (also known as k^2)
66 
67 } // end of namespace
68 
69 
70 /////////////////////////////////////////////////////////////////////
71 /////////////////////////////////////////////////////////////////////
72 /////////////////////////////////////////////////////////////////////
73 
74 //========= start_of_problem_class=====================================
75 /// Problem class to demonstrate use of perfectly matched layers
76 /// for Helmholtz problems.
77 //=====================================================================
78 template<class ELEMENT>
79 class PMLProblem : public Problem
80 {
81 
82 public:
83 
84  /// Constructor
85  PMLProblem();
86 
87  /// Destructor (empty)
89 
90  /// Update the problem specs before solve (empty)
92 
93  /// Update the problem specs after solve (empty)
95 
96  /// \short Doc the solution. DocInfo object stores flags/labels for where the
97  /// output gets written to
98  void doc_solution(DocInfo& doc_info);
99 
100  /// Create PML meshes
101  void create_pml_meshes();
102 
103  // Apply boundary conditions
104  void apply_boundary_conditions();
105 
106 #ifdef ADAPTIVE
107 
108  /// Actions before adapt: Wipe the PML meshes
109  void actions_before_adapt();
110 
111  /// Actions after adapt: Rebuild the PML meshes
112  void actions_after_adapt();
113 
114 #endif
115 
116 
117 #ifdef ADAPTIVE
118 
119 private:
120 
121  /// Pointer to the refineable "bulk" mesh
122  RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
123 
124 #else
125 
126 private:
127 
128  /// Pointer to the "bulk" mesh
129  TriangleMesh<ELEMENT>* Bulk_mesh_pt;
130 
131 #endif
132 
133 
134  /// Pointer to the right PML mesh
136 
137  /// Pointer to the top PML mesh
139 
140  /// Pointer to the left PML mesh
142 
143  /// Pointer to the bottom PML mesh
145 
146  /// Pointer to the top right corner PML mesh
148 
149  /// Pointer to the top left corner PML mesh
151 
152  /// Pointer to the bottom right corner PML mesh
154 
155  /// Pointer to the bottom left corner PML mesh
157 
158  /// Trace file
159  ofstream Trace_file;
160 
161 }; // end of problem class
162 
163 
164 
165 //=======start_of_constructor=============================================
166 /// Constructor for Helmholtz problem
167 //========================================================================
168 template<class ELEMENT>
170 {
171 
172  // Open trace file
173  Trace_file.open("RESLT/trace.dat");
174 
175  // Create circle representing inner boundary
176  double a=0.2;
177  double x_c=0.0;
178  double y_c=0.0;
179  Circle* inner_circle_pt=new Circle(x_c,y_c,a);
180 
181  // Outer boundary
182  //---------------
183  TriangleMeshClosedCurve* outer_boundary_pt=0;
184 
185  unsigned n_segments = 20;
186  Vector<TriangleMeshCurveSection*> outer_boundary_line_pt(4);
187 
188  // Each polyline only has three vertices, provide storage for their
189  // coordinates
190  Vector<Vector<double> > vertex_coord(2);
191  for(unsigned i=0;i<2;i++)
192  {
193  vertex_coord[i].resize(2);
194  }
195 
196  // First polyline
197  vertex_coord[0][0]=-2.0;
198  vertex_coord[0][1]=-2.0;
199  vertex_coord[1][0]=-2.0;
200  vertex_coord[1][1]=2.0;
201 
202  // Build the 1st boundary polyline
203  unsigned boundary_id=2;
204  outer_boundary_line_pt[0] = new TriangleMeshPolyLine(vertex_coord,
205  boundary_id);
206 
207  // Second boundary polyline
208  vertex_coord[0][0]=-2.0;
209  vertex_coord[0][1]=2.0;
210  vertex_coord[1][0]=2.0;
211  vertex_coord[1][1]=2.0;
212 
213  // Build the 2nd boundary polyline
214  boundary_id=3;
215  outer_boundary_line_pt[1] = new TriangleMeshPolyLine(vertex_coord,
216  boundary_id);
217 
218  // Third boundary polyline
219  vertex_coord[0][0]=2.0;
220  vertex_coord[0][1]=2.0;
221  vertex_coord[1][0]=2.0;
222  vertex_coord[1][1]=-2.0;
223 
224  // Build the 3rd boundary polyline
225  boundary_id=4;
226  outer_boundary_line_pt[2] = new TriangleMeshPolyLine(vertex_coord,
227  boundary_id);
228 
229  // Fourth boundary polyline
230  vertex_coord[0][0]=2.0;
231  vertex_coord[0][1]=-2.0;
232  vertex_coord[1][0]=-2.0;
233  vertex_coord[1][1]=-2.0;
234 
235  // Build the 4th boundary polyline
236  boundary_id=5;
237  outer_boundary_line_pt[3] = new TriangleMeshPolyLine(vertex_coord,
238  boundary_id);
239 
240  // Create the triangle mesh polygon for outer boundary
241  outer_boundary_pt = new TriangleMeshPolygon(outer_boundary_line_pt);
242 
243  // Inner circular boundary
244  //------------------------
245  Vector<TriangleMeshCurveSection*> inner_boundary_line_pt(2);
246 
247  // The intrinsic coordinates for the beginning and end of the curve
248  double s_start = 0.0;
249  double s_end = MathematicalConstants::Pi;
250  boundary_id = 0;
251  inner_boundary_line_pt[0]=
252  new TriangleMeshCurviLine(inner_circle_pt,
253  s_start,
254  s_end,
255  n_segments,
256  boundary_id);
257 
258  // The intrinsic coordinates for the beginning and end of the curve
259  s_start = MathematicalConstants::Pi;
260  s_end = 2.0*MathematicalConstants::Pi;
261  boundary_id = 1;
262  inner_boundary_line_pt[1]=
263  new TriangleMeshCurviLine(inner_circle_pt,
264  s_start,
265  s_end,
266  n_segments,
267  boundary_id);
268 
269 
270  // Combine to hole
271  //----------------
272  Vector<TriangleMeshClosedCurve*> hole_pt(1);
273  Vector<double> hole_coords(2);
274  hole_coords[0]=0.0;
275  hole_coords[1]=0.0;
276  hole_pt[0]=new TriangleMeshClosedCurve(inner_boundary_line_pt,
277  hole_coords);
278 
279 
280  // Use the TriangleMeshParameters object for helping on the manage
281  // of the TriangleMesh parameters. The only parameter that needs to take
282  // is the outer boundary.
283  TriangleMeshParameters triangle_mesh_parameters(outer_boundary_pt);
284 
285  // Specify the closed curve using the TriangleMeshParameters object
286  triangle_mesh_parameters.internal_closed_curve_pt() = hole_pt;
287 
288  // Target element size in bulk mesh
289  triangle_mesh_parameters.element_area() = 0.1;
290 
291 #ifdef ADAPTIVE
292 
293  // Build adaptive "bulk" mesh
294  Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
295 
296  // Create/set error estimator
297  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
298 
299  // Choose error tolerances to force some uniform refinement
300  Bulk_mesh_pt->min_permitted_error()=0.00004;
301  Bulk_mesh_pt->max_permitted_error()=0.0001;
302 
303 #else
304 
305  // Build "bulk" mesh
306  Bulk_mesh_pt=new TriangleMesh<ELEMENT>(triangle_mesh_parameters);
307 
308 #endif
309 
310  // Create the main triangular mesh
311  add_sub_mesh(Bulk_mesh_pt);
312 
313  // Create PML meshes and add them to the global mesh
314  create_pml_meshes();
315 
316  // Build the entire mesh from its submeshes
317  build_global_mesh();
318 
319  // Let's have a look where the boundaries are
320  this->mesh_pt()->output("global_mesh.dat");
321  this->mesh_pt()->output_boundaries("global_mesh_boundary.dat");
322 
323  // Complete the build of all elements so they are fully functional
324  unsigned n_element = this->mesh_pt()->nelement();
325  for(unsigned e=0;e<n_element;e++)
326  {
327  // Upcast from GeneralisedElement to Helmholtz bulk element
328  PMLHelmholtzEquations<2> *el_pt =
329  dynamic_cast<PMLHelmholtzEquations<2>*>(mesh_pt()->element_pt(e));
330 
331  //Set the k_squared double pointer
332  el_pt->k_squared_pt() = &GlobalParameters::K_squared;
333  }
334 
335  // Apply boundary conditions
336  apply_boundary_conditions();
337 
338  // Setup equation numbering scheme
339  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
340 
341 } // end of constructor
342 
343 
344 #ifdef ADAPTIVE
345 
346 //=====================start_of_actions_before_adapt======================
347 /// Actions before adapt: Wipe the mesh of face elements
348 //========================================================================
349 template<class ELEMENT>
351 {
352  // Before adapting the added PML meshes must be removed
353  // as they are not refineable and are to be rebuilt from the
354  // newly refined triangular mesh
355  delete PML_right_mesh_pt;
356  PML_right_mesh_pt=0;
357  delete PML_top_mesh_pt;
358  PML_top_mesh_pt=0;
359  delete PML_left_mesh_pt;
360  PML_left_mesh_pt=0;
361  delete PML_bottom_mesh_pt;
362  PML_bottom_mesh_pt=0;
363  delete PML_top_right_mesh_pt;
364  PML_top_right_mesh_pt=0;
365  delete PML_top_left_mesh_pt;
366  PML_top_left_mesh_pt=0;
367  delete PML_bottom_right_mesh_pt;
368  PML_bottom_right_mesh_pt=0;
369  delete PML_bottom_left_mesh_pt;
370  PML_bottom_left_mesh_pt=0;
371 
372  // Rebuild the Problem's global mesh from its various sub-meshes
373  // but first flush all its submeshes
374  flush_sub_meshes();
375 
376  // Then add the triangular mesh back
377  add_sub_mesh(Bulk_mesh_pt);
378 
379  // Rebuild the global mesh such that it now stores
380  // the triangular mesh only
381  rebuild_global_mesh();
382 
383 }// end of actions_before_adapt
384 
385 
386 //=====================start_of_actions_after_adapt=======================
387 /// Actions after adapt: Rebuild the face element meshes
388 //========================================================================
389 template<class ELEMENT>
391 {
392 
393  // Build PML meshes and add them to the global mesh
394  create_pml_meshes();
395 
396  // Build the entire mesh from its submeshes
397  rebuild_global_mesh();
398 
399  // Complete the build of all elements so they are fully functional
400 
401  // Loop over the entire mesh elements to set up element-specific
402  // things that cannot be handled by constructor
403  unsigned n_element = this->mesh_pt()->nelement();
404 
405  for(unsigned e=0;e<n_element;e++)
406  {
407  // Upcast from GeneralisedElement to PMLHelmholtz bulk element
408  PMLHelmholtzEquations<2> *el_pt =
409  dynamic_cast<PMLHelmholtzEquations<2>*>(mesh_pt()->element_pt(e));
410 
411  //Set the frequency function pointer
412  el_pt->k_squared_pt() = &GlobalParameters::K_squared;
413  }
414 
415  // Re-apply boundary conditions
416  apply_boundary_conditions();
417 
418 }// end of actions_after_adapt
419 
420 #endif
421 
422 //==================start_of_apply_boundary_conditions====================
423 /// Apply boundary conditions
424 //========================================================================
425 template<class ELEMENT>
427 {
428 
429  // Boundary conditions are set on the surface of the circle
430  // as a constant nonzero Dirichlet boundary condition
431  unsigned n_bound = Bulk_mesh_pt->nboundary();
432 
433  for(unsigned b=0;b<n_bound;b++)
434  {
435  unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
436  for (unsigned n=0;n<n_node;n++)
437  {
438  if ((b==0) || (b==1))
439  {
440  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,n);
441  nod_pt->pin(0);
442  nod_pt->pin(1);
443 
444  nod_pt->set_value(0,0.1);
445  nod_pt->set_value(1,0.0);
446  }
447  }
448  }
449 
450 }// end of apply_boundary_conditions
451 
452 
453 //=====================start_of_doc=======================================
454 /// Doc the solution: doc_info contains labels/output directory etc.
455 //========================================================================
456 template<class ELEMENT>
457 void PMLProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
458 {
459 
460  ofstream some_file,some_file2;
461  char filename[100];
462 
463  // Number of plot points
464  unsigned npts;
465  npts=5;
466 
467  // Output solution
468  //-----------------
469  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
470  doc_info.number());
471  some_file.open(filename);
472  Bulk_mesh_pt->output(some_file,npts);
473  some_file.close();
474 
475  // Output coarse solution
476  //-----------------------
477  sprintf(filename,"%s/coarse_soln%i.dat",doc_info.directory().c_str(),
478  doc_info.number());
479  some_file.open(filename);
480  unsigned npts_coarse=2;
481  Bulk_mesh_pt->output(some_file,npts_coarse);
482  some_file.close();
483 
484 
485  // Output solution within pml domains
486  //-----------------------------------
487  sprintf(filename,"%s/pml_soln%i.dat",doc_info.directory().c_str(),
488  doc_info.number());
489  some_file.open(filename);
490  PML_top_mesh_pt->output(some_file,npts);
491  PML_right_mesh_pt->output(some_file,npts);
492  PML_bottom_mesh_pt->output(some_file,npts);
493  PML_left_mesh_pt->output(some_file,npts);
494  PML_top_right_mesh_pt->output(some_file,npts);
495  PML_bottom_right_mesh_pt->output(some_file,npts);
496  PML_top_left_mesh_pt->output(some_file,npts);
497  PML_bottom_left_mesh_pt->output(some_file,npts);
498  some_file.close();
499 
500 
501 
502 
503  // Write norm of solution to trace file
504  double norm=0.0;
505  Bulk_mesh_pt->compute_norm(norm);
506  Trace_file << norm << std::endl;
507 
508  // //Do animation of Helmholtz solution
509  // //-----------------------------------
510  // unsigned nstep=40;
511  // for (unsigned i=0;i<nstep;i++)
512  // {
513  // sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",
514  // doc_info.directory().c_str(),
515  // doc_info.number(),i);
516  // some_file.open(filename);
517  // double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);
518  // unsigned nelem=Bulk_mesh_pt->nelement();
519  // for (unsigned e=0;e<nelem;e++)
520  // {
521  // ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
522  // Bulk_mesh_pt->element_pt(e));
523  // el_pt->output_real(some_file,phi,npts);
524  // }
525  // some_file.close();
526  // }
527 
528 } // end of doc
529 
530 //============start_of_create_pml_meshes======================================
531 /// Create PML meshes and add them to the problem's sub-meshes
532 //============================================================================
533 template<class ELEMENT>
535 {
536 
537  // Bulk mesh left boundary id
538  unsigned int left_boundary_id = 2;
539 
540  // Bulk mesh top boundary id
541  unsigned int top_boundary_id = 3;
542 
543  // Bulk mesh right boundary id
544  unsigned int right_boundary_id = 4;
545 
546  // Bulk mesh bottom boundary id
547  unsigned int bottom_boundary_id = 5;
548 
549  // PML width in elements for the right layer
550  unsigned n_x_right_pml = 3;
551 
552  // PML width in elements for the top layer
553  unsigned n_y_top_pml = 3;
554 
555  // PML width in elements for the left layer
556  unsigned n_x_left_pml = 3;
557 
558  // PML width in elements for the left layer
559  unsigned n_y_bottom_pml = 3;
560 
561  // Outer physical length of the PML layers
562  // defaults to 0.2, so 10% of the size of the
563  // physical domain
564  double width_x_right_pml = 0.2;
565  double width_y_top_pml = 0.2;
566  double width_x_left_pml = 0.2;
567  double width_y_bottom_pml = 0.2;
568 
569  // Build the PML meshes based on the new adapted triangular mesh
570  PML_right_mesh_pt =
571  TwoDimensionalPMLHelper::create_right_pml_mesh
572  <PMLLayerElement<ELEMENT> >
573  (Bulk_mesh_pt,right_boundary_id,
574  n_x_right_pml, width_x_right_pml);
575  PML_top_mesh_pt =
576  TwoDimensionalPMLHelper::create_top_pml_mesh
577  <PMLLayerElement<ELEMENT> >
578  (Bulk_mesh_pt, top_boundary_id,
579  n_y_top_pml, width_y_top_pml);
580  PML_left_mesh_pt =
581  TwoDimensionalPMLHelper::create_left_pml_mesh
582  <PMLLayerElement<ELEMENT> >
583  (Bulk_mesh_pt, left_boundary_id,
584  n_x_left_pml, width_x_left_pml);
585  PML_bottom_mesh_pt=
586  TwoDimensionalPMLHelper::create_bottom_pml_mesh
587  <PMLLayerElement<ELEMENT> >
588  (Bulk_mesh_pt, bottom_boundary_id,
589  n_y_bottom_pml, width_y_bottom_pml);
590 
591  // Add submeshes to the global mesh
592  add_sub_mesh(PML_right_mesh_pt);
593  add_sub_mesh(PML_top_mesh_pt);
594  add_sub_mesh(PML_left_mesh_pt);
595  add_sub_mesh(PML_bottom_mesh_pt);
596 
597  // Rebuild corner PML meshes
598  PML_top_right_mesh_pt =
599  TwoDimensionalPMLHelper::create_top_right_pml_mesh
600  <PMLLayerElement<ELEMENT> >
601  (PML_right_mesh_pt, PML_top_mesh_pt,
602  Bulk_mesh_pt, right_boundary_id);
603 
604  PML_bottom_right_mesh_pt =
605  TwoDimensionalPMLHelper::create_bottom_right_pml_mesh
606  <PMLLayerElement<ELEMENT> >
607  (PML_right_mesh_pt, PML_bottom_mesh_pt,
608  Bulk_mesh_pt, right_boundary_id);
609 
610  PML_top_left_mesh_pt =
611  TwoDimensionalPMLHelper::create_top_left_pml_mesh
612  <PMLLayerElement<ELEMENT> >
613  (PML_left_mesh_pt, PML_top_mesh_pt,
614  Bulk_mesh_pt, left_boundary_id);
615 
616  PML_bottom_left_mesh_pt =
617  TwoDimensionalPMLHelper::create_bottom_left_pml_mesh
618  <PMLLayerElement<ELEMENT> >
619  (PML_left_mesh_pt, PML_bottom_mesh_pt,
620  Bulk_mesh_pt, left_boundary_id);
621 
622  // Add submeshes to the global mesh
623  add_sub_mesh(PML_top_right_mesh_pt);
624  add_sub_mesh(PML_bottom_right_mesh_pt);
625  add_sub_mesh(PML_top_left_mesh_pt);
626  add_sub_mesh(PML_bottom_left_mesh_pt);
627 
628 } // end of create_pml_meshes
629 
630 
631 
632 //==========start_of_main=================================================
633 /// Solve 2D Helmholtz problem
634 //========================================================================
635 int main(int argc, char **argv)
636 {
637  //Set up the problem
638  //------------------
639 
640 #ifdef ADAPTIVE
641 
642  // Set up the problem with projectable 2D six-node elements from the
643  // TPMLHelmholtzElement family.
644  PMLProblem<ProjectablePMLHelmholtzElement
645  <TPMLHelmholtzElement<2,3> > > problem;
646 
647  // Set up the problem with 2D ten-node elements from the
648  // TPMLHelmholtzElement family.
649  // PMLProblem<ProjectablePMLHelmholtzElement
650  // <TPMLHelmholtzElement<2,4> > > problem;
651 
652 #else
653 
654  // Set up the problem with 2D six-node elements from the
655  // TPMLHelmholtzElement family.
657 
658  // Set up the problem with 2D ten-node elements from the
659  // TPMLHelmholtzElement family.
660  // PMLProblem<TPMLHelmholtzElement<2,4> > problem;
661 
662 #endif
663 
664  // Create label for output
665  //------------------------
666  DocInfo doc_info;
667 
668  // Set output directory
669  doc_info.set_directory("RESLT");
670 
671 
672 #ifdef ADAPTIVE
673 
674  // Max. number of adaptations
675  unsigned max_adapt=1;
676 
677  // Solve the problem with the adaptive Newton's method, allowing
678  // up to max_adapt mesh adaptations after every solve.
679  problem.newton_solve(max_adapt);
680 
681 #else
682 
683  // Solve the problem with Newton's method
684  problem.newton_solve();
685 
686 #endif
687 
688  //Output solution
689  problem.doc_solution(doc_info);
690 
691 } //end of main
double K_squared
Square of the wavenumber (also known as k^2)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void create_pml_meshes()
Create PML meshes.
Mesh * PML_bottom_mesh_pt
Pointer to the bottom PML mesh.
int main(int argc, char **argv)
Solve 2D Helmholtz problem.
double Wavenumber
Wavenumber (also known as k), k=omega/c.
Mesh * PML_right_mesh_pt
Pointer to the right PML mesh.
RefineableTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the refineable "bulk" mesh.
~PMLProblem()
Destructor (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to...
Mesh * PML_left_mesh_pt
Pointer to the left PML mesh.
Mesh * PML_top_right_mesh_pt
Pointer to the top right corner PML mesh.
void actions_after_adapt()
Actions after adapt: Rebuild the PML meshes.
void apply_boundary_conditions()
Apply boundary conditions.
Mesh * PML_bottom_right_mesh_pt
Pointer to the bottom right corner PML mesh.
ofstream Trace_file
Trace file.
Mesh * PML_top_left_mesh_pt
Pointer to the top left corner PML mesh.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_adapt()
Actions before adapt: Wipe the PML meshes.
Mesh * PML_bottom_left_mesh_pt
Pointer to the bottom left corner PML mesh.
Mesh * PML_top_mesh_pt
Pointer to the top PML mesh.
TriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.