30 #include "helmholtz.h"
31 #include "time_harmonic_linear_elasticity.h"
44 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
46 :
public virtual FaceGeometry<ELASTICITY_BULK_ELEMENT>,
47 public virtual FaceElement,
48 public virtual ElementWithExternalElement
61 Vector<std::complex<unsigned>>
69 Vector<double>& residuals);
75 FiniteElement*
const& element_pt,
const int& face_index)
76 : FaceGeometry<ELASTICITY_BULK_ELEMENT>(),
82 element_pt->build_face_element(face_index,
this);
87 ELASTICITY_BULK_ELEMENT* elem_pt =
88 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
90 if (elem_pt->dim() == 3)
93 RefineableElement* ref_el_pt =
94 dynamic_cast<RefineableElement*
>(elem_pt);
97 if (this->has_hanging_nodes())
99 throw OomphLibError(
"This flux element will not work correctly "
100 "if nodes are hanging\n",
101 OOMPH_CURRENT_FUNCTION,
102 OOMPH_EXCEPTION_LOCATION);
111 this->set_ninteraction(1);
114 unsigned n_dim = element_pt->nodal_dimension();
117 ELASTICITY_BULK_ELEMENT* cast_element_pt =
118 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(element_pt);
121 for (
unsigned i = 0; i < n_dim; i++)
124 cast_element_pt->u_index_time_harmonic_linear_elasticity(i);
138 DenseMatrix<double>& jacobian)
144 fill_in_jacobian_from_external_interaction_by_fd(residuals, jacobian);
152 const double&
q()
const
175 void output(std::ostream& outfile,
const unsigned& n_plot)
178 unsigned n_dim = this->nodal_dimension();
181 Vector<std::complex<double>> traction(n_dim);
184 const double q_value =
q();
189 unsigned n_intpt = integral_pt()->nweight();
192 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
194 Vector<double> s_int(n_dim - 1);
195 for (
unsigned i = 0; i < n_dim - 1; i++)
197 s_int[i] = integral_pt()->knot(ipt, i);
201 Vector<double> interpolated_normal(n_dim);
202 outer_unit_normal(ipt, interpolated_normal);
205 Vector<double> zeta(1);
206 interpolated_zeta(s_int, zeta);
209 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
210 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(external_element_pt(0, ipt));
211 Vector<double> s_ext(external_element_local_coord(0, ipt));
212 std::complex<double> u_helmholtz =
213 ext_el_pt->interpolated_u_helmholtz(s_ext);
216 ext_el_pt->interpolated_u_helmholtz(s_ext);
217 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
218 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
220 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
221 << ext_el_pt->interpolated_x(s_ext, 1) <<
" "
222 << traction[0].real() <<
" " << traction[1].real() <<
" "
223 << traction[0].imag() <<
" " << traction[1].imag() <<
" "
224 << interpolated_normal[0] <<
" " << interpolated_normal[1]
225 <<
" " << u_helmholtz.real() <<
" " << u_helmholtz.imag() <<
" "
226 << interpolated_x(s_int, 0) <<
" " << interpolated_x(s_int, 1)
229 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
230 interpolated_x(s_int, 0),
232 pow(ext_el_pt->interpolated_x(s_ext, 1) -
233 interpolated_x(s_int, 1),
235 <<
" " << zeta[0] << std::endl;
242 FaceGeometry<ELASTICITY_BULK_ELEMENT>::output(file_pt);
246 void output(FILE* file_pt,
const unsigned& n_plot)
248 FaceGeometry<ELASTICITY_BULK_ELEMENT>::output(file_pt, n_plot);
261 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
262 double TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
263 ELASTICITY_BULK_ELEMENT,
264 HELMHOLTZ_BULK_ELEMENT>::Default_Q_Value = 1.0;
270 template<
class ELASTICITY_BULK_ELEMENT,
class HELMHOLTZ_BULK_ELEMENT>
271 void TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement<
272 ELASTICITY_BULK_ELEMENT,
273 HELMHOLTZ_BULK_ELEMENT>::
274 fill_in_contribution_to_residuals_helmholtz_traction(
275 Vector<double>& residuals)
278 unsigned n_node = nnode();
282 unsigned n_position_type = this->nnodal_position_type();
283 if (n_position_type != 1)
285 throw OomphLibError(
"LinearElasticity is not yet implemented for more "
286 "than one position type.",
287 OOMPH_CURRENT_FUNCTION,
288 OOMPH_EXCEPTION_LOCATION);
293 const unsigned n_dim = this->nodal_dimension();
296 std::vector<std::complex<unsigned>> u_nodal_index(n_dim);
297 for (
unsigned i = 0; i < n_dim; i++)
300 this->U_index_time_harmonic_linear_elasticity_helmholtz_traction[i];
308 DShape dpsids(n_node, n_dim - 1);
311 const double q_value = q();
314 Vector<std::complex<double>> traction(2);
317 unsigned n_intpt = integral_pt()->nweight();
320 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
323 double w = integral_pt()->weight(ipt);
326 dshape_local_at_knot(ipt, psi, dpsids);
329 Vector<double> interpolated_x(n_dim, 0.0);
332 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
335 for (
unsigned l = 0; l < n_node; l++)
338 for (
unsigned i = 0; i < n_dim; i++)
341 const double x_local = nodal_position(l, i);
342 interpolated_x[i] += x_local * psi(l);
345 for (
unsigned j = 0; j < n_dim - 1; j++)
347 interpolated_A(j, i) += x_local * dpsids(l, j);
353 DenseMatrix<double> A(n_dim - 1);
354 for (
unsigned i = 0; i < n_dim - 1; i++)
356 for (
unsigned j = 0; j < n_dim - 1; j++)
362 for (
unsigned k = 0; k < n_dim; k++)
364 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
370 Vector<double> interpolated_normal(n_dim);
371 outer_unit_normal(ipt, interpolated_normal);
381 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
385 "Wrong dimension in TimeHarmonicLinElastLoadedByPressureElement",
386 "TimeHarmonicLinElastLoadedByPressureElement::fill_in_contribution_"
388 OOMPH_EXCEPTION_LOCATION);
393 double W = w * sqrt(Adet);
396 HELMHOLTZ_BULK_ELEMENT* ext_el_pt =
397 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(external_element_pt(0, ipt));
398 Vector<double> s_ext(external_element_local_coord(0, ipt));
401 std::complex<double> u_helmholtz =
402 ext_el_pt->interpolated_u_helmholtz(s_ext);
403 traction[0] = -q_value * interpolated_normal[0] * u_helmholtz;
404 traction[1] = -q_value * interpolated_normal[1] * u_helmholtz;
407 for (
unsigned l = 0; l < n_node; l++)
410 for (
unsigned i = 0; i < n_dim; i++)
412 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].real());
417 residuals[local_eqn] -= traction[i].real() * psi(l) * W;
420 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].imag());
425 residuals[local_eqn] -= traction[i].imag() * psi(l) * W;
448 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
450 :
public virtual FaceGeometry<HELMHOLTZ_BULK_ELEMENT>,
451 public virtual FaceElement,
452 public virtual ElementWithExternalElement
458 FiniteElement*
const& bulk_el_pt,
const int& face_index);
480 residuals, GeneralisedElement::Dummy_matrix, 0);
487 DenseMatrix<double>& jacobian)
491 residuals, jacobian, 1);
494 fill_in_jacobian_from_external_interaction_by_fd(residuals, jacobian);
506 void output(std::ostream& outfile,
const unsigned& n_plot)
511 const unsigned n_intpt = integral_pt()->nweight();
514 Vector<double> s_int(
Dim - 1);
517 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
520 for (
unsigned i = 0; i < (
Dim - 1); i++)
522 s_int[i] = integral_pt()->knot(ipt, i);
526 Vector<double> interpolated_normal(
Dim);
527 outer_unit_normal(ipt, interpolated_normal);
529 Vector<double> zeta(1);
530 interpolated_zeta(s_int, zeta);
533 ELASTICITY_BULK_ELEMENT* ext_el_pt =
534 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(external_element_pt(0, ipt));
535 Vector<double> s_ext(external_element_local_coord(0, ipt));
536 Vector<std::complex<double>> displ(2);
537 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
542 std::complex<double> flux = (displ[0] * interpolated_normal[0] +
543 displ[1] * interpolated_normal[1]);
546 outfile << ext_el_pt->interpolated_x(s_ext, 0) <<
" "
547 << ext_el_pt->interpolated_x(s_ext, 1) <<
" "
548 << flux.real() * interpolated_normal[0] <<
" "
549 << flux.real() * interpolated_normal[1] <<
" "
550 << flux.imag() * interpolated_normal[0] <<
" "
551 << flux.imag() * interpolated_normal[1] <<
" "
552 << interpolated_normal[0] <<
" " << interpolated_normal[1]
553 <<
" " << flux.real() <<
" " << flux.imag() <<
" "
554 << interpolated_x(s_int, 0) <<
" " << interpolated_x(s_int, 1)
556 << sqrt(pow(ext_el_pt->interpolated_x(s_ext, 0) -
557 interpolated_x(s_int, 0),
559 pow(ext_el_pt->interpolated_x(s_ext, 1) -
560 interpolated_x(s_int, 1),
562 <<
" " << zeta[0] << std::endl;
570 FaceGeometry<HELMHOLTZ_BULK_ELEMENT>::output(file_pt);
574 void output(FILE* file_pt,
const unsigned& n_plot)
576 FaceGeometry<HELMHOLTZ_BULK_ELEMENT>::output(file_pt, n_plot);
589 unsigned n_node = nnode();
595 for (
unsigned i = 0; i < n_node; i++)
601 return J_eulerian(s);
613 unsigned n_node = nnode();
616 shape_at_knot(ipt, psi);
619 for (
unsigned i = 0; i < n_node; i++)
625 return J_eulerian_at_knot(ipt);
634 Vector<double>& residuals,
635 DenseMatrix<double>& jacobian,
636 const unsigned& flag);
655 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
657 ELASTICITY_BULK_ELEMENT>::
658 HelmholtzFluxFromNormalDisplacementBCElement(
659 FiniteElement*
const& bulk_el_pt,
const int& face_index)
660 : FaceGeometry<HELMHOLTZ_BULK_ELEMENT>(), FaceElement()
665 bulk_el_pt->build_face_element(face_index,
this);
670 HELMHOLTZ_BULK_ELEMENT* elem_pt =
671 dynamic_cast<HELMHOLTZ_BULK_ELEMENT*
>(bulk_el_pt);
673 if (elem_pt->dim() == 3)
676 RefineableElement* ref_el_pt =
677 dynamic_cast<RefineableElement*
>(elem_pt);
680 if (this->has_hanging_nodes())
682 throw OomphLibError(
"This flux element will not work correctly if "
683 "nodes are hanging\n",
684 OOMPH_CURRENT_FUNCTION,
685 OOMPH_EXCEPTION_LOCATION);
695 this->set_ninteraction(1);
699 Dim = this->node_pt(0)->ndim();
716 HelmholtzEquations<1>* eqn_pt =
717 dynamic_cast<HelmholtzEquations<1>*
>(bulk_el_pt);
721 std::string error_string =
722 "Bulk element must inherit from HelmholtzEquations.";
724 "Nodes are one dimensional, but cannot cast the bulk element to\n";
725 error_string +=
"HelmholtzEquations<1>\n.";
726 error_string +=
"If you desire this functionality, you must "
727 "implement it yourself\n";
730 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
744 HelmholtzEquations<2>* eqn_pt =
745 dynamic_cast<HelmholtzEquations<2>*
>(bulk_el_pt);
749 std::string error_string =
750 "Bulk element must inherit from HelmholtzEquations.";
752 "Nodes are two dimensional, but cannot cast the bulk element to\n";
753 error_string +=
"HelmholtzEquations<2>\n.";
754 error_string +=
"If you desire this functionality, you must "
755 "implement it yourself\n";
758 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
771 HelmholtzEquations<3>* eqn_pt =
772 dynamic_cast<HelmholtzEquations<3>*
>(bulk_el_pt);
776 std::string error_string =
777 "Bulk element must inherit from HelmholtzEquations.";
778 error_string +=
"Nodes are three dimensional, but cannot cast the "
780 error_string +=
"HelmholtzEquations<3>\n.";
781 error_string +=
"If you desire this functionality, you must "
782 "implement it yourself\n";
785 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
797 std::ostringstream error_stream;
798 error_stream <<
"Dimension of node is " <<
Dim
799 <<
". It should be 1,2, or 3!" << std::endl;
802 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
812 template<
class HELMHOLTZ_BULK_ELEMENT,
class ELASTICITY_BULK_ELEMENT>
814 ELASTICITY_BULK_ELEMENT>::
815 fill_in_generic_residual_contribution_helmholtz_flux_from_displacement(
816 Vector<double>& residuals,
817 DenseMatrix<double>& jacobian,
818 const unsigned& flag)
821 const unsigned n_node = nnode();
824 Shape psif(n_node), testf(n_node);
827 const unsigned n_intpt = integral_pt()->nweight();
830 Vector<double> s(Dim - 1);
836 const std::complex<unsigned> u_index_helmholtz =
837 U_index_helmholtz_from_displacement;
841 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
844 for (
unsigned i = 0; i < (Dim - 1); i++)
846 s[i] = integral_pt()->knot(ipt, i);
850 double w = integral_pt()->weight(ipt);
854 double J = shape_and_test(s, psif, testf);
860 Vector<double> interpolated_x(Dim, 0.0);
863 for (
unsigned l = 0; l < n_node; l++)
866 for (
unsigned i = 0; i < Dim; i++)
868 interpolated_x[i] += nodal_position(l, i) * psif[l];
873 Vector<double> interpolated_normal(2);
874 outer_unit_normal(ipt, interpolated_normal);
877 ELASTICITY_BULK_ELEMENT* ext_el_pt =
878 dynamic_cast<ELASTICITY_BULK_ELEMENT*
>(external_element_pt(0, ipt));
879 Vector<double> s_ext(external_element_local_coord(0, ipt));
880 Vector<std::complex<double>> displ(2);
881 ext_el_pt->interpolated_u_time_harmonic_linear_elasticity(s_ext, displ);
886 std::complex<double> flux =
887 (displ[0] * interpolated_normal[0] + displ[1] * interpolated_normal[1]);
892 for (
unsigned l = 0; l < n_node; l++)
894 local_eqn = nodal_local_eqn(l, u_index_helmholtz.real());
899 residuals[local_eqn] -= flux.real() * testf[l] * W;
906 local_eqn = nodal_local_eqn(l, u_index_helmholtz.imag());
911 residuals[local_eqn] -= flux.imag() * testf[l] * W;
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
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...
unsigned Dim
The spatial dimension of the problem.
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.
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 ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Broken assignment operator.
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 ...
HelmholtzFluxFromNormalDisplacementBCElement(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...
HelmholtzFluxFromNormalDisplacementBCElement(const HelmholtzFluxFromNormalDisplacementBCElement &dummy)=delete
Broken copy constructor.
void output(FILE *file_pt)
C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: flux etc at Gauss points; n_plot is ignored.
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
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.
void output(std::ostream &outfile)
Output function.
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.
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_helmholtz_traction
Index at which the i-th displacement component is stored.
void output(FILE *file_pt)
C_style output function.
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
double *& q_pt()
Return a pointer the ratio of stress scales used to non-dimensionalise the fluid and solid equations.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
double * Q_pt
Pointer to the ratio, , of the stress used to non-dimensionalise the fluid stresses to the stress us...
TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and elasticity equations....