oscillating_sphere.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: 1259 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-11-08 23:52:03 +0000 (Tue, 08 Nov 2016) $
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 pml Fourier-decomposed Helmholtz problem
31 
32 #include <complex>
33 #include <cmath>
34 
35 //Generic routines
36 #include "generic.h"
37 
38 // The Helmholtz equations
39 #include "pml_fourier_decomposed_helmholtz.h"
40 
41 // The mesh
42 #include "meshes/triangle_mesh.h"
43 
44 // Get the Bessel functions
45 #include "oomph_crbond_bessel.h"
46 
47 using namespace oomph;
48 using namespace std;
49 
50 
51 /////////////////////////////////////////////////////////////////////
52 /////////////////////////////////////////////////////////////////////
53 /////////////////////////////////////////////////////////////////////
54 
55 
56 
57 //===== start_of_namespace=============================================
58 /// Namespace for the Fourier decomposed Helmholtz problem parameters
59 //=====================================================================
60 namespace ProblemParameters
61 {
62  /// Output directory
63  string Directory="RESLT";
64 
65  /// Frequency
66  double K_squared = 10.0;
67 
68  /// Default physical PML thickness
69  double PML_thickness=4.0;
70 
71  /// Default number of elements within PMLs
72  unsigned Nel_pml=15;
73 
74  /// Target area for initial mesh
75  double Element_area = 0.1;
76 
77  /// The default Fourier wave number
78  int N_fourier=0;
79 
80  /// Number of terms in the exact solution
81  unsigned N_terms=6;
82 
83  /// Coefficients in the exact solution
84  Vector<double> Coeff(N_terms,1.0);
85 
86  /// Imaginary unit
87  std::complex<double> I(0.0,1.0);
88 
89  /// Exact solution as a Vector of size 2, containing real and imag parts
90  void get_exact_u(const Vector<double>& x, Vector<double>& u)
91  {
92  double K = sqrt(K_squared);
93 
94  // Switch to spherical coordinates
95  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
96 
97  double theta;
98  theta=atan2(x[0],x[1]);
99 
100  // Argument for Bessel/Hankel functions
101  double kr= K*R;
102 
103  // Need half-order Bessel functions
104  double bessel_offset=0.5;
105 
106  // Evaluate Bessel/Hankel functions
107  Vector<double> jv(N_terms);
108  Vector<double> yv(N_terms);
109  Vector<double> djv(N_terms);
110  Vector<double> dyv(N_terms);
111  double order_max_in=double(N_terms-1)+bessel_offset;
112  double order_max_out=0;
113 
114  // This function returns vectors containing
115  // J_k(x), Y_k(x) and their derivatives
116  // up to k=order_max, with k increasing in
117  // integer increments starting with smallest
118  // positive value. So, e.g. for order_max=3.5
119  // jv[0] contains J_{1/2}(x),
120  // jv[1] contains J_{3/2}(x),
121  // jv[2] contains J_{5/2}(x),
122  // jv[3] contains J_{7/2}(x).
123  CRBond_Bessel::bessjyv(order_max_in,
124  kr,
125  order_max_out,
126  &jv[0],&yv[0],
127  &djv[0],&dyv[0]);
128 
129  // Assemble exact solution (actually no need to add terms
130  // below i=N_fourier as Legendre polynomial would be zero anyway)
131  complex<double> u_ex(0.0,0.0);
132  for(unsigned i=N_fourier;i<N_terms;i++)
133  {
134  //Associated_legendre_functions
135  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
136  cos(theta));
137  // Set exact solution
138  u_ex+=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*(jv[i]+I*yv[i])*p;
139  }
140 
141  // Get the real & imaginary part of the result
142  u[0]=u_ex.real();
143  u[1]=u_ex.imag();
144 
145  }//end of get_exact_u
146 
147 
148 
149  /// \short Get -du/dr (spherical r) for exact solution. Equal to prescribed
150  /// flux on inner boundary.
151  void exact_minus_dudr(const Vector<double>& x, std::complex<double>& flux)
152  {
153  double K = sqrt(K_squared);
154 
155  // Initialise flux
156  flux=std::complex<double>(0.0,0.0);
157 
158  // Switch to spherical coordinates
159  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
160 
161  double theta;
162  theta=atan2(x[0],x[1]);
163 
164  // Argument for Bessel/Hankel functions
165  double kr=K*R;
166 
167  // Need half-order Bessel functions
168  double bessel_offset=0.5;
169 
170  // Evaluate Bessel/Hankel functions
171  Vector<double> jv(N_terms);
172  Vector<double> yv(N_terms);
173  Vector<double> djv(N_terms);
174  Vector<double> dyv(N_terms);
175  double order_max_in=double(N_terms-1)+bessel_offset;
176  double order_max_out=0;
177 
178  // This function returns vectors containing
179  // J_k(x), Y_k(x) and their derivatives
180  // up to k=order_max, with k increasing in
181  // integer increments starting with smallest
182  // positive value. So, e.g. for order_max=3.5
183  // jv[0] contains J_{1/2}(x),
184  // jv[1] contains J_{3/2}(x),
185  // jv[2] contains J_{5/2}(x),
186  // jv[3] contains J_{7/2}(x).
187  CRBond_Bessel::bessjyv(order_max_in,
188  kr,
189  order_max_out,
190  &jv[0],&yv[0],
191  &djv[0],&dyv[0]);
192 
193 
194  // Assemble exact solution (actually no need to add terms
195  // below i=N_fourier as Legendre polynomial would be zero anyway)
196  complex<double> u_ex(0.0,0.0);
197  for(unsigned i=N_fourier;i<N_terms;i++)
198  {
199  //Associated_legendre_functions
200  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
201  cos(theta));
202  // Set flux of exact solution
203  flux-=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*p*
204  ( K*(djv[i]+I*dyv[i]) - (0.5*(jv[i]+I*yv[i])/R) );
205  }
206  } // end of exact_normal_derivative
207 
208 
209  /// Radial position of point source
210  double R_source = 2.0;
211 
212  /// Axial position of point source
213  double Z_source = 2.0;
214 
215  /// Point source magnitude (Complex)
216  std::complex<double> Magnitude(100.0,100.0);
217 
218 } // end of namespace
219 
220 
221 /////////////////////////////////////////////////////////////////////
222 /////////////////////////////////////////////////////////////////////
223 /////////////////////////////////////////////////////////////////////
224 
225 
226 namespace oomph
227 {
228 
229 //========= start_of_point_source_wrapper==============================
230 /// Class to impose point source to (wrapped) Helmholtz element
231 //=====================================================================
232 template<class ELEMENT>
233 class PMLHelmholtzPointSourceElement : public virtual ELEMENT
234 {
235 
236 public:
237 
238  /// Constructor
240  {
241  // Initialise
242  Point_source_magnitude=std::complex<double>(0.0,0.0);
243  }
244 
245  /// Destructor (empty)
247 
248  /// Set local coordinate and magnitude of point source
249  void setup(const Vector<double>& s_point_source,
250  const std::complex<double>& magnitude)
251  {
252  S_point_source=s_point_source;
253  Point_source_magnitude=magnitude;
254  }
255 
256 
257  /// Add the element's contribution to its residual vector (wrapper)
258  void fill_in_contribution_to_residuals(Vector<double> &residuals)
259  {
260  //Call the generic residuals function with flag set to 0
261  //using a dummy matrix argument
262  ELEMENT::fill_in_generic_residual_contribution_pml_fourier_decomposed_helmholtz(
263  residuals,GeneralisedElement::Dummy_matrix,0);
264 
265  // Add point source contribution
266  fill_in_point_source_contribution_to_residuals(residuals);
267  }
268 
269 
270  /// \short Add the element's contribution to its residual vector and
271  /// element Jacobian matrix (wrapper)
272  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
273  DenseMatrix<double> &jacobian)
274  {
275  //Call the generic routine with the flag set to 1
276  ELEMENT::fill_in_generic_residual_contribution_pml_fourier_decomposed_helmholtz(residuals,
277  jacobian,1);
278 
279  // Add point source contribution
280  fill_in_point_source_contribution_to_residuals(residuals);
281  }
282 
283 
284 private:
285 
286 
287  /// Add the point source contribution to the residual vector
288  void fill_in_point_source_contribution_to_residuals(Vector<double> &residuals)
289  {
290  // No further action
291  if (S_point_source.size()==0) return;
292 
293  //Find out how many nodes there are
294  const unsigned n_node = this->nnode();
295 
296  //Set up memory for the shape/test functions
297  Shape psi(n_node);
298 
299  //Integers to store the local equation and unknown numbers
300  int local_eqn_real=0;
301  int local_eqn_imag=0;
302 
303  // Get shape/test fcts
304  this->shape(S_point_source,psi);
305 
306  // Assemble residuals
307  //--------------------------------
308 
309  // Loop over the test functions
310  for(unsigned l=0;l<n_node;l++)
311  {
312  // first, compute the real part contribution
313  //-------------------------------------------
314 
315  //Get the local equation
316  local_eqn_real =
317  this->nodal_local_eqn
318  (l,this->u_index_pml_fourier_decomposed_helmholtz().real());
319 
320  /*IF it's not a boundary condition*/
321  if(local_eqn_real >= 0)
322  {
323  residuals[local_eqn_real] -= 2.0*Point_source_magnitude.real()*psi(l);
324  }
325 
326  // Second, compute the imaginary part contribution
327  //------------------------------------------------
328 
329  //Get the local equation
330  local_eqn_imag =
331  this->nodal_local_eqn
332  (l,this->u_index_pml_fourier_decomposed_helmholtz().imag());
333 
334  /*IF it's not a boundary condition*/
335  if(local_eqn_imag >= 0)
336  {
337  // Add body force/source term and Helmholtz bit
338  residuals[local_eqn_imag] -= 2.0*Point_source_magnitude.imag()*psi(l);
339  }
340  }
341  }
342 
343  /// Local coordinates of point at which point source is applied
344  Vector<double> S_point_source;
345 
346  /// Magnitude of point source (complex!)
347  std::complex<double> Point_source_magnitude;
348 
349  };
350 
351 
352 
353 
354 //=======================================================================
355 /// Face geometry for element is the same as that for the underlying
356 /// wrapped element
357 //=======================================================================
358  template<class ELEMENT>
359  class FaceGeometry<PMLHelmholtzPointSourceElement<ELEMENT> >
360  : public virtual FaceGeometry<ELEMENT>
361  {
362  public:
363  FaceGeometry() : FaceGeometry<ELEMENT>() {}
364  };
365 
366 
367 //=======================================================================
368 /// Face geometry of the Face Geometry for element is the same as
369 /// that for the underlying wrapped element
370 //=======================================================================
371  template<class ELEMENT>
372  class FaceGeometry<FaceGeometry<PMLHelmholtzPointSourceElement<ELEMENT> > >
373  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
374  {
375  public:
376  FaceGeometry() : FaceGeometry<FaceGeometry<ELEMENT> >() {}
377  };
378 
379 
380 //=======================================================================
381 /// Policy class defining the elements to be used in the actual
382 /// PML layers.
383 //=======================================================================
384  template<class ELEMENT>
385 class PMLLayerElement<PMLHelmholtzPointSourceElement<ELEMENT> > :
386  public virtual PMLLayerElement<ELEMENT>
387 {
388 
389  public:
390 
391  /// \short Constructor: Call the constructor for the
392  /// appropriate Element
393  PMLLayerElement() : PMLLayerElement<ELEMENT>()
394  {}
395 
396 };
397 
398 
399 
400 //=======================================================================
401 /// Policy class defining the elements to be used in the actual
402 /// PML layers.
403 //=======================================================================
404  template<class ELEMENT>
405 class PMLLayerElement<
406  ProjectablePMLFourierDecomposedHelmholtzElement<
407  PMLHelmholtzPointSourceElement<ELEMENT> > >:
408  public virtual PMLLayerElement<ELEMENT>
409 {
410 
411  public:
412 
413  /// \short Constructor: Call the constructor for the
414  /// appropriate Element
415  PMLLayerElement() : PMLLayerElement<ELEMENT>()
416  {}
417 
418 };
419 
420 }
421 
422 
423 
424 
425 /////////////////////////////////////////////////////////////////////
426 /////////////////////////////////////////////////////////////////////
427 /////////////////////////////////////////////////////////////////////
428 
429 //========= start_of_problem_class=====================================
430 /// Problem class
431 //=====================================================================
432 template<class ELEMENT>
434 {
435 
436 public:
437 
438  /// Constructor
440 
441  /// Destructor (empty)
443 
444  /// Update the problem specs before solve (empty)
446 
447  /// Update the problem after solve (empty)
449 
450  /// \short Doc the solution. DocInfo object stores flags/labels for where the
451  /// output gets written to
452  void doc_solution(DocInfo& doc_info);
453 
454  /// Create PML meshes
455  void create_pml_meshes();
456 
457  /// Create mesh of face elements that monitor the radiated power
458  void create_power_monitor_mesh();
459 
460  #ifdef ADAPTIVE
461 
462  /// Actions before adapt: Wipe the mesh of prescribed flux elements
463  void actions_before_adapt();
464 
465  /// Actions after adapt: Rebuild the mesh of prescribed flux elements
466  void actions_after_adapt();
467 
468  #endif
469 
470 
471  // Apply boundary condtions for odd Fourier wavenumber
472  void complete_problem_setup();
473 
474 private:
475 
476  /// Create flux elements on inner boundary
477  void create_flux_elements_on_inner_boundary();
478 
479  /// \short Delete boundary face elements and wipe the surface mesh
480  void delete_face_elements(Mesh* const & boundary_mesh_pt)
481  {
482  // Loop over the surface elements
483  unsigned n_element = boundary_mesh_pt->nelement();
484  for(unsigned e=0;e<n_element;e++)
485  {
486  // Kill surface element
487  delete boundary_mesh_pt->element_pt(e);
488  }
489 
490  // Wipe the mesh
491  boundary_mesh_pt->flush_element_and_node_storage();
492  }
493 
494  // Apply boundary condtions for odd Fourier wavenumber
495  void apply_zero_dirichlet_boundary_conditions();
496 
497  /// Pointer to mesh that stores the power monitor elements
499 
500 
501 #ifdef ADAPTIVE
502 
503  // Create point source element (only used in adaptive run)
504  void setup_point_source();
505 
506  /// Pointer to the refineable "bulk" mesh
507  RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
508 
509  /// Mesh as geometric object representation of bulk mesh
510  MeshAsGeomObject* Mesh_as_geom_obj_pt;
511 
512 #else
513 
514  /// Pointer to the "bulk" mesh
515  TriangleMesh<ELEMENT>* Bulk_mesh_pt;
516 
517 #endif
518 
519  /// Mesh of FaceElements that apply the flux bc on the inner boundary
521 
522  /// Pointer to the right PML mesh
524 
525  /// Pointer to the top PML mesh
527 
528  /// Pointer to the bottom PML mesh
530 
531  /// Pointer to the top right corner PML mesh
533 
534  /// Pointer to the bottom right corner PML mesh
536 
537  /// Trace file
538  ofstream Trace_file;
539 
540 }; // end of problem class
541 
542 
543 //===================start_of_create_power_monitor_mesh===================
544 /// Create BC elements on outer boundary
545 //========================================================================
546 template<class ELEMENT>
549 {
550  // Loop over outer boundaries
551  for (unsigned b=1;b<4;b++)
552  {
553  // Loop over the bulk elements adjacent to boundary b?
554  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
555  for(unsigned e=0;e<n_element;e++)
556  {
557  // Get pointer to the bulk element that is adjacent to boundary b
558  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
559  Bulk_mesh_pt->boundary_element_pt(b,e));
560 
561  //Find the index of the face of element e along boundary b
562  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
563 
564  // Build the corresponding element
565  PMLFourierDecomposedHelmholtzPowerMonitorElement<ELEMENT>*
566  flux_element_pt = new
567  PMLFourierDecomposedHelmholtzPowerMonitorElement<ELEMENT>
568  (bulk_elem_pt,face_index);
569 
570  //Add the flux boundary element
571  Power_monitor_mesh_pt->add_element_pt(flux_element_pt);
572  }
573  }
574 } // end of create_power_monitor_mesh
575 
576 
577 
578 #ifdef ADAPTIVE
579 
580 //=====================start_of_actions_before_adapt======================
581 /// Actions before adapt: Wipe the mesh of face elements
582 //========================================================================
583 template<class ELEMENT>
586 {
587  // Before adapting the added PML meshes must be removed
588  // as they are not refineable and are to be rebuilt from the
589  // newly refined triangular mesh
590  delete PML_right_mesh_pt;
591  PML_right_mesh_pt=0;
592  delete PML_top_mesh_pt;
593  PML_top_mesh_pt=0;
594  delete PML_bottom_mesh_pt;
595  PML_bottom_mesh_pt=0;
596  delete PML_top_right_mesh_pt;
597  PML_top_right_mesh_pt=0;
598  delete PML_bottom_right_mesh_pt;
599  PML_bottom_right_mesh_pt=0;
600 
601  // Wipe the power monitor elements
602  delete_face_elements(Power_monitor_mesh_pt);
603 
604  // Rebuild the Problem's global mesh from its various sub-meshes
605  // but first flush all its submeshes
606  flush_sub_meshes();
607 
608  // Then add the triangular mesh back
609  add_sub_mesh(Bulk_mesh_pt);
610 
611  // Rebuild the Problem's global mesh from its various sub-meshes
612  rebuild_global_mesh();
613 
614 }// end of actions_before_adapt
615 
616 
617 
618 //=====================start_of_actions_after_adapt=======================
619 /// Actions after adapt: Rebuild the face element meshes
620 //========================================================================
621 template<class ELEMENT>
624 {
625 
626  // Build PML meshes and add them to the global mesh
627  create_pml_meshes();
628 
629  // Re-attach the power monitor elements
630  create_power_monitor_mesh();
631 
632  // Build the entire mesh from its submeshes
633  rebuild_global_mesh();
634 
635  // Complete the build of all elements
636  complete_problem_setup();
637 
638  // Setup point source
639  setup_point_source();
640 
641 }// end of actions_after_adapt
642 
643 #endif
644 
645 
646 
647 //=================start_of_complete_problem_setup==================
648 // Complete the build of all elements so that they are fully
649 // functional
650 //==================================================================
651 template<class ELEMENT>
654 {
655  // Complete the build of all elements so they are fully functional
656  unsigned n_element = this->mesh_pt()->nelement();
657  for(unsigned i=0;i<n_element;i++)
658  {
659  // Upcast from GeneralsedElement to the present element
660  PMLFourierDecomposedHelmholtzEquations *el_pt = dynamic_cast<
661  PMLFourierDecomposedHelmholtzEquations*>(
662  mesh_pt()->element_pt(i));
663 
664  if (!(el_pt==0))
665  {
666  //Set the frequency pointer
667  el_pt->k_squared_pt()=&ProblemParameters::K_squared;
668 
669  // Set pointer to Fourier wave number
670  el_pt->pml_fourier_wavenumber_pt()=&ProblemParameters::N_fourier;
671  }
672  }
673 
674  // If the Fourier wavenumber is odd, then apply zero dirichlet boundary
675  // conditions on the two straight boundaries on the symmetry line.
676  if (ProblemParameters::N_fourier % 2 == 1)
677  {
678  cout
679  << "Zero Dirichlet boundary condition has been applied on symmetry line\n";
680  cout << "due to an odd Fourier wavenumber\n" << std::endl;
681  apply_zero_dirichlet_boundary_conditions();
682  }
683 
684 } // end of complete_problem_setup
685 
686 
687 #ifdef ADAPTIVE
688 
689 //==================start_of_setup_point_source===========================
690 /// Set point source
691 //========================================================================
692 template<class ELEMENT>
695 {
696  // Create mesh as geometric object
697  delete Mesh_as_geom_obj_pt;
698  Mesh_as_geom_obj_pt= new MeshAsGeomObject(Bulk_mesh_pt);
699 
700  // Position of point source
701  Vector<double> x_point_source;
702  x_point_source.resize(2);
703  x_point_source[0]=ProblemParameters::R_source;
704  x_point_source[1]=ProblemParameters::Z_source;
705 
706  GeomObject* sub_geom_object_pt=0;
707  Vector<double> s_point_source(2);
708  Mesh_as_geom_obj_pt->locate_zeta(x_point_source,sub_geom_object_pt,
709  s_point_source);
710 
711  // Set point force
712  if(x_point_source[0]==0.0)
713  {
715  {
716  // if source on z axis, only contribution to residual comes
717  // from Fourier wavenumber zero
718  ProblemParameters::Magnitude=complex<double>(0.0,0.0);
719  }
720  }
721 
722  // Set point force
723  dynamic_cast<ELEMENT*>(sub_geom_object_pt)->
724  setup(s_point_source,ProblemParameters::Magnitude);
725 
726 }
727 
728 
729 #endif
730 
731 //=========start_of_apply_zero_dirichlet_boundary_conditions========
732 // Apply extra bounday conditions if given an odd Fourier wavenumber
733 //==================================================================
734 template<class ELEMENT>
737 {
738  // Apply zero dirichlet conditions on the bottom straight boundary
739  // and the top straight boundary located on the symmetry line.
740 
741  // Bottom straight boundary on symmetry line:
742  {
743  //Boundary id
744  unsigned b=0;
745 
746  // How many nodes are there?
747  unsigned n_node=Bulk_mesh_pt->nboundary_node(b);
748  for (unsigned n=0;n<n_node;n++)
749  {
750  // Get the node
751  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,n);
752 
753  // Pin the node
754  nod_pt->pin(0);
755  nod_pt->pin(1);
756 
757  // Set the node's value
758  nod_pt->set_value(0, 0.0);
759  nod_pt->set_value(1, 0.0);
760  }
761  }
762 
763 // Top straight boundary on symmetry line:
764  {
765  //Boundary id
766  unsigned b=4;
767 
768  // How many nodes are there?
769  unsigned n_node=Bulk_mesh_pt->nboundary_node(b);
770  for (unsigned n=0;n<n_node;n++)
771  {
772  // Get the node
773  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,n);
774 
775  // Pin the node
776  nod_pt->pin(0);
777  nod_pt->pin(1);
778 
779  // Set the node's value
780  nod_pt->set_value(0, 0.0);
781  nod_pt->set_value(1, 0.0);
782  }
783  }
784 
785 
786 } // end of apply_zero_dirichlet_boundary_conditions
787 
788 
789 //=======start_of_constructor=============================================
790 /// Constructor for Pml Fourier-decomposed Helmholtz problem
791 //========================================================================
792 template<class ELEMENT>
795 {
796  string trace_file_location = ProblemParameters::Directory + "/trace.dat";
797 
798  // Open trace file
799  Trace_file.open(trace_file_location.c_str());
800 
801  /// Setup "bulk" mesh
802 
803  // Create the circle that represents the inner boundary
804  double x_c=0.0;
805  double y_c=0.0;
806  double r_min=1.0;
807  double r_max=3.0;
808  Circle* inner_circle_pt=new Circle(x_c,y_c,r_min);
809 
810 
811  // Edges/boundary segments making up outer boundary
812  //-------------------------------------------------
813  Vector<TriangleMeshCurveSection*> outer_boundary_line_pt(6);
814 
815  // All poly boundaries are defined by two vertices
816  Vector<Vector<double> > boundary_vertices(2);
817 
818 
819  // Bottom straight boundary on symmetry line
820  //------------------------------------------
821  boundary_vertices[0].resize(2);
822  boundary_vertices[0][0]=0.0;
823  boundary_vertices[0][1]=-r_min;
824  boundary_vertices[1].resize(2);
825  boundary_vertices[1][0]=0.0;
826  boundary_vertices[1][1]=-r_max;
827 
828  unsigned boundary_id=0;
829  outer_boundary_line_pt[0]=
830  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
831 
832 
833  // Bottom boundary of bulk mesh
834  //-----------------------------
835  boundary_vertices[0][0]=0.0;
836  boundary_vertices[0][1]=-r_max;
837  boundary_vertices[1][0]=r_max;;
838  boundary_vertices[1][1]=-r_max;
839 
840  boundary_id=1;
841  outer_boundary_line_pt[1]=
842  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
843 
844 
845  // Right boundary of bulk mesh
846  //----------------------------
847  boundary_vertices[0][0]=r_max;
848  boundary_vertices[0][1]=-r_max;
849  boundary_vertices[1][0]=r_max;;
850  boundary_vertices[1][1]=r_max;
851 
852  boundary_id=2;
853  outer_boundary_line_pt[2]=
854  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
855 
856 
857  // Top boundary of bulk mesh
858  //--------------------------
859  boundary_vertices[0][0]=r_max;
860  boundary_vertices[0][1]=r_max;
861  boundary_vertices[1][0]=0.0;
862  boundary_vertices[1][1]=r_max;
863 
864  boundary_id=3;
865  outer_boundary_line_pt[3]=
866  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
867 
868 // Top straight boundary on symmetry line
869  //---------------------------------------
870  boundary_vertices[0][0]=0.0;
871  boundary_vertices[0][1]=r_max;
872  boundary_vertices[1][0]=0.0;
873  boundary_vertices[1][1]=r_min;
874 
875  boundary_id=4;
876  outer_boundary_line_pt[4]=
877  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
878 
879 
880  // Inner circular boundary:
881  //-------------------------
882 
883  // Number of segments used for representing the curvilinear boundary
884  unsigned n_segments = 20;
885 
886  // The intrinsic coordinates for the beginning and end of the curve
887  double s_start = 0.5*MathematicalConstants::Pi;
888  double s_end = -0.5*MathematicalConstants::Pi;
889 
890  boundary_id = 5;
891  outer_boundary_line_pt[5]=
892  new TriangleMeshCurviLine(inner_circle_pt,
893  s_start,
894  s_end,
895  n_segments,
896  boundary_id);
897 
898 
899  // Create closed curve that defines outer boundary
900  //------------------------------------------------
901  TriangleMeshClosedCurve *outer_boundary_pt =
902  new TriangleMeshClosedCurve(outer_boundary_line_pt);
903 
904 
905  // Use the TriangleMeshParameters object for helping on the manage of the
906  // TriangleMesh parameters. The only parameter that needs to take is the
907  // outer boundary.
908  TriangleMeshParameters triangle_mesh_parameters(outer_boundary_pt);
909 
910 
911  // Specify maximum element area
912  double element_area = ProblemParameters::Element_area;
913  triangle_mesh_parameters.element_area() = element_area;
914 
915 #ifdef ADAPTIVE
916 
917  // Build "bulk" mesh
918  Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
919 
920  // Add the bulk mesh to the problem
921  add_sub_mesh(Bulk_mesh_pt);
922 
923  // Initialise mesh as geom object
924  Mesh_as_geom_obj_pt=0;
925 
926  // Create/set error estimator
927  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
928 
929  // Choose error tolerances to force some uniform refinement
930  Bulk_mesh_pt->min_permitted_error()=0.00004;
931  Bulk_mesh_pt->max_permitted_error()=0.0001;
932 
933  // Reduce wavenumber to make effect of singularity more prominent
935 
936 #else
937 
938  // Create the bulk mesh
939  Bulk_mesh_pt= new TriangleMesh<ELEMENT>(triangle_mesh_parameters);
940 
941  // Add the bulk mesh to the problem
942  add_sub_mesh(Bulk_mesh_pt);
943 
944  // Create flux elements on inner boundary
945  Helmholtz_inner_boundary_mesh_pt=new Mesh;
946  create_flux_elements_on_inner_boundary();
947 
948  // ...and add the mesh to the problem
949  add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
950 
951 #endif
952 
953 
954  // Attach the power monitor elements
955  Power_monitor_mesh_pt=new Mesh;
956  create_power_monitor_mesh();
957 
958  // Create the pml meshes
959  create_pml_meshes();
960 
961  // Build the Problem's global mesh from its various sub-meshes
962  build_global_mesh();
963 
964  // Complete the build of all elements
965  complete_problem_setup();
966 
967  // Setup equation numbering scheme
968  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
969 
970 } // end of constructor
971 
972 
973 
974 //===============start_of_doc=============================================
975 /// Doc the solution: doc_info contains labels/output directory etc.
976 //========================================================================
977 template<class ELEMENT>
979 doc_solution(DocInfo& doc_info)
980 {
981 
982  ofstream some_file;
983  char filename[100];
984 
985  // Number of plot points: npts x npts
986  unsigned npts=5;
987 
988 
989  // Total radiated power
990  double power=0.0;
991 
992  // Compute/output the radiated power
993  //----------------------------------
994  sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
995  doc_info.number());
996  some_file.open(filename);
997 
998  // Accumulate contribution from elements
999  unsigned nn_element=Power_monitor_mesh_pt->nelement();
1000  for(unsigned e=0;e<nn_element;e++)
1001  {
1002  PMLFourierDecomposedHelmholtzPowerMonitorElement<ELEMENT> *el_pt =
1003  dynamic_cast<PMLFourierDecomposedHelmholtzPowerMonitorElement
1004  <ELEMENT>*>(Power_monitor_mesh_pt->element_pt(e));
1005  power += el_pt->global_power_contribution(some_file);
1006  }
1007  some_file.close();
1008  oomph_info << "Total radiated power: " << power << std::endl;
1009 
1010 
1011  // Output solution within the bulk mesh
1012  //-------------------------------------
1013  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
1014  doc_info.number());
1015  some_file.open(filename);
1016  Bulk_mesh_pt->output(some_file,npts);
1017  some_file.close();
1018 
1019  // Output solution within pml domains
1020  //-----------------------------------
1021  sprintf(filename,"%s/pml_soln%i.dat",doc_info.directory().c_str(),
1022  doc_info.number());
1023  some_file.open(filename);
1024  PML_top_mesh_pt->output(some_file,npts);
1025  PML_right_mesh_pt->output(some_file,npts);
1026  PML_bottom_mesh_pt->output(some_file,npts);
1027  PML_top_right_mesh_pt->output(some_file,npts);
1028  PML_bottom_right_mesh_pt->output(some_file,npts);
1029  some_file.close();
1030 
1031 
1032  // Output exact solution
1033  //----------------------
1034  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
1035  doc_info.number());
1036  some_file.open(filename);
1037  Bulk_mesh_pt->output_fct(some_file,npts,ProblemParameters::get_exact_u);
1038  some_file.close();
1039 
1040 
1041  // Doc error and return of the square of the L2 error
1042  //---------------------------------------------------
1043  double error,norm;
1044  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
1045  doc_info.number());
1046  some_file.open(filename);
1047  Bulk_mesh_pt->compute_error(some_file,ProblemParameters::get_exact_u,
1048  error,norm);
1049  some_file.close();
1050 
1051  // Doc L2 error and norm of solution
1052  cout << "\nNorm of error : " << sqrt(error) << std::endl;
1053  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
1054 
1055  sprintf(filename,"%s/int_error%i.dat",doc_info.directory().c_str(),
1056  doc_info.number());
1057  some_file.open(filename);
1058  // Doc L2 error and norm of solution
1059  some_file << "\nNorm of error : " << sqrt(error) << std::endl;
1060  some_file << "Norm of solution: " << sqrt(norm) << std::endl <<std::endl;
1061  some_file << "Relative error: " << sqrt(error)/sqrt(norm) << std::endl;
1062  some_file.close();
1063 
1064  // Write norm and radiated power of solution to trace file
1065  Bulk_mesh_pt->compute_norm(norm);
1066  Trace_file << norm << " " << power << std::endl;
1067 
1068 } // end of doc
1069 
1070 
1071 //============start_of_create_flux_elements=================
1072 /// Create flux elements on inner boundary
1073 //==========================================================
1074 template<class ELEMENT>
1077 {
1078  // Apply flux bc on inner boundary (boundary 5)
1079  unsigned b=5;
1080 
1081 // Loop over the bulk elements adjacent to boundary b
1082  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
1083  for(unsigned e=0;e<n_element;e++)
1084  {
1085  // Get pointer to the bulk element that is adjacent to boundary b
1086  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1087  Bulk_mesh_pt->boundary_element_pt(b,e));
1088 
1089  //Find the index of the face of element e along boundary b
1090  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
1091 
1092  // Build the corresponding prescribed incoming-flux element
1093  PMLFourierDecomposedHelmholtzFluxElement<ELEMENT>*
1094  flux_element_pt = new
1095  PMLFourierDecomposedHelmholtzFluxElement<ELEMENT>
1096  (bulk_elem_pt,face_index);
1097 
1098  //Add the prescribed incoming-flux element to the surface mesh
1099  Helmholtz_inner_boundary_mesh_pt->add_element_pt(flux_element_pt);
1100 
1101  // Set the pointer to the prescribed flux function
1102  flux_element_pt->flux_fct_pt() = &ProblemParameters::exact_minus_dudr;
1103 
1104  } //end of loop over bulk elements adjacent to boundary b
1105 
1106 } // end of create flux elements on inner boundary
1107 
1108 
1109 
1110 //============start_of_create_pml_meshes======================================
1111 /// Create PML meshes and add them to the problem's sub-meshes
1112 //============================================================================
1113 template<class ELEMENT>
1115 {
1116  // Bulk mesh bottom boundary id
1117  unsigned int bottom_boundary_id = 1;
1118 
1119  // Bulk mesh right boundary id
1120  unsigned int right_boundary_id = 2;
1121 
1122  // Bulk mesh top boundary id
1123  unsigned int top_boundary_id = 3;
1124 
1125  // PML width in elements for the right layer
1126  unsigned n_x_right_pml = ProblemParameters::Nel_pml;
1127 
1128  // PML width in elements for the top layer
1129  unsigned n_y_top_pml = ProblemParameters::Nel_pml;
1130 
1131  // PML width in elements for the bottom layer
1132  unsigned n_y_bottom_pml = ProblemParameters::Nel_pml;
1133 
1134 
1135  // Outer physical length of the PML layers
1136  // defaults to 4.0
1137  double width_x_right_pml = ProblemParameters::PML_thickness;
1138  double width_y_top_pml = ProblemParameters::PML_thickness;
1139  double width_y_bottom_pml = ProblemParameters::PML_thickness;
1140 
1141  // Build the PML meshes based on the new adapted triangular mesh
1142  PML_right_mesh_pt = TwoDimensionalPMLHelper::create_right_pml_mesh
1143  <PMLLayerElement<ELEMENT> >(Bulk_mesh_pt,
1144  right_boundary_id,
1145  n_x_right_pml,
1146  width_x_right_pml);
1147  PML_top_mesh_pt = TwoDimensionalPMLHelper::create_top_pml_mesh
1148  <PMLLayerElement<ELEMENT> >(Bulk_mesh_pt,
1149  top_boundary_id,
1150  n_y_top_pml,
1151  width_y_top_pml);
1152  PML_bottom_mesh_pt= TwoDimensionalPMLHelper::create_bottom_pml_mesh
1153  <PMLLayerElement<ELEMENT> >(Bulk_mesh_pt,
1154  bottom_boundary_id,
1155  n_y_bottom_pml,
1156  width_y_bottom_pml);
1157 
1158  // Add submeshes to the global mesh
1159  add_sub_mesh(PML_right_mesh_pt);
1160  add_sub_mesh(PML_top_mesh_pt);
1161  add_sub_mesh(PML_bottom_mesh_pt);
1162 
1163  // Rebuild corner PML meshes
1164  PML_top_right_mesh_pt =
1165  TwoDimensionalPMLHelper::create_top_right_pml_mesh
1166  <PMLLayerElement<ELEMENT> >(PML_right_mesh_pt,
1167  PML_top_mesh_pt,
1168  Bulk_mesh_pt,
1169  right_boundary_id);
1170 
1171  PML_bottom_right_mesh_pt =
1172  TwoDimensionalPMLHelper::create_bottom_right_pml_mesh
1173  <PMLLayerElement<ELEMENT> >(PML_right_mesh_pt,
1174  PML_bottom_mesh_pt,
1175  Bulk_mesh_pt,
1176  right_boundary_id);
1177 
1178  // Add submeshes to the global mesh
1179  add_sub_mesh(PML_top_right_mesh_pt);
1180  add_sub_mesh(PML_bottom_right_mesh_pt);
1181 
1182 } // end of create_pml_meshes
1183 
1184 
1185 //===== start_of_main=====================================================
1186 /// Driver code for Pml Fourier decomposed Helmholtz problem
1187 //========================================================================
1188 int main(int argc, char **argv)
1189 {
1190  // Store command line arguments
1191  CommandLineArgs::setup(argc,argv);
1192 
1193  // Define possible command line arguments and parse the ones that
1194  // were actually specified
1195 
1196  // Fourier wavenumber
1197  CommandLineArgs::specify_command_line_flag("--Fourier_wavenumber",
1199 
1200  // Output directory
1201  CommandLineArgs::specify_command_line_flag("--dir",
1203 
1204  // PML thickness
1205  CommandLineArgs::specify_command_line_flag("--pml_thick",
1207 
1208  // Number of elements within PMLs
1209  CommandLineArgs::specify_command_line_flag("--npml_element",
1211 
1212  // Target Element size on first mesh generation
1213  CommandLineArgs::specify_command_line_flag("--element_area",
1215  // k squared (wavenumber squared)
1216  CommandLineArgs::specify_command_line_flag("--k_squared",
1218 
1219  // Validation run?
1220  CommandLineArgs::specify_command_line_flag("--validate");
1221 
1222 
1223  // Demonstrate across a range of omega (or k) values with good mesh for a
1224  // visual test of accuracy (put in by Jonathan Deakin)
1225  CommandLineArgs::specify_command_line_flag("--demonstrate");
1226 
1227  // Parse command line
1228  CommandLineArgs::parse_and_assign();
1229 
1230  // Doc what has actually been specified on the command line
1231  CommandLineArgs::doc_specified_flags();
1232 
1233 
1234  // Create label for output
1235  DocInfo doc_info;
1236 
1237  // Set output directory
1238  doc_info.set_directory(ProblemParameters::Directory);
1239 
1240 #ifdef ADAPTIVE
1241 
1242  // Create the problem with 2D projectable six-node elements from the
1243  // TPMLFourierDecomposedHelmholtzElement family,
1244  // allowing for the imposition of a point source (via another
1245  // templated wrapper)
1247  ProjectablePMLFourierDecomposedHelmholtzElement<
1249  TPMLFourierDecomposedHelmholtzElement<3> > > > problem;
1250 
1251 #else
1252 
1253  // Create the problem with 2D six-node elements from the
1254  // TPMLFourierDecomposedHelmholtzElement family.
1256  <TPMLFourierDecomposedHelmholtzElement<3> >
1257  problem;
1258 
1259 #endif
1260 
1261  // Step number
1262  doc_info.number()=ProblemParameters::N_fourier;
1263 
1264 #ifdef ADAPTIVE
1265 
1266  // Max. number of adaptations
1267  unsigned max_adapt=4;
1268 
1269  // Validation run?
1270  if (CommandLineArgs::command_line_flag_has_been_set("--validate"))
1271  {
1272  max_adapt=1;
1273  }
1274 
1275  // Solve the problem, allowing
1276  // up to max_adapt mesh adaptations after every solve.
1277  problem.newton_solve(max_adapt);
1278 
1279 #else
1280 
1281 
1282 // Solve the problem with Newton's method
1283 problem.newton_solve();
1284 
1285 
1286 
1287 #endif
1288 
1289  //Output the solution
1290  problem.doc_solution(doc_info);
1291 
1292 
1293 } //end of main
Mesh * Helmholtz_inner_boundary_mesh_pt
Mesh of FaceElements that apply the flux bc on the inner boundary.
~PMLHelmholtzPointSourceElement()
Destructor (empty)
unsigned N_terms
Number of terms in the exact solution.
void setup(const Vector< double > &s_point_source, const std::complex< double > &magnitude)
Set local coordinate and magnitude of point source.
Mesh * PML_top_mesh_pt
Pointer to the top PML mesh.
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete boundary face elements and wipe the surface mesh.
double Z_source
Axial position of point source.
MeshAsGeomObject * Mesh_as_geom_obj_pt
Mesh as geometric object representation of bulk mesh.
TriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
double Element_area
Target area for initial mesh.
void fill_in_point_source_contribution_to_residuals(Vector< double > &residuals)
Add the point source contribution to the residual vector.
Class to impose point source to (wrapped) Helmholtz element.
Vector< double > Coeff(N_terms, 1.0)
Coefficients in the exact solution.
double PML_thickness
Default physical PML thickness.
RefineableTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the refineable "bulk" mesh.
Mesh * PML_right_mesh_pt
Pointer to the right PML mesh.
void create_pml_meshes()
Create PML meshes.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and element Jacobian matrix (wrapper) ...
void actions_before_newton_solve()
Update the problem specs before solve (empty)
PMLLayerElement()
Constructor: Call the constructor for the appropriate Element.
unsigned Nel_pml
Default number of elements within PMLs.
double R_source
Radial position of point source.
string Directory
Output directory.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to...
std::complex< double > I(0.0, 1.0)
Imaginary unit.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
Mesh * PML_bottom_mesh_pt
Pointer to the bottom PML mesh.
void actions_after_newton_solve()
Update the problem after solve (empty)
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector of size 2, containing real and imag parts.
int main(int argc, char **argv)
Driver code for Pml Fourier decomposed Helmholtz problem.
void setup_point_source()
Set point source.
double K_squared
Frequency.
Mesh * Power_monitor_mesh_pt
Pointer to mesh that stores the power monitor elements.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
std::complex< double > Point_source_magnitude
Magnitude of point source (complex!)
int N_fourier
The default Fourier wave number.
~PMLFourierDecomposedHelmholtzProblem()
Destructor (empty)
Mesh * PML_bottom_right_mesh_pt
Pointer to the bottom right corner PML mesh.
std::complex< double > Magnitude(100.0, 100.0)
Point source magnitude (Complex)
void create_flux_elements_on_inner_boundary()
Create flux elements on inner boundary.
Vector< double > S_point_source
Local coordinates of point at which point source is applied.
void exact_minus_dudr(const Vector< double > &x, std::complex< double > &flux)
Get -du/dr (spherical r) for exact solution. Equal to prescribed flux on inner boundary.
void create_power_monitor_mesh()
Create mesh of face elements that monitor the radiated power.
Mesh * PML_top_right_mesh_pt
Pointer to the top right corner PML mesh.