time_harmonic_fourier_decomposed_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-2026 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
29#include "time_harmonic_fourier_decomposed_linear_elasticity_elements.h"
30
31
32namespace oomph
33{
34 /// Static default value for square of angular frequency: Zero
35 std::complex<double> TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::
36 Default_omega_sq_value(0.0, 0.0);
37
38 /// Static default value for Young's modulus (1.0 -- for natural
39 /// scaling, i.e. all stresses have been non-dimensionalised by
40 /// the same reference Young's modulus. Setting the "non-dimensional"
41 /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
42 /// to a number larger than one means that the material is stiffer
43 /// than assumed in the non-dimensionalisation.
44 std::complex<double> TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::
45 Default_youngs_modulus_value(1.0, 0.0);
46
47
48 //////////////////////////////////////////////////////////////////
49 //////////////////////////////////////////////////////////////////
50 //////////////////////////////////////////////////////////////////
51
52
53 //=======================================================================
54 /// Compute norm of the solution
55 //=======================================================================
56 void TimeHarmonicFourierDecomposedLinearElasticityEquations::compute_norm(
57 double& norm)
58 {
59 // Initialise
60 norm = 0.0;
61
62 // Vector of local coordinates
63 Vector<double> s(2);
64
65 // Vector for coordintes
66 Vector<double> x(2);
67
68 // Displacement vector
69 Vector<std::complex<double>> disp(3);
70
71 // Find out how many nodes there are in the element
72 unsigned n_node = this->nnode();
73
74 Shape psi(n_node);
75
76 // Set the value of n_intpt
77 unsigned n_intpt = this->integral_pt()->nweight();
78
79 // Loop over the integration points
80 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
81 {
82 // Assign values of s
83 for (unsigned i = 0; i < 2; i++)
84 {
85 s[i] = this->integral_pt()->knot(ipt, i);
86 }
87
88 // Get the integral weight
89 double w = this->integral_pt()->weight(ipt);
90
91 // Get jacobian of mapping
92 double J = this->J_eulerian(s);
93
94 // Premultiply the weights and the Jacobian
95 double W = w * J;
96
97 // Get FE function value
98 this->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
99 s, disp);
100
101 // Add to norm
102 for (unsigned ii = 0; ii < 3; ii++)
103 {
104 norm += (disp[ii].real() * disp[ii].real() +
105 disp[ii].imag() * disp[ii].imag()) *
106 W;
107 }
108 }
109 }
110
111 //=======================================================================
112 /// Get strain tensor
113 //=======================================================================
114 void TimeHarmonicFourierDecomposedLinearElasticityEquations::get_strain(
115 const Vector<double>& s, DenseMatrix<std::complex<double>>& strain)
116 {
117 // Find out how many nodes there are
118 unsigned n_node = this->nnode();
119
120#ifdef PARANOID
121 // Find out how many positional dofs there are
122 unsigned n_position_type = this->nnodal_position_type();
123
124 if (n_position_type != 1)
125 {
126 throw OomphLibError("TimeHarmonicFourierDecomposedLinearElasticity is "
127 "not yet implemented for more than one position type",
128 OOMPH_CURRENT_FUNCTION,
129 OOMPH_EXCEPTION_LOCATION);
130 }
131#endif
132
133 // Find the indices at which the local displacements are stored
134 std::complex<unsigned> u_nodal_index[3];
135 for (unsigned i = 0; i < 3; i++)
136 {
137 u_nodal_index[i] =
138 this->u_index_time_harmonic_fourier_decomposed_linear_elasticity(i);
139 }
140
141 // Fourier wavenumber as double
142 double n = double(this->fourier_wavenumber());
143
144 // Set up memory for the shape functions
145 Shape psi(n_node);
146
147 // Derivs only w.r.t. r [0] and z [1]
148 DShape dpsidx(n_node, 2);
149
150 // Call the derivatives of the shape functions (ignore Jacobian)
151 this->dshape_eulerian(s, psi, dpsidx);
152
153 // Storage for Eulerian coordinates (r,z; initialised to zero)
154 Vector<double> interpolated_x(2, 0.0);
155
156 // Displacements u_r,u_z,u_theta
157 Vector<std::complex<double>> interpolated_u(3,
158 std::complex<double>(0.0, 0.0));
159
160 // Calculate interpolated values of the derivatives w.r.t.
161 // Eulerian coordinates
162 DenseMatrix<std::complex<double>> interpolated_dudx(
163 3, 2, std::complex<double>(0.0, 0.0));
164
165 // Calculate displacements and derivatives
166 for (unsigned l = 0; l < n_node; l++)
167 {
168 // Calculate the coordinates
169 for (unsigned i = 0; i < 2; i++)
170 {
171 interpolated_x[i] += this->nodal_position(l, i) * psi(l);
172 }
173 // Get the nodal displacements
174 for (unsigned i = 0; i < 3; i++)
175 {
176 const std::complex<double> u_value =
177 std::complex<double>(this->nodal_value(l, u_nodal_index[i].real()),
178 this->nodal_value(l, u_nodal_index[i].imag()));
179
180 interpolated_u[i] += u_value * psi(l);
181
182 // Loop over derivative directions
183 for (unsigned j = 0; j < 2; j++)
184 {
185 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
186 }
187 }
188 }
189
190 // define shorthand notation for regularly-occurring terms
191 double r = interpolated_x[0];
192 const std::complex<double> I(0.0, 1.0);
193
194 // r component of displacement
195 std::complex<double> ur = interpolated_u[0];
196
197 // z component of displacement
198 std::complex<double> uz = interpolated_u[1];
199
200 // theta component of displacement
201 std::complex<double> uth = interpolated_u[2];
202
203 // derivatives w.r.t. r and z:
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);
210
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);
220 }
221
222 ///////////////////////////////////////////////////////////////
223 ///////////////////////////////////////////////////////////////
224 ///////////////////////////////////////////////////////////////
225
226 //=======================================================================
227 /// Compute the residuals for the Fourier decomposed (in cyl. polars)
228 /// time harmonic linear elasticity equations in. Flag indicates if
229 /// we want Jacobian too.
230 //=======================================================================
231 void TimeHarmonicFourierDecomposedLinearElasticityEquations::
232 fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity(
233 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
234 {
235 // Find out how many nodes there are
236 unsigned n_node = this->nnode();
237
238#ifdef PARANOID
239 // Find out how many positional dofs there are
240 unsigned n_position_type = this->nnodal_position_type();
241
242 if (n_position_type != 1)
243 {
244 throw OomphLibError("TimeHarmonicFourierDecomposedLinearElasticity is "
245 "not yet implemented for more than one position type",
246 OOMPH_CURRENT_FUNCTION,
247 OOMPH_EXCEPTION_LOCATION);
248 }
249#endif
250
251 // Find the indices at which the local displacements are stored
252 std::complex<unsigned> u_nodal_index[3];
253 for (unsigned i = 0; i < 3; i++)
254 {
255 u_nodal_index[i] =
256 this->u_index_time_harmonic_fourier_decomposed_linear_elasticity(i);
257 }
258
259 // Get (complex) elastic parameters
260 std::complex<double> nu_local = this->nu();
261 std::complex<double> youngs_modulus_local = this->youngs_modulus();
262
263 // Obtain Lame parameters from Young's modulus and Poisson's ratio
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);
267
268
269 // Fourier wavenumber as double
270 double n = double(this->fourier_wavenumber());
271
272 // Square of non-dimensional frequency
273 const std::complex<double> omega_sq = this->omega_sq();
274
275 // Set up memory for the shape functions
276 Shape psi(n_node);
277
278 // Derivs only w.r.t. r [0] and z [1]
279 DShape dpsidx(n_node, 2);
280
281 // Set the value of Nintpt -- the number of integration points
282 unsigned n_intpt = this->integral_pt()->nweight();
283
284 // Set the vector to hold the local coordinates in the element
285 Vector<double> s(2);
286
287 // Integers to store the local equation numbers
288 int local_eqn_real = 0, local_eqn_imag = 0, local_unknown_real = 0,
289 local_unknown_imag = 0;
290
291 // Loop over the integration points
292 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
293 {
294 // Assign the values of s
295 for (unsigned i = 0; i < 2; ++i)
296 {
297 s[i] = this->integral_pt()->knot(ipt, i);
298 }
299
300 // Get the integral weight
301 double w = this->integral_pt()->weight(ipt);
302
303 // Call the derivatives of the shape functions (and get Jacobian)
304 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
305
306 // Storage for Eulerian coordinates (r,z; initialised to zero)
307 Vector<double> interpolated_x(2, 0.0);
308
309 // Displacements u_r,u_z,u_theta
310 Vector<std::complex<double>> interpolated_u(
311 3, std::complex<double>(0.0, 0.0));
312
313 // Calculate interpolated values of the derivatives w.r.t.
314 // Eulerian coordinates
315 DenseMatrix<std::complex<double>> interpolated_dudx(
316 3, 2, std::complex<double>(0.0, 0.0));
317
318 // Calculate displacements and derivatives
319 for (unsigned l = 0; l < n_node; l++)
320 {
321 // Calculate the coordinates
322 for (unsigned i = 0; i < 2; i++)
323 {
324 interpolated_x[i] += this->raw_nodal_position(l, i) * psi(l);
325 }
326 // Get the nodal displacements
327 for (unsigned i = 0; i < 3; i++)
328 {
329 const std::complex<double> u_value = std::complex<double>(
330 this->raw_nodal_value(l, u_nodal_index[i].real()),
331 this->raw_nodal_value(l, u_nodal_index[i].imag()));
332
333 interpolated_u[i] += u_value * psi(l);
334
335 // Loop over derivative directions
336 for (unsigned j = 0; j < 2; j++)
337 {
338 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
339 }
340 }
341 }
342
343 // Get body force
344 Vector<std::complex<double>> b(3);
345 this->body_force(interpolated_x, b);
346
347 // Premultiply the weights and the Jacobian
348 double W = w * J;
349
350 //=====EQUATIONS OF FOURIER DECOMPOSED TIME HARMONIC LINEAR ELASTICITY
351 //========
352
353 // define shorthand notation for regularly-occurring terms
354 double r = interpolated_x[0];
355 double rsq = pow(r, 2);
356 double nsq = pow(n, 2);
357 const std::complex<double> I(0.0, 1.0);
358
359 // r component of displacement
360 std::complex<double> ur = interpolated_u[0];
361
362 // z component of displacement
363 std::complex<double> uz = interpolated_u[1];
364
365 // theta component of displacement
366 std::complex<double> uth = interpolated_u[2];
367
368 // derivatives w.r.t. r and z:
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);
375
376 // storage for (complex) terms required for analytic Jacobian
377 std::complex<double> G_r, G_z, G_theta;
378
379 // Loop over the test functions, nodes of the element
380 for (unsigned l = 0; l < n_node; l++)
381 {
382 // Loop over the displacement components
383 for (unsigned a = 0; a < 3; a++)
384 {
385 // Get the REAL equation number
386 local_eqn_real = this->nodal_local_eqn(l, u_nodal_index[a].real());
387
388 /*IF it's not a boundary condition*/
389 if (local_eqn_real >= 0)
390 {
391 // Acceleration and body force
392 residuals[local_eqn_real] +=
393 (-omega_sq.real() * interpolated_u[a].real() +
394 omega_sq.imag() * interpolated_u[a].imag() - b[a].real()) *
395 psi(l) * r * W;
396
397 // Three components of the stress divergence term:
398 // a=0: r; a=1: z; a=2: theta
399
400 // Real part of r-equation
401 if (a == 0)
402 {
403 residuals[local_eqn_real] +=
404 (mu.real() *
405 (2.0 * durdr.real() * dpsidx(l, 0) +
406 n * psi(l) / r *
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())) +
410 mu.imag() *
411 (-2.0 * durdr.imag() * dpsidx(l, 0) +
412 n * psi(l) / r *
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())) +
416 lambda.real() *
417 (durdr.real() + ur.real() / r - n / r * uth.imag() +
418 duzdz.real()) *
419 (dpsidx(l, 0) + psi(l) / r) -
420 lambda.imag() *
421 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
422 duzdz.imag()) *
423 (dpsidx(l, 0) + psi(l) / r)) *
424 r * W;
425 }
426 // Real part of z-equation
427 else if (a == 1)
428 {
429 residuals[local_eqn_real] +=
430 (mu.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)) +
434 mu.imag() *
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)) +
438 lambda.real() *
439 (durdr.real() + ur.real() / r - n / r * uth.imag() +
440 duzdz.real()) *
441 dpsidx(l, 1) -
442 lambda.imag() *
443 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
444 duzdz.imag()) *
445 dpsidx(l, 1)) *
446 r * W;
447 }
448 // Real part of theta-equation
449 else if (a == 2)
450 {
451 residuals[local_eqn_real] +=
452 (mu.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())) +
458 mu.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())) +
464 lambda.real() *
465 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
466 duzdz.imag()) *
467 n * psi(l) / r +
468 lambda.imag() *
469 (durdr.real() + ur.real() / r - n / r * uth.imag() +
470 duzdz.real()) *
471 n * psi(l) / r) *
472 r * W;
473 }
474 // error: a should be 0, 1 or 2
475 else
476 {
477 throw OomphLibError("a should equal 0, 1 or 2",
478 OOMPH_CURRENT_FUNCTION,
479 OOMPH_EXCEPTION_LOCATION);
480 }
481
482 // Jacobian entries
483 if (flag)
484 {
485 // Loop over the displacement basis functions again
486 for (unsigned l2 = 0; l2 < n_node; l2++)
487 {
488 // define complex terms used to obtain entries of current row in
489 // the Jacobian:
490
491 // terms for rows of Jacobian matrix corresponding to r-equation
492 if (a == 0)
493 {
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) -
499 omega_sq * psi(l2) * psi(l)) *
500 r * W;
501
502 G_z = (mu * dpsidx(l2, 0) * dpsidx(l, 1) +
503 lambda * dpsidx(l2, 1) * (dpsidx(l, 0) + psi(l) / r)) *
504 r * W;
505
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))) *
510 r * W;
511 }
512 // terms for rows of Jacobian matrix corresponding to z-equation
513 else if (a == 1)
514 {
515 G_r =
516 (mu * dpsidx(l2, 1) * dpsidx(l, 0) +
517 lambda * (dpsidx(l2, 0) + psi(l2) / r) * dpsidx(l, 1)) *
518 r * W;
519
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) -
524 omega_sq * psi(l2) * psi(l)) *
525 r * W;
526
527 G_theta = (-I * (mu * n / r * dpsidx(l2, 1) * psi(l) -
528 lambda * n / r * dpsidx(l, 1) * psi(l2))) *
529 r * W;
530 }
531 // terms for rows of Jacobian matrix corresponding to
532 // theta-equation
533 else if (a == 2)
534 {
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) *
538 psi(l))) *
539 r * W;
540
541 G_z = (-I * (-mu * n / r * psi(l2) * dpsidx(l, 1) +
542 lambda * n / r * dpsidx(l2, 1) * psi(l))) *
543 r * W;
544
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) -
550 omega_sq * psi(l2) * psi(l)) *
551 r * W;
552 }
553
554 // Loop over the displacement components
555 for (unsigned c = 0; c < 3; c++)
556 {
557 // Get real and imaginary parts of local unknown
558 local_unknown_real =
559 this->nodal_local_eqn(l2, u_nodal_index[c].real());
560 local_unknown_imag =
561 this->nodal_local_eqn(l2, u_nodal_index[c].imag());
562
563 // If the real part of the local unknown is not pinned
564 if (local_unknown_real >= 0)
565 {
566 if (c == 0)
567 {
568 jacobian(local_eqn_real, local_unknown_real) +=
569 G_r.real();
570 }
571 else if (c == 1)
572 {
573 jacobian(local_eqn_real, local_unknown_real) +=
574 G_z.real();
575 }
576 else if (c == 2)
577 {
578 jacobian(local_eqn_real, local_unknown_real) +=
579 G_theta.real();
580 }
581 }
582
583 // If the imaginary part of the local unknown is not pinned
584 if (local_unknown_imag >= 0)
585 {
586 if (c == 0)
587 {
588 jacobian(local_eqn_real, local_unknown_imag) +=
589 -G_r.imag();
590 }
591 else if (c == 1)
592 {
593 jacobian(local_eqn_real, local_unknown_imag) +=
594 -G_z.imag();
595 }
596 else if (c == 2)
597 {
598 jacobian(local_eqn_real, local_unknown_imag) +=
599 -G_theta.imag();
600 }
601 }
602 }
603 }
604 } // End of jacobian calculation
605
606 } // End of if not boundary condition for real eqn
607
608
609 // Get the IMAG equation number
610 local_eqn_imag = this->nodal_local_eqn(l, u_nodal_index[a].imag());
611
612 /*IF it's not a boundary condition*/
613 if (local_eqn_imag >= 0)
614 {
615 // Acceleration and body force
616 residuals[local_eqn_imag] +=
617 (-omega_sq.real() * interpolated_u[a].imag() -
618 omega_sq.imag() * interpolated_u[a].real() - b[a].imag()) *
619 psi(l) * r * W;
620
621 // Three components of the stress divergence term:
622 // a=0: r; a=1: z; a=2: theta
623
624 // Imag part of r-equation
625 if (a == 0)
626 {
627 residuals[local_eqn_imag] +=
628 (mu.real() *
629 (2 * durdr.imag() * dpsidx(l, 0) +
630 n * psi(l) / r *
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())) +
634 mu.imag() *
635 (2.0 * durdr.real() * dpsidx(l, 0) +
636 n * psi(l) / r *
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())) +
640 lambda.real() *
641 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
642 duzdz.imag()) *
643 (dpsidx(l, 0) + psi(l) / r) +
644 lambda.imag() *
645 (durdr.real() + ur.real() / r - n / r * uth.imag() +
646 duzdz.real()) *
647 (dpsidx(l, 0) + psi(l) / r)) *
648 r * W;
649 }
650 // Imag part of z-equation
651 else if (a == 1)
652 {
653 residuals[local_eqn_imag] +=
654 (mu.real() *
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)) +
658 mu.imag() *
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)) +
662 lambda.real() *
663 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
664 duzdz.imag()) *
665 dpsidx(l, 1) +
666 lambda.imag() *
667 (durdr.real() + ur.real() / r - n / r * uth.imag() +
668 duzdz.real()) *
669 dpsidx(l, 1)) *
670 r * W;
671 }
672 // Imag part of theta-equation
673 else if (a == 2)
674 {
675 residuals[local_eqn_imag] +=
676 (mu.real() *
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())) +
682 mu.imag() *
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())) +
688 lambda.real() *
689 (-durdr.real() - ur.real() / r + n / r * uth.imag() -
690 duzdz.real()) *
691 n * psi(l) / r +
692 lambda.imag() *
693 (durdr.imag() + ur.imag() / r + n / r * uth.real() +
694 duzdz.imag()) *
695 n * psi(l) / r) *
696 r * W;
697 }
698 // error: a should be 0, 1 or 2
699 else
700 {
701 throw OomphLibError("a should equal 0, 1 or 2",
702 OOMPH_CURRENT_FUNCTION,
703 OOMPH_EXCEPTION_LOCATION);
704 }
705
706 // Jacobian entries
707 if (flag)
708 {
709 // Loop over the displacement basis functions again
710 for (unsigned l2 = 0; l2 < n_node; l2++)
711 {
712 // define complex terms used to obtain entries of current row in
713 // the Jacobian:
714
715 // terms for rows of Jacobian matrix corresponding to r-equation
716 if (a == 0)
717 {
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) -
723 omega_sq * psi(l2) * psi(l)) *
724 r * W;
725
726 G_z = (mu * dpsidx(l2, 0) * dpsidx(l, 1) +
727 lambda * dpsidx(l2, 1) * (dpsidx(l, 0) + psi(l) / r)) *
728 r * W;
729
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))) *
734 r * W;
735 }
736 // terms for rows of Jacobian matrix corresponding to z-equation
737 else if (a == 1)
738 {
739 G_r =
740 (mu * dpsidx(l2, 1) * dpsidx(l, 0) +
741 lambda * (dpsidx(l2, 0) + psi(l2) / r) * dpsidx(l, 1)) *
742 r * W;
743
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) -
748 omega_sq * psi(l2) * psi(l)) *
749 r * W;
750
751 G_theta = (-I * (mu * n / r * dpsidx(l2, 1) * psi(l) -
752 lambda * n / r * dpsidx(l, 1) * psi(l2))) *
753 r * W;
754 }
755 // terms for rows of Jacobian matrix corresponding to
756 // theta-equation
757 else if (a == 2)
758 {
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) *
762 psi(l))) *
763 r * W;
764
765 G_z = (-I * (-mu * n / r * psi(l2) * dpsidx(l, 1) +
766 lambda * n / r * dpsidx(l2, 1) * psi(l))) *
767 r * W;
768
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) -
774 omega_sq * psi(l2) * psi(l)) *
775 r * W;
776 }
777
778 // Loop over the displacement components
779 for (unsigned c = 0; c < 3; c++)
780 {
781 // Get real and imaginary parts of local unknown
782 local_unknown_real =
783 this->nodal_local_eqn(l2, u_nodal_index[c].real());
784 local_unknown_imag =
785 this->nodal_local_eqn(l2, u_nodal_index[c].imag());
786
787 // If the real part of the local unknown is not pinned
788 if (local_unknown_real >= 0)
789 {
790 if (c == 0)
791 {
792 jacobian(local_eqn_imag, local_unknown_real) +=
793 G_r.imag();
794 }
795 else if (c == 1)
796 {
797 jacobian(local_eqn_imag, local_unknown_real) +=
798 G_z.imag();
799 }
800 else if (c == 2)
801 {
802 jacobian(local_eqn_imag, local_unknown_real) +=
803 G_theta.imag();
804 }
805 }
806
807 // If the imaginary part of the local unknown is not pinned
808 if (local_unknown_imag >= 0)
809 {
810 if (c == 0)
811 {
812 jacobian(local_eqn_imag, local_unknown_imag) +=
813 G_r.real();
814 }
815 else if (c == 1)
816 {
817 jacobian(local_eqn_imag, local_unknown_imag) +=
818 G_z.real();
819 }
820 else if (c == 2)
821 {
822 jacobian(local_eqn_imag, local_unknown_imag) +=
823 G_theta.real();
824 }
825 }
826 }
827 }
828 } // End of jacobian calculation
829
830 } // End of if not boundary condition for imag eqn
831
832 } // End of loop over coordinate directions
833 } // End of loop over shape functions
834 } // End of loop over integration points
835 }
836
837 //=======================================================================
838 /// Output exact solution r,z, u_r_real, u_z_real, ..., u_theta_imag
839 //=======================================================================
840 void TimeHarmonicFourierDecomposedLinearElasticityEquations::output_fct(
841 std::ostream& outfile,
842 const unsigned& nplot,
843 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
844 {
845 // Vector of local coordinates
846 Vector<double> s(2);
847
848 // Vector for coordintes
849 Vector<double> x(2);
850
851 // Tecplot header info
852 outfile << this->tecplot_zone_string(nplot);
853
854 // Exact solution Vector
855 Vector<double> exact_soln(6);
856
857 // Loop over plot points
858 unsigned num_plot_points = this->nplot_points(nplot);
859 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
860 {
861 // Get local coordinates of plot point
862 this->get_s_plot(iplot, nplot, s);
863
864 // Get x position as Vector
865 this->interpolated_x(s, x);
866
867 // Get exact solution at this point
868 (*exact_soln_pt)(x, exact_soln);
869
870 // Output r,z,...,u_exact,...
871 for (unsigned i = 0; i < 2; i++)
872 {
873 outfile << x[i] << " ";
874 }
875 for (unsigned i = 0; i < 6; i++)
876 {
877 outfile << exact_soln[i] << " ";
878 }
879 outfile << std::endl;
880 }
881
882 // Write tecplot footer (e.g. FE connectivity lists)
883 this->write_tecplot_zone_footer(outfile, nplot);
884 }
885
886 //=======================================================================
887 /// Output: r,z, u_r_real, u_z_real, ..., u_theta_imag
888 //=======================================================================
889 void TimeHarmonicFourierDecomposedLinearElasticityEquations::output(
890 std::ostream& outfile, const unsigned& nplot)
891 {
892 // Set output Vector
893 Vector<double> s(2);
894 Vector<double> x(2);
895 Vector<std::complex<double>> u(3);
896
897 // Tecplot header info
898 outfile << this->tecplot_zone_string(nplot);
899
900 // Loop over plot points
901 unsigned num_plot_points = this->nplot_points(nplot);
902 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
903 {
904 // Get local coordinates of plot point
905 this->get_s_plot(iplot, nplot, s);
906
907 // Get Eulerian coordinates and displacements
908 this->interpolated_x(s, x);
909 this->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
910 s, u);
911
912 // Output the r,z,..
913 for (unsigned i = 0; i < 2; i++)
914 {
915 outfile << x[i] << " ";
916 }
917
918 // Output real part of displacement
919 for (unsigned i = 0; i < 3; i++)
920 {
921 outfile << u[i].real() << " ";
922 }
923
924 // Output imag part of displacement
925 for (unsigned i = 0; i < 3; i++)
926 {
927 outfile << u[i].imag() << " ";
928 }
929
930 outfile << std::endl;
931 }
932
933 // Write tecplot footer (e.g. FE connectivity lists)
934 this->write_tecplot_zone_footer(outfile, nplot);
935 }
936
937
938 //=======================================================================
939 /// C-style output:r,z, u_r_real, u_z_real, ..., u_theta_imag
940 //=======================================================================
941 void TimeHarmonicFourierDecomposedLinearElasticityEquations::output(
942 FILE* file_pt, const unsigned& nplot)
943 {
944 // Vector of local coordinates
945 Vector<double> s(2);
946
947 // Tecplot header info
948 fprintf(file_pt, "%s", this->tecplot_zone_string(nplot).c_str());
949
950 // Loop over plot points
951 unsigned num_plot_points = this->nplot_points(nplot);
952 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
953 {
954 // Get local coordinates of plot point
955 this->get_s_plot(iplot, nplot, s);
956
957 // Coordinates
958 for (unsigned i = 0; i < 2; i++)
959 {
960 fprintf(file_pt, "%g ", this->interpolated_x(s, i));
961 }
962
963 // Displacement
964 for (unsigned i = 0; i < 3; i++)
965 {
966 fprintf(
967 file_pt,
968 "%g ",
969 this
970 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
971 s, i)
972 .real());
973 }
974 for (unsigned i = 0; i < 3; i++)
975 {
976 fprintf(
977 file_pt,
978 "%g ",
979 this
980 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
981 s, i)
982 .imag());
983 }
984 }
985 fprintf(file_pt, "\n");
986
987 // Write tecplot footer (e.g. FE connectivity lists)
988 this->write_tecplot_zone_footer(file_pt, nplot);
989 }
990
991 //======================================================================
992 /// Validate against exact solution
993 /// Solution is provided via function pointer.
994 /// Plot at a given number of plot points and compute L2 error
995 /// and L2 norm of velocity solution over element.
996 //=======================================================================
997 void TimeHarmonicFourierDecomposedLinearElasticityEquations::compute_error(
998 std::ostream& outfile,
999 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
1000 double& error,
1001 double& norm)
1002 {
1003 error = 0.0;
1004 norm = 0.0;
1005
1006 // Vector of local coordinates
1007 Vector<double> s(2);
1008
1009 // Vector for coordinates
1010 Vector<double> x(2);
1011
1012 // Set the value of n_intpt
1013 unsigned n_intpt = this->integral_pt()->nweight();
1014
1015 outfile << "ZONE" << std::endl;
1016
1017 // Exact solution Vector (u_r_real, u_z_real, ..., u_theta_imag)
1018 Vector<double> exact_soln(6);
1019
1020 // Loop over the integration points
1021 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1022 {
1023 // Assign values of s
1024 for (unsigned i = 0; i < 2; i++)
1025 {
1026 s[i] = this->integral_pt()->knot(ipt, i);
1027 }
1028
1029 // Get the integral weight
1030 double w = this->integral_pt()->weight(ipt);
1031
1032 // Get jacobian of mapping
1033 double J = this->J_eulerian(s);
1034
1035 // Premultiply the weights and the Jacobian
1036 double W = w * J;
1037
1038 // Get x position as Vector
1039 this->interpolated_x(s, x);
1040
1041 // Get exact solution at this point
1042 (*exact_soln_pt)(x, exact_soln);
1043
1044 // Displacement error
1045 for (unsigned i = 0; i < 3; i++)
1046 {
1047 norm += (exact_soln[i] * exact_soln[i] +
1048 exact_soln[i + 3] * exact_soln[i + 3]) *
1049 W;
1050 error +=
1051 ((exact_soln[i] -
1052 this
1053 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1054 s, i)
1055 .real()) *
1056 (exact_soln[i] -
1057 this
1058 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1059 s, i)
1060 .real()) +
1061 (exact_soln[i + 3] -
1062 this
1063 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1064 s, i)
1065 .imag()) *
1066 (exact_soln[i + 3] -
1067 this
1068 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1069 s, i)
1070 .imag())) *
1071 W;
1072 }
1073
1074
1075 // Output r,z coordinates
1076 for (unsigned i = 0; i < 2; i++)
1077 {
1078 outfile << x[i] << " ";
1079 }
1080
1081 // Output ur_error, uz_error, uth_error
1082 for (unsigned i = 0; i < 3; i++)
1083 {
1084 outfile
1085 << exact_soln[i] -
1086 this
1087 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1088 s, i)
1089 .real()
1090 << " ";
1091 }
1092 for (unsigned i = 0; i < 3; i++)
1093 {
1094 outfile
1095 << exact_soln[i + 3] -
1096 this
1097 ->interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity(
1098 s, i)
1099 .imag()
1100 << " ";
1101 }
1102 outfile << std::endl;
1103 }
1104 }
1105
1106 // Instantiate required elements
1107 template class QTimeHarmonicFourierDecomposedLinearElasticityElement<2>;
1108 template class QTimeHarmonicFourierDecomposedLinearElasticityElement<3>;
1109 template class QTimeHarmonicFourierDecomposedLinearElasticityElement<4>;
1110
1111
1112} // namespace oomph
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition shape.h:278
oomph::DenseMatrix
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition matrices.h:386
oomph::FiniteElement::J_eulerian
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition elements.cc:4133
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition elements.h:1967
oomph::FiniteElement::nodal_value
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...
Definition elements.h:2597
oomph::FiniteElement::tecplot_zone_string
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition elements.h:3165
oomph::FiniteElement::dshape_eulerian_at_knot
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-...
Definition elements.cc:3355
oomph::FiniteElement::nnodal_position_type
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition elements.h:2467
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition elements.cc:3992
oomph::FiniteElement::nodal_local_eqn
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.
Definition elements.h:1436
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
oomph::FiniteElement::SteadyExactSolutionFctPt
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition elements.h:1763
oomph::FiniteElement::get_s_plot
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...
Definition elements.h:3152
oomph::FiniteElement::nplot_points
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition elements.h:3190
oomph::FiniteElement::nodal_position
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 ...
Definition elements.h:2321
oomph::FiniteElement::dshape_eulerian
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...
Definition elements.cc:3328
oomph::FiniteElement::raw_nodal_value
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.
Definition elements.h:2580
oomph::FiniteElement::raw_nodal_position
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....
Definition elements.cc:1714
oomph::FiniteElement::write_tecplot_zone_footer
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"...
Definition elements.h:3178
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition shape.h:76
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::interpolated_u_time_harmonic_fourier_decomposed_linear_elasticity
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.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:72
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::omega_sq
const std::complex< double > & omega_sq() const
Access function for square of non-dim frequency.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:159
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::Default_youngs_modulus_value
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-...
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:337
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::Default_omega_sq_value
static std::complex< double > Default_omega_sq_value
Static default value for squared frequency.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:329
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::body_force
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...
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:241
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquationsBase::u_index_time_harmonic_fourier_decomposed_linear_elasticity
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...
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:65
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquations::output
void output(std::ostream &outfile)
Output: r,z, u_r_real, u_z_real, ..., u_theta_imag.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.h:398
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquations::compute_error
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...
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.cc:997
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquations::output_fct
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.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.cc:840
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquations::fill_in_generic_contribution_to_residuals_fourier_decomp_time_harmonic_linear_elasticity
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.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.cc:232
oomph::TimeHarmonicFourierDecomposedLinearElasticityEquations::get_strain
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain)
Get strain tensor.
Definition time_harmonic_fourier_decomposed_linear_elasticity_elements.cc:114
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30
time_harmonic_fourier_decomposed_linear_elasticity_elements.h