unstructured_two_d_helmholtz_scattering.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 
36 #include "math.h"
37 #include <complex>
38 
39 // Generic routines
40 #include "generic.h"
41 
42 // The Helmholtz equations
43 #include "helmholtz.h"
44 
45 // The pml Helmholtz equations
46 #include "pml_helmholtz.h"
47 
48 // The meshes needed in the PML constructions
49 #include "meshes/triangle_mesh.h"
50 #include "meshes/rectangular_quadmesh.h"
51 
52 // Get the Bessel functions
53 #include "oomph_crbond_bessel.h"
54 
55 using namespace oomph;
56 using namespace std;
57 
58 
59 /////////////////////////////////////////////////////////////////////
60 /////////////////////////////////////////////////////////////////////
61 /////////////////////////////////////////////////////////////////////
62 
63 //===== start_of_namespace=============================================
64 /// Namespace for the Helmholtz problem parameters
65 //=====================================================================
66 namespace GlobalParameters
67 {
68 
69  /// Wavenumber (also known as k), k=omega/c
70  double Wavenumber = sqrt(10.0);
71 
72  /// Square of the wavenumber (also known as k^2)
73  double K_squared = Wavenumber * Wavenumber;
74 
75  /// \short Number of terms used in the computation
76  /// of the exact solution
77  unsigned N_fourier=100;
78 
79  /// Imaginary unit
80  std::complex<double> I(0.0,1.0);
81 
82  /// \short Exact solution for scattered field
83  /// (vector returns real and impaginary parts).
84  void get_exact_u(const Vector<double>& x, Vector<double>& u)
85  {
86  // Switch to polar coordinates
87  double r;
88  r=sqrt(x[0]*x[0]+x[1]*x[1]);
89  double theta;
90  theta=atan2(x[1],x[0]);
91 
92  // Argument for Bessel/Hankel functions
93  double rr=Wavenumber*r;
94 
95  // Evaluate Bessel/Hankel functions
96  complex <double > u_ex(0.0,0.0);
97  Vector<double> jn(N_fourier+1), yn(N_fourier+1),
98  jnp(N_fourier+1), ynp(N_fourier+1);
99  Vector<double> jn_a(N_fourier+1),yn_a(N_fourier+1),
100  jnp_a(N_fourier+1), ynp_a(N_fourier+1);
101  Vector<complex<double> > h(N_fourier+1),h_a(N_fourier+1),
102  hp(N_fourier+1), hp_a(N_fourier+1);
103 
104  // We want to compute N_fourier terms but the function
105  // may return fewer than that.
106  int n_actual=0;
107  CRBond_Bessel::bessjyna(N_fourier,Wavenumber,n_actual,
108  &jn_a[0],&yn_a[0],
109  &jnp_a[0],&ynp_a[0]);
110 
111  // Shout if things went wrong
112 #ifdef PARANOID
113  if (n_actual!=int(N_fourier))
114  {
115  std::ostringstream error_stream;
116  error_stream << "CRBond_Bessel::bessjyna() only computed "
117  << n_actual << " rather than " << N_fourier
118  << " Bessel functions.\n";
119  throw OomphLibError(error_stream.str(),
120  OOMPH_CURRENT_FUNCTION,
121  OOMPH_EXCEPTION_LOCATION);
122  }
123 #endif
124 
125  // Evaluate Hankel at actual radius
126  Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier,rr,h,hp);
127 
128  // Evaluate Hankel at inner (unit) radius
129  Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier
130  ,Wavenumber,
131  h_a,hp_a);
132 
133  // Compute the sum: Separate the computation of the negative
134  // and positive terms
135  // ALH: The construction with the static cast to a double is to avoid
136  // a floating point exception when running unoptimised on my machine
137  for (unsigned i=0;i<N_fourier;i++)
138  {
139  u_ex-=pow(I,i)*h[i]*jnp_a[i]*pow(exp(I*theta),static_cast<double>(i))/hp_a[i];
140  }
141  for (unsigned i=1;i<N_fourier;i++)
142  {
143  u_ex-=pow(I,i)*h[i]*jnp_a[i]*pow(exp(-I*theta),static_cast<double>(i))/hp_a[i];
144  }
145 
146  // Get the real & imaginary part of the result
147  u[0]=real(u_ex);
148  u[1]=imag(u_ex);
149 
150  }// end of get_exact_u
151 
152 
153 
154  /// \short Flux (normal derivative) on the unit disk
155  /// for a planar incoming wave
156  void prescribed_incoming_flux(const Vector<double>& x,
157  complex<double>& flux)
158  {
159  // Switch to polar coordinates
160  double r;
161  r=sqrt(x[0]*x[0]+x[1]*x[1]);
162  double theta;
163  theta=atan2(x[1],x[0]);
164 
165  // Argument of the Bessel/Hankel fcts
166  double rr=Wavenumber*r;
167 
168  // Compute Bessel/Hankel functions
169  Vector<double> jn(N_fourier+1), yn(N_fourier+1),
170  jnp(N_fourier+1), ynp(N_fourier+1);
171 
172  // We want to compute N_fourier terms but the function
173  // may return fewer than that.
174  int n_actual=0;
175  CRBond_Bessel::bessjyna(N_fourier,rr,n_actual,&jn[0],&yn[0],
176  &jnp[0],&ynp[0]);
177 
178  // Shout if things went wrong...
179 #ifdef PARANOID
180  if (n_actual!=int(N_fourier))
181  {
182  std::ostringstream error_stream;
183  error_stream << "CRBond_Bessel::bessjyna() only computed "
184  << n_actual << " rather than " << N_fourier
185  << " Bessel functions.\n";
186  throw OomphLibError(error_stream.str(),
187  OOMPH_CURRENT_FUNCTION,
188  OOMPH_EXCEPTION_LOCATION);
189  }
190 #endif
191 
192  // Compute the sum: Separate the computation of the negative and
193  // positive terms
194  flux=std::complex<double>(0.0,0.0);
195  for (unsigned i=0;i<N_fourier;i++)
196  {
197  flux+=pow(I,i)*(Wavenumber)*pow(exp(I*theta),i)*jnp[i];
198  }
199  for (unsigned i=1;i<N_fourier;i++)
200  {
201  flux+=pow(I,i)*(Wavenumber)*pow(exp(-I*theta),i)*jnp[i];
202  }
203 
204  }// end of prescribed_incoming_flux
205 
206  class TestPMLMapping : public PMLMapping
207  {
208 
209  public:
210 
211  /// Default constructor (empty)
213 
214  /// \short Overwrite the pure PML mapping coefficient function to return the
215  /// coeffcients proposed by Bermudez et al
216  std::complex<double> gamma(const double& nu_i,
217  const double& pml_width_i,
218  const double& k_squared_local,
219  const double& alpha_shift=0.0)
220  {
221  // (return) gamma = 1 + (1/k) * (i/|outer_boundary - x|)
222  /*return 1.0 + (1.0 / std::complex<double> (sqrt(k_squared_local), 0))
223  * std::complex<double>
224  (0.0, 1.0/(std::fabs(pml_width_i - nu_i)));*/
225  return 1.0 + (1.0 / std::complex<double> (sqrt(k_squared_local), 0))
226  * std::complex<double>
227  (0.0, 2.0/(std::fabs(pml_width_i - nu_i)));
228  }
229 
230  };
231 
233 
234 } // end of namespace
235 
236 
237 /////////////////////////////////////////////////////////////////////
238 /////////////////////////////////////////////////////////////////////
239 /////////////////////////////////////////////////////////////////////
240 
241 //========= start_of_problem_class=====================================
242 /// Problem class to demonstrate use of perfectly matched layers
243 /// for Helmholtz problems.
244 //=====================================================================
245 template<class ELEMENT>
246 class PMLProblem : public Problem
247 {
248 
249 public:
250 
251  /// Constructor
252  PMLProblem();
253 
254  /// Destructor (empty)
256 
257  /// Update the problem specs before solve (empty)
259 
260  /// Update the problem specs after solve (empty)
262 
263  /// \short Doc the solution. DocInfo object stores flags/labels for where the
264  /// output gets written to
265  void doc_solution(DocInfo& doc_info);
266 
267  /// Create PML meshes
268  void create_pml_meshes();
269 
270  /// \short Create Helmholtz flux elements on boundary b of the Mesh pointed
271  /// to by bulk_mesh_pt and add them to the specified surface Mesh
272  void create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
273  Mesh* const & helmholtz_inner_boundary_mesh_pt);
274 
275  /// \short Create Helmholtz power elements on boundary b of the Mesh pointed
276  /// to by bulk_mesh_pt and add them to the specified surface Mesh
277  void create_power_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
278  Mesh* const & helmholtz_power_boundary_mesh_pt);
279 
280 #ifdef ADAPTIVE
281 
282  /// Actions before adapt: Wipe the PML meshes
283  void actions_before_adapt();
284 
285  /// Actions after adapt: Rebuild the PML meshes
286  void actions_after_adapt();
287 
288 
289 
290 #endif
291 
292 
293 #ifdef ADAPTIVE
294 
295 private:
296 
297  /// Pointer to the refineable "bulk" mesh
298  RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
299 
300 #else
301 
302 private:
303 
304  /// Pointer to the "bulk" mesh
305  TriangleMesh<ELEMENT>* Bulk_mesh_pt;
306 
307 #endif
308 
309 
310  /// Pointer to the right PML mesh
311  Mesh* PML_right_mesh_pt;
312 
313  /// Pointer to the top PML mesh
314  Mesh* PML_top_mesh_pt;
315 
316  /// Pointer to the left PML mesh
317  Mesh* PML_left_mesh_pt;
318 
319  /// Pointer to the bottom PML mesh
320  Mesh* PML_bottom_mesh_pt;
321 
322  /// Pointer to the top right corner PML mesh
323  Mesh* PML_top_right_mesh_pt;
324 
325  /// Pointer to the top left corner PML mesh
326  Mesh* PML_top_left_mesh_pt;
327 
328  /// Pointer to the bottom right corner PML mesh
329  Mesh* PML_bottom_right_mesh_pt;
330 
331  /// Pointer to the bottom left corner PML mesh
332  Mesh* PML_bottom_left_mesh_pt;
333 
334  /// \short Pointer to the mesh containing
335  /// the Helmholtz inner boundary condition elements
337 
338  /// Pointer to mesh of elements that compute the radiated power
340 
341  /// Trace file
342  ofstream Trace_file;
343 
344 }; // end of problem class
345 
346 
347 
348 //=======start_of_constructor=============================================
349 /// Constructor for Helmholtz problem
350 //========================================================================
351 template<class ELEMENT>
353 {
354 
355  // Open trace file
356  Trace_file.open("RESLT/trace.dat");
357 
358  // Create circle representing inner boundary
359  double a=1.0;
360  double x_c=0.0;
361  double y_c=0.0;
362  Circle* inner_circle_pt=new Circle(x_c,y_c,a);
363 
364  // Outer boundary
365  //---------------
366  TriangleMeshClosedCurve* outer_boundary_pt=0;
367 
368  unsigned n_segments = 20;
369  Vector<TriangleMeshCurveSection*> outer_boundary_line_pt(4);
370 
371  // Each polyline only has three vertices, provide storage for their
372  // coordinates
373  Vector<Vector<double> > vertex_coord(2);
374  for(unsigned i=0;i<2;i++)
375  {
376  vertex_coord[i].resize(2);
377  }
378 
379  // First polyline
380  vertex_coord[0][0]=-2.0;
381  vertex_coord[0][1]=-2.0;
382  vertex_coord[1][0]=-2.0;
383  vertex_coord[1][1]=2.0;
384 
385  // Build the 1st boundary polyline
386  unsigned boundary_id=2;
387  outer_boundary_line_pt[0] = new TriangleMeshPolyLine(vertex_coord,
388  boundary_id);
389 
390  // Second boundary polyline
391  vertex_coord[0][0]=-2.0;
392  vertex_coord[0][1]=2.0;
393  vertex_coord[1][0]=2.0;
394  vertex_coord[1][1]=2.0;
395 
396  // Build the 2nd boundary polyline
397  boundary_id=3;
398  outer_boundary_line_pt[1] = new TriangleMeshPolyLine(vertex_coord,
399  boundary_id);
400 
401  // Third boundary polyline
402  vertex_coord[0][0]=2.0;
403  vertex_coord[0][1]=2.0;
404  vertex_coord[1][0]=2.0;
405  vertex_coord[1][1]=-2.0;
406 
407  // Build the 3rd boundary polyline
408  boundary_id=4;
409  outer_boundary_line_pt[2] = new TriangleMeshPolyLine(vertex_coord,
410  boundary_id);
411 
412  // Fourth boundary polyline
413  vertex_coord[0][0]=2.0;
414  vertex_coord[0][1]=-2.0;
415  vertex_coord[1][0]=-2.0;
416  vertex_coord[1][1]=-2.0;
417 
418  // Build the 4th boundary polyline
419  boundary_id=5;
420  outer_boundary_line_pt[3] = new TriangleMeshPolyLine(vertex_coord,
421  boundary_id);
422 
423  // Create the triangle mesh polygon for outer boundary
424  outer_boundary_pt = new TriangleMeshPolygon(outer_boundary_line_pt);
425 
426  // Inner circular boundary
427  //------------------------
428  Vector<TriangleMeshCurveSection*> inner_boundary_line_pt(2);
429 
430  // The intrinsic coordinates for the beginning and end of the curve
431  double s_start = 0.0;
432  double s_end = MathematicalConstants::Pi;
433  boundary_id = 0;
434  inner_boundary_line_pt[0]=
435  new TriangleMeshCurviLine(inner_circle_pt,
436  s_start,
437  s_end,
438  n_segments,
439  boundary_id);
440 
441  // The intrinsic coordinates for the beginning and end of the curve
442  s_start = MathematicalConstants::Pi;
443  s_end = 2.0*MathematicalConstants::Pi;
444  boundary_id = 1;
445  inner_boundary_line_pt[1]=
446  new TriangleMeshCurviLine(inner_circle_pt,
447  s_start,
448  s_end,
449  n_segments,
450  boundary_id);
451 
452 
453  // Combine to hole
454  //----------------
455  Vector<TriangleMeshClosedCurve*> hole_pt(1);
456  Vector<double> hole_coords(2);
457  hole_coords[0]=0.0;
458  hole_coords[1]=0.0;
459  hole_pt[0]=new TriangleMeshClosedCurve(inner_boundary_line_pt,
460  hole_coords);
461 
462 
463  // Use the TriangleMeshParameters object for helping on the manage
464  // of the TriangleMesh parameters. The only parameter that needs to take
465  // is the outer boundary.
466  TriangleMeshParameters triangle_mesh_parameters(outer_boundary_pt);
467 
468  // Specify the closed curve using the TriangleMeshParameters object
469  triangle_mesh_parameters.internal_closed_curve_pt() = hole_pt;
470 
471  // Target element size in bulk mesh
472  triangle_mesh_parameters.element_area() = 0.05;
473 
474 #ifdef ADAPTIVE
475 
476  // Build adaptive "bulk" mesh
477  Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
478 
479  // Create/set error estimator
480  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
481 
482  // Choose error tolerances to force some uniform refinement
483  Bulk_mesh_pt->min_permitted_error()=0.00004;
484  Bulk_mesh_pt->max_permitted_error()=0.0001;
485 
486 #else
487 
488  // Build "bulk" mesh
489  Bulk_mesh_pt=new TriangleMesh<ELEMENT>(triangle_mesh_parameters);
490 
491 #endif
492 
493  // Pointer to mesh containing the Helmholtz inner boundary condition
494  // elements. Specify outer radius
495  Helmholtz_inner_boundary_mesh_pt = new Mesh;
496 
497  // Create prescribed-flux elements from all elements that are
498  // adjacent to the inner boundary , but add them to a separate mesh.
499  create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
500  create_flux_elements(1,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
501 
502  // Pointer to mesh containing the Helmholtz power condition
503  // elements.
504  Helmholtz_power_boundary_mesh_pt = new Mesh;
505 
506  // Create power elements from all elements that are
507  // adjacent to the inner boundary , but add them to a separate mesh.
508  create_power_elements(0,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
509  create_power_elements(1,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
510 
511  // Create the main triangular mesh
512  add_sub_mesh(Bulk_mesh_pt);
513 
514  // Create PML meshes and add them to the global mesh
515  create_pml_meshes();
516 
517  // Add the flux mesh
518  add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
519 
520  // Build the entire mesh from its submeshes
521  build_global_mesh();
522 
523  // Complete the build of all elements so they are fully functional
524  unsigned n_element = this->mesh_pt()->nelement();
525  for(unsigned e=0;e<n_element;e++)
526  {
527  // Upcast from GeneralisedElement to Pml Helmholtz bulk element
528  PMLHelmholtzEquations<2> *el_pt =
529  dynamic_cast<PMLHelmholtzEquations<2>*>(mesh_pt()->element_pt(e));
530 
531  if (el_pt!=0)
532  {
533  //Set the k_squared function pointer
534  el_pt->k_squared_pt() = &GlobalParameters::K_squared;
535  // Pass in a pointer to the class containing the PML mapping function
536  //el_pt->pml_mapping_pt() = GlobalParameters::Test_pml_mapping_pt;
537  }
538 
539  }
540 
541  // Setup equation numbering scheme
542  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
543 
544 } // end of constructor
545 
546 
547 #ifdef ADAPTIVE
548 
549 //=====================start_of_actions_before_adapt======================
550 /// Actions before adapt: Wipe the mesh of face elements
551 //========================================================================
552 template<class ELEMENT>
554 {
555  // Before adapting the added PML meshes must be removed
556  // as they are not refineable and are to be rebuilt from the
557  // newly refined triangular mesh
558  delete PML_right_mesh_pt;
559  PML_right_mesh_pt=0;
560  delete PML_top_mesh_pt;
561  PML_top_mesh_pt=0;
562  delete PML_left_mesh_pt;
563  PML_left_mesh_pt=0;
564  delete PML_bottom_mesh_pt;
565  PML_bottom_mesh_pt=0;
566  delete PML_top_right_mesh_pt;
567  PML_top_right_mesh_pt=0;
568  delete PML_top_left_mesh_pt;
569  PML_top_left_mesh_pt=0;
570  delete PML_bottom_right_mesh_pt;
571  PML_bottom_right_mesh_pt=0;
572  delete PML_bottom_left_mesh_pt;
573  PML_bottom_left_mesh_pt=0;
574 
575  // Rebuild the Problem's global mesh from its various sub-meshes
576  // but first flush all its submeshes
577  flush_sub_meshes();
578 
579  // Then add the triangular mesh back
580  add_sub_mesh(Bulk_mesh_pt);
581 
582  // Rebuild the global mesh such that it now stores
583  // the triangular mesh only
584  rebuild_global_mesh();
585 
586 }// end of actions_before_adapt
587 
588 
589 //=====================start_of_actions_after_adapt=======================
590 /// Actions after adapt: Rebuild the face element meshes
591 //========================================================================
592 template<class ELEMENT>
594 {
595  create_inner_bc_elements(0,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
596  create_inner_bc_elements(1,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
597  // Build PML meshes and add them to the global mesh
598  create_pml_meshes();
599 
600  // Complete the build of all elements so they are fully functional
601 
602  // Loop over the entire mesh elements to set up element-specific
603  // things that cannot be handled by constructor
604  unsigned n_element = this->mesh_pt()->nelement();
605 
606  for(unsigned e=0;e<n_element;e++)
607  {
608  // Upcast from GeneralisedElement to PMLHelmholtz bulk element
609  PMLHelmholtzEquations<2> *el_pt =
610  dynamic_cast<PMLHelmholtzEquations<2>*>(mesh_pt()->element_pt(e));
611 
612  if (el_pt!=0)
613  {
614  //Set the k_squared double pointer
615  el_pt->k_squared_pt() = &GlobalParameters::K_squared;
616  }
617  }
618 
619  // Create prescribed-flux elements and BC elements
620  // from all elements that are adjacent to the boundaries and add them to
621  // Helmholtz_boundary_meshes
622  create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
623  create_flux_elements(1,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
624 er_eleme
625  // Create power elements from all elements that are
626  // adjacent to the inner boundary , but add them to a separate mesh.
627  create_power_elements(0,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
628  create_power_elements(1,Bulk_mesh_pt,Helmholtz_power_boundary_mesh_pt);
629 
630  // Add the flux mesh
631  add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
632 
633  // Rebuild the Problem's global mesh from its various sub-meshes
634  rebuild_global_mesh();
635 
636 }// end of actions_after_adapt
637 
638 #endif
639 
640 //=====================start_of_doc=======================================
641 /// Doc the solution: doc_info contains labels/output directory etc.
642 //========================================================================
643 template<class ELEMENT>
644 void PMLProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
645 {
646 
647  ofstream some_file,some_file2;
648  char filename[100];
649 
650  // Number of plot points
651  unsigned npts;
652  npts=5;
653 
654  // Compute/output the radiated power
655  //----------------------------------
656  sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
657  doc_info.number());
658  some_file.open(filename);
659 
660  // Accumulate contribution from elements
661  double power=0.0;
662  unsigned nn_element=Helmholtz_power_boundary_mesh_pt->nelement();
663  for(unsigned e=0;e<nn_element;e++)
664  {
665  PMLHelmholtzPowerElement<ELEMENT> *el_pt =
666  dynamic_cast<PMLHelmholtzPowerElement<ELEMENT>*>(
667  Helmholtz_power_boundary_mesh_pt->element_pt(e));
668  power += el_pt->global_power_contribution(some_file);
669  }
670  some_file.close();
671  oomph_info << "Total radiated power: " << power << std::endl;
672 
673  // Output solution
674  //-----------------
675  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
676  doc_info.number());
677  some_file.open(filename);
678  Bulk_mesh_pt->output(some_file,npts);
679  some_file.close();
680 
681 // Output exact solution
682 //----------------------
683  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
684  doc_info.number());
685  some_file.open(filename);
686  Bulk_mesh_pt->output_fct(some_file,npts,GlobalParameters::get_exact_u);
687  some_file.close();
688 
689  double error,norm;
690  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
691  doc_info.number());
692  some_file.open(filename);
693  Bulk_mesh_pt->compute_error(some_file,GlobalParameters::get_exact_u,
694  error,norm);
695  some_file.close();
696 
697  // Doc L2 error and norm of solution
698  oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
699  oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
700 
701  // Output PML layers
702  //-----------------
703  sprintf(filename,"%s/pml_soln%i.dat",doc_info.directory().c_str(),
704  doc_info.number());
705  some_file.open(filename);
706  PML_right_mesh_pt->output(some_file,npts);
707  PML_left_mesh_pt->output(some_file,npts);
708  PML_top_mesh_pt->output(some_file,npts);
709  PML_bottom_mesh_pt->output(some_file,npts);
710  PML_bottom_right_mesh_pt->output(some_file,npts);
711  PML_top_right_mesh_pt->output(some_file,npts);
712  PML_bottom_left_mesh_pt->output(some_file,npts);
713  PML_top_left_mesh_pt->output(some_file,npts);
714  some_file.close();
715 
716 
717 
718  // Write norm of solution to trace file
719  double norm2=0.0;
720  Bulk_mesh_pt->compute_norm(norm2);
721  Trace_file << norm2 << std::endl;
722 
723  // //Do animation of Helmholtz solution
724  // //-----------------------------------
725  // unsigned nstep=40;
726  // for (unsigned i=0;i<nstep;i++)
727  // {
728  // sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",
729  // doc_info.directory().c_str(),
730  // doc_info.number(),i);
731  // some_file.open(filename);
732  // double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);
733  // unsigned nelem=Bulk_mesh_pt->nelement();
734  // for (unsigned e=0;e<nelem;e++)
735  // {
736  // ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
737  // Bulk_mesh_pt->element_pt(e));
738  // el_pt->output_real(some_file,phi,npts);
739  // }
740  // some_file.close();
741  // }
742 
743 } // end of doc
744 
745 //============start_of_create_flux_elements==================================
746 /// Create Helmholtz inner Flux Elements on the b-th boundary of
747 /// the Mesh object pointed to by bulk_mesh_pt and add the elements
748 /// to the Mesh object pointed to by helmholtz_inner_boundary_mesh_pt
749 //============================================================================
750 template<class ELEMENT>
752 create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
753  Mesh* const & helmholtz_inner_boundary_mesh_pt)
754 {
755  // Loop over the bulk elements adjacent to boundary b
756  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
757  for(unsigned e=0;e<n_element;e++)
758  {
759  // Get pointer to the bulk element that is adjacent to boundary b
760  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
761  bulk_mesh_pt->boundary_element_pt(b,e));
762 
763  //Find the index of the face of element e along boundary b
764  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
765 
766  // Build the corresponding prescribed incoming-flux element
767  PMLHelmholtzFluxElement<ELEMENT>* flux_element_pt = new
768  PMLHelmholtzFluxElement<ELEMENT>(bulk_elem_pt,face_index);
769 
770  //Add the prescribed incoming-flux element to the surface mesh
771  helmholtz_inner_boundary_mesh_pt->add_element_pt(flux_element_pt);
772 
773  // Set the pointer to the prescribed flux function
774  flux_element_pt->flux_fct_pt() =
776 
777  } //end of loop over bulk elements adjacent to boundary b
778 
779 } // end of create_flux_elements
780 
781 //============start_of_create_power_elements==================================
782 /// Create Helmholtz inner Flux Elements on the b-th boundary of
783 /// the Mesh object pointed to by bulk_mesh_pt and add the elements
784 /// to the Mesh object pointed to by helmholtz_power_boundary_mesh_pt
785 //============================================================================
786 template<class ELEMENT>
788 create_power_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
789  Mesh* const & helmholtz_power_boundary_mesh_pt)
790 {
791  // Loop over the bulk elements adjacent to boundary b
792  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
793  for(unsigned e=0;e<n_element;e++)
794  {
795  // Get pointer to the bulk element that is adjacent to boundary b
796  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
797  bulk_mesh_pt->boundary_element_pt(b,e));
798 
799  //Find the index of the face of element e along boundary b
800  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
801 
802  // Build the corresponding prescribed power element
803  PMLHelmholtzPowerElement<ELEMENT>* power_element_pt = new
804  PMLHelmholtzPowerElement<ELEMENT>(bulk_elem_pt,face_index);
805 
806  //Add the prescribed power element to the surface mesh
807  helmholtz_power_boundary_mesh_pt->add_element_pt(power_element_pt);
808 
809  } //end of loop over bulk elements adjacent to boundary b
810 
811 } // end of create_power_elements
812 
813 
814 //============start_of_create_pml_meshes======================================
815 /// Create PML meshes and add them to the problem's sub-meshes
816 //============================================================================
817 template<class ELEMENT>
819 {
820 
821  // Bulk mesh left boundary id
822  unsigned int left_boundary_id = 2;
823 
824  // Bulk mesh top boundary id
825  unsigned int top_boundary_id = 3;
826 
827  // Bulk mesh right boundary id
828  unsigned int right_boundary_id = 4;
829 
830  // Bulk mesh bottom boundary id
831  unsigned int bottom_boundary_id = 5;
832 
833  // PML width in elements for the right layer
834  unsigned n_x_right_pml = 3;
835 
836  // PML width in elements for the top layer
837  unsigned n_y_top_pml = 3;
838 
839  // PML width in elements for the left layer
840  unsigned n_x_left_pml = 3;
841 
842  // PML width in elements for the bottom layer
843  unsigned n_y_bottom_pml = 3;
844 
845  // Outer physical length of the PML layers
846  // defaults to 0.2, so 10% of the size of the
847  // physical domain
848  double width_x_right_pml = 0.2;
849  double width_y_top_pml = 0.2;
850  double width_x_left_pml = 0.2;
851  double width_y_bottom_pml = 0.2;
852 
853  // Build the PML meshes based on the new adapted triangular mesh
854  PML_right_mesh_pt =
855  TwoDimensionalPMLHelper::create_right_pml_mesh
856  <PMLLayerElement<ELEMENT> >
857  (Bulk_mesh_pt,right_boundary_id,
858  n_x_right_pml, width_x_right_pml);
859  PML_top_mesh_pt =
860  TwoDimensionalPMLHelper::create_top_pml_mesh
861  <PMLLayerElement<ELEMENT> >
862  (Bulk_mesh_pt, top_boundary_id,
863  n_y_top_pml, width_y_top_pml);
864  PML_left_mesh_pt =
865  TwoDimensionalPMLHelper::create_left_pml_mesh
866  <PMLLayerElement<ELEMENT> >
867  (Bulk_mesh_pt, left_boundary_id,
868  n_x_left_pml, width_x_left_pml);
869  PML_bottom_mesh_pt=
870  TwoDimensionalPMLHelper::create_bottom_pml_mesh
871  <PMLLayerElement<ELEMENT> >
872  (Bulk_mesh_pt, bottom_boundary_id,
873  n_y_bottom_pml, width_y_bottom_pml);
874 
875  // Add submeshes to the global mesh
876  add_sub_mesh(PML_right_mesh_pt);
877  add_sub_mesh(PML_top_mesh_pt);
878  add_sub_mesh(PML_left_mesh_pt);
879  add_sub_mesh(PML_bottom_mesh_pt);
880 
881  // Rebuild corner PML meshes
882  PML_top_right_mesh_pt =
883  TwoDimensionalPMLHelper::create_top_right_pml_mesh
884  <PMLLayerElement<ELEMENT> >
885  (PML_right_mesh_pt, PML_top_mesh_pt,
886  Bulk_mesh_pt, right_boundary_id);
887 
888  PML_bottom_right_mesh_pt =
889  TwoDimensionalPMLHelper::create_bottom_right_pml_mesh
890  <PMLLayerElement<ELEMENT> >
891  (PML_right_mesh_pt, PML_bottom_mesh_pt,
892  Bulk_mesh_pt, right_boundary_id);
893 
894  PML_top_left_mesh_pt =
895  TwoDimensionalPMLHelper::create_top_left_pml_mesh
896  <PMLLayerElement<ELEMENT> >
897  (PML_left_mesh_pt, PML_top_mesh_pt,
898  Bulk_mesh_pt, left_boundary_id);
899 
900  PML_bottom_left_mesh_pt =
901  TwoDimensionalPMLHelper::create_bottom_left_pml_mesh
902  <PMLLayerElement<ELEMENT> >
903  (PML_left_mesh_pt, PML_bottom_mesh_pt,
904  Bulk_mesh_pt, left_boundary_id);
905 
906  // Add submeshes to the global mesh
907  add_sub_mesh(PML_top_right_mesh_pt);
908  add_sub_mesh(PML_bottom_right_mesh_pt);
909  add_sub_mesh(PML_top_left_mesh_pt);
910  add_sub_mesh(PML_bottom_left_mesh_pt);
911 
912 } // end of create_pml_meshes
913 
914 
915 
916 //==========start_of_main=================================================
917 /// Solve 2D Helmholtz problem
918 //========================================================================
919 int main(int argc, char **argv)
920 {
921  //Set up the problem
922  //------------------
923 
924 #ifdef ADAPTIVE
925 
926  // Set up the problem with projectable 2D six-node elements from the
927  // TPMLHelmholtzElement family.
928  PMLProblem<ProjectablePMLHelmholtzElement
929  <TPMLHelmholtzElement<2,3> > > problem;
930 
931 #else
932 
933  // Set up the problem with 2D six-node elements from the
934  // TPMLHelmholtzElement family.
936 #endif
937 
938  // Create label for output
939  //------------------------
940  DocInfo doc_info;
941 
942  // Set output directory
943  doc_info.set_directory("RESLT");
944 
945 
946 #ifdef ADAPTIVE
947 
948  // Max. number of adaptations
949  unsigned max_adapt=1;
950 
951  // Solve the problem with the adaptive Newton's method, allowing
952  // up to max_adapt mesh adaptations after every solve.
953  problem.newton_solve(max_adapt);
954 
955 #else
956 
957  // Solve the problem with Newton's method
958  problem.newton_solve();
959 
960 #endif
961 
962  //Output solution
963  problem.doc_solution(doc_info);
964 
965 } //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_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_inner_boundary_mesh_pt)
Create Helmholtz flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to t...
void create_pml_meshes()
Create PML meshes.
std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &k_squared_local, const double &alpha_shift=0.0)
Overwrite the pure PML mapping coefficient function to return the coeffcients proposed by Bermudez et...
double Wavenumber
Wavenumber (also known as k), k=omega/c.
void create_power_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_power_boundary_mesh_pt)
Create Helmholtz power elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to ...
std::complex< double > I(0.0, 1.0)
Imaginary unit.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to...
void actions_after_adapt()
Actions after adapt: Rebuild the PML meshes.
unsigned N_fourier
Number of terms used in the computation of the exact solution.
int main(int argc, char **argv)
Solve 2D Helmholtz problem.
Mesh * Helmholtz_inner_boundary_mesh_pt
Pointer to the mesh containing the Helmholtz inner boundary condition elements.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_adapt()
Actions before adapt: Wipe the PML meshes.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution for scattered field (vector returns real and impaginary parts).
Mesh * Helmholtz_power_boundary_mesh_pt
Pointer to mesh of elements that compute the radiated power.
void prescribed_incoming_flux(const Vector< double > &x, complex< double > &flux)
Flux (normal derivative) on the unit disk for a planar incoming wave.