pml_time_harmonic_linear_elasticity_elements.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for elements that solve the equations of linear
27// elasticity in cartesian coordinates
28
30
31
32namespace oomph
33{
34 /// Static default value for square of frequency
35 template<unsigned DIM>
36 double
38 1.0;
39
40
41 /// ///////////////////////////////////////////////////////////////
42 /// ///////////////////////////////////////////////////////////////
43 /// ///////////////////////////////////////////////////////////////
44
45 //======================================================================
46 /// Compute the strain tensor at local coordinate s
47 //======================================================================
48 template<unsigned DIM>
50 const Vector<double>& s, DenseMatrix<std::complex<double>>& strain) const
51 {
52#ifdef PARANOID
53 if ((strain.ncol() != DIM) || (strain.nrow() != DIM))
54 {
55 std::ostringstream error_message;
56 error_message << "Strain matrix is " << strain.ncol() << " x "
57 << strain.nrow() << ", but dimension of the equations is "
58 << DIM << std::endl;
59 throw OomphLibError(
60 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
61 }
62
63 // Find out how many position types there are
64 unsigned n_position_type = this->nnodal_position_type();
65
66 if (n_position_type != 1)
67 {
68 std::ostringstream error_message;
69 error_message << "PMLTimeHarmonicLinearElasticity is not yet "
70 << "implemented for more than one position type"
71 << std::endl;
72 throw OomphLibError(
73 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
74 }
75#endif
76
77
78 // Find out how many nodes there are in the element
79 unsigned n_node = nnode();
80
81 // Find the indices at which the local velocities are stored
82 std::complex<unsigned> u_nodal_index[DIM];
83 for (unsigned i = 0; i < DIM; i++)
84 {
85 u_nodal_index[i] = u_index_time_harmonic_linear_elasticity(i);
86 }
87
88 // Set up memory for the shape and derivative functions
89 Shape psi(n_node);
90 DShape dpsidx(n_node, DIM);
91
92 // Call the derivatives of the shape functions
93 (void)dshape_eulerian(s, psi, dpsidx);
94
95 // Calculate interpolated values of the derivative of global position
96 DenseMatrix<std::complex<double>> interpolated_dudx(
97 DIM, DIM, std::complex<double>(0.0, 0.0));
98
99 // Loop over nodes
100 for (unsigned l = 0; l < n_node; l++)
101 {
102 // Loop over velocity components
103 for (unsigned i = 0; i < DIM; i++)
104 {
105 // Get the nodal value
106 const std::complex<double> u_value =
107 std::complex<double>(this->nodal_value(l, u_nodal_index[i].real()),
108 this->nodal_value(l, u_nodal_index[i].imag()));
109
110 // Loop over derivative directions
111 for (unsigned j = 0; j < DIM; j++)
112 {
113 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
114 }
115 }
116 }
117
118 /// Now fill in the entries of the strain tensor
119 for (unsigned i = 0; i < DIM; i++)
120 {
121 // Do upper half of matrix
122 // Note that j must be signed here for the comparison test to work
123 // Also i must be cast to an int
124 for (int j = (DIM - 1); j >= static_cast<int>(i); j--)
125 {
126 // Off diagonal terms
127 if (static_cast<int>(i) != j)
128 {
129 strain(i, j) =
130 0.5 * (interpolated_dudx(i, j) + interpolated_dudx(j, i));
131 }
132 // Diagonal terms will including growth factor when it comes back in
133 else
134 {
135 strain(i, i) = interpolated_dudx(i, i);
136 }
137 }
138 // Matrix is symmetric so just copy lower half
139 for (int j = (i - 1); j >= 0; j--)
140 {
141 strain(i, j) = strain(j, i);
142 }
143 }
144 }
145
146
147 //======================================================================
148 /// Compute the Cauchy stress tensor at local coordinate s for
149 /// displacement formulation.
150 //======================================================================
151 template<unsigned DIM>
153 const Vector<double>& s, DenseMatrix<std::complex<double>>& stress) const
154 {
155#ifdef PARANOID
156 if ((stress.ncol() != DIM) || (stress.nrow() != DIM))
157 {
158 std::ostringstream error_message;
159 error_message << "Stress matrix is " << stress.ncol() << " x "
160 << stress.nrow() << ", but dimension of the equations is "
161 << DIM << std::endl;
162 throw OomphLibError(
163 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
164 }
165#endif
166
167 // Get strain
168 DenseMatrix<std::complex<double>> strain(DIM, DIM);
169 this->get_strain(s, strain);
170
171 // Now fill in the entries of the stress tensor without exploiting
172 // symmetry -- sorry too lazy. This fct is only used for
173 // postprocessing anyway...
174 for (unsigned i = 0; i < DIM; i++)
175 {
176 for (unsigned j = 0; j < DIM; j++)
177 {
178 stress(i, j) = 0.0;
179 for (unsigned k = 0; k < DIM; k++)
180 {
181 for (unsigned l = 0; l < DIM; l++)
182 {
183 stress(i, j) += this->E(i, j, k, l) * strain(k, l);
184 }
185 }
186 }
187 }
188 }
189
190
191 //=======================================================================
192 /// Compute the residuals for the linear elasticity equations in
193 /// cartesian coordinates. Flag indicates if we want Jacobian too.
194 //=======================================================================
195 template<unsigned DIM>
198 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
199 {
200 // Find out how many nodes there are
201 unsigned n_node = this->nnode();
202
203#ifdef PARANOID
204 // Find out how many positional dofs there are
205 unsigned n_position_type = this->nnodal_position_type();
206
207 if (n_position_type != 1)
208 {
209 std::ostringstream error_message;
210 error_message << "PMLTimeHarmonicLinearElasticity is not yet "
211 << "implemented for more than one position type"
212 << std::endl;
213 throw OomphLibError(
214 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
215 }
216
217 // Throw and error if an elasticity tensor has not been set
218 if (this->Elasticity_tensor_pt == 0)
219 {
220 throw OomphLibError("No elasticity tensor set.",
221 OOMPH_CURRENT_FUNCTION,
222 OOMPH_EXCEPTION_LOCATION);
223 }
224#endif
225
226 // Find the indices at which the local velocities are stored
227 std::complex<unsigned> u_nodal_index[DIM];
228 for (unsigned i = 0; i < DIM; i++)
229 {
230 u_nodal_index[i] = this->u_index_time_harmonic_linear_elasticity(i);
231 }
232
233
234 // Square of non-dimensional frequency
235 const double omega_sq_local = this->omega_sq();
236
237 // Set up memory for the shape functions
238 Shape psi(n_node);
239 DShape dpsidx(n_node, DIM);
240
241 // Set the value of Nintpt -- the number of integration points
242 unsigned n_intpt = this->integral_pt()->nweight();
243
244 // Set the vector to hold the local coordinates in the element
245 Vector<double> s(DIM);
246
247 // Integers to store the local equation numbers
248 int local_eqn_real = 0, local_unknown_real = 0;
249 int local_eqn_imag = 0, local_unknown_imag = 0;
250
251 // Loop over the integration points
252 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
253 {
254 // Assign the values of s
255 for (unsigned i = 0; i < DIM; ++i)
256 {
257 s[i] = this->integral_pt()->knot(ipt, i);
258 }
259
260 // Get the integral weight
261 double w = this->integral_pt()->weight(ipt);
262
263 // Call the derivatives of the shape functions (and get Jacobian)
264 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
265
266 // Storage for Eulerian coordinates (initialised to zero)
267 Vector<double> interpolated_x(DIM, 0.0);
268
269 // Displacement
270 Vector<std::complex<double>> interpolated_u(
271 DIM, std::complex<double>(0.0, 0.0));
272
273 // Calculate interpolated values of the derivative of global position
274 // wrt lagrangian coordinates
275 DenseMatrix<std::complex<double>> interpolated_dudx(
276 DIM, DIM, std::complex<double>(0.0, 0.0));
277
278 // Calculate displacements and derivatives
279 for (unsigned l = 0; l < n_node; l++)
280 {
281 // Loop over displacement components (deformed position)
282 for (unsigned i = 0; i < DIM; i++)
283 {
284 // Calculate the Lagrangian coordinates and the accelerations
285 interpolated_x[i] += this->raw_nodal_position(l, i) * psi(l);
286
287
288 // Get the nodal displacements
289 const std::complex<double> u_value = std::complex<double>(
290 this->raw_nodal_value(l, u_nodal_index[i].real()),
291 this->raw_nodal_value(l, u_nodal_index[i].imag()));
292
293 interpolated_u[i] += u_value * psi(l);
294
295 // Loop over derivative directions
296 for (unsigned j = 0; j < DIM; j++)
297 {
298 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
299 }
300 }
301 }
302
303 // Get body force
304 Vector<std::complex<double>> body_force_vec(DIM);
305 this->body_force(interpolated_x, body_force_vec);
306
307 // Premultiply the weights and the Jacobian
308 double W = w * J;
309
310
311 /// All the PML weights that participate in the assemby process
312 /// are computed here. pml_inverse_jacobian_diagonals are are used to
313 /// transform derivatives in real x to derivatives in transformed space
314 /// \f$\tilde x \f$.
315 /// pml_jacobian_det allows us to transform volume integrals in
316 /// transformed space to real space.
317 /// If the PML is not enabled via enable_pml, both default to 1.0.
318 Vector<std::complex<double>> pml_inverse_jacobian_diagonal(DIM);
319 std::complex<double> pml_jacobian_det;
320 this->compute_pml_coefficients(
321 ipt, interpolated_x, pml_inverse_jacobian_diagonal, pml_jacobian_det);
322
323 // Loop over the test functions, nodes of the element
324 for (unsigned l = 0; l < n_node; l++)
325 {
326 // Loop over the displacement components
327 for (unsigned a = 0; a < DIM; a++)
328 {
329 // Get real and imaginary equation numbers
330 local_eqn_real = this->nodal_local_eqn(l, u_nodal_index[a].real());
331 local_eqn_imag = this->nodal_local_eqn(l, u_nodal_index[a].imag());
332
333 //===== EQUATIONS OF PML TIME HARMONIC LINEAR ELASTICITY ========
334 // (This is where the maths happens)
335
336 // Calculate jacobian and residual contributions from acceleration and
337 // body force
338 std::complex<double> mass_jacobian_contribution =
339 -omega_sq_local * pml_jacobian_det;
340 std::complex<double> mass_residual_contribution =
341 (-omega_sq_local * interpolated_u[a] - body_force_vec[a]) *
342 pml_jacobian_det;
343
344 // Calculate jacobian and residual contributions from stress term
345 std::complex<double> stress_jacobian_contributions[DIM][DIM][DIM];
346 std::complex<double> stress_residual_contributions[DIM];
347 for (unsigned b = 0; b < DIM; b++)
348 {
349 stress_residual_contributions[b] = 0.0;
350 for (unsigned c = 0; c < DIM; c++)
351 {
352 for (unsigned d = 0; d < DIM; d++)
353 {
354 stress_jacobian_contributions[b][c][d] =
355 this->E(a, b, c, d) * pml_jacobian_det *
356 pml_inverse_jacobian_diagonal[b] *
357 pml_inverse_jacobian_diagonal[d];
358
359 stress_residual_contributions[b] +=
360 stress_jacobian_contributions[b][c][d] *
361 interpolated_dudx(c, d);
362 }
363 }
364 }
365
366 //===== ADD CONTRIBUTIONS TO GLOBAL RESIDUALS AND JACOBIAN ========
367
368 /*IF it's not a boundary condition*/
369 if (local_eqn_real >= 0)
370 {
371 // Acceleration and body force
372 residuals[local_eqn_real] +=
373 mass_residual_contribution.real() * psi(l) * W;
374
375 // Stress term
376 for (unsigned b = 0; b < DIM; b++)
377 {
378 // Add the stress terms to the residuals
379 residuals[local_eqn_real] +=
380 stress_residual_contributions[b].real() * dpsidx(l, b) * W;
381 }
382
383 // Jacobian entries
384 if (flag)
385 {
386 // Loop over the displacement basis functions again
387 for (unsigned l2 = 0; l2 < n_node; l2++)
388 {
389 // Loop over the displacement components again
390 for (unsigned c = 0; c < DIM; c++)
391 {
392 local_unknown_real =
393 this->nodal_local_eqn(l2, u_nodal_index[c].real());
394 local_unknown_imag =
395 this->nodal_local_eqn(l2, u_nodal_index[c].imag());
396 // If it's not pinned
397 if (local_unknown_real >= 0)
398 {
399 // Inertial term
400 if (a == c)
401 {
402 jacobian(local_eqn_real, local_unknown_real) +=
403 mass_jacobian_contribution.real() * psi(l) * psi(l2) *
404 W;
405 }
406
407 // Stress term
408 for (unsigned b = 0; b < DIM; b++)
409 {
410 for (unsigned d = 0; d < DIM; d++)
411 {
412 // Add the contribution to the Jacobian matrix
413 jacobian(local_eqn_real, local_unknown_real) +=
414 stress_jacobian_contributions[b][c][d].real() *
415 dpsidx(l2, d) * dpsidx(l, b) * W;
416 }
417 }
418 } // End of if not boundary condition
419
420 if (local_unknown_imag >= 0)
421 {
422 // Inertial term
423 if (a == c)
424 {
425 jacobian(local_eqn_real, local_unknown_imag) -=
426 mass_jacobian_contribution.imag() * psi(l) * psi(l2) *
427 W;
428 }
429
430 // Stress term
431 for (unsigned b = 0; b < DIM; b++)
432 {
433 for (unsigned d = 0; d < DIM; d++)
434 {
435 // Add the contribution to the Jacobian matrix
436 jacobian(local_eqn_real, local_unknown_imag) -=
437 stress_jacobian_contributions[b][c][d].imag() *
438 dpsidx(l2, d) * dpsidx(l, b) * W;
439 }
440 }
441 } // End of if not boundary condition
442 }
443 }
444 } // End of jacobian calculation
445
446 } // End of if not boundary condition for real eqn
447
448
449 /*IF it's not a boundary condition*/
450 if (local_eqn_imag >= 0)
451 {
452 // Acceleration and body force
453 residuals[local_eqn_imag] +=
454 mass_residual_contribution.imag() * psi(l) * W;
455
456 // Stress term
457 for (unsigned b = 0; b < DIM; b++)
458 {
459 // Add the stress terms to the residuals
460 residuals[local_eqn_imag] +=
461 stress_residual_contributions[b].imag() * dpsidx(l, b) * W;
462 }
463
464 // Jacobian entries
465 if (flag)
466 {
467 // Loop over the displacement basis functions again
468 for (unsigned l2 = 0; l2 < n_node; l2++)
469 {
470 // Loop over the displacement components again
471 for (unsigned c = 0; c < DIM; c++)
472 {
473 local_unknown_imag =
474 this->nodal_local_eqn(l2, u_nodal_index[c].imag());
475 local_unknown_real =
476 this->nodal_local_eqn(l2, u_nodal_index[c].real());
477 // If it's not pinned
478 if (local_unknown_imag >= 0)
479 {
480 // Inertial term
481 if (a == c)
482 {
483 jacobian(local_eqn_imag, local_unknown_imag) +=
484 mass_jacobian_contribution.real() * psi(l) * psi(l2) *
485 W;
486 }
487
488 // Stress term
489 for (unsigned b = 0; b < DIM; b++)
490 {
491 for (unsigned d = 0; d < DIM; d++)
492 {
493 // Add the contribution to the Jacobian matrix
494 jacobian(local_eqn_imag, local_unknown_imag) +=
495 stress_jacobian_contributions[b][c][d].real() *
496 dpsidx(l2, d) * dpsidx(l, b) * W;
497 }
498 }
499 } // End of if not boundary condition
500
501 if (local_unknown_real >= 0)
502 {
503 // Inertial term
504 if (a == c)
505 {
506 jacobian(local_eqn_imag, local_unknown_real) +=
507 mass_jacobian_contribution.imag() * psi(l) * psi(l2) *
508 W;
509 }
510
511 // Stress term
512 for (unsigned b = 0; b < DIM; b++)
513 {
514 for (unsigned d = 0; d < DIM; d++)
515 {
516 // Add the contribution to the Jacobian matrix
517 jacobian(local_eqn_imag, local_unknown_real) +=
518 stress_jacobian_contributions[b][c][d].imag() *
519 dpsidx(l2, d) * dpsidx(l, b) * W;
520 }
521 }
522 } // End of if not boundary condition
523 }
524 }
525 } // End of jacobian calculation
526
527 } // End of if not boundary condition for imag eqn
528
529 } // End of loop over coordinate directions
530 } // End of loop over shape functions
531 } // End of loop over integration points
532 }
533
534 //=======================================================================
535 /// Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
536 //=======================================================================
537 template<unsigned DIM>
539 std::ostream& outfile,
540 const unsigned& nplot,
542 {
543 // Vector of local coordinates
544 Vector<double> s(DIM);
545
546 // Vector for coordintes
547 Vector<double> x(DIM);
548
549 // Tecplot header info
550 outfile << this->tecplot_zone_string(nplot);
551
552 // Exact solution Vector
553 Vector<double> exact_soln(2 * DIM);
554
555 // Loop over plot points
556 unsigned num_plot_points = this->nplot_points(nplot);
557 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
558 {
559 // Get local coordinates of plot point
560 this->get_s_plot(iplot, nplot, s);
561
562 // Get x position as Vector
563 this->interpolated_x(s, x);
564
565 // Get exact solution at this point
566 (*exact_soln_pt)(x, exact_soln);
567
568 // Output x,y,...,u_exact,...
569 for (unsigned i = 0; i < DIM; i++)
570 {
571 outfile << x[i] << " ";
572 }
573 for (unsigned i = 0; i < 2 * DIM; i++)
574 {
575 outfile << exact_soln[i] << " ";
576 }
577 outfile << std::endl;
578 }
579
580 // Write tecplot footer (e.g. FE connectivity lists)
581 this->write_tecplot_zone_footer(outfile, nplot);
582 }
583
584
585 //=======================================================================
586 /// Output: x,y,[z],u,v,[w]
587 //=======================================================================
588 template<unsigned DIM>
590 std::ostream& outfile, const unsigned& nplot)
591 {
592 // Initialise local coord, global coord and solution vectors
593 Vector<double> s(DIM);
594 Vector<double> x(DIM);
596
597 // Tecplot header info
598 outfile << this->tecplot_zone_string(nplot);
599
600 // Loop over plot points
601 unsigned num_plot_points = this->nplot_points(nplot);
602 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
603 {
604 // Get local coordinates of plot point
605 this->get_s_plot(iplot, nplot, s);
606
607 // Get Eulerian coordinates and displacements
608 this->interpolated_x(s, x);
609 this->interpolated_u_time_harmonic_linear_elasticity(s, u);
610
611 // Output the x,y,..
612 for (unsigned i = 0; i < DIM; i++)
613 {
614 outfile << x[i] << " ";
615 }
616
617 // Output u,v,..
618 for (unsigned i = 0; i < DIM; i++)
619 {
620 outfile << u[i].real() << " ";
621 }
622
623 // Output u,v,..
624 for (unsigned i = 0; i < DIM; i++)
625 {
626 outfile << u[i].imag() << " ";
627 }
628
629 outfile << std::endl;
630 }
631
632 // Write tecplot footer (e.g. FE connectivity lists)
633 this->write_tecplot_zone_footer(outfile, nplot);
634 }
635
636
637 //======================================================================
638 /// Output function for real part of full time-dependent solution
639 /// constructed by adding the scattered field
640 ///
641 /// u = Re( (u_r +i u_i) exp(-i omega t)
642 ///
643 /// at phase angle omega t = phi computed here, to the corresponding
644 /// incoming wave specified via the function pointer.
645 ///
646 /// x,y,u or x,y,z,u
647 ///
648 /// Output at nplot points in each coordinate direction
649 //======================================================================
650 template<unsigned DIM>
652 std::ostream& outfile,
653 FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt,
654 const double& phi,
655 const unsigned& nplot)
656 {
657 // Initialise local coord, global coord and solution vectors
658 Vector<double> s(DIM);
659 Vector<double> x(DIM);
661
662 // Real and imag part of incoming wave
663 Vector<double> incoming_soln(2 * DIM);
664
665 // Tecplot header info
666 outfile << this->tecplot_zone_string(nplot);
667
668 // Loop over plot points
669 unsigned num_plot_points = this->nplot_points(nplot);
670 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
671 {
672 // Get local coordinates of plot point
673 this->get_s_plot(iplot, nplot, s);
674
675 // Get Eulerian coordinates and displacements
676 this->interpolated_x(s, x);
677 this->interpolated_u_time_harmonic_linear_elasticity(s, u);
678
679 // Get exact solution at this point
680 (*incoming_wave_fct_pt)(x, incoming_soln);
681
682 // Output the x,y,..
683 for (unsigned i = 0; i < DIM; i++)
684 {
685 outfile << x[i] << " ";
686 }
687
688 // Output u,v,..
689 for (unsigned i = 0; i < DIM; i++)
690 {
691 outfile << (u[i].real() + incoming_soln[2 * i]) * cos(phi) +
692 (u[i].imag() + incoming_soln[2 * i + 1]) * sin(phi)
693 << " ";
694 }
695
696 outfile << std::endl;
697 }
698
699 // Write tecplot footer (e.g. FE connectivity lists)
700 this->write_tecplot_zone_footer(outfile, nplot);
701 }
702
703 //======================================================================
704 /// Output function for imaginary part of full time-dependent solution
705 ///
706 /// u = Im( (u_r +i u_i) exp(-i omega t))
707 ///
708 /// at phase angle omega t = phi.
709 ///
710 /// x,y,u or x,y,z,u
711 ///
712 /// Output at nplot points in each coordinate direction
713 //======================================================================
714 template<unsigned DIM>
716 std::ostream& outfile, const double& phi, const unsigned& nplot)
717 {
718 // Initialise local coord, global coord and solution vectors
719 Vector<double> s(DIM);
720 Vector<double> x(DIM);
722
723 // Tecplot header info
724 outfile << this->tecplot_zone_string(nplot);
725
726 // Loop over plot points
727 unsigned num_plot_points = this->nplot_points(nplot);
728 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
729 {
730 // Get local coordinates of plot point
731 this->get_s_plot(iplot, nplot, s);
732
733 // Get Eulerian coordinates and displacements
734 this->interpolated_x(s, x);
735 this->interpolated_u_time_harmonic_linear_elasticity(s, u);
736
737 // Output the x,y,..
738 for (unsigned i = 0; i < DIM; i++)
739 {
740 outfile << x[i] << " ";
741 }
742
743 // Output u,v,..
744 for (unsigned i = 0; i < DIM; i++)
745 {
746 outfile << u[i].real() * cos(phi) + u[i].imag() * sin(phi) << " ";
747 }
748
749 outfile << std::endl;
750 }
751
752 // Write tecplot footer (e.g. FE connectivity lists)
753 this->write_tecplot_zone_footer(outfile, nplot);
754 }
755
756 //======================================================================
757 /// Output function for imaginary part of full time-dependent solution
758 ///
759 /// u = Im( (u_r +i u_i) exp(-i omega t))
760 ///
761 /// at phase angle omega t = phi.
762 ///
763 /// x,y,u or x,y,z,u
764 ///
765 /// Output at nplot points in each coordinate direction
766 //======================================================================
767 template<unsigned DIM>
769 std::ostream& outfile, const double& phi, const unsigned& nplot)
770 {
771 // Initialise local coord, global coord and solution vectors
772 Vector<double> s(DIM);
773 Vector<double> x(DIM);
775
776 // Tecplot header info
777 outfile << this->tecplot_zone_string(nplot);
778
779 // Loop over plot points
780 unsigned num_plot_points = this->nplot_points(nplot);
781 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
782 {
783 // Get local coordinates of plot point
784 this->get_s_plot(iplot, nplot, s);
785
786 // Get Eulerian coordinates and displacements
787 this->interpolated_x(s, x);
788 this->interpolated_u_time_harmonic_linear_elasticity(s, u);
789
790 // Output the x,y,..
791 for (unsigned i = 0; i < DIM; i++)
792 {
793 outfile << x[i] << " ";
794 }
795
796 // Output u,v,..
797 for (unsigned i = 0; i < DIM; i++)
798 {
799 outfile << u[i].imag() * cos(phi) - u[i].real() * sin(phi) << " ";
800 }
801
802 outfile << std::endl;
803 }
804
805 // Write tecplot footer (e.g. FE connectivity lists)
806 this->write_tecplot_zone_footer(outfile, nplot);
807 }
808
809
810 //=======================================================================
811 /// C-style output: x,y,[z],u,v,[w]
812 //=======================================================================
813 template<unsigned DIM>
815 FILE* file_pt, const unsigned& nplot)
816 {
817 // Vector of local coordinates
818 Vector<double> s(DIM);
819
820 // Tecplot header info
821 fprintf(file_pt, "%s", this->tecplot_zone_string(nplot).c_str());
822
823 // Loop over plot points
824 unsigned num_plot_points = this->nplot_points(nplot);
825 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
826 {
827 // Get local coordinates of plot point
828 this->get_s_plot(iplot, nplot, s);
829
830 // Coordinates
831 for (unsigned i = 0; i < DIM; i++)
832 {
833 fprintf(file_pt, "%g ", this->interpolated_x(s, i));
834 }
835
836 // Displacement
837 for (unsigned i = 0; i < DIM; i++)
838 {
839 fprintf(
840 file_pt,
841 "%g ",
842 this->interpolated_u_time_harmonic_linear_elasticity(s, i).real());
843 }
844 for (unsigned i = 0; i < DIM; i++)
845 {
846 fprintf(
847 file_pt,
848 "%g ",
849 this->interpolated_u_time_harmonic_linear_elasticity(s, i).imag());
850 }
851 }
852 fprintf(file_pt, "\n");
853
854 // Write tecplot footer (e.g. FE connectivity lists)
855 this->write_tecplot_zone_footer(file_pt, nplot);
856 }
857
858
859 //=======================================================================
860 /// Compute norm of the solution
861 //=======================================================================
862 template<unsigned DIM>
864 {
865 // Initialise
866 norm = 0.0;
867
868 // Vector of local coordinates
869 Vector<double> s(DIM);
870
871 // Vector for coordintes
872 Vector<double> x(DIM);
873
874 // Displacement vector
876
877 // Find out how many nodes there are in the element
878 unsigned n_node = this->nnode();
879
880 Shape psi(n_node);
881
882 // Set the value of n_intpt
883 unsigned n_intpt = this->integral_pt()->nweight();
884
885 // Loop over the integration points
886 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
887 {
888 // Assign values of s
889 for (unsigned i = 0; i < DIM; i++)
890 {
891 s[i] = this->integral_pt()->knot(ipt, i);
892 }
893
894 // Get the integral weight
895 double w = this->integral_pt()->weight(ipt);
896
897 // Get jacobian of mapping
898 double J = this->J_eulerian(s);
899
900 // Premultiply the weights and the Jacobian
901 double W = w * J;
902
903 // Get FE function value
904 this->interpolated_u_time_harmonic_linear_elasticity(s, disp);
905
906 // Add to norm
907 for (unsigned ii = 0; ii < DIM; ii++)
908 {
909 norm += (disp[ii].real() * disp[ii].real() +
910 disp[ii].imag() * disp[ii].imag()) *
911 W;
912 }
913 }
914 }
915
916 //======================================================================
917 /// Validate against exact solution
918 ///
919 /// Solution is provided via function pointer.
920 ///
921 //======================================================================
922 template<unsigned DIM>
924 std::ostream& outfile,
926 double& error,
927 double& norm)
928 {
929 // Initialise
930 error = 0.0;
931 norm = 0.0;
932
933 // Vector of local coordinates
934 Vector<double> s(DIM);
935
936 // Vector for coordintes
937 Vector<double> x(DIM);
938
939 // Displacement vector
941
942 // Find out how many nodes there are in the element
943 unsigned n_node = this->nnode();
944
945 Shape psi(n_node);
946
947 // Set the value of n_intpt
948 unsigned n_intpt = this->integral_pt()->nweight();
949
950 // Exact solution Vector
951 Vector<double> exact_soln(2 * DIM);
952
953 // Loop over the integration points
954 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
955 {
956 // Assign values of s
957 for (unsigned i = 0; i < DIM; i++)
958 {
959 s[i] = this->integral_pt()->knot(ipt, i);
960 }
961
962 // Get the integral weight
963 double w = this->integral_pt()->weight(ipt);
964
965 // Get jacobian of mapping
966 double J = this->J_eulerian(s);
967
968 // Premultiply the weights and the Jacobian
969 double W = w * J;
970
971 // Get x position as Vector
972 this->interpolated_x(s, x);
973
974 // Get exact solution at this point
975 (*exact_soln_pt)(x, exact_soln);
976
977 // Get FE function value
978 this->interpolated_u_time_harmonic_linear_elasticity(s, disp);
979
980 // Add to norm
981 for (unsigned ii = 0; ii < DIM; ii++)
982 {
983 // Add to error and norm
984 error += ((exact_soln[ii] - disp[ii].real()) *
985 (exact_soln[ii] - disp[ii].real()) +
986 (exact_soln[ii + DIM] - disp[ii].imag()) *
987 (exact_soln[ii + DIM] - disp[ii].imag())) *
988 W;
989 norm += (disp[ii].real() * disp[ii].real() +
990 disp[ii].imag() * disp[ii].imag()) *
991 W;
992 }
993 }
994 }
995
996
997 // Instantiate the required elements
1000
1004
1005 template<unsigned DIM>
1008
1009
1010} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A mapping function proposed by Bermudez et al, similar to the one above but is continuous across the ...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
An OomphLibError object which should be thrown when an run-time error is encountered....
A base class for elements that solve the equations of time-harmonic linear elasticity in Cartesian co...
static double Default_omega_sq_value
Static default value for square of frequency.
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain) const
Return the strain tensor.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void output_imag(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for imaginary part of full time-dependent solution u = Im( (u_r +i u_i) exp(-i omega ...
void get_stress(const Vector< double > &s, DenseMatrix< std::complex< double > > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
void output(std::ostream &outfile)
Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
virtual void fill_in_generic_contribution_to_residuals_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.
void output_total_real(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt, const double &phi, const unsigned &nplot)
Output function for real part of full time-dependent solution constructed by adding the scattered fie...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t)) a...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...