sphere_scattering.cc File Reference

#include <complex>
#include <cmath>
#include "generic.h"
#include "fourier_decomposed_helmholtz.h"
#include "meshes/simple_rectangular_quadmesh.h"
#include "oomph_crbond_bessel.h"

Classes
class	AnnularQuadMesh< ELEMENT >
class	FourierDecomposedHelmholtzProblem< ELEMENT >
	Problem class. More...

Namespaces
	PlanarWave
	ProblemParameters
	Namespace for exact solution and problem parameters.

Functions
std::complex< double >	PlanarWave::I (0.0, 1.0)
	Imaginary unit. More...
void	PlanarWave::get_exact_u (const Vector< double > &x, Vector< double > &u)
	Exact solution as a Vector of size 2, containing real and imag parts. More...
void	PlanarWave::plot ()
	Plot. More...
std::complex< double >	ProblemParameters::I (0.0, 1.0)
	Imaginary unit. More...
void	ProblemParameters::get_exact_u (const Vector< double > &x, Vector< double > &u)
	Exact solution as a Vector of size 2, containing real and imag parts. More...
void	ProblemParameters::exact_minus_dudr (const Vector< double > &x, std::complex< double > &flux)
int	main (int argc, char **argv)
	Driver code for Fourier decomposed Helmholtz problem. More...

Variables
unsigned	PlanarWave::N_terms =100
	Number of terms in series. More...
double	PlanarWave::K =3.0*MathematicalConstants::Pi
	Wave number. More...
double	ProblemParameters::K_squared =10.0
	Square of the wavenumber. More...
int	ProblemParameters::N_fourier =3
	Fourier wave number. More...
unsigned	ProblemParameters::Nterms_for_DtN =6
	Number of terms in computation of DtN boundary condition. More...
unsigned	ProblemParameters::N_terms =6
	Number of terms in the exact solution. More...
Vector< double >	ProblemParameters::Coeff (N_terms, 1.0)
	Coefficients in the exact solution. More...
unsigned	ProblemParameters::El_multiplier =1
	Multiplier for number of elements. More...

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

Driver code for Fourier decomposed Helmholtz problem.

{
 
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Define possible command line arguments and parse the ones that
 // were actually specified
 
 // Multiplier for number of elements
 CommandLineArgs::specify_command_line_flag("--el_multiplier",
                                            &ProblemParameters::El_multiplier);
 
 // Parse command line
 CommandLineArgs::parse_and_assign(); 
 
 // Doc what has actually been specified on the command line
 CommandLineArgs::doc_specified_flags();
 
 
 // Check if the claimed representation of a planar wave in
 // the tutorial is correct -- of course it is!
 //PlanarWave::plot();
 
 // Test Bessel/Hankel functions
 //-----------------------------
 {
  // Number of Bessel functions to be computed
  unsigned n=3;
  
  // Offset of Bessel function order (less than 1!)
  double bessel_offset=0.5;
  
  ofstream bessely_file("besselY.dat");
  ofstream bessely_deriv_file("dbesselY.dat");
  
  ofstream besselj_file("besselJ.dat");
  ofstream besselj_deriv_file("dbesselJ.dat");
  
  // Evaluate Bessel/Hankel functions
  Vector<double> jv(n+1);
  Vector<double> yv(n+1);
  Vector<double> djv(n+1);
  Vector<double> dyv(n+1);
  double x_min=0.5;
  double x_max=5.0;
  unsigned nplot=100;
  for (unsigned i=0;i<nplot;i++)
   {
    double x=x_min+(x_max-x_min)*double(i)/double(nplot-1);
    double order_max_in=double(n)+bessel_offset;
    double order_max_out=0;
    
    // This function returns vectors containing 
    // J_k(x), Y_k(x) and their derivatives
    // up to k=order_max, with k increasing in
    // integer increments starting with smallest
    // positive value. So, e.g. for order_max=3.5
    // jv[0] contains J_{1/2}(x),
    // jv[1] contains J_{3/2}(x),
    // jv[2] contains J_{5/2}(x),
    // jv[3] contains J_{7/2}(x).
    CRBond_Bessel::bessjyv(order_max_in,x,
                           order_max_out,
                           &jv[0],&yv[0],
                           &djv[0],&dyv[0]);
    bessely_file << x << " ";
    for (unsigned j=0;j<=n;j++)
     {
      bessely_file << yv[j] << " ";
     }
    bessely_file << std::endl;
    
    besselj_file << x << " ";
    for (unsigned j=0;j<=n;j++)
     {
      besselj_file << jv[j] << " ";
     }
    besselj_file << std::endl;
    
    bessely_deriv_file << x << " ";
    for (unsigned j=0;j<=n;j++)
     {
      bessely_deriv_file << dyv[j] << " ";
     }
    bessely_deriv_file << std::endl;
    
    besselj_deriv_file << x << " ";
    for (unsigned j=0;j<=n;j++)
     {
      besselj_deriv_file << djv[j] << " ";
     }
    besselj_deriv_file << std::endl;
    
   }
  bessely_file.close();
  besselj_file.close();
  bessely_deriv_file.close();
  besselj_deriv_file.close();
 }
 
 
 // Test Legendre Polynomials
 //--------------------------
 {
  // Fourier wavenumber
  unsigned n=3;
    
  ofstream some_file("legendre3.dat");
  unsigned nplot=100;
  for (unsigned i=0;i<nplot;i++)
   {
    double x=double(i)/double(nplot-1)*2.0*MathematicalConstants::Pi;
 
    some_file << x << " ";
    for (unsigned j=n;j<=5;j++)
     {
      some_file <<  Legendre_functions_helper::plgndr2(j,n,cos(x)) << " ";
     }
    some_file << std::endl;
   }
  some_file.close();
 }
 
 
 {
  ofstream some_file("legendre.dat");
  unsigned nplot=100;
  for (unsigned i=0;i<nplot;i++)
   {
    double x=double(i)/double(nplot-1);
 
    some_file << x << " ";
    for (unsigned j=0;j<=3;j++)
     {
      some_file <<  Legendre_functions_helper::plgndr2(j,0,x) << " ";
     }
    some_file << std::endl;
   }
  some_file.close();
 }
 
 
 
 // Create the problem with 2D nine-node elements from the
 // QFourierDecomposedHelmholtzElement family. 
 FourierDecomposedHelmholtzProblem<QFourierDecomposedHelmholtzElement<3> > 
  problem;
 
 // Create label for output
 DocInfo doc_info;
 
 // Set output directory
 doc_info.set_directory("RESLT");
 
 // Solve for a few Fourier wavenumbers
 for (ProblemParameters::N_fourier=0;ProblemParameters::N_fourier<4;
      ProblemParameters::N_fourier++)
  {
   // Step number
   doc_info.number()=ProblemParameters::N_fourier;
   
   // Solve the problem
   problem.newton_solve();
   
   //Output the solution
   problem.doc_solution(doc_info);
  }
 
} //end of main

References CRBond_Bessel::bessjyv(), cos(), oomph::CommandLineArgs::doc_specified_flags(), ProblemParameters::El_multiplier, i, j, n, ProblemParameters::N_fourier, oomph::DocInfo::number(), oomph::CommandLineArgs::parse_and_assign(), oomph::MathematicalConstants::Pi, oomph::Legendre_functions_helper::plgndr2(), problem, oomph::DocInfo::set_directory(), oomph::CommandLineArgs::setup(), oomph::CommandLineArgs::specify_command_line_flag(), and plotDoE::x.