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