26 #ifndef OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
27 #define OOMPH_FD_HH_TIME_HARMONIC_LIN_ELAST_HEADER
32 #include "fourier_decomposed_helmholtz.h"
33 #include "time_harmonic_fourier_decomposed_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);
119 ELASTICITY_BULK_ELEMENT* cast_element_pt =
120 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
122 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
124 for (
unsigned i = 0;
i < n_dim + 1;
i++)
127 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
130 ->u_index_time_harmonic_fourier_decomposed_linear_elasticity(
i);
158 const double&
q()
const
181 void output(std::ostream& outfile,
const unsigned& n_plot)
190 const double q_value =
q();
198 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
213 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
216 std::complex<double> u_helmholtz =
217 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
220 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
221 for (
unsigned i = 0;
i < n_dim;
i++)
223 traction[
i] = -q_value * interpolated_normal[
i] * u_helmholtz;
226 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
227 << ext_el_pt->interpolated_x(s_ext, 1) <<
" " << traction[0]
228 <<
" " << traction[1] <<
" " << interpolated_normal[0] <<
" "
229 << interpolated_normal[1] <<
" " << u_helmholtz <<
" "
230 << zeta[0] << std::endl;
241 void output(FILE* file_pt,
const unsigned& n_plot)
256 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
258 FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
259 ELASTICITY_BULK_ELEMENT,
260 HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value = 1.0;
266 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
267 void FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
268 ELASTICITY_BULK_ELEMENT,
269 HELMHOLTZ_BULK_ELEMENT>::
270 fill_in_contribution_to_residuals_helmholtz_traction(
274 unsigned n_node = nnode();
278 unsigned n_position_type = this->nnodal_position_type();
279 if (n_position_type != 1)
281 throw OomphLibError(
"LinearElasticity is not yet implemented for more "
282 "than one position type.",
283 OOMPH_CURRENT_FUNCTION,
284 OOMPH_EXCEPTION_LOCATION);
289 const unsigned n_dim = this->nodal_dimension();
292 std::vector<std::complex<unsigned>> u_nodal_index(n_dim + 1);
293 for (
unsigned i = 0;
i < n_dim + 1;
i++)
297 ->U_index_fourier_decomposed_time_harmonic_linear_elasticity_helmholtz_traction
306 DShape dpsids(n_node, n_dim - 1);
309 const double q_value = q();
315 unsigned n_intpt = integral_pt()->nweight();
318 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
321 double w = integral_pt()->weight(ipt);
324 dshape_local_at_knot(ipt, psi, dpsids);
333 for (
unsigned l = 0; l < n_node; l++)
336 for (
unsigned i = 0;
i < n_dim;
i++)
339 const double x_local = nodal_position(l,
i);
340 interpolated_x[
i] += x_local * psi(l);
343 for (
unsigned j = 0; j < n_dim - 1; j++)
345 interpolated_A(j,
i) += x_local * dpsids(l, j);
351 double r = interpolated_x[0];
355 for (
unsigned i = 0;
i < n_dim - 1;
i++)
357 for (
unsigned j = 0; j < n_dim - 1; j++)
363 for (
unsigned k = 0; k < n_dim; k++)
365 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
372 outer_unit_normal(ipt, interpolated_normal);
382 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
386 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
387 "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_"
389 OOMPH_EXCEPTION_LOCATION);
394 double W = w * sqrt(Adet);
397 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
398 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(external_element_pt(0, ipt));
402 std::complex<double> u_helmholtz =
403 ext_el_pt->interpolated_u_fourier_decomposed_helmholtz(s_ext);
404 for (
unsigned i = 0;
i < n_dim;
i++)
406 traction[
i] = -q_value * interpolated_normal[
i] * u_helmholtz;
409 for (
unsigned l = 0; l < n_node; l++)
412 for (
unsigned i = 0;
i < n_dim + 1;
i++)
414 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i].real());
419 residuals[local_eqn] -= traction[
i].real() * psi(l) * r *
W;
422 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i].imag());
427 residuals[local_eqn] -= traction[
i].imag() * psi(l) * r *
W;
450 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
495 residuals, jacobian, 1);
510 void output(std::ostream& outfile,
const unsigned& n_plot)
521 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
524 for (
unsigned i = 0;
i < (
Dim - 1);
i++)
536 ELASTICITY_BULK_ELEMENT* ext_el_pt =
541 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
547 std::complex<double> flux = (displ[0] * interpolated_normal[0] +
548 displ[1] * interpolated_normal[1]);
551 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
552 << ext_el_pt->interpolated_x(s_ext, 1) <<
" "
553 << flux * interpolated_normal[0] <<
" "
554 << flux * interpolated_normal[1] <<
" "
555 << interpolated_normal[0] <<
" " << interpolated_normal[1]
556 <<
" " << flux <<
" " << zeta[0] << std::endl;
568 void output(FILE* file_pt,
const unsigned& n_plot)
583 unsigned n_node =
nnode();
589 for (
unsigned i = 0;
i < n_node;
i++)
607 unsigned n_node =
nnode();
613 for (
unsigned i = 0;
i < n_node;
i++)
630 const unsigned& flag);
649 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
651 HELMHOLTZ_BULK_ELEMENT,
652 ELASTICITY_BULK_ELEMENT>::
653 FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(
665 HELMHOLTZ_BULK_ELEMENT* elem_pt =
666 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(bulk_el_pt);
668 if (elem_pt->dim() == 3)
677 throw OomphLibError(
"This flux element will not work correctly if "
678 "nodes are hanging\n",
679 OOMPH_CURRENT_FUNCTION,
680 OOMPH_EXCEPTION_LOCATION);
714 "Bulk element must inherit from HelmholtzEquations.";
716 "Nodes are three dimensional, but cannot cast the bulk element to\n";
717 error_string +=
"HelmholtzEquations<3>\n.";
719 "If you desire this functionality, you must implement it yourself\n";
722 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
737 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
739 HELMHOLTZ_BULK_ELEMENT,
740 ELASTICITY_BULK_ELEMENT>::
741 fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
744 const unsigned& flag)
747 const unsigned n_node = nnode();
750 Shape psif(n_node), testf(n_node);
753 const unsigned n_intpt = integral_pt()->nweight();
762 const std::complex<unsigned> u_index_helmholtz =
763 U_index_helmholtz_from_displacement;
767 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
770 for (
unsigned i = 0;
i < (
Dim - 1);
i++)
772 s[
i] = integral_pt()->knot(ipt,
i);
776 double w = integral_pt()->weight(ipt);
780 double J = shape_and_test(
s, psif, testf);
789 for (
unsigned l = 0; l < n_node; l++)
792 for (
unsigned i = 0;
i <
Dim;
i++)
794 interpolated_x[
i] += nodal_position(l,
i) * psif[l];
799 double r = interpolated_x[0];
802 outer_unit_normal(ipt, interpolated_normal);
805 ELASTICITY_BULK_ELEMENT* ext_el_pt =
806 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(external_element_pt(0, ipt));
810 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
816 std::complex<double> flux =
817 (displ[0] * interpolated_normal[0] + displ[1] * interpolated_normal[1]);
822 for (
unsigned l = 0; l < n_node; l++)
824 local_eqn = nodal_local_eqn(l, u_index_helmholtz.real());
829 residuals[local_eqn] -= flux.real() * testf[l] * r *
W;
836 local_eqn = nodal_local_eqn(l, u_index_helmholtz.imag());
841 residuals[local_eqn] -= flux.imag() * testf[l] * r *
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.
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.
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 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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Broken assignment operator.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(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...
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...
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(const FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
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_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Broken assignment operator.
void output(FILE *file_pt)
C-style output function.
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 ...
unsigned Dim
The spatial dimension of the problem.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
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 ...
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
void output(FILE *file_pt)
C_style output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Plot traction etc at Gauss points nplot is ignored.
Vector< std::complex< unsigned > > U_index_fourier_decomposed_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...
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations....
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
void fill_in_contribution_to_residuals_helmholtz_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void output(std::ostream &outfile)
Output function.
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.
An OomphLibError object which should be thrown when an run-time error is encountered....
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 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 ...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...