axisym_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 Young's modulus (1.0 -- for natural
35 /// scaling, i.e. all stresses have been non-dimensionalised by
36 /// the same reference Young's modulus. Setting the "non-dimensional"
37 /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
38 /// to a number larger than one means that the material is stiffer
39 /// than assumed in the non-dimensionalisation.
40 double
42 1.0;
43
44 /// Static default value for timescale ratio (1.0 -- for natural scaling)
46 1.0;
47
48 /// ///////////////////////////////////////////////////////////////
49 /// ///////////////////////////////////////////////////////////////
50 /// ///////////////////////////////////////////////////////////////
51
52 //=======================================================================
53 /// Get strain (3x3 entries; r, z, phi)
54 //=======================================================================
56 const Vector<double>& s, DenseMatrix<double>& strain)
57 {
58 // Find out how many nodes there are
59 unsigned n_node = this->nnode();
60
61
62#ifdef PARANOID
63 // Find out how many positional dofs there are
64 unsigned n_position_type = this->nnodal_position_type();
65
66 if (n_position_type != 1)
67 {
68 throw OomphLibError("AxisymmetricLinearElasticity is not yet implemented "
69 "for more than one position type",
70 OOMPH_CURRENT_FUNCTION,
71 OOMPH_EXCEPTION_LOCATION);
72 }
73#endif
74
75 // Find the indices at which the local displacements are stored
76 unsigned u_nodal_index[3];
77 for (unsigned i = 0; i < 3; i++)
78 {
79 u_nodal_index[i] = this->u_index_axisymmetric_linear_elasticity(i);
80 }
81
82 // Set up memory for the shape functions
83 Shape psi(n_node);
84
85 // Derivs only w.r.t. r [0] and z [1]
86 DShape dpsidx(n_node, 2);
87
88 // Storage for Eulerian coordinates (r,z; initialised to zero)
90
91 // Displacements u_r,u_z,u_theta
92 Vector<double> interpolated_u(3, 0.0);
93
94 // Calculate interpolated values of the derivatives w.r.t.
95 // Eulerian coordinates
96 DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
97
98
99 // Calculate displacements and derivatives
100 for (unsigned l = 0; l < n_node; l++)
101 {
102 // Calculate the coordinates
103 for (unsigned i = 0; i < 2; i++)
104 {
105 interpolated_x[i] += this->raw_nodal_position(l, i) * psi(l);
106 }
107
108 // Get the nodal displacements
109 for (unsigned i = 0; i < 3; i++)
110 {
111 const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
112 interpolated_u[i] += u_value * psi(l);
113
114 // Loop over derivative directions
115 for (unsigned j = 0; j < 2; j++)
116 {
117 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
118 }
119 }
120 }
121
122
123 // define shorthand notation for regularly-occurring terms
124 double r = interpolated_x[0];
125
126 // r component of displacement
127 double ur = interpolated_u[0];
128
129 // theta component of displacement
130 double uth = interpolated_u[2];
131
132 // derivatives w.r.t. r and z:
133 double durdr = interpolated_dudx(0, 0);
134 double durdz = interpolated_dudx(0, 1);
135 double duzdr = interpolated_dudx(1, 0);
136 double duzdz = interpolated_dudx(1, 1);
137 double duthdr = interpolated_dudx(2, 0);
138 double duthdz = interpolated_dudx(2, 1);
139
140
141 // e_rr
142 strain(0, 0) = durdr;
143 // e_rz
144 strain(0, 1) = 0.5 * (durdz + duzdr);
145 strain(1, 0) = 0.5 * (durdz + duzdr);
146 // e_rphi
147 strain(0, 2) = 0.5 * (duthdr - uth / r);
148 strain(2, 0) = 0.5 * (duthdr - uth / r);
149 // e_zz
150 strain(1, 1) = duzdz;
151 // e_zphi
152 strain(1, 2) = 0.5 * duthdz;
153 strain(2, 1) = 0.5 * duthdz;
154 // e_phiphi
155 strain(2, 2) = ur / r;
156 }
157
158 //=======================================================================
159 /// Compute the residuals for the axisymmetric (in cyl. polars)
160 /// linear elasticity equations in. Flag indicates if we want Jacobian too.
161 //=======================================================================
164 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
165 {
166 // Find out how many nodes there are
167 unsigned n_node = this->nnode();
168
169 // Get continuous time from timestepper of first node
170 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
171
172#ifdef PARANOID
173 // Find out how many positional dofs there are
174 unsigned n_position_type = this->nnodal_position_type();
175
176 if (n_position_type != 1)
177 {
178 throw OomphLibError("AxisymmetricLinearElasticity is not yet implemented "
179 "for more than one position type",
180 OOMPH_CURRENT_FUNCTION,
181 OOMPH_EXCEPTION_LOCATION);
182 }
183#endif
184
185 // Find the indices at which the local displacements are stored
186 unsigned u_nodal_index[3];
187 for (unsigned i = 0; i < 3; i++)
188 {
189 u_nodal_index[i] = this->u_index_axisymmetric_linear_elasticity(i);
190 }
191
192 // Get elastic parameters
193 double nu_local = this->nu();
194 double youngs_modulus_local = this->youngs_modulus();
195
196 // Obtain Lame parameters from Young's modulus and Poisson's ratio
197 double lambda = youngs_modulus_local * nu_local / (1.0 + nu_local) /
198 (1.0 - 2.0 * nu_local);
199 double mu = youngs_modulus_local / 2.0 / (1.0 + nu_local);
200
201
202 // Lambda squared --- time scaling, NOT sqaure of Lame parameter lambda
203 const double lambda_sq = this->lambda_sq();
204
205 // Set up memory for the shape functions
206 Shape psi(n_node);
207
208 // Derivs only w.r.t. r [0] and z [1]
209 DShape dpsidx(n_node, 2);
210
211 // Set the value of Nintpt -- the number of integration points
212 unsigned n_intpt = this->integral_pt()->nweight();
213
214 // Set the vector to hold the local coordinates in the element
215 Vector<double> s(2);
216
217 // Integers to store the local equation numbers
218 int local_eqn = 0, local_unknown = 0;
219
220 // Loop over the integration points
221 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
222 {
223 // Assign the values of s
224 for (unsigned i = 0; i < 2; ++i)
225 {
226 s[i] = this->integral_pt()->knot(ipt, i);
227 }
228
229 // Get the integral weight
230 double w = this->integral_pt()->weight(ipt);
231
232 // Call the derivatives of the shape functions (and get Jacobian)
233 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
234
235 // Storage for Eulerian coordinates (r,z; initialised to zero)
237
238 // Displacements u_r,u_z,u_theta
239 Vector<double> interpolated_u(3, 0.0);
240
241 // Calculate interpolated values of the derivatives w.r.t.
242 // Eulerian coordinates
243 DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
244
245 Vector<double> d2u_dt2(3, 0.0);
246
247 // Calculate displacements and derivatives
248 for (unsigned l = 0; l < n_node; l++)
249 {
250 // Calculate the coordinates
251 for (unsigned i = 0; i < 2; i++)
252 {
253 interpolated_x[i] += this->raw_nodal_position(l, i) * psi(l);
254 }
255 // Get the nodal displacements
256 for (unsigned i = 0; i < 3; i++)
257 {
258 const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
259
260 interpolated_u[i] += u_value * psi(l);
261
262 d2u_dt2[i] += d2u_dt2_axisymmetric_linear_elasticity(l, i) * psi(l);
263
264 // Loop over derivative directions
265 for (unsigned j = 0; j < 2; j++)
266 {
267 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
268 }
269 }
270 }
271
272 // Get body force
273 Vector<double> b(3);
274 this->body_force(time, interpolated_x, b);
275
276 // Premultiply the weights and the Jacobian
277 double W = w * J;
278
279 //=====EQUATIONS OF AXISYMMETRIC LINEAR ELASTICITY ========
280
281 // define shorthand notation for regularly-occurring terms
282 double r = interpolated_x[0];
283 double rsq = pow(r, 2);
284
285 // r component of displacement
286 double ur = interpolated_u[0];
287
288 // theta component of displacement
289 double uth = interpolated_u[2];
290
291 // derivatives w.r.t. r and z:
292 double durdr = interpolated_dudx(0, 0);
293 double durdz = interpolated_dudx(0, 1);
294 double duzdr = interpolated_dudx(1, 0);
295 double duzdz = interpolated_dudx(1, 1);
296 double duthdr = interpolated_dudx(2, 0);
297 double duthdz = interpolated_dudx(2, 1);
298
299 // storage for terms required for analytic Jacobian
300 double G_r, G_z, G_theta;
301
302 // Loop over the test functions, nodes of the element
303 for (unsigned l = 0; l < n_node; l++)
304 {
305 // Loop over the displacement components
306 for (unsigned a = 0; a < 3; a++)
307 {
308 // Get the equation number
309 local_eqn = this->nodal_local_eqn(l, u_nodal_index[a]);
310
311 /*IF it's not a boundary condition*/
312 if (local_eqn >= 0)
313 {
314 // Acceleration and body force
315 residuals[local_eqn] +=
316 (lambda_sq * d2u_dt2[a] - b[a]) * psi(l) * r * W;
317
318 // Three components of the stress divergence term:
319 // a=0: r; a=1: z; a=2: theta
320
321 // r-equation
322 if (a == 0)
323 {
324 residuals[local_eqn] += (mu * (2.0 * durdr * dpsidx(l, 0) +
325 +dpsidx(l, 1) * (durdz + duzdr) +
326 2.0 * psi(l) / pow(r, 2) * (ur)) +
327 lambda * (durdr + ur / r + duzdz) *
328 (dpsidx(l, 0) + psi(l) / r)) *
329 r * W;
330 }
331 // z-equation
332 else if (a == 1)
333 {
334 residuals[local_eqn] +=
335 (mu * (dpsidx(l, 0) * (durdz + duzdr) +
336 2.0 * duzdz * dpsidx(l, 1)) +
337 lambda * (durdr + ur / r + duzdz) * dpsidx(l, 1)) *
338 r * W;
339 }
340 // theta-equation
341 else if (a == 2)
342 {
343 residuals[local_eqn] +=
344 (mu * ((duthdr - uth / r) * (dpsidx(l, 0) - psi(l) / r) +
345 dpsidx(l, 1) * (duthdz))) *
346 r * W;
347 }
348 // error: a should be 0, 1 or 2
349 else
350 {
351 throw OomphLibError("a should equal 0, 1 or 2",
352 OOMPH_CURRENT_FUNCTION,
353 OOMPH_EXCEPTION_LOCATION);
354 }
355
356 // Jacobian entries
357 if (flag)
358 {
359 // Loop over the displacement basis functions again
360 for (unsigned l2 = 0; l2 < n_node; l2++)
361 {
362 // define terms used to obtain entries of current row in the
363 // Jacobian:
364
365 // terms for rows of Jacobian matrix corresponding to r-equation
366 if (a == 0)
367 {
368 G_r =
369 (mu * (2.0 * dpsidx(l2, 0) * dpsidx(l, 0) +
370 2.0 / rsq * psi(l2) * psi(l) +
371 dpsidx(l2, 1) * dpsidx(l, 1)) +
372 lambda * (dpsidx(l2, 0) + psi(l2) / r) *
373 (dpsidx(l, 0) + psi(l) / r) +
374 lambda_sq * node_pt(l2)->time_stepper_pt()->weight(2, 0) *
375 psi(l2) * psi(l)) *
376 r * W;
377
378 G_z = (mu * dpsidx(l2, 0) * dpsidx(l, 1) +
379 lambda * dpsidx(l2, 1) * (dpsidx(l, 0) + psi(l) / r)) *
380 r * W;
381
382 G_theta = 0;
383 }
384
385 // terms for rows of Jacobian matrix corresponding to z-equation
386 else if (a == 1)
387 {
388 G_r =
389 (mu * dpsidx(l2, 1) * dpsidx(l, 0) +
390 lambda * (dpsidx(l2, 0) + psi(l2) / r) * dpsidx(l, 1)) *
391 r * W;
392
393 G_z =
394 (mu * (dpsidx(l2, 0) * dpsidx(l, 0) +
395 2.0 * dpsidx(l2, 1) * dpsidx(l, 1)) +
396 lambda * dpsidx(l2, 1) * dpsidx(l, 1) +
397 lambda_sq * node_pt(l2)->time_stepper_pt()->weight(2, 0) *
398 psi(l2) * psi(l)) *
399 r * W;
400
401 G_theta = 0;
402 }
403
404 // terms for rows of Jacobian matrix corresponding to
405 // theta-equation
406 else if (a == 2)
407 {
408 G_r = 0;
409
410 G_z = 0;
411
412 G_theta =
413 (mu * ((dpsidx(l2, 0) - psi(l2) / r) *
414 (dpsidx(l, 0) - psi(l) / r) +
415 dpsidx(l2, 1) * dpsidx(l, 1)) +
416 lambda_sq * node_pt(l2)->time_stepper_pt()->weight(2, 0) *
417 psi(l2) * psi(l)) *
418 r * W;
419 }
420
421 // Loop over the displacement components
422 for (unsigned c = 0; c < 3; c++)
423 {
424 // Get local unknown
425 local_unknown = this->nodal_local_eqn(l2, u_nodal_index[c]);
426
427 // If the local unknown is not pinned
428 if (local_unknown >= 0)
429 {
430 if (c == 0)
431 {
432 jacobian(local_eqn, local_unknown) += G_r;
433 }
434 else if (c == 1)
435 {
436 jacobian(local_eqn, local_unknown) += G_z;
437 }
438 else if (c == 2)
439 {
440 jacobian(local_eqn, local_unknown) += G_theta;
441 }
442 }
443 }
444 }
445 } // End of jacobian calculation
446
447 } // End of if not boundary condition
448 } // End of loop over coordinate directions
449 } // End of loop over shape functions
450 } // End of loop over integration points
451 }
452
453 //=======================================================================
454 /// Output exact solution r,z, u_r, u_z, u_theta
455 //=======================================================================
457 std::ostream& outfile,
458 const unsigned& nplot,
460 {
461 // Vector of local coordinates
462 Vector<double> s(2);
463
464 // Vector for coordintes
465 Vector<double> x(2);
466
467 // Tecplot header info
468 outfile << this->tecplot_zone_string(nplot);
469
470 // Exact solution Vector
471 Vector<double> exact_soln(9);
472
473 // Loop over plot points
474 unsigned num_plot_points = this->nplot_points(nplot);
475 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
476 {
477 // Get local coordinates of plot point
478 this->get_s_plot(iplot, nplot, s);
479
480 // Get x position as Vector
481 this->interpolated_x(s, x);
482
483 // Get exact solution at this point
484 (*exact_soln_pt)(x, exact_soln);
485
486 // Output r,z,...,u_exact,...
487 for (unsigned i = 0; i < 2; i++)
488 {
489 outfile << x[i] << " ";
490 }
491
492 for (unsigned i = 0; i < 3; i++)
493 {
494 outfile << exact_soln[i] << " ";
495 }
496 outfile << std::endl;
497 }
498
499 // Write tecplot footer (e.g. FE connectivity lists)
500 this->write_tecplot_zone_footer(outfile, nplot);
501 }
502
503 //=======================================================================
504 /// Output exact solution r,z, u_r, u_z, u_theta
505 /// Time dependent version
506 //=======================================================================
508 std::ostream& outfile,
509 const unsigned& nplot,
510 const double& time,
512 {
513 // Vector of local coordinates
514 Vector<double> s(2);
515
516 // Vector for coordintes
517 Vector<double> x(2);
518
519 // Tecplot header info
520 outfile << this->tecplot_zone_string(nplot);
521
522 // Exact solution Vector
523 Vector<double> exact_soln(9);
524
525 // Loop over plot points
526 unsigned num_plot_points = this->nplot_points(nplot);
527 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
528 {
529 // Get local coordinates of plot point
530 this->get_s_plot(iplot, nplot, s);
531
532 // Get x position as Vector
533 this->interpolated_x(s, x);
534
535 // Get exact solution at this point
536 (*exact_soln_pt)(time, x, exact_soln);
537
538 // Output r,z,...,u_exact,...
539 for (unsigned i = 0; i < 2; i++)
540 {
541 outfile << x[i] << " ";
542 }
543
544 for (unsigned i = 0; i < 9; i++)
545 {
546 outfile << exact_soln[i] << " ";
547 }
548 outfile << std::endl;
549 }
550
551 // Write tecplot footer (e.g. FE connectivity lists)
552 this->write_tecplot_zone_footer(outfile, nplot);
553 }
554
555 //=======================================================================
556 /// Output: r,z, u_r, u_z, u_theta
557 //=======================================================================
559 const unsigned& nplot)
560 {
561 // Set output Vector
562 Vector<double> s(2);
563 Vector<double> x(2);
564 Vector<double> u(3);
565 Vector<double> du_dt(3);
566 Vector<double> d2u_dt2(3);
567
568
569 // Tecplot header info
570 outfile << this->tecplot_zone_string(nplot);
571
572 // Loop over plot points
573 unsigned num_plot_points = this->nplot_points(nplot);
574 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
575 {
576 // Get local coordinates of plot point
577 this->get_s_plot(iplot, nplot, s);
578
579 // Get Eulerian coordinates and displacements
580 this->interpolated_x(s, x);
584
585 // Output the r,z,..
586 for (unsigned i = 0; i < 2; i++)
587 {
588 outfile << x[i] << " ";
589 }
590
591 // Output displacement
592 for (unsigned i = 0; i < 3; i++)
593 {
594 outfile << u[i] << " ";
595 }
596
597 // Output veloc
598 for (unsigned i = 0; i < 3; i++)
599 {
600 outfile << du_dt[i] << " ";
601 }
602
603 // Output accel
604 for (unsigned i = 0; i < 3; i++)
605 {
606 outfile << d2u_dt2[i] << " ";
607 }
608
609 outfile << std::endl;
610 }
611
612 // Write tecplot footer (e.g. FE connectivity lists)
613 this->write_tecplot_zone_footer(outfile, nplot);
614 }
615
616
617 //=======================================================================
618 /// C-style output:r,z, u_r, u_z, u_theta
619 //=======================================================================
621 const unsigned& nplot)
622 {
623 // Vector of local coordinates
624 Vector<double> s(2);
625
626 // Tecplot header info
627 fprintf(file_pt, "%s", this->tecplot_zone_string(nplot).c_str());
628
629 // Loop over plot points
630 unsigned num_plot_points = this->nplot_points(nplot);
631 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
632 {
633 // Get local coordinates of plot point
634 this->get_s_plot(iplot, nplot, s);
635
636 // Coordinates
637 for (unsigned i = 0; i < 2; i++)
638 {
639 fprintf(file_pt, "%g ", this->interpolated_x(s, i));
640 }
641
642 // Displacement
643 for (unsigned i = 0; i < 3; i++)
644 {
645 fprintf(file_pt,
646 "%g ",
648 }
649 }
650 fprintf(file_pt, "\n");
651
652 // Write tecplot footer (e.g. FE connectivity lists)
653 this->write_tecplot_zone_footer(file_pt, nplot);
654 }
655
656 //======================================================================
657 /// Validate against exact solution
658 /// Solution is provided via function pointer.
659 /// Plot at a given number of plot points and compute L2 error
660 /// and L2 norm of velocity solution over element.
661 //=======================================================================
663 std::ostream& outfile,
665 double& error,
666 double& norm)
667 {
668 error = 0.0;
669 norm = 0.0;
670
671 // Vector of local coordinates
672 Vector<double> s(2);
673
674 // Vector for coordinates
675 Vector<double> x(2);
676
677 // Set the value of n_intpt
678 unsigned n_intpt = this->integral_pt()->nweight();
679
680 outfile << "ZONE" << std::endl;
681
682 // Exact solution Vector (u_r, u_z, u_theta)
683 Vector<double> exact_soln(9);
684
685 // Loop over the integration points
686 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
687 {
688 // Assign values of s
689 for (unsigned i = 0; i < 2; i++)
690 {
691 s[i] = this->integral_pt()->knot(ipt, i);
692 }
693
694 // Get the integral weight
695 double w = this->integral_pt()->weight(ipt);
696
697 // Get jacobian of mapping
698 double J = this->J_eulerian(s);
699
700 // Premultiply the weights and the Jacobian
701 double W = w * J;
702
703 // Get x position as Vector
704 this->interpolated_x(s, x);
705
706 // Get exact solution at this point
707 (*exact_soln_pt)(x, exact_soln);
708
709 // Displacement error
710 for (unsigned i = 0; i < 3; i++)
711 {
712 norm += (exact_soln[i] * exact_soln[i]) * W;
713 error += ((exact_soln[i] -
715 (exact_soln[i] -
717 W;
718 }
719
720
721 // Output r,z coordinates
722 for (unsigned i = 0; i < 2; i++)
723 {
724 outfile << x[i] << " ";
725 }
726
727 // Output ur_error, uz_error, uth_error
728 for (unsigned i = 0; i < 3; i++)
729 {
730 outfile << exact_soln[i] -
732 << " ";
733 }
734 outfile << std::endl;
735 }
736 }
737
738 //======================================================================
739 /// Validate against exact solution
740 /// Solution is provided via function pointer.
741 /// Plot at a given number of plot points and compute L2 error
742 /// and L2 norm of velocity solution over element.
743 //=======================================================================
745 std::ostream& outfile,
747 const double& time,
748 double& error,
749 double& norm)
750 {
751 error = 0.0;
752 norm = 0.0;
753
754 // Vector of local coordinates
755 Vector<double> s(2);
756
757 // Vector for coordinates
758 Vector<double> x(2);
759
760 // Set the value of n_intpt
761 unsigned n_intpt = this->integral_pt()->nweight();
762
763 outfile << "ZONE" << std::endl;
764
765 // Exact solution Vector (u_r, u_z, u_theta)
766 Vector<double> exact_soln(9);
767
768 // Loop over the integration points
769 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
770 {
771 // Assign values of s
772 for (unsigned i = 0; i < 2; i++)
773 {
774 s[i] = this->integral_pt()->knot(ipt, i);
775 }
776
777 // Get the integral weight
778 double w = this->integral_pt()->weight(ipt);
779
780 // Get jacobian of mapping
781 double J = this->J_eulerian(s);
782
783 // Premultiply the weights and the Jacobian
784 double W = w * J;
785
786 // Get x position as Vector
787 this->interpolated_x(s, x);
788
789 // Get exact solution at this point
790 (*exact_soln_pt)(time, x, exact_soln);
791
792 // Displacement error
793 for (unsigned i = 0; i < 3; i++)
794 {
795 norm += (exact_soln[i] * exact_soln[i]) * W;
796 error += ((exact_soln[i] -
798 (exact_soln[i] -
800 W;
801 }
802
803
804 // Output r,z coordinates
805 for (unsigned i = 0; i < 2; i++)
806 {
807 outfile << x[i] << " ";
808 }
809
810 // Output ur_error, uz_error, uth_error
811 for (unsigned i = 0; i < 3; i++)
812 {
813 outfile << exact_soln[i] -
815 << " ";
816 }
817 outfile << std::endl;
818 }
819 }
820
821 // Instantiate required elements
825
826
827} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
double d2u_dt2_axisymmetric_linear_elasticity(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
void interpolated_d2u_dt2_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &d2u_dt2) const
Compute vector of FE interpolated accel d2u/dt2 at local coordinate s.
virtual unsigned u_index_axisymmetric_linear_elasticity(const unsigned &i) const
Return the index at which the i-th (i=0: r, i=1: z; i=2: theta) unknown displacement component is sto...
double youngs_modulus() const
Access function to Young's modulus.
double & nu() const
Access function for Poisson's ratio.
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 for natural scaling)
void interpolated_du_dt_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &du_dt) const
Compute vector of FE interpolated velocity du/dt at local coordinate s.
static double Default_youngs_modulus_value
Static default value for Young's modulus (1.0 – for natural scaling, i.e. all stresses have been non-...
void body_force(const double &time, const Vector< double > &x, Vector< double > &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
void interpolated_u_axisymmetric_linear_elasticity(const Vector< double > &s, Vector< double > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain)
Get strain (3x3 entries; r, z, phi)
virtual void fill_in_generic_contribution_to_residuals_axisymmetric_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(std::ostream &outfile)
Output: r,z, u_r, u_z, u_theta.
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, u_z, u_theta.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
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
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
//////////////////////////////////////////////////////////////////// ////////////////////////////////...