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/simple_rectangular_quadmesh.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 // AnnularQuadMesh, derived from SimpleRectangularQuadMesh.
57 //========================================================================
58 template<class ELEMENT>
59 class AnnularQuadMesh : public SimpleRectangularQuadMesh<ELEMENT>
60 {
61 
62  public:
63 
64  // \short Constructor for angular mesh with n_r x n_phi
65  // 2D quad elements. Calls constructor for the underlying
66  // SimpleRectangularQuadMesh; then deforms the mesh so that it fits
67  // into the annular region bounded by the radii r_min and r_max
68  // and angles (in degree) of phi_min and phi_max.
69  AnnularQuadMesh(const unsigned& n_r, const unsigned& n_phi,
70  const double& r_min, const double& r_max,
71  const double& phi_min, const double& phi_max) :
72  SimpleRectangularQuadMesh<ELEMENT>(n_r,n_phi,1.0,1.0)
73  {
74 
75  // The constructor for the SimpleRectangularQuadMesh has
76  // built the mesh with n_x x n_y = n_r x n_phi elements in the unit
77  // square. Let's reposition the nodal points so that the mesh
78  // gets mapped into the required annular region:
79 
80  // Find out how many nodes there are
81  unsigned n_node=this->nnode();
82 
83  // Loop over all nodes
84  for (unsigned n=0;n<n_node;n++)
85  {
86  // Pointer to node:
87  Node* nod_pt=this->node_pt(n);
88 
89  // Get the x/y coordinates
90  double x_old=nod_pt->x(0);
91  double y_old=nod_pt->x(1);
92 
93  // Map from the old x/y to the new r/phi:
94  double r=r_min+(r_max-r_min)*x_old;
95  double phi=(phi_min+(phi_max-phi_min)*y_old)*
96  MathematicalConstants::Pi/180.0;
97 
98  // Set new nodal coordinates
99  nod_pt->x(0)=r*cos(phi);
100  nod_pt->x(1)=r*sin(phi);
101  }
102  }
103 
104 };
105 
106 
107 
108 //===== start_of_namespace_planar_wave=================================
109 /// Namespace to test representation of planar wave in spherical
110 /// polars
111 //=====================================================================
112 namespace PlanarWave
113 {
114 
115  /// Number of terms in series
116  unsigned N_terms=100;
117 
118  /// Wave number
119  double K=3.0*MathematicalConstants::Pi;
120 
121  /// Imaginary unit
122  std::complex<double> I(0.0,1.0);
123 
124  /// Exact solution as a Vector of size 2, containing real and imag parts
125  void get_exact_u(const Vector<double>& x, Vector<double>& u)
126  {
127  // Switch to spherical coordinates
128  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
129 
130  double theta;
131  theta=atan2(x[0],x[1]);
132 
133  // Argument for Bessel/Hankel functions
134  double kr = K*R;
135 
136  // Need half-order Bessel functions
137  double bessel_offset=0.5;
138 
139  // Evaluate Bessel/Hankel functions
140  Vector<double> jv(N_terms);
141  Vector<double> yv(N_terms);
142  Vector<double> djv(N_terms);
143  Vector<double> dyv(N_terms);
144  double order_max_in=double(N_terms-1)+bessel_offset;
145  double order_max_out=0;
146 
147  // This function returns vectors containing
148  // J_k(x), Y_k(x) and their derivatives
149  // up to k=order_max, with k increasing in
150  // integer increments starting with smallest
151  // positive value. So, e.g. for order_max=3.5
152  // jv[0] contains J_{1/2}(x),
153  // jv[1] contains J_{3/2}(x),
154  // jv[2] contains J_{5/2}(x),
155  // jv[3] contains J_{7/2}(x).
156  CRBond_Bessel::bessjyv(order_max_in,
157  kr,
158  order_max_out,
159  &jv[0],&yv[0],
160  &djv[0],&dyv[0]);
161 
162 
163  // Assemble exact solution (actually no need to add terms
164  // below i=N_fourier as Legendre polynomial would be zero anyway)
165  complex<double> u_ex(0.0,0.0);
166  for(unsigned i=0;i<N_terms;i++)
167  {
168  //Associated_legendre_functions
169  double p=Legendre_functions_helper::plgndr2(i,0,cos(theta));
170 
171  // Set exact solution
172  u_ex+=(2.0*i+1.0)*pow(I,i)*
173  sqrt(MathematicalConstants::Pi/(2.0*kr))*jv[i]*p;
174  }
175 
176  // Get the real & imaginary part of the result
177  u[0]=u_ex.real();
178  u[1]=u_ex.imag();
179 
180  }//end of get_exact_u
181 
182 
183  /// Plot
184  void plot()
185  {
186  unsigned nr=20;
187  unsigned nz=100;
188  unsigned nt=40;
189 
190  ofstream some_file("planar_wave.dat");
191 
192  for (unsigned i_t=0;i_t<nt;i_t++)
193  {
194  double t=2.0*MathematicalConstants::Pi*double(i_t)/double(nt-1);
195 
196  some_file << "ZONE I="<< nz << ", J="<< nr << std::endl;
197 
198  Vector<double> x(2);
199  Vector<double> u(2);
200  for (unsigned i=0;i<nr;i++)
201  {
202  x[0]=0.001+double(i)/double(nr-1);
203  for (unsigned j=0;j<nz;j++)
204  {
205  x[1]=double(j)/double(nz-1);
206  get_exact_u(x,u);
207  complex<double> uu=complex<double>(u[0],u[1])*exp(-I*t);
208  some_file << x[0] << " " << x[1] << " "
209  << uu.real() << " " << uu.imag() << "\n";
210  }
211  }
212  }
213  }
214 
215 }
216 
217 
218 /////////////////////////////////////////////////////////////////////
219 /////////////////////////////////////////////////////////////////////
220 /////////////////////////////////////////////////////////////////////
221 
222 
223 //===== start_of_namespace=============================================
224 /// Namespace for the Fourier decomposed Helmholtz problem parameters
225 //=====================================================================
226 namespace ProblemParameters
227 {
228  /// \short Square of the wavenumber
229  double K_squared=10.0;
230 
231  /// Fourier wave number
232  int N_fourier=3;
233 
234  /// Number of terms in computation of DtN boundary condition
235  unsigned Nterms_for_DtN=6;
236 
237  /// Number of terms in the exact solution
238  unsigned N_terms=6;
239 
240  /// Coefficients in the exact solution
241  Vector<double> Coeff(N_terms,1.0);
242 
243  /// Imaginary unit
244  std::complex<double> I(0.0,1.0);
245 
246  /// Exact solution as a Vector of size 2, containing real and imag parts
247  void get_exact_u(const Vector<double>& x, Vector<double>& u)
248  {
249  // Switch to spherical coordinates
250  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
251 
252  double theta;
253  theta=atan2(x[0],x[1]);
254 
255  // Argument for Bessel/Hankel functions
256  double kr = sqrt(K_squared)*R;
257 
258  // Need half-order Bessel functions
259  double bessel_offset=0.5;
260 
261  // Evaluate Bessel/Hankel functions
262  Vector<double> jv(N_terms);
263  Vector<double> yv(N_terms);
264  Vector<double> djv(N_terms);
265  Vector<double> dyv(N_terms);
266  double order_max_in=double(N_terms-1)+bessel_offset;
267  double order_max_out=0;
268 
269  // This function returns vectors containing
270  // J_k(x), Y_k(x) and their derivatives
271  // up to k=order_max, with k increasing in
272  // integer increments starting with smallest
273  // positive value. So, e.g. for order_max=3.5
274  // jv[0] contains J_{1/2}(x),
275  // jv[1] contains J_{3/2}(x),
276  // jv[2] contains J_{5/2}(x),
277  // jv[3] contains J_{7/2}(x).
278  CRBond_Bessel::bessjyv(order_max_in,
279  kr,
280  order_max_out,
281  &jv[0],&yv[0],
282  &djv[0],&dyv[0]);
283 
284 
285  // Assemble exact solution (actually no need to add terms
286  // below i=N_fourier as Legendre polynomial would be zero anyway)
287  complex<double> u_ex(0.0,0.0);
288  for(unsigned i=N_fourier;i<N_terms;i++)
289  {
290  //Associated_legendre_functions
291  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
292  cos(theta));
293  // Set exact solution
294  u_ex+=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*(jv[i]+I*yv[i])*p;
295  }
296 
297  // Get the real & imaginary part of the result
298  u[0]=u_ex.real();
299  u[1]=u_ex.imag();
300 
301  }//end of get_exact_u
302 
303 
304  /// \short Get -du/dr (spherical r) for exact solution. Equal to prescribed
305  /// flux on inner boundary.
306  void exact_minus_dudr(const Vector<double>& x, std::complex<double>& flux)
307  {
308  // Initialise flux
309  flux=std::complex<double>(0.0,0.0);
310 
311  // Switch to spherical coordinates
312  double R=sqrt(x[0]*x[0]+x[1]*x[1]);
313 
314  double theta;
315  theta=atan2(x[0],x[1]);
316 
317  // Argument for Bessel/Hankel functions
318  double kr=sqrt(K_squared)*R;
319 
320  // Helmholtz wavenumber
321  double k=sqrt(K_squared);
322 
323  // Need half-order Bessel functions
324  double bessel_offset=0.5;
325 
326  // Evaluate Bessel/Hankel functions
327  Vector<double> jv(N_terms);
328  Vector<double> yv(N_terms);
329  Vector<double> djv(N_terms);
330  Vector<double> dyv(N_terms);
331  double order_max_in=double(N_terms-1)+bessel_offset;
332  double order_max_out=0;
333 
334  // This function returns vectors containing
335  // J_k(x), Y_k(x) and their derivatives
336  // up to k=order_max, with k increasing in
337  // integer increments starting with smallest
338  // positive value. So, e.g. for order_max=3.5
339  // jv[0] contains J_{1/2}(x),
340  // jv[1] contains J_{3/2}(x),
341  // jv[2] contains J_{5/2}(x),
342  // jv[3] contains J_{7/2}(x).
343  CRBond_Bessel::bessjyv(order_max_in,
344  kr,
345  order_max_out,
346  &jv[0],&yv[0],
347  &djv[0],&dyv[0]);
348 
349 
350  // Assemble exact solution (actually no need to add terms
351  // below i=N_fourier as Legendre polynomial would be zero anyway)
352  complex<double> u_ex(0.0,0.0);
353  for(unsigned i=N_fourier;i<N_terms;i++)
354  {
355  //Associated_legendre_functions
356  double p=Legendre_functions_helper::plgndr2(i,N_fourier,
357  cos(theta));
358  // Set flux of exact solution
359  flux-=Coeff[i]*sqrt(MathematicalConstants::Pi/(2.0*kr))*p*
360  ( k*(djv[i]+I*dyv[i]) - (0.5*(jv[i]+I*yv[i])/R) );
361  }
362 
363  }// end of exact_normal_derivative
364 
365 
366  /// Multiplier for number of elements
367  unsigned El_multiplier=1;
368 
369 
370 } // end of namespace
371 
372 
373 
374 /////////////////////////////////////////////////////////////////////
375 /////////////////////////////////////////////////////////////////////
376 /////////////////////////////////////////////////////////////////////
377 
378 
379 //========= start_of_problem_class=====================================
380 /// Problem class
381 //=====================================================================
382 template<class ELEMENT>
383 class FourierDecomposedHelmholtzProblem : public Problem
384 {
385 
386 public:
387 
388  /// Constructor
390 
391  /// Destructor (empty)
393 
394  /// Update the problem specs before solve (empty)
396 
397  /// Update the problem after solve (empty)
399 
400  /// \short Doc the solution. DocInfo object stores flags/labels for where the
401  /// output gets written to
402  void doc_solution(DocInfo& doc_info);
403 
404  /// Recompute gamma integral before checking Newton residuals
406  {
407  Helmholtz_outer_boundary_mesh_pt->setup_gamma();
408  }
409 
410  /// Check gamma computation
411  void check_gamma(DocInfo& doc_info);
412 
413 private:
414 
415  /// \short Create BC elements on outer boundary
416  void create_outer_bc_elements();
417 
418  /// Create flux elements on inner boundary
419  void create_flux_elements_on_inner_boundary();
420 
421  /// Pointer to bulk mesh
423 
424  /// \short Pointer to mesh containing the DtN boundary
425  /// condition elements
426  FourierDecomposedHelmholtzDtNMesh<ELEMENT>* Helmholtz_outer_boundary_mesh_pt;
427 
428  /// \short Mesh of face elements that apply the prescribed flux
429  /// on the inner boundary
431 
432 }; // end of problem class
433 
434 
435 
436 //=======start_of_constructor=============================================
437 /// Constructor for Fourier-decomposed Helmholtz problem
438 //========================================================================
439 template<class ELEMENT>
442 {
443 
444 // Build annular mesh
445 
446 // # of elements in r-direction
447  unsigned n_r=10*ProblemParameters::El_multiplier;
448 
449  // # of elements in theta-direction
450  unsigned n_theta=10*ProblemParameters::El_multiplier;
451 
452  // Domain boundaries in theta-direction
453  double theta_min=-90.0;
454  double theta_max=90.0;
455 
456  // Domain boundaries in r-direction
457  double r_min=1.0;
458  double r_max=3.0;
459 
460  // Build and assign mesh
461  Bulk_mesh_pt =
462  new AnnularQuadMesh<ELEMENT>(n_r,n_theta,r_min,r_max,theta_min,theta_max);
463 
464  // Create mesh for DtN elements on outer boundary
465  Helmholtz_outer_boundary_mesh_pt=
466  new FourierDecomposedHelmholtzDtNMesh<ELEMENT>(
468 
469  // Populate it with elements
470  create_outer_bc_elements();
471 
472  // Create flux elements on inner boundary
473  Helmholtz_inner_boundary_mesh_pt=new Mesh;
474  create_flux_elements_on_inner_boundary();
475 
476  // Add the several sub meshes to the problem
477  add_sub_mesh(Bulk_mesh_pt);
478  add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
479  add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);
480 
481  // Build the Problem's global mesh from its various sub-meshes
482  build_global_mesh();
483 
484  // Complete the build of all elements so they are fully functional
485  unsigned n_element = Bulk_mesh_pt->nelement();
486  for(unsigned i=0;i<n_element;i++)
487  {
488  // Upcast from GeneralsedElement to the present element
489  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
490 
491  //Set the k_squared pointer
492  el_pt->k_squared_pt()=&ProblemParameters::K_squared;
493 
494  // Set pointer to Fourier wave number
495  el_pt->fourier_wavenumber_pt()=&ProblemParameters::N_fourier;
496  }
497 
498  // Setup equation numbering scheme
499  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
500 
501 } // end of constructor
502 
503 
504 
505 //=================================start_of_check_gamma===================
506 /// Check gamma computation: \f$ \gamma = -du/dn \f$
507 //========================================================================
508 template<class ELEMENT>
510 {
511 
512 
513  // Compute gamma stuff
514  Helmholtz_outer_boundary_mesh_pt->setup_gamma();
515 
516  ofstream some_file;
517  char filename[100];
518 
519  sprintf(filename,"%s/gamma_test%i.dat",doc_info.directory().c_str(),
520  doc_info.number());
521  some_file.open(filename);
522 
523  //first loop over elements e
524  unsigned nel=Helmholtz_outer_boundary_mesh_pt->nelement();
525  for (unsigned e=0;e<nel;e++)
526  {
527  // Get a pointer to element
528  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* el_pt=
529  dynamic_cast<FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>*>
530  (Helmholtz_outer_boundary_mesh_pt->element_pt(e));
531 
532  //Set the value of n_intpt
533  const unsigned n_intpt =el_pt->integral_pt()->nweight();
534 
535  // Get gamma at all gauss points in element
536  Vector<std::complex<double> > gamma(
537  Helmholtz_outer_boundary_mesh_pt->gamma_at_gauss_point(el_pt));
538 
539  //Loop over the integration points
540  for(unsigned ipt=0;ipt<n_intpt;ipt++)
541  {
542  //Allocate and initialise coordiante
543  Vector<double> x(el_pt->dim()+1,0.0);
544 
545  //Set the Vector to hold local coordinates
546  unsigned n=el_pt->dim();
547  Vector<double> s(n,0.0);
548  for(unsigned i=0;i<n;i++)
549  {
550  s[i]=el_pt->integral_pt()->knot(ipt,i);
551  }
552 
553  //Get the coordinates of the integration point
554  el_pt->interpolated_x(s,x);
555 
556  complex<double> flux;
558  some_file << atan2(x[0],x[1]) << " "
559  << gamma[ipt].real() << " "
560  << gamma[ipt].imag() << " "
561  << flux.real() << " "
562  << flux.imag() << " "
563  << std::endl;
564 
565  }// end of loop over integration points
566 
567  }// end of loop over elements
568 
569  some_file.close();
570 
571 }//end of output_gamma
572 
573 
574 //===============start_of_doc=============================================
575 /// Doc the solution: doc_info contains labels/output directory etc.
576 //========================================================================
577 template<class ELEMENT>
579 {
580 
581  ofstream some_file;
582  char filename[100];
583 
584  // Number of plot points: npts x npts
585  unsigned npts=5;
586 
587  // Output solution
588  //-----------------
589  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
590  doc_info.number());
591  some_file.open(filename);
592  Bulk_mesh_pt->output(some_file,npts);
593  some_file.close();
594 
595 
596  // Output exact solution
597  //----------------------
598  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
599  doc_info.number());
600  some_file.open(filename);
601  Bulk_mesh_pt->output_fct(some_file,npts,ProblemParameters::get_exact_u);
602  some_file.close();
603 
604 
605  // Doc error and return of the square of the L2 error
606  //---------------------------------------------------
607  double error,norm;
608  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
609  doc_info.number());
610  some_file.open(filename);
611  Bulk_mesh_pt->compute_error(some_file,ProblemParameters::get_exact_u,
612  error,norm);
613  some_file.close();
614 
615  // Doc L2 error and norm of solution
616  cout << "\nNorm of error : " << sqrt(error) << std::endl;
617  cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
618 
619  // Check gamma computation
620  check_gamma(doc_info);
621 
622 
623  // Compute/output the radiated power
624  //----------------------------------
625  sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
626  doc_info.number());
627  some_file.open(filename);
628  sprintf(filename,"%s/total_power%i.dat",doc_info.directory().c_str(),
629  doc_info.number());
630 
631  // Accumulate contribution from elements
632  double power=0.0;
633  unsigned nn_element=Helmholtz_outer_boundary_mesh_pt->nelement();
634  for(unsigned e=0;e<nn_element;e++)
635  {
636  FourierDecomposedHelmholtzBCElementBase<ELEMENT> *el_pt =
637  dynamic_cast<FourierDecomposedHelmholtzBCElementBase<ELEMENT>*>(
638  Helmholtz_outer_boundary_mesh_pt->element_pt(e));
639  power += el_pt->global_power_contribution(some_file);
640  }
641  some_file.close();
642 
643  // Output total power
644  oomph_info << "Radiated power: "
647  << power << std::endl;
648  some_file.open(filename);
649  some_file << ProblemParameters::N_fourier << " "
651  << power << std::endl;
652  some_file.close();
653 
654 } // end of doc
655 
656 
657 
658 //============start_of_create_outer_bc_elements==============================
659 /// Create BC elements on outer boundary
660 //========================================================================
661 template<class ELEMENT>
663 {
664  // Outer boundary is boundary 1:
665  unsigned b=1;
666 
667  // Loop over the bulk elements adjacent to boundary b?
668  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
669  for(unsigned e=0;e<n_element;e++)
670  {
671  // Get pointer to the bulk element that is adjacent to boundary b
672  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
673  Bulk_mesh_pt->boundary_element_pt(b,e));
674 
675  //Find the index of the face of element e along boundary b
676  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
677 
678  // Build the corresponding DtN element
679  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>* flux_element_pt = new
680  FourierDecomposedHelmholtzDtNBoundaryElement<ELEMENT>(bulk_elem_pt,
681  face_index);
682 
683  //Add the flux boundary element to the helmholtz_outer_boundary_mesh
684  Helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
685 
686  // Set pointer to the mesh that contains all the boundary condition
687  // elements on this boundary
688  flux_element_pt->
689  set_outer_boundary_mesh_pt(Helmholtz_outer_boundary_mesh_pt);
690  }
691 
692 } // end of create_outer_bc_elements
693 
694 
695 
696 //============start_of_create_flux_elements=================
697 /// Create flux elements on inner boundary
698 //==========================================================
699 template<class ELEMENT>
702 {
703  // Apply flux bc on inner boundary (boundary 3)
704  unsigned b=3;
705 
706 // Loop over the bulk elements adjacent to boundary b
707  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
708  for(unsigned e=0;e<n_element;e++)
709  {
710  // Get pointer to the bulk element that is adjacent to boundary b
711  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
712  Bulk_mesh_pt->boundary_element_pt(b,e));
713 
714  //Find the index of the face of element e along boundary b
715  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
716 
717  // Build the corresponding prescribed incoming-flux element
718  FourierDecomposedHelmholtzFluxElement<ELEMENT>* flux_element_pt = new
719  FourierDecomposedHelmholtzFluxElement<ELEMENT>(bulk_elem_pt,face_index);
720 
721  //Add the prescribed incoming-flux element to the surface mesh
722  Helmholtz_inner_boundary_mesh_pt->add_element_pt(flux_element_pt);
723 
724  // Set the pointer to the prescribed flux function
725  flux_element_pt->flux_fct_pt() = &ProblemParameters::exact_minus_dudr;
726 
727  } //end of loop over bulk elements adjacent to boundary b
728 
729 } // end of create flux elements on inner boundary
730 
731 
732 
733 //===== start_of_main=====================================================
734 /// Driver code for Fourier decomposed Helmholtz problem
735 //========================================================================
736 int main(int argc, char **argv)
737 {
738 
739  // Store command line arguments
740  CommandLineArgs::setup(argc,argv);
741 
742  // Define possible command line arguments and parse the ones that
743  // were actually specified
744 
745  // Multiplier for number of elements
746  CommandLineArgs::specify_command_line_flag("--el_multiplier",
748 
749  // Parse command line
750  CommandLineArgs::parse_and_assign();
751 
752  // Doc what has actually been specified on the command line
753  CommandLineArgs::doc_specified_flags();
754 
755 
756  // Check if the claimed representation of a planar wave in
757  // the tutorial is correct -- of course it is!
758  //PlanarWave::plot();
759 
760  // Test Bessel/Hankel functions
761  //-----------------------------
762  {
763  // Number of Bessel functions to be computed
764  unsigned n=3;
765 
766  // Offset of Bessel function order (less than 1!)
767  double bessel_offset=0.5;
768 
769  ofstream bessely_file("besselY.dat");
770  ofstream bessely_deriv_file("dbesselY.dat");
771 
772  ofstream besselj_file("besselJ.dat");
773  ofstream besselj_deriv_file("dbesselJ.dat");
774 
775  // Evaluate Bessel/Hankel functions
776  Vector<double> jv(n+1);
777  Vector<double> yv(n+1);
778  Vector<double> djv(n+1);
779  Vector<double> dyv(n+1);
780  double x_min=0.5;
781  double x_max=5.0;
782  unsigned nplot=100;
783  for (unsigned i=0;i<nplot;i++)
784  {
785  double x=x_min+(x_max-x_min)*double(i)/double(nplot-1);
786  double order_max_in=double(n)+bessel_offset;
787  double order_max_out=0;
788 
789  // This function returns vectors containing
790  // J_k(x), Y_k(x) and their derivatives
791  // up to k=order_max, with k increasing in
792  // integer increments starting with smallest
793  // positive value. So, e.g. for order_max=3.5
794  // jv[0] contains J_{1/2}(x),
795  // jv[1] contains J_{3/2}(x),
796  // jv[2] contains J_{5/2}(x),
797  // jv[3] contains J_{7/2}(x).
798  CRBond_Bessel::bessjyv(order_max_in,x,
799  order_max_out,
800  &jv[0],&yv[0],
801  &djv[0],&dyv[0]);
802  bessely_file << x << " ";
803  for (unsigned j=0;j<=n;j++)
804  {
805  bessely_file << yv[j] << " ";
806  }
807  bessely_file << std::endl;
808 
809  besselj_file << x << " ";
810  for (unsigned j=0;j<=n;j++)
811  {
812  besselj_file << jv[j] << " ";
813  }
814  besselj_file << std::endl;
815 
816  bessely_deriv_file << x << " ";
817  for (unsigned j=0;j<=n;j++)
818  {
819  bessely_deriv_file << dyv[j] << " ";
820  }
821  bessely_deriv_file << std::endl;
822 
823  besselj_deriv_file << x << " ";
824  for (unsigned j=0;j<=n;j++)
825  {
826  besselj_deriv_file << djv[j] << " ";
827  }
828  besselj_deriv_file << std::endl;
829 
830  }
831  bessely_file.close();
832  besselj_file.close();
833  bessely_deriv_file.close();
834  besselj_deriv_file.close();
835  }
836 
837 
838  // Test Legendre Polynomials
839  //--------------------------
840  {
841  // Fourier wavenumber
842  unsigned n=3;
843 
844  ofstream some_file("legendre3.dat");
845  unsigned nplot=100;
846  for (unsigned i=0;i<nplot;i++)
847  {
848  double x=double(i)/double(nplot-1)*2.0*MathematicalConstants::Pi;
849 
850  some_file << x << " ";
851  for (unsigned j=n;j<=5;j++)
852  {
853  some_file << Legendre_functions_helper::plgndr2(j,n,cos(x)) << " ";
854  }
855  some_file << std::endl;
856  }
857  some_file.close();
858  }
859 
860 
861  {
862  ofstream some_file("legendre.dat");
863  unsigned nplot=100;
864  for (unsigned i=0;i<nplot;i++)
865  {
866  double x=double(i)/double(nplot-1);
867 
868  some_file << x << " ";
869  for (unsigned j=0;j<=3;j++)
870  {
871  some_file << Legendre_functions_helper::plgndr2(j,0,x) << " ";
872  }
873  some_file << std::endl;
874  }
875  some_file.close();
876  }
877 
878 
879 
880  // Create the problem with 2D nine-node elements from the
881  // QFourierDecomposedHelmholtzElement family.
883  problem;
884 
885  // Create label for output
886  DocInfo doc_info;
887 
888  // Set output directory
889  doc_info.set_directory("RESLT");
890 
891  // Solve for a few Fourier wavenumbers
894  {
895  // Step number
896  doc_info.number()=ProblemParameters::N_fourier;
897 
898  // Solve the problem
899  problem.newton_solve();
900 
901  //Output the solution
902  problem.doc_solution(doc_info);
903  }
904 
905 } //end of main
906 
907 
908 
909 
910 
911 
912 
std::complex< double > I(0.0, 1.0)
Imaginary unit.
unsigned N_terms
Number of terms in the exact solution.
void plot()
Plot.
~FourierDecomposedHelmholtzProblem()
Destructor (empty)
Vector< double > Coeff(N_terms, 1.0)
Coefficients in the exact solution.
void check_gamma(DocInfo &doc_info)
Check gamma computation.
unsigned El_multiplier
Multiplier for number of elements.
void actions_after_newton_solve()
Update the problem after solve (empty)
int main(int argc, char **argv)
Driver code for Fourier decomposed Helmholtz problem.
AnnularQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to bulk mesh.
FourierDecomposedHelmholtzDtNMesh< ELEMENT > * Helmholtz_outer_boundary_mesh_pt
Pointer to mesh containing the DtN boundary condition elements.
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.
int N_fourier
Fourier wave number.
Mesh * Helmholtz_inner_boundary_mesh_pt
Mesh of face elements that apply the prescribed flux on the inner boundary.
unsigned Nterms_for_DtN
Number of terms in computation of DtN boundary condition.
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.
AnnularQuadMesh(const unsigned &n_r, const unsigned &n_phi, const double &r_min, const double &r_max, const double &phi_min, const double &phi_max)
double K
Wave number.