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>
48 :
public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
49 public virtual FaceElement,
50 public virtual ElementWithExternalElement
63 Vector<std::complex<unsigned>>
71 Vector<double>& residuals);
77 FiniteElement*
const& element_pt,
const int& face_index)
78 : FaceGeometry<ELASTICITY_BULK_ELEMENT>(),
84 element_pt->build_face_element(face_index,
this);
89 ELASTICITY_BULK_ELEMENT* elem_pt =
90 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
92 if (elem_pt->dim() == 3)
95 RefineableElement* ref_el_pt =
96 dynamic_cast<RefineableElement*
>(elem_pt);
99 if (this->has_hanging_nodes())
101 throw OomphLibError(
"This flux element will not work correctly "
102 "if nodes are hanging\n",
103 OOMPH_CURRENT_FUNCTION,
104 OOMPH_EXCEPTION_LOCATION);
113 this->set_ninteraction(1);
116 unsigned n_dim = element_pt->nodal_dimension();
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);
144 DenseMatrix<double>& jacobian)
150 fill_in_jacobian_from_external_interaction_by_fd(residuals, jacobian);
158 const double&
q()
const
181 void output(std::ostream& outfile,
const unsigned& n_plot)
184 unsigned n_dim = this->nodal_dimension();
187 Vector<std::complex<double>> traction(n_dim + 1);
190 const double q_value =
q();
195 unsigned n_intpt = integral_pt()->nweight();
198 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
200 Vector<double> s_int(n_dim - 1);
201 s_int[0] = integral_pt()->knot(ipt, 0);
204 Vector<double> interpolated_normal(n_dim);
205 outer_unit_normal(ipt, interpolated_normal);
208 Vector<double> zeta(1);
209 interpolated_zeta(s_int, zeta);
213 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
214 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(external_element_pt(0, ipt));
215 Vector<double> s_ext(external_element_local_coord(0, ipt));
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;
237 FaceGeometry<ELASTICITY_BULK_ELEMENT>::output(file_pt);
241 void output(FILE* file_pt,
const unsigned& n_plot)
243 FaceGeometry<ELASTICITY_BULK_ELEMENT>::output(file_pt, 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(
271 Vector<double>& residuals)
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();
312 Vector<std::complex<double>> traction(n_dim + 1);
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);
327 Vector<double> interpolated_x(n_dim, 0.0);
330 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
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];
354 DenseMatrix<double> A(n_dim - 1);
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);
371 Vector<double> interpolated_normal(n_dim);
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));
399 Vector<double> s_ext(external_element_local_coord(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>
452 :
public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
453 public virtual FaceElement,
454 public virtual ElementWithExternalElement
460 FiniteElement*
const& bulk_el_pt,
const int& face_index);
484 residuals, GeneralisedElement::Dummy_matrix, 0);
491 DenseMatrix<double>& jacobian)
495 residuals, jacobian, 1);
498 fill_in_jacobian_from_external_interaction_by_fd(residuals, jacobian);
510 void output(std::ostream& outfile,
const unsigned& n_plot)
515 const unsigned n_intpt = integral_pt()->nweight();
518 Vector<double> s_int(
Dim - 1);
521 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
524 for (
unsigned i = 0; i < (
Dim - 1); i++)
526 s_int[i] = integral_pt()->knot(ipt, i);
530 Vector<double> interpolated_normal(
Dim);
531 outer_unit_normal(ipt, interpolated_normal);
533 Vector<double> zeta(1);
534 interpolated_zeta(s_int, zeta);
536 ELASTICITY_BULK_ELEMENT* ext_el_pt =
537 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(external_element_pt(0, ipt));
538 Vector<double> s_ext(external_element_local_coord(0, ipt));
539 Vector<std::complex<double>> displ(
Dim + 1);
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;
564 FaceGeometry<HELMHOLTZ_BULK_ELEMENT>::output(file_pt);
568 void output(FILE* file_pt,
const unsigned& n_plot)
570 FaceGeometry<HELMHOLTZ_BULK_ELEMENT>::output(file_pt, n_plot);
583 unsigned n_node = nnode();
589 for (
unsigned i = 0; i < n_node; i++)
595 return J_eulerian(s);
607 unsigned n_node = nnode();
610 shape_at_knot(ipt, psi);
613 for (
unsigned i = 0; i < n_node; i++)
619 return J_eulerian_at_knot(ipt);
628 Vector<double>& residuals,
629 DenseMatrix<double>& jacobian,
630 const unsigned& flag);
649 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
651 HELMHOLTZ_BULK_ELEMENT,
652 ELASTICITY_BULK_ELEMENT>::
653 FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement(
654 FiniteElement*
const& bulk_el_pt,
const int& face_index)
655 : FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
660 bulk_el_pt->build_face_element(face_index,
this);
665 HELMHOLTZ_BULK_ELEMENT* elem_pt =
666 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(bulk_el_pt);
668 if (elem_pt->dim() == 3)
671 RefineableElement* ref_el_pt =
672 dynamic_cast<RefineableElement*
>(elem_pt);
675 if (this->has_hanging_nodes())
677 throw OomphLibError(
"This flux element will not work correctly if "
678 "nodes are hanging\n",
679 OOMPH_CURRENT_FUNCTION,
680 OOMPH_EXCEPTION_LOCATION);
690 this->set_ninteraction(1);
695 Dim = bulk_el_pt->nodal_dimension();
708 FourierDecomposedHelmholtzEquations* eqn_pt =
709 dynamic_cast<FourierDecomposedHelmholtzEquations*
>(bulk_el_pt);
713 std::string error_string =
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);
728 eqn_pt->u_index_fourier_decomposed_helmholtz();
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(
742 Vector<double>& residuals,
743 DenseMatrix<double>& jacobian,
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();
756 Vector<double> s(Dim - 1);
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);
786 Vector<double> interpolated_x(Dim, 0.0);
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];
801 Vector<double> interpolated_normal(Dim);
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));
807 Vector<double> s_ext(external_element_local_coord(0, ipt));
808 Vector<std::complex<double>> displ(Dim + 1);
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;
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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.