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 a specific 2D Helmholtz problem with flux boundary conditions
27//uses two separate meshes for bulk and surface mesh
28
29#include<fenv.h>
30
31
32#include "math.h"
33#include <complex>
34
35
36//Generic routines
37#include "generic.h"
38
39// The Helmholtz equations
40#include "helmholtz.h"
41
42// The mesh
43#include "meshes/annular_mesh.h"
44
45// Get the Bessel functions
46#include "oomph_crbond_bessel.h"
47
48using namespace oomph;
49using namespace std;
50
51
52
53
54/// //////////////////////////////////////////////////////////////////
55/// //////////////////////////////////////////////////////////////////
56/// //////////////////////////////////////////////////////////////////
57
58//===== start_of_namespace=============================================
59/// Namespace for the Helmholtz problem parameters
60//=====================================================================
62{
63 /// Square of the wavenumber
64 double K_squared=10.0;
65
66 /// Number of terms used in the computation
67 /// of the exact solution
68 unsigned N_fourier=10;
69
70 /// Flag to choose the Dirichlet to Neumann BC
71 /// or ABC BC
72 bool DtN_BC=false;
73
74 /// Flag to choose wich order to use
75 // in the ABCs BC: 1 for ABC 1st order...
76 unsigned ABC_order=3;
77
78 /// Radius of outer boundary (must be a circle!)
79 double Outer_radius=1.5;
80
81 /// Imaginary unit
82 std::complex<double> I(0.0,1.0);
83
84 /// Exact solution for scattered field
85 /// (vector returns real and impaginary parts).
86 void get_exact_u(const Vector<double>& x, Vector<double>& u)
87 {
88 // Switch to polar coordinates
89 double r;
90 r=sqrt(x[0]*x[0]+x[1]*x[1]);
91 double theta;
92 theta=atan2(x[1],x[0]);
93
94 // Argument for Bessel/Hankel functions
95 double rr=sqrt(K_squared)*r;
96
97 // Evaluate Bessel/Hankel functions
98 complex <double > u_ex(0.0,0.0);
99 Vector<double> jn(N_fourier+1), yn(N_fourier+1),
100 jnp(N_fourier+1), ynp(N_fourier+1);
101 Vector<double> jn_a(N_fourier+1),yn_a(N_fourier+1),
102 jnp_a(N_fourier+1), ynp_a(N_fourier+1);
103 Vector<complex<double> > h(N_fourier+1),h_a(N_fourier+1),
104 hp(N_fourier+1), hp_a(N_fourier+1);
105
106 // We want to compute N_fourier terms but the function
107 // may return fewer than that.
108 int n_actual=0;
109 CRBond_Bessel::bessjyna(N_fourier,sqrt(K_squared),n_actual,
110 &jn_a[0],&yn_a[0],
111 &jnp_a[0],&ynp_a[0]);
112
113 // Shout if things went wrong
114#ifdef PARANOID
115 if (n_actual!=int(N_fourier))
116 {
117 std::ostringstream error_stream;
118 error_stream << "CRBond_Bessel::bessjyna() only computed "
119 << n_actual << " rather than " << N_fourier
120 << " Bessel functions.\n";
121 throw OomphLibError(error_stream.str(),
122 OOMPH_CURRENT_FUNCTION,
123 OOMPH_EXCEPTION_LOCATION);
124 }
125#endif
126
127 // Evaluate Hankel at actual radius
128 Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier,rr,h,hp);
129
130 // Evaluate Hankel at inner (unit) radius
131 Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier
132 ,sqrt(K_squared),
133 h_a,hp_a);
134
135 // Compute the sum: Separate the computation of the negative
136 // and positive terms
137 for (unsigned i=0;i<N_fourier;i++)
138 {
139 u_ex-=pow(I,i)*h[i]*((jnp_a[i])/hp_a[i])*pow(exp(I*theta),i);
140 }
141 for (unsigned i=1;i<N_fourier;i++)
142 {
143 u_ex-=pow(I,i)*h[i]*((jnp_a[i])/hp_a[i])*pow(exp(-I*theta),i);
144 }
145
146 // Get the real & imaginary part of the result
147 u[0]=real(u_ex);
148 u[1]=imag(u_ex);
149
150 }// end of get_exact_u
151
152
153
154 /// Flux (normal derivative) on the unit disk
155 /// for a planar incoming wave
156 void prescribed_incoming_flux(const Vector<double>& x,
157 complex<double>& flux)
158 {
159 // Switch to polar coordinates
160 double r;
161 r=sqrt(x[0]*x[0]+x[1]*x[1]);
162 double theta;
163 theta=atan2(x[1],x[0]);
164
165 // Argument of the Bessel/Hankel fcts
166 double rr=sqrt(K_squared)*r;
167
168 // Compute Bessel/Hankel functions
169 Vector<double> jn(N_fourier+1), yn(N_fourier+1),
170 jnp(N_fourier+1), ynp(N_fourier+1);
171
172 // We want to compute N_fourier terms but the function
173 // may return fewer than that.
174 int n_actual=0;
175 CRBond_Bessel::bessjyna(N_fourier,rr,n_actual,&jn[0],&yn[0],
176 &jnp[0],&ynp[0]);
177
178 // Shout if things went wrong...
179#ifdef PARANOID
180 if (n_actual!=int(N_fourier))
181 {
182 std::ostringstream error_stream;
183 error_stream << "CRBond_Bessel::bessjyna() only computed "
184 << n_actual << " rather than " << N_fourier
185 << " Bessel functions.\n";
186 throw OomphLibError(error_stream.str(),
187 OOMPH_CURRENT_FUNCTION,
188 OOMPH_EXCEPTION_LOCATION);
189 }
190#endif
191
192 // Compute the sum: Separate the computation of the negative and
193 // positive terms
194 flux=std::complex<double>(0.0,0.0);
195 for (unsigned i=0;i<N_fourier;i++)
196 {
197 flux+=pow(I,i)*(sqrt(K_squared))*pow(exp(I*theta),i)*jnp[i];
198 }
199 for (unsigned i=1;i<N_fourier;i++)
200 {
201 flux+=pow(I,i)*(sqrt(K_squared))*pow(exp(-I*theta),i)*jnp[i];
202 }
203
204
205 }// end of prescribed_incoming_flux
206
207} // end of namespace
208
209
210
211
212/// //////////////////////////////////////////////////////////////////
213/// //////////////////////////////////////////////////////////////////
214/// //////////////////////////////////////////////////////////////////
215
216
217
218//========= start_of_problem_class=====================================
219/// Problem class to compute scattering of planar wave from unit disk
220//=====================================================================
221template<class ELEMENT>
222class ScatteringProblem : public Problem
223{
224
225public:
226
227 /// Constructor
229
230 /// Destructor (empty)
232
233 /// Doc the solution. DocInfo object stores flags/labels for where the
234 /// output gets written to
235 void doc_solution(DocInfo& doc_info);
236
237 /// Update the problem specs before solve (empty)
239
240 /// Update the problem specs after solve (empty)
242
243 /// Recompute gamma integral before checking Newton residuals
245 {
247 {
249 }
250 }
251
252 /// Actions before adapt: Wipe the mesh of prescribed flux elements
254
255 /// Actions after adapt: Rebuild the mesh of prescribed flux elements
256 void actions_after_adapt();
257
258 /// Create BC elements on boundary b of the Mesh pointed
259 /// to by bulk_mesh_pt and add them to the specified survace Mesh
261 const unsigned &b, Mesh* const &bulk_mesh_pt,
262 Mesh* const & helmholtz_outer_boundary_mesh_pt);
263
264 /// Create Helmholtz flux elements on boundary b of the Mesh pointed
265 /// to by bulk_mesh_pt and add them to the specified surface Mesh
266 void create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
267 Mesh* const & helmholtz_inner_boundary_mesh_pt);
268
269 /// Delete boundary face elements and wipe the surface mesh
270 void delete_face_elements( Mesh* const & boundary_mesh_pt);
271
272 /// Set pointer to prescribed-flux function for all
273 /// elements in the surface mesh on the surface of the unit disk
275
276 /// Set up boundary condition elements on outer boundary
278
279#ifdef ADAPTIVE
280
281/// Pointer to the "bulk" mesh
282 RefineableTwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;
283
284#else
285
286 /// Pointer to the "bulk" mesh
287 TwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;
288
289#endif
290
291 /// Pointer to mesh containing the DtN (or ABC) boundary
292 /// condition elements
293 HelmholtzDtNMesh<ELEMENT>* Helmholtz_outer_boundary_mesh_pt;
294
295 /// Pointer to the mesh containing
296 /// the Helmholtz inner boundary condition elements
298
299}; // end of problem class
300
301
302
303//=======start_of_constructor=============================================
304/// Constructor for Helmholtz problem
305//========================================================================
306template<class ELEMENT>
309{
310
311 // Setup "bulk" mesh
312
313 // # of elements in theta
314 unsigned n_theta=15;
315
316 // # of elements in radius
317 unsigned n_r=5;
318
319 // Inner radius
320 double a=1.0;
321
322 // Thickness of annular computational domain
323 double h=0.5;
324
325 // Set outer radius
327
328 // Mesh is periodic
329 bool periodic=true;
330
331 // Full circle
332 double azimuthal_fraction=1.0;
333
334#ifdef ADAPTIVE
335
336 // Build "bulk" mesh
337 Bulk_mesh_pt=
338 new RefineableTwoDAnnularMesh<ELEMENT>(periodic,
339 azimuthal_fraction,n_theta,n_r,a,h);
340
341 // Create/set error estimator
342 Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
343
344 // Choose error tolerances to force some uniform refinement
345 Bulk_mesh_pt->min_permitted_error()=0.004;
346 Bulk_mesh_pt->max_permitted_error()=0.01;
347
348#else
349
350 // Build "bulk" mesh
351 Bulk_mesh_pt=
352 new TwoDAnnularMesh<ELEMENT>(periodic,
353 azimuthal_fraction,n_theta,n_r,a,h);
354
355#endif
356
357 // Pointer to mesh containing the Helmholtz outer boundary condition
358 // elements. Specify outer radius and number of Fourier terms to be
359 // used in gamma integral
360 Helmholtz_outer_boundary_mesh_pt =
361 new HelmholtzDtNMesh<ELEMENT>(a+h,GlobalParameters::N_fourier);
362
363 // Pointer to mesh containing the Helmholtz inner boundary condition
364 // elements. Specify outer radius
365 Helmholtz_inner_boundary_mesh_pt = new Mesh;
366
367 // Create prescribed-flux elements from all elements that are
368 // adjacent to the inner boundary , but add them to a separate mesh.
369 create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
370
371 // Create outer boundary elements from all elements that are
372 // adjacent to the outer boundary , but add them to a separate mesh.
373 create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);
374
375 // Add the several sub meshes to the problem
376 add_sub_mesh(Bulk_mesh_pt);
377 add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);
378 add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);
379
380 // Build the Problem's global mesh from its various sub-meshes
381 build_global_mesh();
382
383 // Complete the build of all elements so they are fully functional
384
385 // Loop over the Helmholtz bulk elements to set up element-specific
386 // things that cannot be handled by constructor: Pass pointer to
387 // wave number squared
388 unsigned n_element = Bulk_mesh_pt->nelement();
389 for(unsigned e=0;e<n_element;e++)
390 {
391 // Upcast from GeneralisedElement to Helmholtz bulk element
392 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
393
394 //Set the k_squared pointer
395 el_pt->k_squared_pt() = &GlobalParameters::K_squared;
396 }
397
398 // Set up elements on outer boundary
399 setup_outer_boundary();
400
401 // Set pointer to prescribed flux function for flux elements
402 set_prescribed_incoming_flux_pt();
403
404 // Setup equation numbering scheme
405 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
406
407} // end of constructor
408
409//=====================start_of_actions_before_adapt======================
410/// Actions before adapt: Wipe the mesh of face elements
411//========================================================================
412template<class ELEMENT>
414{
415 // Kill the flux elements and wipe the boundary meshs
416 delete_face_elements(Helmholtz_outer_boundary_mesh_pt);
417 delete_face_elements(Helmholtz_inner_boundary_mesh_pt);
418
419 // Rebuild the Problem's global mesh from its various sub-meshes
420 rebuild_global_mesh();
421
422}// end of actions_before_adapt
423
424
425//=====================start_of_actions_after_adapt=======================
426/// Actions after adapt: Rebuild the face element meshes
427//========================================================================
428template<class ELEMENT>
430{
431 // Create prescribed-flux elements and BC elements
432 // from all elements that are adjacent to the boundaries and add them to
433 // Helmholtz_boundary_meshes
434 create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);
435 create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
436
437 // Rebuild the Problem's global mesh from its various sub-meshes
438 rebuild_global_mesh();
439
440 // Set pointer to prescribed flux function and DtN mesh
441 setup_outer_boundary();
442 set_prescribed_incoming_flux_pt();
443
444}// end of actions_after_adapt
445
446
447//==================start_of_setup_outer_boundary=========================
448/// Set pointers for elements on outer boundary
449//========================================================================
450template<class ELEMENT>
452{
453 // Loop over the flux elements to pass pointer to DtN
454 // BC for the outer boundary
455 unsigned n_element=Helmholtz_outer_boundary_mesh_pt->nelement();
456 for(unsigned e=0;e<n_element;e++)
457 {
458 // Dirichlet to Neumann BC
460 {
461 // Upcast from GeneralisedElement to Helmholtz flux element
462 HelmholtzDtNBoundaryElement<ELEMENT> *el_pt =
463 dynamic_cast< HelmholtzDtNBoundaryElement<ELEMENT>*>(
464 Helmholtz_outer_boundary_mesh_pt->element_pt(e));
465
466 // Set pointer to the mesh that contains all the boundary condition
467 // elements on this boundary
468 el_pt->set_outer_boundary_mesh_pt(Helmholtz_outer_boundary_mesh_pt);
469 }
470 // ABCs BC
471 else
472 {
473 // Upcast from GeneralisedElement to appropriate type
474 HelmholtzAbsorbingBCElement<ELEMENT> *el_pt =
475 dynamic_cast< HelmholtzAbsorbingBCElement<ELEMENT>*>(
476 Helmholtz_outer_boundary_mesh_pt->element_pt(e));
477
478 // Set pointer to outer radius of artificial boundary
479 el_pt->outer_radius_pt()=&GlobalParameters::Outer_radius;
480
481 // Set order of absorbing boundary condition
482 el_pt->abc_order_pt()=&GlobalParameters::ABC_order;
483 }
484 }
485}
486
487
488//==================start_of_set_prescribed_incoming_flux_pt==============
489/// Set pointer to prescribed incoming-flux function for all
490/// elements in the inner boundary
491//========================================================================
492template<class ELEMENT>
494{
495 // Loop over the flux elements to pass pointer to prescribed
496 // flux function for the inner boundary
497 unsigned n_element=Helmholtz_inner_boundary_mesh_pt->nelement();
498 for(unsigned e=0;e<n_element;e++)
499 {
500 // Upcast from GeneralisedElement to Helmholtz flux element
501 HelmholtzFluxElement<ELEMENT> *el_pt =
502 dynamic_cast< HelmholtzFluxElement<ELEMENT>*>(
503 Helmholtz_inner_boundary_mesh_pt->element_pt(e));
504
505 // Set the pointer to the prescribed flux function
506 el_pt->flux_fct_pt() =
508 }
509
510}// end of set prescribed_incoming_flux pt
511
512//=====================start_of_doc=======================================
513/// Doc the solution: doc_info contains labels/output directory etc.
514//========================================================================
515template<class ELEMENT>
517 doc_info)
518{
519
520 ofstream some_file,some_file2;
521 char filename[100];
522
523 // Number of plot points
524 unsigned npts;
525 npts=5;
526
527 // Compute/output the radiated power
528 //----------------------------------
529 sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
530 doc_info.number());
531 some_file.open(filename);
532
533 // Accumulate contribution from elements
534 double power=0.0;
535 unsigned nn_element=Helmholtz_outer_boundary_mesh_pt->nelement();
536 for(unsigned e=0;e<nn_element;e++)
537 {
538 HelmholtzBCElementBase<ELEMENT> *el_pt =
539 dynamic_cast< HelmholtzBCElementBase<ELEMENT>*>(
540 Helmholtz_outer_boundary_mesh_pt->element_pt(e));
541 power += el_pt->global_power_contribution(some_file);
542 }
543 some_file.close();
544 oomph_info << "Total radiated power: " << power << std::endl;
545
546 // Output solution
547 //-----------------
548 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
549 doc_info.number());
550 some_file.open(filename);
551 Bulk_mesh_pt->output(some_file,npts);
552 some_file.close();
553
554 // Output exact solution
555 //----------------------
556 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
557 doc_info.number());
558 some_file.open(filename);
559 Bulk_mesh_pt->output_fct(some_file,npts,GlobalParameters::get_exact_u);
560 some_file.close();
561
562 double error,norm;
563 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
564 doc_info.number());
565 some_file.open(filename);
566 Bulk_mesh_pt->compute_error(some_file,GlobalParameters::get_exact_u,
567 error,norm);
568 some_file.close();
569
570 // Doc L2 error and norm of solution
571 oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
572 oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;
573
574
575 // Do animation of Helmholtz solution
576 //-----------------------------------
577 unsigned nstep=40;
578 for (unsigned i=0;i<nstep;i++)
579 {
580 sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",
581 doc_info.directory().c_str(),
582 doc_info.number(),i);
583 some_file.open(filename);
584 sprintf(filename,"%s/exact_helmholtz_animation%i_frame%i.dat",
585 doc_info.directory().c_str(),
586 doc_info.number(),i);
587 some_file2.open(filename);
588 double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);
589 unsigned nelem=Bulk_mesh_pt->nelement();
590 for (unsigned e=0;e<nelem;e++)
591 {
592 ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
593 Bulk_mesh_pt->element_pt(e));
594 el_pt->output_real(some_file,phi,npts);
595 el_pt->output_real_fct(some_file2,phi,npts,
597 }
598 some_file.close();
599 some_file2.close();
600 }
601
602} // end of doc
603
604//============start_of_create_flux_elements==================================
605/// Create Helmholtz inner Flux Elements on the b-th boundary of
606/// the Mesh object pointed to by bulk_mesh_pt and add the elements
607/// to the Mesh object pointed to by helmholtz_inner_boundary_mesh_pt
608//============================================================================
609template<class ELEMENT>
611create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
612 Mesh* const & helmholtz_inner_boundary_mesh_pt)
613{
614 // Loop over the bulk elements adjacent to boundary b
615 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
616 for(unsigned e=0;e<n_element;e++)
617 {
618 // Get pointer to the bulk element that is adjacent to boundary b
619 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
620 bulk_mesh_pt->boundary_element_pt(b,e));
621
622 //Find the index of the face of element e along boundary b
623 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
624
625 // Build the corresponding prescribed incoming-flux element
626 HelmholtzFluxElement<ELEMENT>* flux_element_pt = new
627 HelmholtzFluxElement<ELEMENT>(bulk_elem_pt,face_index);
628
629 //Add the prescribed incoming-flux element to the surface mesh
630 helmholtz_inner_boundary_mesh_pt->add_element_pt(flux_element_pt);
631
632 } //end of loop over bulk elements adjacent to boundary b
633
634} // end of create_flux_elements
635
636
637
638//============start_of_create_outer_bc_elements==============================
639/// Create outer BC elements on the b-th boundary of
640/// the Mesh object pointed to by bulk_mesh_pt and add the elements
641/// to the Mesh object pointed to by helmholtz_outer_boundary_mesh_pt.
642//===========================================================================
643template<class ELEMENT>
645create_outer_bc_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
646 Mesh* const & helmholtz_outer_boundary_mesh_pt)
647{
648 // Loop over the bulk elements adjacent to boundary b?
649 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
650 for(unsigned e=0;e<n_element;e++)
651 {
652 // Get pointer to the bulk element that is adjacent to boundary b
653 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
654 bulk_mesh_pt->boundary_element_pt(b,e));
655
656 //Find the index of the face of element e along boundary b
657 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
658
659 // Build the corresponding outer flux element
660
661 // Dirichlet to Neumann boundary conditon
663 {
664 HelmholtzDtNBoundaryElement<ELEMENT>* flux_element_pt = new
665 HelmholtzDtNBoundaryElement<ELEMENT>(bulk_elem_pt,face_index);
666
667 //Add the flux boundary element to the helmholtz_outer_boundary_mesh
668 helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
669 }
670 // ABCs BC
671 else
672 {
673 HelmholtzAbsorbingBCElement<ELEMENT>* flux_element_pt = new
674 HelmholtzAbsorbingBCElement<ELEMENT>(bulk_elem_pt,face_index);
675
676 //Add the flux boundary element to the helmholtz_outer_boundary_mesh
677 helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
678 }
679 } //end of loop over bulk elements adjacent to boundary b
680} // end of create_outer_bc_elements
681
682
683//============start_of_delete_face_elements================
684/// Delete face elements and wipe the boundary mesh
685//==========================================================
686template<class ELEMENT>
688delete_face_elements(Mesh* const & boundary_mesh_pt)
689{
690 // Loop over the surface elements
691 unsigned n_element = boundary_mesh_pt->nelement();
692 for(unsigned e=0;e<n_element;e++)
693 {
694 // Kill surface element
695 delete boundary_mesh_pt->element_pt(e);
696 }
697
698 // Wipe the mesh
699 boundary_mesh_pt->flush_element_and_node_storage();
700
701} // end of delete_outer_face_elements
702
703
704
705//==========start_of_main=================================================
706/// Solve 2D Helmholtz problem for scattering of a planar wave from a
707/// unit disk
708//========================================================================
709int main(int argc, char **argv)
710{
711
712 // Store command line arguments
713 CommandLineArgs::setup(argc,argv);
714
715 // Define case to be run
716 unsigned i_case=0;
717 CommandLineArgs::specify_command_line_flag("--case",&i_case);
718
719 // Parse command line
720 CommandLineArgs::parse_and_assign();
721
722 // Doc what has actually been specified on the command line
723 CommandLineArgs::doc_specified_flags();
724
725 // Now set flags accordingly
726 switch(i_case)
727 {
728 case 0:
730 break;
731
732 case 1:
735 break;
736
737 case 2:
740 break;
741
742 case 3:
745 break;
746 }
747
748
749 //Set up the problem
750 //------------------
751
752#ifdef ADAPTIVE
753
754 //Set up the problem with 2D nine-node elements from the
755 //QHelmholtzElement family.
757 problem;
758
759#else
760
761 //Set up the problem with 2D nine-node elements from the
762 //QHelmholtzElement family.
764 problem;
765
766
767#endif
768
769 // Create label for output
770 //------------------------
771 DocInfo doc_info;
772
773 // Set output directory
774 doc_info.set_directory("RESLT");
775
776
777#ifdef ADAPTIVE
778
779 // Max. number of adaptations
780 unsigned max_adapt=1;
781
782 // Solve the problem with Newton's method, allowing
783 // up to max_adapt mesh adaptations after every solve.
784 problem.newton_solve(max_adapt);
785
786#else
787
788 // Solve the problem with Newton's method
789 problem.newton_solve();
790
791#endif
792
793 //Output solution
794 problem.doc_solution(doc_info);
795
796} //end of main
797
798
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: scattering.cc:223
~ScatteringProblem()
Destructor (empty)
Definition: scattering.cc:231
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
Definition: scattering.cc:429
void create_outer_bc_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_outer_boundary_mesh_pt)
Create BC elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to the specified...
Definition: scattering.cc:645
Mesh * Helmholtz_inner_boundary_mesh_pt
Pointer to the mesh containing the Helmholtz inner boundary condition elements.
Definition: scattering.cc:297
void setup_outer_boundary()
Set up boundary condition elements on outer boundary.
Definition: scattering.cc:451
RefineableTwoDAnnularMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
Definition: scattering.cc:282
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Definition: scattering.cc:238
ScatteringProblem()
Constructor.
Definition: scattering.cc:308
HelmholtzDtNMesh< ELEMENT > * Helmholtz_outer_boundary_mesh_pt
Pointer to mesh containing the DtN (or ABC) boundary condition elements.
Definition: scattering.cc:293
TwoDAnnularMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
Definition: scattering.cc:287
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete boundary face elements and wipe the surface mesh.
Definition: scattering.cc:688
void actions_after_newton_solve()
Update the problem specs after solve (empty)
Definition: scattering.cc:241
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
Definition: scattering.cc:413
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.
Definition: scattering.cc:244
void set_prescribed_incoming_flux_pt()
Set pointer to prescribed-flux function for all elements in the surface mesh on the surface of the un...
Definition: scattering.cc:493
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_inner_boundary_mesh_pt)
Create Helmholtz flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to t...
Definition: scattering.cc:611
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to.
Definition: scattering.cc:516
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: scattering.cc:62
void prescribed_incoming_flux(const Vector< double > &x, complex< double > &flux)
Flux (normal derivative) on the unit disk for a planar incoming wave.
Definition: scattering.cc:156
unsigned ABC_order
Flag to choose wich order to use.
Definition: scattering.cc:76
double Outer_radius
Radius of outer boundary (must be a circle!)
Definition: scattering.cc:79
bool DtN_BC
Flag to choose the Dirichlet to Neumann BC or ABC BC.
Definition: scattering.cc:72
std::complex< double > I(0.0, 1.0)
Imaginary unit.
double K_squared
Square of the wavenumber.
Definition: scattering.cc:64
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution for scattered field (vector returns real and impaginary parts).
Definition: scattering.cc:86
unsigned N_fourier
Number of terms used in the computation of the exact solution.
Definition: scattering.cc:68
int main(int argc, char **argv)
Solve 2D Helmholtz problem for scattering of a planar wave from a unit disk.
Definition: scattering.cc:709