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