26 #ifndef OOMPH_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
27 #define OOMPH_GEN_HH_TIME_HARMONIC_LIN_ELAST_HEADER
32 #include "pml_helmholtz.h"
33 #include "time_harmonic_linear_elasticity.h"
46 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
89 ELASTICITY_BULK_ELEMENT* elem_pt =
90 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
92 if (elem_pt->dim() == 3)
101 throw OomphLibError(
"This flux element will not work correctly "
102 "if nodes are hanging\n",
103 OOMPH_CURRENT_FUNCTION,
104 OOMPH_EXCEPTION_LOCATION);
120 ELASTICITY_BULK_ELEMENT* cast_element_pt =
121 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
124 for (
unsigned i = 0;
i < n_dim;
i++)
127 cast_element_pt->u_index_time_harmonic_linear_elasticity(
i);
155 const double&
q()
const
178 void output(std::ostream& outfile,
const unsigned& n_plot)
187 const double q_value =
q();
195 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
198 for (
unsigned i = 0;
i < n_dim - 1;
i++)
212 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
215 std::complex<double> u_helmholtz =
216 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
219 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
220 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
221 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
223 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
224 << ext_el_pt->interpolated_x(s_ext, 1) <<
" "
225 << traction[0].real() <<
" " << traction[1].real() <<
" "
226 << traction[0].imag() <<
" " << traction[1].imag() <<
" "
227 << interpolated_normal[0] <<
" " << interpolated_normal[1]
228 <<
" " << u_helmholtz.real() <<
" " << u_helmholtz.imag() <<
" "
231 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
234 pow(ext_el_pt->interpolated_x(s_ext, 1) -
237 <<
" " << zeta[0] << std::endl;
249 void output(FILE* file_pt,
const unsigned& n_plot)
259 std::ofstream outfile;
267 std::ofstream& outfile)
270 ELASTICITY_BULK_ELEMENT* bulk_elem_pt =
274 const unsigned bulk_dim = bulk_elem_pt->dim();
275 const unsigned local_dim = this->
dim();
285 std::complex<double> flux(0.0, 0.0);
289 if (outfile.is_open())
296 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
299 for (
unsigned i = 0;
i < local_dim;
i++)
319 bulk_elem_pt->interpolated_u_time_harmonic_linear_elasticity(s_bulk, u);
322 std::complex<double> u_dot_n(0.0, 0.0);
323 for (
unsigned i = 0;
i < bulk_dim;
i++)
325 u_dot_n += u[
i] * unit_normal[
i];
329 if (outfile.is_open())
333 outfile << x[0] <<
" " << x[1] <<
" " << u_dot_n.real() <<
" "
334 << u_dot_n.imag() <<
"\n";
350 std::ofstream outfile;
358 std::ofstream& outfile)
361 const unsigned local_dim = this->
dim();
368 std::complex<double> flux(0.0, 0.0);
371 if (outfile.is_open())
378 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
381 for (
unsigned i = 0;
i < local_dim;
i++)
390 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
394 unsigned nnode_ext = ext_el_pt->nnode();
397 const unsigned ext_dim = ext_el_pt->dim();
400 Shape psi_ext(nnode_ext);
401 DShape dpsi_ext_dx(nnode_ext, ext_dim);
406 (void)ext_el_pt->dshape_eulerian(s_bulk, psi_ext, dpsi_ext_dx);
425 ext_dim, std::complex<double>(0.0, 0.0));
426 for (
unsigned l = 0; l < nnode_ext; l++)
429 const std::complex<double> phi_value(
430 ext_el_pt->nodal_value(l, ext_el_pt->u_index_helmholtz().real()),
431 ext_el_pt->nodal_value(l, ext_el_pt->u_index_helmholtz().imag()));
434 for (
unsigned i = 0;
i < ext_dim;
i++)
436 interpolated_dphidx[
i] += phi_value * dpsi_ext_dx(l,
i);
441 std::complex<double> dphi_dn(0.0, 0.0);
442 for (
unsigned i = 0;
i < ext_dim;
i++)
444 dphi_dn += interpolated_dphidx[
i] * unit_normal[
i];
449 double max_permitted_discrepancy = 1.0e-10;
453 if (diff > max_permitted_discrepancy)
455 std::ostringstream error_stream;
457 <<
"Position computed by external el and face element differ by "
458 << diff <<
" which is more than the acceptable tolerance "
459 << max_permitted_discrepancy << std::endl;
461 OOMPH_CURRENT_FUNCTION,
462 OOMPH_EXCEPTION_LOCATION);
467 if (outfile.is_open())
471 outfile << x[0] <<
" " << x[1] <<
" " << dphi_dn.real() <<
" "
472 << dphi_dn.imag() <<
" "
494 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
495 double TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
496 ELASTICITY_BULK_ELEMENT,
497 HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value = 1.0;
503 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
504 void TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement<
505 ELASTICITY_BULK_ELEMENT,
506 HELMHOLTZ_BULK_ELEMENT>::
507 fill_in_contribution_to_residuals_helmholtz_traction(
511 unsigned n_node = nnode();
515 unsigned n_position_type = this->nnodal_position_type();
516 if (n_position_type != 1)
518 throw OomphLibError(
"LinearElasticity is not yet implemented for more "
519 "than one position type.",
520 OOMPH_CURRENT_FUNCTION,
521 OOMPH_EXCEPTION_LOCATION);
526 const unsigned n_dim = this->nodal_dimension();
529 std::complex<unsigned> u_nodal_index[n_dim];
530 for (
unsigned i = 0;
i < n_dim;
i++)
533 this->U_index_time_harmonic_linear_elasticity_helmholtz_traction[
i];
541 DShape dpsids(n_node, n_dim - 1);
544 const double q_value = q();
550 unsigned n_intpt = integral_pt()->nweight();
553 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
556 double w = integral_pt()->weight(ipt);
559 dshape_local_at_knot(ipt, psi, dpsids);
568 for (
unsigned l = 0; l < n_node; l++)
571 for (
unsigned i = 0;
i < n_dim;
i++)
574 const double x_local = nodal_position(l,
i);
575 interpolated_x[
i] += x_local * psi(l);
578 for (
unsigned j = 0; j < n_dim - 1; j++)
580 interpolated_A(j,
i) += x_local * dpsids(l, j);
587 for (
unsigned i = 0;
i < n_dim - 1;
i++)
589 for (
unsigned j = 0; j < n_dim - 1; j++)
595 for (
unsigned k = 0; k < n_dim; k++)
597 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
604 outer_unit_normal(ipt, interpolated_normal);
614 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
618 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
619 OOMPH_CURRENT_FUNCTION,
620 OOMPH_EXCEPTION_LOCATION);
625 double W = w * sqrt(Adet);
628 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
629 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(external_element_pt(0, ipt));
633 std::complex<double> u_helmholtz =
634 ext_el_pt->interpolated_u_pml_helmholtz(s_ext);
635 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
636 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
639 for (
unsigned l = 0; l < n_node; l++)
642 for (
unsigned i = 0;
i < n_dim;
i++)
644 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i].real());
649 residuals[local_eqn] -= traction[
i].real() * psi(l) *
W;
652 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i].imag());
657 residuals[local_eqn] -= traction[
i].imag() * psi(l) *
W;
680 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
723 residuals, jacobian, 1);
738 void output(std::ostream& outfile,
const unsigned& n_plot)
749 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
752 for (
unsigned i = 0;
i < (
Dim - 1);
i++)
765 ELASTICITY_BULK_ELEMENT* ext_el_pt =
769 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
774 std::complex<double> flux = (displ[0] * interpolated_normal[0] +
775 displ[1] * interpolated_normal[1]);
778 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
779 << ext_el_pt->interpolated_x(s_ext, 1) <<
" "
780 << flux.real() * interpolated_normal[0] <<
" "
781 << flux.real() * interpolated_normal[1] <<
" "
782 << flux.imag() * interpolated_normal[0] <<
" "
783 << flux.imag() * interpolated_normal[1] <<
" "
784 << interpolated_normal[0] <<
" " << interpolated_normal[1]
785 <<
" " << flux.real() <<
" " << flux.imag() <<
" "
788 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
791 pow(ext_el_pt->interpolated_x(s_ext, 1) -
794 <<
" " << zeta[0] << std::endl;
806 void output(FILE* file_pt,
const unsigned& n_plot)
821 unsigned n_node =
nnode();
827 for (
unsigned i = 0;
i < n_node;
i++)
845 unsigned n_node =
nnode();
851 for (
unsigned i = 0;
i < n_node;
i++)
868 const unsigned& flag);
887 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
889 ELASTICITY_BULK_ELEMENT>::
890 PMLHelmholtzFluxFromNormalDisplacementBCElement(
902 HELMHOLTZ_BULK_ELEMENT* elem_pt =
903 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(bulk_el_pt);
905 if (elem_pt->dim() == 3)
914 throw OomphLibError(
"This flux element will not work correctly if "
915 "nodes are hanging\n",
916 OOMPH_CURRENT_FUNCTION,
917 OOMPH_EXCEPTION_LOCATION);
954 "Bulk element must inherit from PMLHelmholtzEquations.";
956 "Nodes are one dimensional, but cannot cast the bulk element to\n";
957 error_string +=
"PMLHelmholtzEquations<1>\n.";
958 error_string +=
"If you desire this functionality, you must "
959 "implement it yourself\n";
962 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
982 "Bulk element must inherit from PMLHelmholtzEquations.";
984 "Nodes are two dimensional, but cannot cast the bulk element to\n";
985 error_string +=
"PMLHelmholtzEquations<2>\n.";
986 error_string +=
"If you desire this functionality, you must "
987 "implement it yourself\n";
990 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1009 "Bulk element must inherit from PMLHelmholtzEquations.";
1010 error_string +=
"Nodes are three dimensional, but cannot cast the "
1011 "bulk element to\n";
1012 error_string +=
"PMLHelmholtzEquations<3>\n.";
1013 error_string +=
"If you desire this functionality, you must "
1014 "implement it yourself\n";
1017 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1029 std::ostringstream error_stream;
1030 error_stream <<
"Dimension of node is " <<
Dim
1031 <<
". It should be 1,2, or 3!" << std::endl;
1034 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1044 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
1046 HELMHOLTZ_BULK_ELEMENT,
1047 ELASTICITY_BULK_ELEMENT>::
1048 fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
1051 const unsigned& flag)
1054 const unsigned n_node = nnode();
1057 Shape psif(n_node), testf(n_node);
1060 const unsigned n_intpt = integral_pt()->nweight();
1069 const std::complex<unsigned> u_index_helmholtz =
1070 U_index_helmholtz_from_displacement;
1074 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1077 for (
unsigned i = 0;
i < (
Dim - 1);
i++)
1079 s[
i] = integral_pt()->knot(ipt,
i);
1083 double w = integral_pt()->weight(ipt);
1087 double J = shape_and_test(
s, psif, testf);
1096 for (
unsigned l = 0; l < n_node; l++)
1099 for (
unsigned i = 0;
i <
Dim;
i++)
1101 interpolated_x[
i] += nodal_position(l,
i) * psif[l];
1107 outer_unit_normal(ipt, interpolated_normal);
1110 ELASTICITY_BULK_ELEMENT* ext_el_pt =
1111 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(external_element_pt(0, ipt));
1114 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
1119 std::complex<double> flux =
1120 (displ[0] * interpolated_normal[0] + displ[1] * interpolated_normal[1]);
1125 for (
unsigned l = 0; l < n_node; l++)
1127 local_eqn = nodal_local_eqn(l, u_index_helmholtz.real());
1132 residuals[local_eqn] -= flux.real() * testf[l] *
W;
1139 local_eqn = nodal_local_eqn(l, u_index_helmholtz.imag());
1144 residuals[local_eqn] -= flux.imag() * testf[l] *
W;
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
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.
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.
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
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...
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
A general Finite Element class.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
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...
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.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
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....
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
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.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
An OomphLibError object which should be thrown when an run-time error is encountered....
A class for all isoparametric elements that solve the Helmholtz equations with pml capabilities....
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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.
void fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
PMLHelmholtzFluxFromNormalDisplacementBCElement(const PMLHelmholtzFluxFromNormalDisplacementBCElement &dummy)=delete
Broken copy constructor.
void output(FILE *file_pt)
C-style output function.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
unsigned Dim
The spatial dimension of the problem.
void output(std::ostream &outfile)
Output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
std::complex< unsigned > U_index_helmholtz_from_displacement
The index at which the unknown is stored at the nodes.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Broken assignment operator.
PMLHelmholtzFluxFromNormalDisplacementBCElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Constructor, takes the pointer to the "bulk" element and the face index identifying the face to which...
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_helmholtz_traction
Index at which the i-th displacement component is stored.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
void output(FILE *file_pt)
C_style output function.
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
std::complex< double > global_flux_contribution_from_helmholtz()
Compute the global_flux_contribution through the helmholtz elements : we compute dphidn....
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
void fill_in_contribution_to_residuals_helmholtz_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
std::complex< double > global_flux_contribution_from_solid(std::ofstream &outfile)
Compute the element's contribution to the integral of the flux over the artificial boundary....
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
std::complex< double > global_flux_contribution_from_helmholtz(std::ofstream &outfile)
Compute the element's contribution to the integral of the flux over the artificial boundary....
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations....
TimeHarmonicLinElastLoadedByPMLHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
std::complex< double > global_flux_contribution_from_solid()
Compute the global_flux_contribution through the template elasticity elements : we compute u....
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void output(std::ostream &outfile)
Output function.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
void output()
Doc the command line arguments.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...