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-2023 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 
32 namespace oomph
33 {
34  /// Static default value for square of angular frequency: Zero
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.
46 
47 
48  /// ///////////////////////////////////////////////////////////////
49  /// ///////////////////////////////////////////////////////////////
50  /// ///////////////////////////////////////////////////////////////
51 
52 
53  //=======================================================================
54  /// Compute norm of the solution
55  //=======================================================================
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
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
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  //=======================================================================
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] =
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)
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  //=======================================================================
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] =
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)
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
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  //=======================================================================
841  std::ostream& outfile,
842  const unsigned& nplot,
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  //=======================================================================
890  std::ostream& outfile, const unsigned& nplot)
891  {
892  // Set output Vector
893  Vector<double> s(2);
894  Vector<double> x(2);
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);
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  //=======================================================================
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
971  s, i)
972  .real());
973  }
974  for (unsigned i = 0; i < 3; i++)
975  {
976  fprintf(
977  file_pt,
978  "%g ",
979  this
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  //=======================================================================
998  std::ostream& outfile,
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
1059  s, i)
1060  .real()) +
1061  (exact_soln[i + 3] -
1062  this
1064  s, i)
1065  .imag()) *
1066  (exact_soln[i + 3] -
1067  this
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
1110 
1111 
1112 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
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
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:4103
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:2593
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:3161
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:3325
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
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:3962
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:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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:3148
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:3186
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:2317
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:3298
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:2576
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:1686
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:3174
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...
Definition: shape.h:76
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...
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.
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...
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 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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...