72 unsigned n_node = this->
nnode();
80 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
83 for (
unsigned i = 0;
i < 2;
i++)
102 for (
unsigned ii = 0; ii < 3; ii++)
104 norm += (disp[ii].real() * disp[ii].real() +
105 disp[ii].imag() * disp[ii].imag()) *
118 unsigned n_node = this->
nnode();
124 if (n_position_type != 1)
126 throw OomphLibError(
"TimeHarmonicFourierDecomposedLinearElasticity is "
127 "not yet implemented for more than one position type",
128 OOMPH_CURRENT_FUNCTION,
129 OOMPH_EXCEPTION_LOCATION);
134 std::complex<unsigned> u_nodal_index[3];
135 for (
unsigned i = 0;
i < 3;
i++)
158 std::complex<double>(0.0, 0.0));
163 3, 2, std::complex<double>(0.0, 0.0));
166 for (
unsigned l = 0; l < n_node; l++)
169 for (
unsigned i = 0;
i < 2;
i++)
174 for (
unsigned i = 0;
i < 3;
i++)
176 const std::complex<double> u_value =
177 std::complex<double>(this->
nodal_value(l, u_nodal_index[
i].real()),
180 interpolated_u[
i] += u_value * psi(l);
183 for (
unsigned j = 0; j < 2; j++)
185 interpolated_dudx(
i, j) += u_value * dpsidx(l, j);
192 const std::complex<double>
I(0.0, 1.0);
195 std::complex<double> ur = interpolated_u[0];
198 std::complex<double> uz = interpolated_u[1];
201 std::complex<double> uth = interpolated_u[2];
204 std::complex<double> durdr = interpolated_dudx(0, 0);
205 std::complex<double> durdz = interpolated_dudx(0, 1);
206 std::complex<double> duzdr = interpolated_dudx(1, 0);
207 std::complex<double> duzdz = interpolated_dudx(1, 1);
208 std::complex<double> duthdr = interpolated_dudx(2, 0);
209 std::complex<double> duthdz = interpolated_dudx(2, 1);
211 strain(0, 0) = durdr;
212 strain(2, 2) =
I * double(n) * uth / r + ur / r;
213 strain(1, 1) = duzdz;
214 strain(0, 2) = 0.5 * (
I * double(n) * ur / r - uth / r + duthdr);
215 strain(2, 0) = 0.5 * (
I * double(n) * ur / r - uth / r + duthdr);
216 strain(0, 1) = 0.5 * (durdz + duzdr);
217 strain(1, 0) = 0.5 * (durdz + duzdr);
218 strain(2, 1) = 0.5 * (duthdz +
I * double(n) * uz / r);
219 strain(1, 2) = 0.5 * (duthdz +
I * double(n) * uz / r);
236 unsigned n_node = this->
nnode();
242 if (n_position_type != 1)
244 throw OomphLibError(
"TimeHarmonicFourierDecomposedLinearElasticity is "
245 "not yet implemented for more than one position type",
246 OOMPH_CURRENT_FUNCTION,
247 OOMPH_EXCEPTION_LOCATION);
252 std::complex<unsigned> u_nodal_index[3];
253 for (
unsigned i = 0;
i < 3;
i++)
260 std::complex<double> nu_local = this->
nu();
261 std::complex<double> youngs_modulus_local = this->
youngs_modulus();
264 std::complex<double> lambda = youngs_modulus_local * nu_local /
265 (1.0 + nu_local) / (1.0 - 2.0 * nu_local);
266 std::complex<double> mu = youngs_modulus_local / 2.0 / (1.0 + nu_local);
288 int local_eqn_real = 0, local_eqn_imag = 0, local_unknown_real = 0,
289 local_unknown_imag = 0;
292 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
295 for (
unsigned i = 0;
i < 2; ++
i)
311 3, std::complex<double>(0.0, 0.0));
316 3, 2, std::complex<double>(0.0, 0.0));
319 for (
unsigned l = 0; l < n_node; l++)
322 for (
unsigned i = 0;
i < 2;
i++)
327 for (
unsigned i = 0;
i < 3;
i++)
329 const std::complex<double> u_value = std::complex<double>(
333 interpolated_u[
i] += u_value * psi(l);
336 for (
unsigned j = 0; j < 2; j++)
338 interpolated_dudx(
i, j) += u_value * dpsidx(l, j);
355 double rsq = pow(r, 2);
356 double nsq = pow(n, 2);
357 const std::complex<double>
I(0.0, 1.0);
360 std::complex<double> ur = interpolated_u[0];
363 std::complex<double> uz = interpolated_u[1];
366 std::complex<double> uth = interpolated_u[2];
369 std::complex<double> durdr = interpolated_dudx(0, 0);
370 std::complex<double> durdz = interpolated_dudx(0, 1);
371 std::complex<double> duzdr = interpolated_dudx(1, 0);
372 std::complex<double> duzdz = interpolated_dudx(1, 1);
373 std::complex<double> duthdr = interpolated_dudx(2, 0);
374 std::complex<double> duthdz = interpolated_dudx(2, 1);
377 std::complex<double> G_r, G_z, G_theta;
380 for (
unsigned l = 0; l < n_node; l++)
383 for (
unsigned a = 0; a < 3; a++)
389 if (local_eqn_real >= 0)
392 residuals[local_eqn_real] +=
393 (-
omega_sq.real() * interpolated_u[a].real() +
394 omega_sq.imag() * interpolated_u[a].imag() - b[a].real()) *
403 residuals[local_eqn_real] +=
405 (2.0 * durdr.real() * dpsidx(l, 0) +
407 (n * ur.real() / r + duthdr.imag() - uth.imag() / r) +
408 dpsidx(l, 1) * (durdz.real() + duzdr.real()) +
409 2.0 * psi(l) / pow(r, 2) * (ur.real() - n * uth.imag())) +
411 (-2.0 * durdr.imag() * dpsidx(l, 0) +
413 (-n * ur.imag() / r + duthdr.real() - uth.real() / r) -
414 dpsidx(l, 1) * (durdz.imag() + duzdr.imag()) -
415 2.0 * psi(l) / pow(r, 2) * (ur.imag() + n * uth.real())) +
417 (durdr.real() + ur.real() / r - n / r * uth.imag() +
419 (dpsidx(l, 0) + psi(l) / r) -
421 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
423 (dpsidx(l, 0) + psi(l) / r)) *
429 residuals[local_eqn_real] +=
431 (dpsidx(l, 0) * (durdz.real() + duzdr.real()) +
432 n * psi(l) / r * (n * uz.real() / r + duthdz.imag()) +
433 2.0 * duzdz.real() * dpsidx(l, 1)) +
435 (-dpsidx(l, 0) * (durdz.imag() + duzdr.imag()) +
436 n * psi(l) / r * (-n * uz.imag() / r + duthdz.real()) -
437 2.0 * duzdz.imag() * dpsidx(l, 1)) +
439 (durdr.real() + ur.real() / r - n / r * uth.imag() +
443 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
451 residuals[local_eqn_real] +=
453 ((duthdr.real() - uth.real() / r - n * ur.imag() / r) *
454 (dpsidx(l, 0) - psi(l) / r) +
455 2.0 * n * psi(l) / pow(r, 2) *
456 (n * uth.real() + ur.imag()) +
457 dpsidx(l, 1) * (duthdz.real() - n / r * uz.imag())) +
459 ((-duthdr.imag() + uth.imag() / r - n * ur.real() / r) *
460 (dpsidx(l, 0) - psi(l) / r) +
461 2.0 * n * psi(l) / pow(r, 2) *
462 (-n * uth.imag() + ur.real()) +
463 dpsidx(l, 1) * (-duthdz.imag() - n / r * uz.real())) +
465 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
469 (durdr.real() + ur.real() / r - n / r * uth.imag() +
478 OOMPH_CURRENT_FUNCTION,
479 OOMPH_EXCEPTION_LOCATION);
486 for (
unsigned l2 = 0; l2 < n_node; l2++)
494 G_r = (mu * (2.0 * dpsidx(l2, 0) * dpsidx(l, 0) +
495 (nsq + 2.0) / rsq * psi(l2) * psi(l) +
496 dpsidx(l2, 1) * dpsidx(l, 1)) +
497 lambda * (dpsidx(l2, 0) + psi(l2) / r) *
498 (dpsidx(l, 0) + psi(l) / r) -
502 G_z = (mu * dpsidx(l2, 0) * dpsidx(l, 1) +
503 lambda * dpsidx(l2, 1) * (dpsidx(l, 0) + psi(l) / r)) *
506 G_theta = (-
I * (mu * (n / r * dpsidx(l2, 0) * psi(l) -
507 3.0 * n / rsq * psi(l2) * psi(l)) -
508 lambda * n / r * psi(l2) *
509 (dpsidx(l, 0) + psi(l) / r))) *
516 (mu * dpsidx(l2, 1) * dpsidx(l, 0) +
517 lambda * (dpsidx(l2, 0) + psi(l2) / r) * dpsidx(l, 1)) *
520 G_z = (mu * (dpsidx(l2, 0) * dpsidx(l, 0) +
521 nsq / rsq * psi(l2) * psi(l) +
522 2.0 * dpsidx(l2, 1) * dpsidx(l, 1)) +
523 lambda * dpsidx(l2, 1) * dpsidx(l, 1) -
527 G_theta = (-
I * (mu * n / r * dpsidx(l2, 1) * psi(l) -
528 lambda * n / r * dpsidx(l, 1) * psi(l2))) *
535 G_r = (-
I * (mu * (-n / r * psi(l2) * dpsidx(l, 0) +
536 3.0 * n / rsq * psi(l2) * psi(l)) +
537 lambda * n / r * (dpsidx(l2, 0) + psi(l2) / r) *
541 G_z = (-
I * (-mu * n / r * psi(l2) * dpsidx(l, 1) +
542 lambda * n / r * dpsidx(l2, 1) * psi(l))) *
545 G_theta = (mu * ((dpsidx(l2, 0) - psi(l2) / r) *
546 (dpsidx(l, 0) - psi(l) / r) +
547 2.0 * nsq / rsq * psi(l2) * psi(l) +
548 dpsidx(l2, 1) * dpsidx(l, 1)) +
549 lambda * nsq / rsq * psi(l2) * psi(l) -
555 for (
unsigned c = 0; c < 3; c++)
564 if (local_unknown_real >= 0)
568 jacobian(local_eqn_real, local_unknown_real) +=
573 jacobian(local_eqn_real, local_unknown_real) +=
578 jacobian(local_eqn_real, local_unknown_real) +=
584 if (local_unknown_imag >= 0)
588 jacobian(local_eqn_real, local_unknown_imag) +=
593 jacobian(local_eqn_real, local_unknown_imag) +=
598 jacobian(local_eqn_real, local_unknown_imag) +=
613 if (local_eqn_imag >= 0)
616 residuals[local_eqn_imag] +=
617 (-
omega_sq.real() * interpolated_u[a].imag() -
618 omega_sq.imag() * interpolated_u[a].real() - b[a].imag()) *
627 residuals[local_eqn_imag] +=
629 (2 * durdr.imag() * dpsidx(l, 0) +
631 (n * ur.imag() / r - duthdr.real() + uth.real() / r) +
632 dpsidx(l, 1) * (durdz.imag() + duzdr.imag()) +
633 2 * psi(l) / pow(r, 2) * (ur.imag() + n * uth.real())) +
635 (2.0 * durdr.real() * dpsidx(l, 0) +
637 (n * ur.real() / r + duthdr.imag() - uth.imag() / r) +
638 dpsidx(l, 1) * (durdz.real() + duzdr.real()) +
639 2.0 * psi(l) / pow(r, 2) * (ur.real() - n * uth.imag())) +
641 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
643 (dpsidx(l, 0) + psi(l) / r) +
645 (durdr.real() + ur.real() / r - n / r * uth.imag() +
647 (dpsidx(l, 0) + psi(l) / r)) *
653 residuals[local_eqn_imag] +=
655 (dpsidx(l, 0) * (durdz.imag() + duzdr.imag()) +
656 n * psi(l) / r * (n * uz.imag() / r - duthdz.real()) +
657 2 * duzdz.imag() * dpsidx(l, 1)) +
659 (dpsidx(l, 0) * (durdz.real() + duzdr.real()) +
660 n * psi(l) / r * (n * uz.real() / r + duthdz.imag()) +
661 2.0 * duzdz.real() * dpsidx(l, 1)) +
663 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
667 (durdr.real() + ur.real() / r - n / r * uth.imag() +
675 residuals[local_eqn_imag] +=
677 ((duthdr.imag() - uth.imag() / r + n * ur.real() / r) *
678 (dpsidx(l, 0) - psi(l) / r) +
679 2.0 * n * psi(l) / pow(r, 2.) *
680 (n * uth.imag() - ur.real()) +
681 dpsidx(l, 1) * (duthdz.imag() + n / r * uz.real())) +
683 ((duthdr.real() - uth.real() / r - n * ur.imag() / r) *
684 (dpsidx(l, 0) - psi(l) / r) +
685 2.0 * n * psi(l) / pow(r, 2) *
686 (n * uth.real() + ur.imag()) +
687 dpsidx(l, 1) * (duthdz.real() - n / r * uz.imag())) +
689 (-durdr.real() - ur.real() / r + n / r * uth.imag() -
693 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
702 OOMPH_CURRENT_FUNCTION,
703 OOMPH_EXCEPTION_LOCATION);
710 for (
unsigned l2 = 0; l2 < n_node; l2++)
718 G_r = (mu * (2.0 * dpsidx(l2, 0) * dpsidx(l, 0) +
719 (nsq + 2.0) / rsq * psi(l2) * psi(l) +
720 dpsidx(l2, 1) * dpsidx(l, 1)) +
721 lambda * (dpsidx(l2, 0) + psi(l2) / r) *
722 (dpsidx(l, 0) + psi(l) / r) -
726 G_z = (mu * dpsidx(l2, 0) * dpsidx(l, 1) +
727 lambda * dpsidx(l2, 1) * (dpsidx(l, 0) + psi(l) / r)) *
730 G_theta = (-
I * (mu * (n / r * dpsidx(l2, 0) * psi(l) -
731 3.0 * n / rsq * psi(l2) * psi(l)) -
732 lambda * n / r * psi(l2) *
733 (dpsidx(l, 0) + psi(l) / r))) *
740 (mu * dpsidx(l2, 1) * dpsidx(l, 0) +
741 lambda * (dpsidx(l2, 0) + psi(l2) / r) * dpsidx(l, 1)) *
744 G_z = (mu * (dpsidx(l2, 0) * dpsidx(l, 0) +
745 nsq / rsq * psi(l2) * psi(l) +
746 2.0 * dpsidx(l2, 1) * dpsidx(l, 1)) +
747 lambda * dpsidx(l2, 1) * dpsidx(l, 1) -
751 G_theta = (-
I * (mu * n / r * dpsidx(l2, 1) * psi(l) -
752 lambda * n / r * dpsidx(l, 1) * psi(l2))) *
759 G_r = (-
I * (mu * (-n / r * psi(l2) * dpsidx(l, 0) +
760 3.0 * n / rsq * psi(l2) * psi(l)) +
761 lambda * n / r * (dpsidx(l2, 0) + psi(l2) / r) *
765 G_z = (-
I * (-mu * n / r * psi(l2) * dpsidx(l, 1) +
766 lambda * n / r * dpsidx(l2, 1) * psi(l))) *
769 G_theta = (mu * ((dpsidx(l2, 0) - psi(l2) / r) *
770 (dpsidx(l, 0) - psi(l) / r) +
771 2.0 * nsq / rsq * psi(l2) * psi(l) +
772 dpsidx(l2, 1) * dpsidx(l, 1)) +
773 lambda * nsq / rsq * psi(l2) * psi(l) -
779 for (
unsigned c = 0; c < 3; c++)
788 if (local_unknown_real >= 0)
792 jacobian(local_eqn_imag, local_unknown_real) +=
797 jacobian(local_eqn_imag, local_unknown_real) +=
802 jacobian(local_eqn_imag, local_unknown_real) +=
808 if (local_unknown_imag >= 0)
812 jacobian(local_eqn_imag, local_unknown_imag) +=
817 jacobian(local_eqn_imag, local_unknown_imag) +=
822 jacobian(local_eqn_imag, local_unknown_imag) +=
841 std::ostream& outfile,
842 const unsigned& nplot,
859 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
868 (*exact_soln_pt)(x, exact_soln);
871 for (
unsigned i = 0;
i < 2;
i++)
873 outfile << x[
i] <<
" ";
875 for (
unsigned i = 0;
i < 6;
i++)
877 outfile << exact_soln[
i] <<
" ";
879 outfile << std::endl;
890 std::ostream& outfile,
const unsigned& nplot)
902 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
913 for (
unsigned i = 0;
i < 2;
i++)
915 outfile << x[
i] <<
" ";
919 for (
unsigned i = 0;
i < 3;
i++)
921 outfile << u[
i].real() <<
" ";
925 for (
unsigned i = 0;
i < 3;
i++)
927 outfile << u[
i].imag() <<
" ";
930 outfile << std::endl;
942 FILE* file_pt,
const unsigned& nplot)
952 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
958 for (
unsigned i = 0;
i < 2;
i++)
964 for (
unsigned i = 0;
i < 3;
i++)
974 for (
unsigned i = 0;
i < 3;
i++)
985 fprintf(file_pt,
"\n");
998 std::ostream& outfile,
1015 outfile <<
"ZONE" << std::endl;
1021 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1024 for (
unsigned i = 0;
i < 2;
i++)
1042 (*exact_soln_pt)(x, exact_soln);
1045 for (
unsigned i = 0;
i < 3;
i++)
1047 norm += (exact_soln[
i] * exact_soln[
i] +
1048 exact_soln[
i + 3] * exact_soln[
i + 3]) *
1053 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1061 (exact_soln[
i + 3] -
1066 (exact_soln[
i + 3] -
1076 for (
unsigned i = 0;
i < 2;
i++)
1078 outfile << x[
i] <<
" ";
1082 for (
unsigned i = 0;
i < 3;
i++)
1087 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1092 for (
unsigned i = 0;
i < 3;
i++)
1095 << exact_soln[
i + 3] -
1097 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1102 outfile << std::endl;
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
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 std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
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 nnode() const
Return the number of nodes.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
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...
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
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.
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...
static std::complex< double > Default_youngs_modulus_value
Static default value for Young's modulus (1.0 – for natural scaling, i.e. all stresses have been non-...
static std::complex< double > Default_omega_sq_value
Static default value for squared frequency.
virtual std::complex< unsigned > u_index_time_harmonic_fourier_decomposed_linear_elasticity(const unsigned i) const
Return the index at which the i-th (i=0: r, i=1: z; i=2: theta) real or imag unknown displacement com...
int & fourier_wavenumber() const
Access function for Fourier wavenumber.
void interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(const Vector< double > &s, Vector< std::complex< double >> &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
std::complex< double > & nu() const
Access function for Poisson's ratio.
void body_force(const Vector< double > &x, Vector< std::complex< double >> &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
std::complex< double > youngs_modulus() const
Access function to Young's modulus.
const std::complex< double > & omega_sq() const
Access function for square of non-dim frequency.
void output(std::ostream &outfile)
Output: r,z, u_r_real, u_z_real, ..., u_theta_imag.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution: r,z, u_r_real, u_z_real, ..., u_theta_imag.
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double >> &strain)
Get strain tensor.
virtual void fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Private helper function to compute residuals and (if requested via flag) also the Jacobian matrix.
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...