unstructured_sphere_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 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 "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 //===== start_of_namespace_planar_wave=================================
57 /// Namespace to test representation of planar wave in spherical
58 /// polars
59 //=====================================================================
60 namespace PlanarWave
61 {
62 
63  /// Number of terms in series
64  unsigned N_terms=100;
65 
66  /// Wave number
67  double K=3.0*MathematicalConstants::Pi;
68 
69  /// Imaginary unit
70  std::complex<double> I(0.0,1.0);
71 
72  /// Exact solution as a Vector of size 2, containing real and imag parts
73  void get_exact_u(const Vector<double>& x, Vector<double>& u)
74  {
75  // Switch to spherical coordinates
76  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
77 
78  double theta;
79  theta=atan2(x[0],x[1]);
80 
81  // Argument for Bessel/Hankel functions
82  double kr = K*R;
83 
84  // Need half-order Bessel functions
85  double bessel_offset=0.5;
86 
87  // Evaluate Bessel/Hankel functions
88  Vector<double> jv(N_terms);
89  Vector<double> yv(N_terms);
90  Vector<double> djv(N_terms);
91  Vector<double> dyv(N_terms);
92  double order_max_in=double(N_terms-1)+bessel_offset;
93  double order_max_out=0;
94 
95  // This function returns vectors containing
96  // J_k(x), Y_k(x) and their derivatives
97  // up to k=order_max, with k increasing in
98  // integer increments starting with smallest
99  // positive value. So, e.g. for order_max=3.5
100  // jv[0] contains J_{1/2}(x),
101  // jv[1] contains J_{3/2}(x),
102  // jv[2] contains J_{5/2}(x),
103  // jv[3] contains J_{7/2}(x).
104  CRBond_Bessel::bessjyv(order_max_in,
105  kr,
106  order_max_out,
107  &jv[0],&yv[0],
108  &djv[0],&dyv[0]);
109 
110 
111  // Assemble exact solution (actually no need to add terms
112  // below i=N_fourier as Legendre polynomial would be zero anyway)
113  complex<double> u_ex(0.0,0.0);
114  for(unsigned i=0;i<N_terms;i++)
115  {
116  //Associated_legendre_functions
117  double p=Legendre_functions_helper::plgndr2(i,0,cos(theta));
118 
119  // Set exact solution
120  u_ex+=(2.0*i+1.0)*pow(I,i)*
121  sqrt(MathematicalConstants::Pi/(2.0*kr))*jv[i]*p;
122  }
123 
124  // Get the real & imaginary part of the result
125  u[0]=u_ex.real();
126  u[1]=u_ex.imag();
127 
128  }//end of get_exact_u
129 
130 
131  /// Plot
132  void plot()
133  {
134  unsigned nr=20;
135  unsigned nz=100;
136  unsigned nt=40;
137 
138  ofstream some_file("planar_wave.dat");
139 
140  for (unsigned i_t=0;i_t<nt;i_t++)
141  {
142  double t=2.0*MathematicalConstants::Pi*double(i_t)/double(nt-1);
143 
144  some_file << "ZONE I="<< nz << ", J="<< nr << std::endl;
145 
146  Vector<double> x(2);
147  Vector<double> u(2);
148  for (unsigned i=0;i<nr;i++)
149  {
150  x[0]=0.001+double(i)/double(nr-1);
151  for (unsigned j=0;j<nz;j++)
152  {
153  x[1]=double(j)/double(nz-1);
154  get_exact_u(x,u);
155  complex<double> uu=complex<double>(u[0],u[1])*exp(-I*t);
156  some_file << x[0] << " " << x[1] << " "
157  << uu.real() << " " << uu.imag() << "\n";
158  }
159  }
160  }
161  }
162 
163 }
164 
165 
166 /////////////////////////////////////////////////////////////////////
167 /////////////////////////////////////////////////////////////////////
168 /////////////////////////////////////////////////////////////////////
169 
170 
171 //===== start_of_namespace=============================================
172 /// Namespace for the Fourier decomposed Helmholtz problem parameters
173 //=====================================================================
174 namespace ProblemParameters
175 {
176  /// \short Square of the wavenumber
177  double K_squared=10.0;
178 
179  /// Fourier wave number
180  int N_fourier=3;
181 
182  /// Number of terms in computation of DtN boundary condition
183  unsigned Nterms_for_DtN=6;
184 
185  /// Number of terms in the exact solution
186  unsigned N_terms=6;
187 
188  /// Coefficients in the exact solution
189  Vector<double> Coeff(N_terms,1.0);
190 
191  /// Imaginary unit
192  std::complex<double> I(0.0,1.0);
193 
194  /// Exact solution as a Vector of size 2, containing real and imag parts
195  void get_exact_u(const Vector<double>& x, Vector<double>& u)
196  {
197  // Switch to spherical coordinates
198  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
199 
200  double theta;
201  theta=atan2(x[0],x[1]);
202 
203  // Argument for Bessel/Hankel functions
204  double kr = sqrt(K_squared)*R;
205 
206  // Need half-order Bessel functions
207  double bessel_offset=0.5;
208 
209  // Evaluate Bessel/Hankel functions
210  Vector<double> jv(N_terms);
211  Vector<double> yv(N_terms);
212  Vector<double> djv(N_terms);
213  Vector<double> dyv(N_terms);
214  double order_max_in=double(N_terms-1)+bessel_offset;
215  double order_max_out=0;
216 
217  // This function returns vectors containing
218  // J_k(x), Y_k(x) and their derivatives
219  // up to k=order_max, with k increasing in
220  // integer increments starting with smallest
221  // positive value. So, e.g. for order_max=3.5
222  // jv[0] contains J_{1/2}(x),
223  // jv[1] contains J_{3/2}(x),
224  // jv[2] contains J_{5/2}(x),
225  // jv[3] contains J_{7/2}(x).
226  CRBond_Bessel::bessjyv(order_max_in,
227  kr,
228  order_max_out,
229  &jv[0],&yv[0],
230  &djv[0],&dyv[0]);
231 
232 
233  // Assemble exact solution (actually no need to add terms
234  // below i=N_fourier as Legendre polynomial would be zero anyway)
235  complex<double> u_ex(0.0,0.0);
236  for(unsigned i=N_fourier;i<N_terms;i++)
237  {
238  //Associated_legendre_functions
239  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
240  cos(theta));
241  // Set exact solution
242  u_ex+=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*(jv[i]+I*yv[i])*p;
243  }
244 
245  // Get the real & imaginary part of the result
246  u[0]=u_ex.real();
247  u[1]=u_ex.imag();
248 
249  }//end of get_exact_u
250 
251 
252  /// \short Get -du/dr (spherical r) for exact solution. Equal to prescribed
253  /// flux on inner boundary.
254  void exact_minus_dudr(const Vector<double>& x, std::complex<double>& flux)
255  {
256  // Initialise flux
257  flux=std::complex<double>(0.0,0.0);
258 
259  // Switch to spherical coordinates
260  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
261 
262  double theta;
263  theta=atan2(x[0],x[1]);
264 
265  // Argument for Bessel/Hankel functions
266  double kr=sqrt(K_squared)*R;
267 
268  // Helmholtz wavenumber
269  double k=sqrt(K_squared);
270 
271  // Need half-order Bessel functions
272  double bessel_offset=0.5;
273 
274  // Evaluate Bessel/Hankel functions
275  Vector<double> jv(N_terms);
276  Vector<double> yv(N_terms);
277  Vector<double> djv(N_terms);
278  Vector<double> dyv(N_terms);
279  double order_max_in=double(N_terms-1)+bessel_offset;
280  double order_max_out=0;
281 
282  // This function returns vectors containing
283  // J_k(x), Y_k(x) and their derivatives
284  // up to k=order_max, with k increasing in
285  // integer increments starting with smallest
286  // positive value. So, e.g. for order_max=3.5
287  // jv[0] contains J_{1/2}(x),
288  // jv[1] contains J_{3/2}(x),
289  // jv[2] contains J_{5/2}(x),
290  // jv[3] contains J_{7/2}(x).
291  CRBond_Bessel::bessjyv(order_max_in,
292  kr,
293  order_max_out,
294  &jv[0],&yv[0],
295  &djv[0],&dyv[0]);
296 
297 
298  // Assemble exact solution (actually no need to add terms
299  // below i=N_fourier as Legendre polynomial would be zero anyway)
300  complex<double> u_ex(0.0,0.0);
301  for(unsigned i=N_fourier;i<N_terms;i++)
302  {
303  //Associated_legendre_functions
304  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
305  cos(theta));
306  // Set flux of exact solution
307  flux-=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*p*
308  ( k*(djv[i]+I*dyv[i]) - (0.5*(jv[i]+I*yv[i])/R) );
309  }
310 
311  }// end of exact_normal_derivative
312 
313 
314 
315 } // end of namespace
316 
317 
318 
319 /////////////////////////////////////////////////////////////////////
320 /////////////////////////////////////////////////////////////////////
321 /////////////////////////////////////////////////////////////////////
322 
323 
324 //========= start_of_problem_class=====================================
325 /// Problem class
326 //=====================================================================
327 template<class ELEMENT>
328 class FourierDecomposedHelmholtzProblem : public Problem
329 {
330 
331 public:
332 
333  /// Constructor
335 
336  /// Destructor (empty)
338 
339  /// Update the problem specs before solve (empty)
341 
342  /// Update the problem after solve (empty)
344 
345  /// \short Doc the solution. DocInfo object stores flags/labels for where the
346  /// output gets written to
347  void doc_solution(DocInfo& doc_info);
348 
349  /// Recompute gamma integral before checking Newton residuals
351  {
352  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
353  {
354  Helmholtz_outer_boundary_mesh_pt->setup_gamma();
355  }
356  }
357 
358  /// Actions before adapt: Wipe the mesh of prescribed flux elements
359  void actions_before_adapt();
360 
361  /// Actions after adapt: Rebuild the mesh of prescribed flux elements
362  void actions_after_adapt();
363 
364  /// Check gamma computation
365  void check_gamma(DocInfo& doc_info);
366 
367 private:
368 
369  /// \short Create BC elements on outer boundary
370  void create_outer_bc_elements();
371 
372  /// Create flux elements on inner boundary
373  void create_flux_elements_on_inner_boundary();
374 
375 
376  /// \short Delete boundary face elements and wipe the surface mesh
377  void delete_face_elements( Mesh* const & boundary_mesh_pt)
378  {
379  // Loop over the surface elements
380  unsigned n_element = boundary_mesh_pt->nelement();
381  for(unsigned e=0;e<n_element;e++)
382  {
383  // Kill surface element
384  delete boundary_mesh_pt->element_pt(e);
385  }
386 
387  // Wipe the mesh
388  boundary_mesh_pt->flush_element_and_node_storage();
389 
390  }
391 
392 #ifdef ADAPTIVE
393 
394  /// Pointer to the "bulk" mesh
395  RefineableTriangleMesh<ELEMENT>* Bulk_mesh_pt;
396 
397 #else
398 
399  /// Pointer to the "bulk" mesh
400  TriangleMesh<ELEMENT>* Bulk_mesh_pt;
401 
402 #endif
403 
404  /// \short Pointer to mesh containing the DtN boundary
405  /// condition elements
406  FourierDecomposedHelmholtzDtNMesh<ELEMENT>* Helmholtz_outer_boundary_mesh_pt;
407 
408  /// on the inner boundary
409  Mesh* Helmholtz_inner_boundary_mesh_pt;
410 
411  /// Trace file
412  ofstream Trace_file;
413 
414 }; // end of problem class
415 
416 
417 
418 //=====================start_of_actions_before_adapt======================
419 /// Actions before adapt: Wipe the mesh of face elements
420 //========================================================================
421 template<class ELEMENT>
423 {
424  // Kill the flux elements and wipe the boundary meshs
425  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
426  {
427  delete_face_elements(Helmholtz_outer_boundary_mesh_pt);
428  }
429  delete_face_elements(Helmholtz_inner_boundary_mesh_pt);
430 
431  // Rebuild the Problem's global mesh from its various sub-meshes
432  rebuild_global_mesh();
433 
434 }// end of actions_before_adapt
435 
436 
437 //=====================start_of_actions_after_adapt=======================
438 /// Actions after adapt: Rebuild the face element meshes
439 //========================================================================
440 template<class ELEMENT>
442 {
443 
444 
445  // Complete the build of all elements so they are fully functional
446 
447  // Loop over the Helmholtz bulk elements to set up element-specific
448  // things that cannot be handled by constructor: Pass pointer to
449  // wave number squared
450  unsigned n_element = Bulk_mesh_pt->nelement();
451  for(unsigned e=0;e<n_element;e++)
452  {
453  // Upcast from GeneralisedElement to Helmholtz bulk element
454  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
455 
456  //Set the k_squared pointer
457  el_pt->k_squared_pt() = &ProblemParameters::K_squared;
458 
459  // Set pointer to Fourier wave number
460  el_pt->fourier_wavenumber_pt()=&ProblemParameters::N_fourier;
461  }
462 
463  // Create prescribed-flux elements and BC elements
464  // from all elements that are adjacent to the boundaries and add them to
465  // Helmholtz_boundary_meshes
466  create_flux_elements_on_inner_boundary();
467  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
468  {
469  create_outer_bc_elements();
470  }
471 
472  // Rebuild the Problem's global mesh from its various sub-meshes
473  rebuild_global_mesh();
474 
475 }// end of actions_after_adapt
476 
477 
478 //=======start_of_constructor=============================================
479 /// Constructor for Fourier-decomposed Helmholtz problem
480 //========================================================================
481 template<class ELEMENT>
484 {
485 
486  // Open trace file
487  Trace_file.open("RESLT/trace.dat");
488 
489  // Create circles representing inner and outer boundary
490  double x_c=0.0;
491  double y_c=0.0;
492  double r_min=1.0;
493  double r_max=3.0;
494  Circle* inner_circle_pt=new Circle(x_c,y_c,r_min);
495  Circle* outer_circle_pt=new Circle(x_c,y_c,r_max);
496 
497  // Edges/boundary segments making up outer boundary
498  //-------------------------------------------------
499  Vector<TriangleMeshCurveSection*> outer_boundary_line_pt(4);
500 
501  // Number of segments used for representing the curvilinear boundaries
502  unsigned n_segments = 20;
503 
504  // All poly boundaries are defined by two vertices
505  Vector<Vector<double> > boundary_vertices(2);
506 
507 
508  // Bottom straight boundary on symmetry line
509  //------------------------------------------
510  boundary_vertices[0].resize(2);
511  boundary_vertices[0][0]=0.0;
512  boundary_vertices[0][1]=-r_min;
513  boundary_vertices[1].resize(2);
514  boundary_vertices[1][0]=0.0;
515  boundary_vertices[1][1]=-r_max;
516 
517  unsigned boundary_id=0;
518  outer_boundary_line_pt[0]=
519  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
520 
521 
522  if (CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
523  {
524  // Square outer boundary:
525  //-----------------------
526 
527  Vector<Vector<double> > boundary_vertices(4);
528  boundary_vertices[0].resize(2);
529  boundary_vertices[0][0]=0.0;
530  boundary_vertices[0][1]=-r_max;
531  boundary_vertices[1].resize(2);
532  boundary_vertices[1][0]=r_max;
533  boundary_vertices[1][1]=-r_max;
534  boundary_vertices[2].resize(2);
535  boundary_vertices[2][0]=r_max;
536  boundary_vertices[2][1]=r_max;
537  boundary_vertices[3].resize(2);
538  boundary_vertices[3][0]=0.0;
539  boundary_vertices[3][1]=r_max;
540 
541  boundary_id=1;
542  outer_boundary_line_pt[1]=
543  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
544  }
545  else
546  {
547  // Outer circular boundary:
548  //-------------------------
549  // The intrinsic coordinates for the beginning and end of the curve
550  double s_start = -0.5*MathematicalConstants::Pi;
551  double s_end = 0.5*MathematicalConstants::Pi;
552 
553  boundary_id = 1;
554  outer_boundary_line_pt[1]=
555  new TriangleMeshCurviLine(outer_circle_pt,
556  s_start,
557  s_end,
558  n_segments,
559  boundary_id);
560  }
561 
562 
563  // Top straight boundary on symmetry line
564  //---------------------------------------
565  boundary_vertices[0][0]=0.0;
566  boundary_vertices[0][1]=r_max;
567  boundary_vertices[1][0]=0.0;
568  boundary_vertices[1][1]=r_min;
569 
570  boundary_id=2;
571  outer_boundary_line_pt[2]=
572  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
573 
574 
575  // Inner circular boundary:
576  //-------------------------
577 
578  // The intrinsic coordinates for the beginning and end of the curve
579  double s_start = 0.5*MathematicalConstants::Pi;
580  double s_end = -0.5*MathematicalConstants::Pi;
581 
582  boundary_id = 3;
583  outer_boundary_line_pt[3]=
584  new TriangleMeshCurviLine(inner_circle_pt,
585  s_start,
586  s_end,
587  n_segments,
588  boundary_id);
589 
590 
591  // Create closed curve that defines outer boundary
592  //------------------------------------------------
593  TriangleMeshClosedCurve *outer_boundary_pt =
594  new TriangleMeshClosedCurve(outer_boundary_line_pt);
595 
596 
597  // Use the TriangleMeshParameters object for helping on the manage of the
598  // TriangleMesh parameters. The only parameter that needs to take is the
599  // outer boundary.
600  TriangleMeshParameters triangle_mesh_parameters(outer_boundary_pt);
601 
602  // Specify maximum element area
603  double element_area = 0.1;
604  triangle_mesh_parameters.element_area() = element_area;
605 
606 #ifdef ADAPTIVE
607 
608  // Build "bulk" mesh
609  Bulk_mesh_pt=new RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
610 
611  // Create/set error estimator
612  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
613 
614  // Choose error tolerances to force some uniform refinement
615  Bulk_mesh_pt->min_permitted_error()=0.00004;
616  Bulk_mesh_pt->max_permitted_error()=0.0001;
617 
618 #else
619 
620  // Pass the TriangleMeshParameters object to the TriangleMesh one
621  Bulk_mesh_pt= new TriangleMesh<ELEMENT>(triangle_mesh_parameters);
622 
623 #endif
624 
625  // Check what we've built so far...
626  Bulk_mesh_pt->output("mesh.dat");
627  Bulk_mesh_pt->output_boundaries("boundaries.dat");
628 
629 
630  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
631  {
632  // Create mesh for DtN elements on outer boundary
633  Helmholtz_outer_boundary_mesh_pt=
634  new FourierDecomposedHelmholtzDtNMesh<ELEMENT>(
636 
637  // Populate it with elements
638  create_outer_bc_elements();
639  }
640 
641  // Create flux elements on inner boundary
642  Helmholtz_inner_boundary_mesh_pt=new Mesh;
643  create_flux_elements_on_inner_boundary();
644 
645  // Add the several sub meshes to the problem
646  add_sub_mesh(Bulk_mesh_pt);
647  add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
648  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
649  {
650  add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);
651  }
652 
653  // Build the Problem's global mesh from its various sub-meshes
654  build_global_mesh();
655 
656  // Complete the build of all elements so they are fully functional
657  unsigned n_element = Bulk_mesh_pt->nelement();
658  for(unsigned i=0;i<n_element;i++)
659  {
660  // Upcast from GeneralsedElement to the present element
661  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
662 
663  //Set the k_squared pointer
664  el_pt->k_squared_pt()=&ProblemParameters::K_squared;
665 
666  // Set pointer to Fourier wave number
667  el_pt->fourier_wavenumber_pt()=&ProblemParameters::N_fourier;
668  }
669 
670  // Setup equation numbering scheme
671  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
672 
673 } // end of constructor
674 
675 
676 
677 //=================================start_of_check_gamma===================
678 /// Check gamma computation: \f$ \gamma = -du/dn \f$
679 //========================================================================
680 template<class ELEMENT>
682 {
683 
684  // Compute gamma stuff
685  Helmholtz_outer_boundary_mesh_pt->setup_gamma();
686 
687  ofstream some_file;
688  char filename[100];
689 
690  sprintf(filename,"%s/gamma_test%i.dat",doc_info.directory().c_str(),
691  doc_info.number());
692  some_file.open(filename);
693 
694  //first loop over elements e
695  unsigned nel=Helmholtz_outer_boundary_mesh_pt->nelement();
696  for (unsigned e=0;e<nel;e++)
697  {
698  // Get a pointer to element
699  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* el_pt=
700  dynamic_cast<FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>*>
701  (Helmholtz_outer_boundary_mesh_pt->element_pt(e));
702 
703  //Set the value of n_intpt
704  const unsigned n_intpt =el_pt->integral_pt()->nweight();
705 
706  // Get gamma at all gauss points in element
707  Vector<std::complex<double> > gamma(
708  Helmholtz_outer_boundary_mesh_pt->gamma_at_gauss_point(el_pt));
709 
710  //Loop over the integration points
711  for(unsigned ipt=0;ipt<n_intpt;ipt++)
712  {
713  //Allocate and initialise coordiante
714  Vector<double> x(el_pt->dim()+1,0.0);
715 
716  //Set the Vector to hold local coordinates
717  unsigned n=el_pt->dim();
718  Vector<double> s(n,0.0);
719  for(unsigned i=0;i<n;i++)
720  {
721  s[i]=el_pt->integral_pt()->knot(ipt,i);
722  }
723 
724  //Get the coordinates of the integration point
725  el_pt->interpolated_x(s,x);
726 
727  complex<double> flux;
729  some_file << atan2(x[0],x[1]) << " "
730  << gamma[ipt].real() << " "
731  << gamma[ipt].imag() << " "
732  << flux.real() << " "
733  << flux.imag() << " "
734  << std::endl;
735 
736  }// end of loop over integration points
737 
738  }// end of loop over elements
739 
740  some_file.close();
741 
742 }//end of output_gamma
743 
744 
745 //===============start_of_doc=============================================
746 /// Doc the solution: doc_info contains labels/output directory etc.
747 //========================================================================
748 template<class ELEMENT>
750 {
751 
752  ofstream some_file;
753  char filename[100];
754 
755  // Number of plot points: npts x npts
756  unsigned npts=5;
757 
758  // Output solution
759  //-----------------
760  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
761  doc_info.number());
762  some_file.open(filename);
763  Bulk_mesh_pt->output(some_file,npts);
764  some_file.close();
765 
766 
767  // Output exact solution
768  //----------------------
769  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
770  doc_info.number());
771  some_file.open(filename);
772  Bulk_mesh_pt->output_fct(some_file,npts,ProblemParameters::get_exact_u);
773  some_file.close();
774 
775 
776  // Doc error and return of the square of the L2 error
777  //---------------------------------------------------
778  double error,norm;
779  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
780  doc_info.number());
781  some_file.open(filename);
782  Bulk_mesh_pt->compute_error(some_file,ProblemParameters::get_exact_u,
783  error,norm);
784  some_file.close();
785 
786  // Doc L2 error and norm of solution
787  cout << "\nNorm of error : " << sqrt(error) << std::endl;
788  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
789 
790 
791  // Write norm of solution to trace file
792  Bulk_mesh_pt->compute_norm(norm);
793  Trace_file << norm << std::endl;
794 
795 
796  if (!CommandLineArgs::command_line_flag_has_been_set("--square_domain"))
797  {
798  // Check gamma computation
799  check_gamma(doc_info);
800  }
801 
802 } // end of doc
803 
804 
805 
806 //============start_of_create_outer_bc_elements==============================
807 /// Create BC elements on outer boundary
808 //========================================================================
809 template<class ELEMENT>
811 {
812 
813  // Outer boundary is boundary 1:
814  unsigned b=1;
815 
816  // Loop over the bulk elements adjacent to boundary b?
817  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
818  for(unsigned e=0;e<n_element;e++)
819  {
820  // Get pointer to the bulk element that is adjacent to boundary b
821  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
822  Bulk_mesh_pt->boundary_element_pt(b,e));
823 
824  //Find the index of the face of element e along boundary b
825  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
826 
827  // Build the corresponding DtN element
828  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* flux_element_pt = new
829  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>(bulk_elem_pt,
830  face_index);
831 
832  //Add the flux boundary element to the helmholtz_outer_boundary_mesh
833  Helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
834 
835  // Set pointer to the mesh that contains all the boundary condition
836  // elements on this boundary
837  flux_element_pt->
838  set_outer_boundary_mesh_pt(Helmholtz_outer_boundary_mesh_pt);
839  }
840 
841 } // end of create_outer_bc_elements
842 
843 
844 
845 //============start_of_create_flux_elements=================
846 /// Create flux elements on inner boundary
847 //==========================================================
848 template<class ELEMENT>
851 {
852  // Apply flux bc on inner boundary (boundary 3)
853  unsigned b=3;
854 
855 // Loop over the bulk elements adjacent to boundary b
856  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
857  for(unsigned e=0;e<n_element;e++)
858  {
859  // Get pointer to the bulk element that is adjacent to boundary b
860  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
861  Bulk_mesh_pt->boundary_element_pt(b,e));
862 
863  //Find the index of the face of element e along boundary b
864  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
865 
866  // Build the corresponding prescribed incoming-flux element
867  FourierDecomposedHelmholtzFluxElement<ELEMENT>* flux_element_pt = new
868  FourierDecomposedHelmholtzFluxElement<ELEMENT>(bulk_elem_pt,face_index);
869 
870  //Add the prescribed incoming-flux element to the surface mesh
871  Helmholtz_inner_boundary_mesh_pt->add_element_pt(flux_element_pt);
872 
873  // Set the pointer to the prescribed flux function
874  flux_element_pt->flux_fct_pt() = &ProblemParameters::exact_minus_dudr;
875 
876  } //end of loop over bulk elements adjacent to boundary b
877 
878 } // end of create flux elements on inner boundary
879 
880 
881 
882 //===== start_of_main=====================================================
883 /// Driver code for Fourier decomposed Helmholtz problem
884 //========================================================================
885 int main(int argc, char **argv)
886 {
887  // Store command line arguments
888  CommandLineArgs::setup(argc,argv);
889 
890  // Define possible command line arguments and parse the ones that
891  // were actually specified
892 
893  // Square domain without DtN
894  CommandLineArgs::specify_command_line_flag("--square_domain");
895 
896  // Parse command line
897  CommandLineArgs::parse_and_assign();
898 
899  // Doc what has actually been specified on the command line
900  CommandLineArgs::doc_specified_flags();
901 
902  // Check if the claimed representation of a planar wave in
903  // the tutorial is correct -- of course it is!
904  //PlanarWave::plot();
905 
906  // Test Bessel/Hankel functions
907  //-----------------------------
908  {
909  // Number of Bessel functions to be computed
910  unsigned n=3;
911 
912  // Offset of Bessel function order (less than 1!)
913  double bessel_offset=0.5;
914 
915  ofstream bessely_file("besselY.dat");
916  ofstream bessely_deriv_file("dbesselY.dat");
917 
918  ofstream besselj_file("besselJ.dat");
919  ofstream besselj_deriv_file("dbesselJ.dat");
920 
921  // Evaluate Bessel/Hankel functions
922  Vector<double> jv(n+1);
923  Vector<double> yv(n+1);
924  Vector<double> djv(n+1);
925  Vector<double> dyv(n+1);
926  double x_min=0.5;
927  double x_max=5.0;
928  unsigned nplot=100;
929  for (unsigned i=0;i<nplot;i++)
930  {
931  double x=x_min+(x_max-x_min)*double(i)/double(nplot-1);
932  double order_max_in=double(n)+bessel_offset;
933  double order_max_out=0;
934 
935  // This function returns vectors containing
936  // J_k(x), Y_k(x) and their derivatives
937  // up to k=order_max, with k increasing in
938  // integer increments starting with smallest
939  // positive value. So, e.g. for order_max=3.5
940  // jv[0] contains J_{1/2}(x),
941  // jv[1] contains J_{3/2}(x),
942  // jv[2] contains J_{5/2}(x),
943  // jv[3] contains J_{7/2}(x).
944  CRBond_Bessel::bessjyv(order_max_in,x,
945  order_max_out,
946  &jv[0],&yv[0],
947  &djv[0],&dyv[0]);
948  bessely_file << x << " ";
949  for (unsigned j=0;j<=n;j++)
950  {
951  bessely_file << yv[j] << " ";
952  }
953  bessely_file << std::endl;
954 
955  besselj_file << x << " ";
956  for (unsigned j=0;j<=n;j++)
957  {
958  besselj_file << jv[j] << " ";
959  }
960  besselj_file << std::endl;
961 
962  bessely_deriv_file << x << " ";
963  for (unsigned j=0;j<=n;j++)
964  {
965  bessely_deriv_file << dyv[j] << " ";
966  }
967  bessely_deriv_file << std::endl;
968 
969  besselj_deriv_file << x << " ";
970  for (unsigned j=0;j<=n;j++)
971  {
972  besselj_deriv_file << djv[j] << " ";
973  }
974  besselj_deriv_file << std::endl;
975 
976  }
977  bessely_file.close();
978  besselj_file.close();
979  bessely_deriv_file.close();
980  besselj_deriv_file.close();
981  }
982 
983 
984  // Test Legrendre Polynomials
985  //---------------------------
986  {
987  // Number of lower indices
988  unsigned n=3;
989 
990  ofstream some_file("legendre3.dat");
991  unsigned nplot=100;
992  for (unsigned i=0;i<nplot;i++)
993  {
994  double x=double(i)/double(nplot-1);
995 
996  some_file << x << " ";
997  for (unsigned j=0;j<=n;j++)
998  {
999  some_file << Legendre_functions_helper::plgndr2(n,j,x) << " ";
1000  }
1001  some_file << std::endl;
1002  }
1003  some_file.close();
1004  }
1005 
1006 
1007 #ifdef ADAPTIVE
1008 
1009  // Create the problem with 2D six-node elements from the
1010  // TFourierDecomposedHelmholtzElement family.
1011  FourierDecomposedHelmholtzProblem<ProjectableFourierDecomposedHelmholtzElement<
1012  TFourierDecomposedHelmholtzElement<3> > > problem;
1013 
1014 #else
1015 
1016  // Create the problem with 2D six-node elements from the
1017  // TFourierDecomposedHelmholtzElement family.
1019  problem;
1020 
1021 #endif
1022 
1023  // Create label for output
1024  DocInfo doc_info;
1025 
1026  // Set output directory
1027  doc_info.set_directory("RESLT");
1028 
1029  // Solve for a few Fourier wavenumbers
1032  {
1033  // Step number
1034  doc_info.number()=ProblemParameters::N_fourier;
1035 
1036 
1037 
1038 #ifdef ADAPTIVE
1039 
1040  // Max. number of adaptations
1041  unsigned max_adapt=1;
1042 
1043  // Solve the problem with Newton's method, allowing
1044  // up to max_adapt mesh adaptations after every solve.
1045  problem.newton_solve(max_adapt);
1046 
1047 #else
1048 
1049  // Solve the problem
1050  problem.newton_solve();
1051 
1052 #endif
1053 
1054  //Output the solution
1055  problem.doc_solution(doc_info);
1056  }
1057 
1058 } //end of main
1059 
1060 
1061 
1062 
1063 
1064 
1065 
void plot()
Plot.
RefineableTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
Vector< double > Coeff(N_terms, 1.0)
Coefficients in the exact solution.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector of size 2, containing real and imag parts.
void check_gamma(DocInfo &doc_info)
Check gamma computation.
void actions_after_newton_solve()
Update the problem after solve (empty)
int main(int argc, char **argv)
Driver code for Fourier decomposed Helmholtz problem.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector of size 2, containing real and imag parts.
void create_flux_elements_on_inner_boundary()
Create flux elements on inner boundary.
double K_squared
Square of the wavenumber.
std::complex< double > I(0.0, 1.0)
Imaginary unit.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
int N_fourier
Fourier wave number.
unsigned Nterms_for_DtN
Number of terms in computation of DtN boundary condition.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to...
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
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_outer_bc_elements()
Create BC elements on outer boundary.
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete boundary face elements and wipe the surface mesh.
TriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
unsigned N_terms
Number of terms in series.
double K
Wave number.