28#ifndef OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER
29#define OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_BC_ELEMENTS_HEADER
34#include <oomph-lib-config.h>
41#include "oomph_crbond_bessel.h"
58 template<
class ELEMENT>
73 "FourierDecomposedHelmholtzBCElementBase",
98 const unsigned&
i)
const
138 return std::complex<unsigned>(
223 std::complex<double>
dphi_dn(0.0, 0.0);
236 l,
bulk_elem_pt->u_index_fourier_decomposed_helmholtz().real()),
238 l,
bulk_elem_pt->u_index_fourier_decomposed_helmholtz().imag()));
295 for (
unsigned i = 0;
i <
nnod;
i++)
323 for (
unsigned j = 0;
j < 1;
j++)
349 template<
class ELEMENT>
408 std::map<FiniteElement*, Vector<std::map<unsigned, std::complex<double>>>>
419 template<
class ELEMENT>
461 std::map<
unsigned, std::complex<double>>&
d_gamma_con);
487 for (
unsigned e = 0;
e <
nel;
e++)
491 for (
unsigned j = 0;
j <
nnod;
j++)
496 if (!(
nod_pt->is_a_copy()))
504 for (
unsigned j = 0;
j <
nnod;
j++)
519 for (std::set<Node*>::iterator
it =
node_set.begin();
535 const unsigned&
flag)
570 for (
unsigned i = 0;
i < 1;
i++)
604 if (local_eqn_real >= 0)
770 template<
class ELEMENT>
776 std::map<
unsigned, std::complex<double>>&
d_gamma_con)
780 dynamic_cast<ELEMENT*
>(this->bulk_element_pt())->fourier_wavenumber();
783 const std::complex<double> I(0.0, 1.0);
803 gamma_con = std::complex<double>(0.0, 0.0);
811 for (
unsigned i = 0;
i < 1;
i++)
827 std::complex<double> interpolated_u(0.0, 0.0);
833 for (
unsigned i = 0;
i < 2;
i++)
842 this->U_index_fourier_decomposed_helmholtz.real()),
844 this->U_index_fourier_decomposed_helmholtz.imag()));
873 l, this->U_index_fourier_decomposed_helmholtz.real());
883 l, this->U_index_fourier_decomposed_helmholtz.imag());
904 namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
917 namespace ToleranceForFourierDecomposedHelmholtzOuterBoundary
935 template<
class ELEMENT>
941 unsigned nel = this->nelement();
942 for (
unsigned e = 0;
e <
nel;
e++)
946 for (
unsigned j = 0;
j <
nnod;
j++)
956 double r =
sqrt(x[0] * x[0] + x[1] * x[1]);
959 if (Outer_radius == 0.0)
961 throw OomphLibError(
"Outer radius for DtN BC must not be zero!",
966 if (std::fabs((
r - this->Outer_radius) / Outer_radius) >
971 <<
"Node at " << x[0] <<
" " << x[1] <<
" has radius " <<
r
972 <<
" which does not "
973 <<
" agree with \nspecified outer radius " << this->Outer_radius
974 <<
" within relative tolerance "
976 <<
".\nYou can adjust the tolerance via\n"
977 <<
"ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol\n"
978 <<
"or recompile without PARANOID.\n";
992 this->element_pt(0));
994 sqrt(
dynamic_cast<ELEMENT*
>(
el_pt->bulk_element_pt())->k_squared());
996 dynamic_cast<ELEMENT*
>(
el_pt->bulk_element_pt())->fourier_wavenumber();
1001 std::complex<double> I(0.0, 1.0);
1005 q(N_terms + 1, std::complex<double>(0.0, 0.0));
1019 for (
unsigned j = 0;
j < N_terms + 1;
j++)
1034 (
hp_a[
i] -
h_a[
i] / (2.0 *
k * Outer_radius)) *
1042 unsigned nel = this->nelement();
1043 for (
unsigned e = 0;
e <
nel;
e++)
1048 this->element_pt(
e));
1055 std::complex<double>(0.0, 0.0));
1078 std::complex<double>
gamma_con(0.0, 0.0);
1079 std::map<unsigned, std::complex<double>>
d_gamma_con;
1091 this->element_pt(
ee));
1094 eel_pt->compute_gamma_contribution(
1104 l,
eel_pt->u_index_fourier_decomposed_helmholtz().real());
1115 l,
eel_pt->u_index_fourier_decomposed_helmholtz().imag());
1138 template<
class ELEMENT>
1141 const int& face_index)
1162 "Bulk element must inherit from FourierDecomposedHelmholtzEquations.";
1171 eqn_pt->u_index_fourier_decomposed_helmholtz();
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement.
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
FaceGeometry class definition: This policy class is used to allow construction of face elements that ...
A general Finite Element class.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
unsigned nnode() const
Return the number of nodes.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
A class for elements that allow the approximation of the Sommerfeld radiation BC for Fourier decompos...
FourierDecomposedHelmholtzBCElementBase()
Broken empty constructor.
void output(FILE *file_pt)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
double d_shape_and_test_local(const Vector< double > &s, Shape &psi, Shape &test, DShape &dpsi_ds, DShape &dtest_ds) const
Function to compute the shape, test functions and derivates and to return the Jacobian of mapping bet...
double global_power_contribution(std::ofstream &outfile)
Compute the element's contribution to the time-averaged radiated power over the artificial boundary....
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the test functions and to return the Jacobian of mapping between local and global...
FourierDecomposedHelmholtzBCElementBase(const FourierDecomposedHelmholtzBCElementBase &dummy)=delete
Broken copy constructor.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Broken assignment operator.
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Return the index at which the real/imag unknown value is stored.
std::complex< unsigned > U_index_fourier_decomposed_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
double global_power_contribution()
Compute the element's contribution to the time-averaged radiated power over the artificial boundary.
FaceElement used to apply Sommerfeld radiation conditon via Dirichlet to Neumann map.
FourierDecomposedHelmholtzDtNMesh< ELEMENT > * outer_boundary_mesh_pt() const
Access function to mesh of all DtN boundary condition elements (needed to get access to gamma values)
void fill_in_generic_residual_contribution_fourier_decomposed_helmholtz_DtN_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the element's residual vector Jacobian matrix. Overloaded version, using the gamma computed i...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
FourierDecomposedHelmholtzDtNBoundaryElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Construct element from specification of bulk element and face index.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void complete_setup_of_dependencies()
Complete the setup of additional dependencies arising through the far-away interaction with other nod...
FourierDecomposedHelmholtzDtNMesh< ELEMENT > * Outer_boundary_mesh_pt
Pointer to mesh of all DtN boundary condition elements (needed to get access to gamma values)
void set_outer_boundary_mesh_pt(FourierDecomposedHelmholtzDtNMesh< ELEMENT > *mesh_pt)
Set mesh of all DtN boundary condition elements.
void compute_gamma_contribution(const double &theta, const unsigned &n, std::complex< double > &gamma_con, std::map< unsigned, std::complex< double > > &d_gamma_con)
Compute the contribution of the element to the Gamma integral and its derivates w....
================================================================= Mesh for DtN boundary condition ele...
std::map< FiniteElement *, Vector< std::complex< double > > > Gamma_at_gauss_point
Container to store the gamma integral for given Gauss point and element.
unsigned & n_terms()
Number of terms used in the computation of the gamma integral.
FourierDecomposedHelmholtzDtNMesh(const double &outer_radius, const unsigned &n_terms)
Constructor: Specify radius of outer boundary and number of terms used in the computation of the gamm...
Vector< std::complex< double > > & gamma_at_gauss_point(FiniteElement *el_pt)
Gamma integral evaluated at Gauss points for specified element.
double & outer_radius()
The outer radius.
void setup_gamma()
Compute and store the gamma integral at all integration points of the constituent elements.
unsigned N_terms
Number of terms used in the Gamma computation.
std::map< FiniteElement *, Vector< std::map< unsigned, std::complex< double > > > > D_Gamma_at_gauss_point
Container to store the derivate of Gamma integral w.r.t global unknown evaluated at Gauss points for ...
double Outer_radius
Outer radius.
Vector< std::map< unsigned, std::complex< double > > > & d_gamma_at_gauss_point(FiniteElement *el_pt)
Derivative of Gamma integral w.r.t global unknown, evaluated at Gauss points for specified element.
A class for all isoparametric elements that solve the Helmholtz equations.
unsigned nexternal_data() const
Return the number of external data objects.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data.
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
double factorial(const unsigned &l)
Factorial.
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
const double Pi
50 digits from maple
double Tol
Relative tolerance to within radius of points on DtN boundary are allowed to deviate from specified v...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).