axisym_solid_elements.h
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 // Header file for axisymmetric solid mechanics elements
27 #ifndef OOMPH_AXISYMM_ELASTICITY_ELEMENTS_HEADER
28 #define OOMPH_AXISYMM_ELASTICITY_ELEMENTS_HEADER
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 // OOMPH-LIB headers
36 #include "../generic/Qelements.h"
37 #include "../generic/hermite_elements.h"
38 #include "../constitutive/constitutive_laws.h"
39 
40 namespace oomph
41 {
42  //=====================================================================
43  /// A class for elements that solve the equations of solid mechanics,
44  /// based on the principle of virtual displacements in
45  /// an axisymmetric formulation. In this case x[0] is the component of
46  /// displacement in the radial direction and x[1] is that in the theta
47  /// direction.
48  //=====================================================================
50  {
51  private:
52  /// Pointer to constitutive law
54 
55  public:
56  /// Constructor
58 
59  /// Return the constitutive law pointer
61  {
62  return Constitutive_law_pt;
63  }
64 
65  /// Return the stress tensor, as calculated from the constitutive law
67  const DenseMatrix<double>& G,
68  DenseMatrix<double>& sigma)
69  {
70 #ifdef PARANOID
71  // If the pointer to the constitutive law hasn't been set, issue an error
72  if (Constitutive_law_pt == 0)
73  {
74  std::string error_message =
75  "Elements derived from AxisymmetricPVDEquations";
76  error_message += " must have a constitutive law :\n ";
77  error_message +=
78  "set one using the constitutive_law_pt() member function\n";
79 
80  throw OomphLibError(
81  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
82  }
83 #endif
85  }
86 
87  /// Fill in the residuals by calling the generic function
89  {
91  }
92 
93  /// Return the residuals for the equations of solid mechanics
95  Vector<double>& residuals)
96  {
97  // Set the number of Lagrangian coordinates
98  unsigned n_lagrangian = 2;
99  // Find out how many nodes there are
100  unsigned n_node = nnode();
101  // Find out how many positional dofs there are
102  unsigned n_position_type = nnodal_position_type();
103 
104  // Integer to store local equation number
105  int local_eqn = 0;
106 
107  // Set up memory for the shape functions
108  Shape psi(n_node, n_position_type);
109  DShape dpsidxi(n_node, n_position_type, n_lagrangian);
110 
111  // Set the value of n_intpt
112  unsigned n_intpt = integral_pt()->nweight();
113 
114  // Loop over the integration points
115  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
116  {
117  // Get the integral weight
118  double w = integral_pt()->weight(ipt);
119  // Call the derivatives of the shape functions
120  double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
121  // Premultiply the weights and the Jacobian
122  double W = w * J;
123 
124  // Calculate the local Lagrangian coordinates, position components
125  // and the derivatives of global position components
126  // wrt lagrangian coordinates
127  double interpolated_xi[2] = {0.0, 0.0};
128  double interpolated_X[2] = {0.0, 0.0};
129  double interpolated_dXdxi[2][2];
130 
131  // Initialise interpolated_dXdxi to zero
132  for (unsigned i = 0; i < 2; i++)
133  {
134  for (unsigned j = 0; j < 2; j++)
135  {
136  interpolated_dXdxi[i][j] = 0.0;
137  }
138  }
139 
140  // Calculate displacements and derivatives
141  for (unsigned l = 0; l < n_node; l++)
142  {
143  // Loop over positional dofs
144  for (unsigned k = 0; k < n_position_type; k++)
145  {
146  // Loop over displacement components (deformed position)
147  for (unsigned i = 0; i < 2; i++)
148  {
149  // Set the value of the lagrangian coordinate
150  interpolated_xi[i] +=
151  lagrangian_position_gen(l, k, i) * psi(l, k);
152  // Set the value of the position component
153  interpolated_X[i] += nodal_position_gen(l, k, i) * psi(l, k);
154  // Loop over Lagrangian derivative directions
155  for (unsigned j = 0; j < 2; j++)
156  {
157  // Calculate dX[i]/dxi_{j}
158  interpolated_dXdxi[i][j] +=
159  nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
160  }
161  }
162  }
163  }
164 
165  // We are now in a position to calculate the undeformed metric tensor
166  DenseMatrix<double> g(3);
167  // r row
168  g(0, 0) = 1.0;
169  g(0, 1) = 0.0;
170  g(0, 2) = 0.0;
171  // theta row
172  g(1, 0) = 0.0;
173  g(1, 1) = interpolated_xi[0] * interpolated_xi[0];
174  g(1, 2) = 0.0;
175  // phi row
176  g(2, 0) = 0.0;
177  g(2, 1) = 0.0;
178  g(2, 2) = interpolated_xi[0] * interpolated_xi[0] *
179  sin(interpolated_xi[1]) * sin(interpolated_xi[1]);
180 
181  // Now multiply the weight by the square-root of the undeformed metric
182  // tensor r^2 sin(theta)
183  W *= sqrt(g(0, 0) * g(1, 1) * g(2, 2));
184 
185  // Now calculate the deformed metric tensor
186  DenseMatrix<double> G(3);
187  // r row
188  G(0, 0) = interpolated_dXdxi[0][0] * interpolated_dXdxi[0][0] +
189  interpolated_dXdxi[1][0] * interpolated_dXdxi[1][0];
190  G(0, 1) = interpolated_dXdxi[0][0] *
191  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
192  interpolated_dXdxi[1][0] *
193  (interpolated_dXdxi[1][1] + interpolated_X[0]);
194  G(0, 2) = 0.0;
195  // theta row
196  G(1, 0) = G(0, 1);
197  G(1, 1) = (interpolated_dXdxi[0][1] - interpolated_X[1]) *
198  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
199  (interpolated_dXdxi[1][1] + interpolated_X[0]) *
200  (interpolated_dXdxi[1][1] + interpolated_X[0]);
201  G(1, 2) = 0.0;
202  // phi row
203  G(2, 0) = 0.0;
204  G(2, 1) = 0.0;
205  G(2, 2) = (interpolated_X[0] * sin(interpolated_xi[1]) +
206  interpolated_X[1] * cos(interpolated_xi[1])) *
207  (interpolated_X[0] * sin(interpolated_xi[1]) +
208  interpolated_X[1] * cos(interpolated_xi[1]));
209 
210 
211  // Now calculate the stress tensor from the constitutive law
212  DenseMatrix<double> sigma(3);
213  get_stress(g, G, sigma);
214 
215  //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
216  // DISPLACEMENTS========
217 
218  // Loop over the test functions, nodes of the element
219  for (unsigned l = 0; l < n_node; l++)
220  {
221  // Loop of types of dofs
222  for (unsigned k = 0; k < n_position_type; k++)
223  {
224  // Radial displacemenet component
225  unsigned i = 0;
226  local_eqn = position_local_eqn(l, k, i);
227  /*IF it's not a boundary condition*/
228  if (local_eqn >= 0)
229  {
230  // Add the term for variations in radial position
231  residuals[local_eqn] +=
232  (sigma(0, 1) * interpolated_dXdxi[1][0] +
233  sigma(1, 1) * (interpolated_dXdxi[1][1] + interpolated_X[0]) +
234  sigma(2, 2) * sin(interpolated_xi[1]) *
235  (interpolated_X[0] * sin(interpolated_xi[1]) +
236  interpolated_X[1] * cos(interpolated_xi[1]))) *
237  psi(l, k) * W;
238 
239  // Add the terms for the variations in dX_{r}/dxi_{j}
240  for (unsigned j = 0; j < 2; j++)
241  {
242  residuals[local_eqn] +=
243  (sigma(j, 0) * interpolated_dXdxi[0][0] +
244  sigma(j, 1) *
245  (interpolated_dXdxi[0][1] - interpolated_X[1])) *
246  dpsidxi(l, k, j) * W;
247  }
248  }
249 
250  // Theta displacement component
251  i = 1;
252  local_eqn = position_local_eqn(l, k, i);
253  /*IF it's not a boundary condition*/
254  if (local_eqn >= 0)
255  {
256  // Add the term for variations in azimuthal position
257  residuals[local_eqn] +=
258  (-sigma(0, 1) * interpolated_dXdxi[0][0] -
259  sigma(1, 1) * (interpolated_dXdxi[0][1] - interpolated_X[1]) +
260  sigma(2, 2) * cos(interpolated_xi[1]) *
261  (interpolated_X[0] * sin(interpolated_xi[1]) +
262  interpolated_X[1] * cos(interpolated_xi[1]))) *
263  psi(l, k) * W;
264 
265  // Add the terms for the variations in dX_{theta}/dxi_{j}
266  for (unsigned j = 0; j < 2; j++)
267  {
268  residuals[local_eqn] +=
269  (sigma(j, 0) * interpolated_dXdxi[1][0] +
270  sigma(j, 1) *
271  (interpolated_dXdxi[1][1] + interpolated_X[0])) *
272  dpsidxi(l, k, j) * W;
273  }
274  }
275  } // End of loop over type of dof
276  } // End of loop over shape functions
277  } // End of loop over integration points
278  }
279 
280  // The jacobian is calculated by finite differences by default,
281  // could overload the get_jacobian function here if desired
282 
283  /// Overload/implement the function to calculate the volume of the element
284  double compute_physical_size() const
285  {
286  unsigned n_node = nnode();
287  unsigned n_position_type = 1;
288 
289  // Set up memory for the shape functions
290  Shape psi(n_node, n_position_type);
291  DShape dpsidxi(n_node, n_position_type, 2);
292 
293  // Set sum to zero
294  double sum = 0.0;
295  // Set the value of n_intpt
296  unsigned n_intpt = integral_pt()->nweight();
297 
298  // Loop over the integration points
299  // Loop in s1 direction*
300  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
301  {
302  // Get the integral weight
303  double w = integral_pt()->weight(ipt);
304  // Call the derivatives of the shape function wrt lagrangian coordinates
305  double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
306  // Premultiply the weights and the Jacobian
307  double W = w * J;
308 
309  // Calculate the local Lagrangian coordinates, position components
310  // and the derivatives of global position components
311  // wrt lagrangian coordinates
312  double interpolated_xi[2] = {0.0, 0.0};
313  double interpolated_X[2] = {0.0, 0.0};
314  double interpolated_dXdxi[2][2];
315 
316  // Initialise interpolated_dXdxi to zero
317  for (unsigned i = 0; i < 2; i++)
318  {
319  for (unsigned j = 0; j < 2; j++)
320  {
321  interpolated_dXdxi[i][j] = 0.0;
322  }
323  }
324 
325  // Calculate displacements and derivatives
326  for (unsigned l = 0; l < n_node; l++)
327  {
328  // Loop over positional dofs
329  for (unsigned k = 0; k < n_position_type; k++)
330  {
331  // Loop over displacement components (deformed position)
332  for (unsigned i = 0; i < 2; i++)
333  {
334  // Set the value of the lagrangian coordinate
335  interpolated_xi[i] +=
336  lagrangian_position_gen(l, k, i) * psi(l, k);
337  // Set the value of the position component
338  interpolated_X[i] += nodal_position_gen(l, k, i) * psi(l, k);
339  // Loop over Lagrangian derivative directions
340  for (unsigned j = 0; j < 2; j++)
341  {
342  // Calculate dX[i]/dxi_{j}
343  interpolated_dXdxi[i][j] +=
344  nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
345  }
346  }
347  }
348  }
349 
350  // Now calculate the deformed metric tensor
351  DenseMatrix<double> G(3);
352  // r row
353  G(0, 0) = interpolated_dXdxi[0][0] * interpolated_dXdxi[0][0] +
354  interpolated_dXdxi[1][0] * interpolated_dXdxi[1][0];
355  G(0, 1) = interpolated_dXdxi[0][0] *
356  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
357  interpolated_dXdxi[1][0] *
358  (interpolated_dXdxi[1][1] + interpolated_X[0]);
359  G(0, 2) = 0.0;
360  // theta row
361  G(1, 0) = G(0, 1);
362  G(1, 1) = (interpolated_dXdxi[0][1] - interpolated_X[1]) *
363  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
364  (interpolated_dXdxi[1][1] + interpolated_X[0]) *
365  (interpolated_dXdxi[1][1] + interpolated_X[0]);
366  G(1, 2) = 0.0;
367  // phi row
368  G(2, 0) = 0.0;
369  G(2, 1) = 0.0;
370  G(2, 2) = (interpolated_X[0] * sin(interpolated_xi[1]) +
371  interpolated_X[1] * cos(interpolated_xi[1])) *
372  (interpolated_X[0] * sin(interpolated_xi[1]) +
373  interpolated_X[1] * cos(interpolated_xi[1]));
374 
375  // Calculate the determinant of the metric tensor
376  double detG = G(0, 0) * G(1, 1) * G(2, 2) - G(0, 1) * G(1, 0) * G(2, 2);
377 
378  // Add the appropriate weight to the integral, i.e. sqrt of metric
379  // tensor
380  sum += W * sqrt(detG);
381  }
382 
383  // Return the volume, need to multiply by 2pi
384  return (2.0 * MathematicalConstants::Pi * sum);
385  }
386 
387  /// Assign the contribution to the residual using only finite differences
389  DenseMatrix<double>& jacobian)
390  {
391  // Add the solid contribution to the residuals
393 
394  // Solve for the consistent acceleration in Newmark scheme?
396  {
398  return;
399  }
400 
401  // Get the solid entries in the jacobian using finite differences
403  }
404 
405 
406  /// Overload the output function
407  void output(std::ostream& outfile)
408  {
409  FiniteElement::output(outfile);
410  }
411 
412  /// Output function
413  void output(std::ostream& outfile, const unsigned& n_plot)
414  {
415  // Set the output Vector
416  Vector<double> s(2);
417 
418  // Tecplot header info
419  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
420 
421  // Loop over element nodes
422  for (unsigned l2 = 0; l2 < n_plot; l2++)
423  {
424  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
425  for (unsigned l1 = 0; l1 < n_plot; l1++)
426  {
427  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
428 
429  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
430  double theta = interpolated_xi(s, 1);
431  // Output the x,y,u,v
432  // First output x and y assuming phi = 0
433  outfile << x_r * sin(theta) + x_theta * cos(theta) << " "
434  << x_r * cos(theta) - x_theta * sin(theta) << " ";
435  // Now output the true variables
436  for (unsigned i = 0; i < 2; i++)
437  outfile << interpolated_x(s, i) << " ";
438  for (unsigned i = 0; i < 2; i++)
439  outfile << interpolated_xi(s, i) << " ";
440  outfile << std::endl;
441  }
442  }
443  outfile << std::endl;
444  }
445 
446 
447  /// Overload the output function
448  void output(FILE* file_pt)
449  {
450  FiniteElement::output(file_pt);
451  }
452 
453 
454  /// Output function
455  void output(FILE* file_pt, const unsigned& n_plot)
456  {
457  // Set the output Vector
458  Vector<double> s(2);
459 
460  // Tecplot header info
461  fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
462 
463  // Loop over element nodes
464  for (unsigned l2 = 0; l2 < n_plot; l2++)
465  {
466  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
467  for (unsigned l1 = 0; l1 < n_plot; l1++)
468  {
469  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
470 
471  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
472  double theta = interpolated_xi(s, 1);
473  // Output the x,y,u,v
474  // First output x and y assuming phi = 0
475  // outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
476  // x_r*cos(theta) - x_theta*sin(theta) << " ";
477  fprintf(file_pt,
478  "%g %g ",
479  x_r * sin(theta) + x_theta * cos(theta),
480  x_r * cos(theta) - x_theta * sin(theta));
481 
482  // Now output the true variables
483  for (unsigned i = 0; i < 2; i++)
484  fprintf(file_pt, "%g ", interpolated_x(s, i));
485  for (unsigned i = 0; i < 2; i++)
486  fprintf(file_pt, "%g ", interpolated_xi(s, i));
487  fprintf(file_pt, "\n");
488  }
489  }
490  fprintf(file_pt, "\n");
491  }
492  };
493 
494 
495  //===========================================================================
496  /// An element that solved the AxisymmetricPVDEquations with
497  /// quadratic interpolation for the positions.
498  //===========================================================================
499  class AxisymQPVDElement : public SolidQElement<2, 3>,
501  {
502  public:
503  /// Constructor, there are no internal data points
505 
506  /// Overload the output function
507  void output(std::ostream& outfile)
508  {
509  FiniteElement::output(outfile);
510  }
511 
512  /// Output function
513  void output(std::ostream& outfile, const unsigned& n_plot)
514  {
515  AxisymmetricPVDEquations::output(outfile, n_plot);
516  }
517 
518  /// Overload the output function
519  void output(FILE* file_pt)
520  {
521  FiniteElement::output(file_pt);
522  }
523 
524  /// Output function
525  void output(FILE* file_pt, const unsigned& n_plot)
526  {
527  AxisymmetricPVDEquations::output(file_pt, n_plot);
528  }
529  };
530 
531  // Explicit definition of the face geometry for the AxisymQPVDElement element
532  template<>
533  class FaceGeometry<AxisymQPVDElement> : public virtual SolidQElement<1, 3>
534  {
535  public:
537  };
538 
539  //===========================================================================
540  /// An element that solved the AxisymmetricPVDEquations with
541  /// (diagonal) Hermite interpolation for the positions -- the
542  /// local and global (Lagrangian) coordinates are assumed to be aligned!
543  //===========================================================================
546  {
547  public:
548  /// Constructor, there are no internal data points
551  {
552  }
553 
554  /// Overload the output function
555  void output(std::ostream& outfile)
556  {
557  FiniteElement::output(outfile);
558  }
559 
560  /// Output function
561  void output(std::ostream& outfile, const unsigned& n_plot)
562  {
563  // Set the output Vector
564  Vector<double> s(2);
565 
566  // Tecplot header info
567  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
568 
569  // Loop over element nodes
570  for (unsigned l2 = 0; l2 < n_plot; l2++)
571  {
572  s[1] = 0.0 + l2 * 1.0 / (n_plot - 1);
573  for (unsigned l1 = 0; l1 < n_plot; l1++)
574  {
575  s[0] = 0.0 + l1 * 1.0 / (n_plot - 1);
576 
577  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
578  double theta = interpolated_xi(s, 1);
579  // Output the x,y,u,v
580  // First output x and y assuming phi = 0
581  outfile << x_r * sin(theta) + x_theta * cos(theta) << " "
582  << x_r * cos(theta) - x_theta * sin(theta) << " ";
583  // Now output the true variables
584  for (unsigned i = 0; i < 2; i++)
585  outfile << interpolated_x(s, i) << " ";
586  for (unsigned i = 0; i < 2; i++)
587  outfile << interpolated_xi(s, i) << " ";
588  outfile << std::endl;
589  }
590  }
591  outfile << std::endl;
592  }
593 
594 
595  /// Overload the output function
596  void output(FILE* file_pt)
597  {
598  FiniteElement::output(file_pt);
599  }
600 
601 
602  /// Output function
603  void output(FILE* file_pt, const unsigned& n_plot)
604  {
605  // Set the output Vector
606  Vector<double> s(2);
607 
608  // Tecplot header info
609  fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
610 
611  // Loop over element nodes
612  for (unsigned l2 = 0; l2 < n_plot; l2++)
613  {
614  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
615  for (unsigned l1 = 0; l1 < n_plot; l1++)
616  {
617  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
618 
619  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
620  double theta = interpolated_xi(s, 1);
621  // Output the x,y,u,v
622  // First output x and y assuming phi = 0
623  // outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
624  // x_r*cos(theta) - x_theta*sin(theta) << " ";
625  fprintf(file_pt,
626  "%g %g ",
627  x_r * sin(theta) + x_theta * cos(theta),
628  x_r * cos(theta) - x_theta * sin(theta));
629 
630  // Now output the true variables
631  for (unsigned i = 0; i < 2; i++)
632  fprintf(file_pt, "%g ", interpolated_x(s, i));
633  for (unsigned i = 0; i < 2; i++)
634  fprintf(file_pt, "%g ", interpolated_xi(s, i));
635  fprintf(file_pt, "\n");
636  }
637  }
638  fprintf(file_pt, "\n");
639  }
640  };
641 
642 
643  /// Explicit definition of the face geometry for the
644  // AxisymDiagHermitePVDElement element
645  template<>
647  : public virtual SolidDiagQHermiteElement<1>
648  {
649  public:
651  };
652 
653 
654  //=========================================================================
655  /// A class for elements that solve the equations of solid mechanics,
656  /// based on the principle of virtual displacements in
657  /// axisymmetric coordinates in a formulation that allows for
658  /// incompressibility or near incompressibility.
659  //==========================================================================
661  {
662  private:
663  /// Pointer to constitutive law
665 
666  /// Boolean to determine whether the solid is incompressible or not
668 
669  protected:
670  /// Access function that returns the local equation number for
671  /// the n-th solid pressure value.
672  virtual int solid_p_local_eqn(const unsigned& i) = 0;
673 
674  /// Return the solid pressure shape functions
675  virtual void solid_pshape(const Vector<double>& s, Shape& psi) const = 0;
676 
677  /// Return the stored solid shape functions at the knots
678  void solid_pshape_at_knot(const unsigned& ipt, Shape& psi) const;
679 
680  public:
681  /// Constructor, by default the element is not incompressible
684  {
685  }
686 
687  /// Return the constitutive law pointer
689  {
690  return Constitutive_law_pt;
691  }
692 
693  /// Return the stress tensor, as calculated from the constitutive law
694  /// in the Near-incompresible formulation
696  const DenseMatrix<double>& G,
697  DenseMatrix<double>& sigma,
698  DenseMatrix<double>& Gup,
699  double& pressure_stress,
700  double& kappa)
701  {
702 #ifdef PARANOID
703  // If the pointer to the constitutive law hasn't been set, issue an error
704  if (Constitutive_law_pt == 0)
705  {
706  std::string error_message =
707  "Elements derived from AxisymmetricPVDEquationsWithPressure";
708  error_message += " must have a constitutive law :\n ";
709  error_message +=
710  "set one using the constitutive_law_pt() member function\n";
711 
712  throw OomphLibError(
713  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
714  }
715 #endif
717  g, G, sigma, Gup, pressure_stress, kappa);
718  }
719 
720  /// Return the stress tensor, as calculated from the constitutive law
721  /// in the "true" incompresible formulation
723  const DenseMatrix<double>& G,
724  DenseMatrix<double>& sigma,
725  DenseMatrix<double>& Gup,
726  double& detG)
727  {
728 #ifdef PARANOID
729  // If the pointer to the constitutive law hasn't been set, issue an error
730  if (Constitutive_law_pt == 0)
731  {
732  std::string error_message =
733  "Elements derived from AxisymmetricPVDEquationsWithPressure";
734  error_message += " must have a constitutive law :\n ";
735  error_message +=
736  "set one using the constitutive_law_pt() member function\n";
737 
738  throw OomphLibError(
739  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
740  }
741 #endif
743  g, G, sigma, Gup, detG);
744  }
745 
746  /// Return whether the material is incompressible
747  bool is_incompressible() const
748  {
749  return Incompressible;
750  }
751 
752  /// Set the material to be incompressible
754  {
755  Incompressible = true;
756  }
757 
758  /// Set the material to be compressible
760  {
761  Incompressible = false;
762  }
763 
764  /// Return the number of solid pressure degrees of freedom
765  virtual unsigned nsolid_pres() const = 0;
766 
767  /// Return the lth solid pressures
768  virtual double solid_p(const unsigned& l) = 0;
769 
770  /// Return the residuals
772  {
773  // Call the generic residuals function with flag set to 0
774  // using a dummy matrix argument
776  residuals, GeneralisedElement::Dummy_matrix, 0);
777  }
778 
779  /// Return the residuals and the jacobian
781  DenseMatrix<double>& jacobian)
782  {
783  // Call the generic routine with the flag set to 1
785  residuals, jacobian, 1);
786  // Call the finite difference routine for the displacements
788  }
789 
790  /// Return the residuals for the equations of solid mechanics
791  /// formulated in the incompressible case!
793  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
794  {
795  // Set the number of Lagrangian coordinates
796  unsigned n_lagrangian = 2;
797  // Find out how many nodes there are
798  unsigned n_node = nnode();
799  // Find out how many positional dofs there are
800  unsigned n_position_type = nnodal_position_type();
801  // Find out how many pressure dofs there are
802  unsigned n_solid_pres = nsolid_pres();
803 
804  // Integers to store the local equation and unknown numbers
805  int local_eqn = 0, local_unknown = 0;
806 
807  // Set up memory for the shape functions
808  Shape psi(n_node, n_position_type);
809  DShape dpsidxi(n_node, n_position_type, n_lagrangian);
810 
811  // Set up memory for the pressure shape functions
812  Shape psisp(n_solid_pres);
813 
814  // Set the value of n_intpt
815  unsigned n_intpt = integral_pt()->nweight();
816 
817  // Loop over the integration points
818  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
819  {
820  // Get the integral weight
821  double w = integral_pt()->weight(ipt);
822  // Call the derivatives of the shape functions
823  double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
824  // Call the pressure shape functions
825  solid_pshape_at_knot(ipt, psisp);
826 
827  // Premultiply the weights and the Jacobian
828  double W = w * J;
829 
830  // Calculate the local Lagrangian coordinates, position components
831  // and the derivatives of global position components
832  // wrt lagrangian coordinates
833  double interpolated_xi[2] = {0.0, 0.0};
834  double interpolated_X[2] = {0.0, 0.0};
835  double interpolated_dXdxi[2][2];
836  double interpolated_solid_p = 0.0;
837 
838  // Initialise interpolated_dXdxi to zero
839  for (unsigned i = 0; i < 2; i++)
840  {
841  for (unsigned j = 0; j < 2; j++)
842  {
843  interpolated_dXdxi[i][j] = 0.0;
844  }
845  }
846 
847  // Calculate displacements and derivatives
848  for (unsigned l = 0; l < n_node; l++)
849  {
850  // Loop over positional dofs
851  for (unsigned k = 0; k < n_position_type; k++)
852  {
853  // Loop over displacement components (deformed position)
854  for (unsigned i = 0; i < 2; i++)
855  {
856  // Set the value of the lagrangian coordinate
857  interpolated_xi[i] +=
858  lagrangian_position_gen(l, k, i) * psi(l, k);
859  // Set the value of the position component
860  interpolated_X[i] += nodal_position_gen(l, k, i) * psi(l, k);
861  // Loop over Lagrangian derivative directions
862  for (unsigned j = 0; j < 2; j++)
863  {
864  // Calculate dX[i]/dxi_{j}
865  interpolated_dXdxi[i][j] +=
866  nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
867  }
868  }
869  }
870  }
871 
872  // Calculate the local internal pressure
873  for (unsigned l = 0; l < n_solid_pres; l++)
874  {
875  interpolated_solid_p += solid_p(l) * psisp[l];
876  }
877 
878  // We are now in a position to calculate the undeformed metric tensor
879  DenseMatrix<double> g(3);
880  // r row
881  g(0, 0) = 1.0;
882  g(0, 1) = 0.0;
883  g(0, 2) = 0.0;
884  // theta row
885  g(1, 0) = 0.0;
886  g(1, 1) = interpolated_xi[0] * interpolated_xi[0];
887  g(1, 2) = 0.0;
888  // phi row
889  g(2, 0) = 0.0;
890  g(2, 1) = 0.0;
891  g(2, 2) = interpolated_xi[0] * interpolated_xi[0] *
892  sin(interpolated_xi[1]) * sin(interpolated_xi[1]);
893 
894  // Find the determinant of the undeformed metric tensor
895  double detg = g(0, 0) * g(1, 1) * g(2, 2);
896 
897  // Now multiply the weight by the square-root of the determinant of the
898  // undeformed metric tensor r^2 sin(theta)
899  W *= sqrt(detg);
900 
901  // Now calculate the deformed metric tensor
902  DenseMatrix<double> G(3);
903  // r row
904  G(0, 0) = interpolated_dXdxi[0][0] * interpolated_dXdxi[0][0] +
905  interpolated_dXdxi[1][0] * interpolated_dXdxi[1][0];
906  G(0, 1) = interpolated_dXdxi[0][0] *
907  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
908  interpolated_dXdxi[1][0] *
909  (interpolated_dXdxi[1][1] + interpolated_X[0]);
910  G(0, 2) = 0.0;
911  // theta row
912  G(1, 0) = G(0, 1);
913  G(1, 1) = (interpolated_dXdxi[0][1] - interpolated_X[1]) *
914  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
915  (interpolated_dXdxi[1][1] + interpolated_X[0]) *
916  (interpolated_dXdxi[1][1] + interpolated_X[0]);
917  G(1, 2) = 0.0;
918  // phi row
919  G(2, 0) = 0.0;
920  G(2, 1) = 0.0;
921  G(2, 2) = (interpolated_X[0] * sin(interpolated_xi[1]) +
922  interpolated_X[1] * cos(interpolated_xi[1])) *
923  (interpolated_X[0] * sin(interpolated_xi[1]) +
924  interpolated_X[1] * cos(interpolated_xi[1]));
925 
926 
927  // Now calculate the stress tensor from the constitutive law
928  DenseMatrix<double> sigma(3), Gup(3);
929  double detG = 0.0, pressure_stress = 0.0, kappa = 0.0;
930  // If it's incompressible call one form of the constitutive law
931  if (Incompressible)
932  {
933  get_stress(g, G, sigma, Gup, detG);
934  }
935  // Otherwise call another form
936  else
937  {
938  get_stress(g, G, sigma, Gup, pressure_stress, kappa);
939  }
940 
941 
942  //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
943  // DISPLACEMENTS========
944 
945  // Loop over the test functions, nodes of the element
946  for (unsigned l = 0; l < n_node; l++)
947  {
948  // Loop of types of dofs
949  for (unsigned k = 0; k < n_position_type; k++)
950  {
951  // Radial displacemenet component
952  unsigned i = 0;
953  local_eqn = position_local_eqn(l, k, i);
954  /*IF it's not a boundary condition*/
955  if (local_eqn >= 0)
956  {
957  // Add the term for variations in radial position
958  residuals[local_eqn] +=
959  ((sigma(0, 1) + interpolated_solid_p * Gup(0, 1)) *
960  interpolated_dXdxi[1][0] +
961  (sigma(1, 1) + interpolated_solid_p * Gup(1, 1)) *
962  (interpolated_dXdxi[1][1] + interpolated_X[0]) +
963  (sigma(2, 2) + interpolated_solid_p * Gup(2, 2)) *
964  sin(interpolated_xi[1]) *
965  (interpolated_X[0] * sin(interpolated_xi[1]) +
966  interpolated_X[1] * cos(interpolated_xi[1]))) *
967  psi(l, k) * W;
968 
969  // Add the terms for the variations in dX_{r}/dxi_{j}
970  for (unsigned j = 0; j < 2; j++)
971  {
972  residuals[local_eqn] +=
973  ((sigma(j, 0) + interpolated_solid_p * Gup(j, 0)) *
974  interpolated_dXdxi[0][0] +
975  (sigma(j, 1) + interpolated_solid_p * Gup(j, 1)) *
976  (interpolated_dXdxi[0][1] - interpolated_X[1])) *
977  dpsidxi(l, k, j) * W;
978  }
979 
980  // Can add in the pressure jacobian terms
981  if (flag)
982  {
983  // Loop over the pressure nodes
984  for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
985  {
986  local_unknown = solid_p_local_eqn(l2);
987  // If it's not a boundary condition
988  if (local_unknown >= 0)
989  {
990  jacobian(local_eqn, local_unknown) +=
991  (psisp[l2] * Gup(0, 1) * interpolated_dXdxi[1][0] +
992  psisp[l2] * Gup(1, 1) *
993  (interpolated_dXdxi[1][1] + interpolated_X[0]) +
994  psisp[l2] * Gup(2, 2) * sin(interpolated_xi[1]) *
995  (interpolated_X[0] * sin(interpolated_xi[1]) +
996  interpolated_X[1] * cos(interpolated_xi[1]))) *
997  psi(l, k) * W;
998 
999  // Add the terms for the variations in dX_{r}/dxi_{j}
1000  for (unsigned j = 0; j < 2; j++)
1001  {
1002  jacobian(local_eqn, local_unknown) +=
1003  (psisp[l2] * Gup(j, 0) * interpolated_dXdxi[0][0] +
1004  psisp[l2] * Gup(j, 1) *
1005  (interpolated_dXdxi[0][1] - interpolated_X[1])) *
1006  dpsidxi(l, k, j) * W;
1007  }
1008  } // End of if not boundary condition
1009  }
1010  } // End of if(flag)
1011  } // End of if Position_local_eqn
1012 
1013  // Theta displacement component
1014  i = 1;
1015  local_eqn = position_local_eqn(l, k, i);
1016  /*IF it's not a boundary condition*/
1017  if (local_eqn >= 0)
1018  {
1019  // Add the term for variations in azimuthal position
1020  residuals[local_eqn] +=
1021  (-(sigma(0, 1) + interpolated_solid_p * Gup(0, 1)) *
1022  interpolated_dXdxi[0][0] -
1023  (sigma(1, 1) + interpolated_solid_p * Gup(1, 1)) *
1024  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
1025  (sigma(2, 2) + interpolated_solid_p * Gup(2, 2)) *
1026  cos(interpolated_xi[1]) *
1027  (interpolated_X[0] * sin(interpolated_xi[1]) +
1028  interpolated_X[1] * cos(interpolated_xi[1]))) *
1029  psi(l, k) * W;
1030 
1031  // Add the terms for the variations in dX_{theta}/dxi_{j}
1032  for (unsigned j = 0; j < 2; j++)
1033  {
1034  residuals[local_eqn] +=
1035  ((sigma(j, 0) + interpolated_solid_p * Gup(j, 0)) *
1036  interpolated_dXdxi[1][0] +
1037  (sigma(j, 1) + interpolated_solid_p * Gup(j, 1)) *
1038  (interpolated_dXdxi[1][1] + interpolated_X[0])) *
1039  dpsidxi(l, k, j) * W;
1040  }
1041 
1042  // Can add in the pressure jacobian terms
1043  if (flag)
1044  {
1045  // Loop over the pressure nodes
1046  for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1047  {
1048  local_unknown = solid_p_local_eqn(l2);
1049  // If it's not a boundary condition
1050  if (local_unknown >= 0)
1051  {
1052  // Add the term for variations in azimuthal position
1053  jacobian(local_eqn, local_unknown) +=
1054  (-psisp[l2] * Gup(0, 1) * interpolated_dXdxi[0][0] -
1055  psisp[l2] * Gup(1, 1) *
1056  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
1057  psisp[l2] * Gup(2, 2) * cos(interpolated_xi[1]) *
1058  (interpolated_X[0] * sin(interpolated_xi[1]) +
1059  interpolated_X[1] * cos(interpolated_xi[1]))) *
1060  psi(l, k) * W;
1061 
1062  // Add the terms for the variations in dX_{theta}/dxi_{j}
1063  for (unsigned j = 0; j < 2; j++)
1064  {
1065  jacobian(local_eqn, local_unknown) +=
1066  (psisp[l2] * Gup(j, 0) * interpolated_dXdxi[1][0] +
1067  psisp[l2] * Gup(j, 1) *
1068  (interpolated_dXdxi[1][1] + interpolated_X[0])) *
1069  dpsidxi(l, k, j) * W;
1070  }
1071  }
1072  }
1073  } // End of if(flag)
1074  } // End of Position_local_eqn
1075  } // End of loop over type of dof
1076  } // End of loop over shape functions
1077 
1078 
1079  //======================CONSTRAINT EQUATIONS FOR THE PRESSURE===========
1080 
1081 
1082  // Now loop over the pressure degrees of freedom
1083  for (unsigned l = 0; l < n_solid_pres; l++)
1084  {
1085  local_eqn = solid_p_local_eqn(l);
1086  // If it's not a bondary condition
1087  if (local_eqn >= 0)
1088  {
1089  // For true incompressibility we need the ratio of determinants of
1090  // the metric tensors to be exactly 1.0
1091  if (Incompressible)
1092  {
1093  residuals[local_eqn] += (detG / detg - 1.0) * psisp[l] * W;
1094 
1095  // No Jacobian terms since the pressure does not feature
1096  // in the incompressibility constraint
1097  }
1098  else
1099  {
1100  // Otherwise the pressure must be that calculated by the
1101  // constitutive law
1102  residuals[local_eqn] +=
1103  (kappa * interpolated_solid_p - pressure_stress) * psisp[l] * W;
1104 
1105  // Add in the jacobian terms
1106  if (flag)
1107  {
1108  // Loop over the pressure nodes again
1109  for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1110  {
1111  local_unknown = solid_p_local_eqn(l2);
1112  // If not a boundary condition
1113  if (local_unknown >= 0)
1114  {
1115  jacobian(local_eqn, local_unknown) +=
1116  kappa * psisp[l2] * psisp[l] * W;
1117  }
1118  }
1119  } // End of jacobian terms
1120  } // End of else
1121 
1122  } // End of if not boundary condition
1123  }
1124 
1125  } // End of loop over integration points
1126  }
1127 
1128  // The jacobian is calculated by finite differences by default,
1129  // could overload the get_jacobian function here if desired
1130 
1131  /// Overload/implement the size function
1132  double compute_physical_size() const
1133  {
1134  unsigned n_node = nnode();
1135  unsigned n_position_type = 1;
1136 
1137  // Set up memory for the shape functions
1138  Shape psi(n_node, n_position_type);
1139  DShape dpsidxi(n_node, n_position_type, 2);
1140 
1141  // Set sum to zero
1142  double sum = 0.0;
1143  // Set the value of n_intpt
1144  unsigned n_intpt = integral_pt()->nweight();
1145 
1146  // Loop over the integration points
1147  // Loop in s1 direction*
1148  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1149  {
1150  // Get the integral weight
1151  double w = integral_pt()->weight(ipt);
1152  // Call the derivatives of the shape function wrt lagrangian coordinates
1153  double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
1154  // Premultiply the weights and the Jacobian
1155  double W = w * J;
1156 
1157  // Calculate the local Lagrangian coordinates, position components
1158  // and the derivatives of global position components
1159  // wrt lagrangian coordinates
1160  double interpolated_xi[2] = {0.0, 0.0};
1161  double interpolated_X[2] = {0.0, 0.0};
1162  double interpolated_dXdxi[2][2];
1163 
1164  // Initialise interpolated_dXdxi to zero
1165  for (unsigned i = 0; i < 2; i++)
1166  {
1167  for (unsigned j = 0; j < 2; j++)
1168  {
1169  interpolated_dXdxi[i][j] = 0.0;
1170  }
1171  }
1172 
1173  // Calculate displacements and derivatives
1174  for (unsigned l = 0; l < n_node; l++)
1175  {
1176  // Loop over positional dofs
1177  for (unsigned k = 0; k < n_position_type; k++)
1178  {
1179  // Loop over displacement components (deformed position)
1180  for (unsigned i = 0; i < 2; i++)
1181  {
1182  // Set the value of the lagrangian coordinate
1183  interpolated_xi[i] +=
1184  lagrangian_position_gen(l, k, i) * psi(l, k);
1185  // Set the value of the position component
1186  interpolated_X[i] += nodal_position_gen(l, k, i) * psi(l, k);
1187  // Loop over Lagrangian derivative directions
1188  for (unsigned j = 0; j < 2; j++)
1189  {
1190  // Calculate dX[i]/dxi_{j}
1191  interpolated_dXdxi[i][j] +=
1192  nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
1193  }
1194  }
1195  }
1196  }
1197 
1198  // Now calculate the deformed metric tensor
1199  DenseMatrix<double> G(3);
1200  // r row
1201  G(0, 0) = interpolated_dXdxi[0][0] * interpolated_dXdxi[0][0] +
1202  interpolated_dXdxi[1][0] * interpolated_dXdxi[1][0];
1203  G(0, 1) = interpolated_dXdxi[0][0] *
1204  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
1205  interpolated_dXdxi[1][0] *
1206  (interpolated_dXdxi[1][1] + interpolated_X[0]);
1207  G(0, 2) = 0.0;
1208  // theta row
1209  G(1, 0) = G(0, 1);
1210  G(1, 1) = (interpolated_dXdxi[0][1] - interpolated_X[1]) *
1211  (interpolated_dXdxi[0][1] - interpolated_X[1]) +
1212  (interpolated_dXdxi[1][1] + interpolated_X[0]) *
1213  (interpolated_dXdxi[1][1] + interpolated_X[0]);
1214  G(1, 2) = 0.0;
1215  // phi row
1216  G(2, 0) = 0.0;
1217  G(2, 1) = 0.0;
1218  G(2, 2) = (interpolated_X[0] * sin(interpolated_xi[1]) +
1219  interpolated_X[1] * cos(interpolated_xi[1])) *
1220  (interpolated_X[0] * sin(interpolated_xi[1]) +
1221  interpolated_X[1] * cos(interpolated_xi[1]));
1222 
1223  // Calculate the determinant of the metric tensor
1224  double detG = G(0, 0) * G(1, 1) * G(2, 2) - G(0, 1) * G(1, 0) * G(2, 2);
1225 
1226  // Add the appropriate weight to the integral, i.e. sqrt of metric
1227  // tensor
1228  sum += W * sqrt(detG);
1229  }
1230 
1231  // Return the volume
1232  return (2.0 * MathematicalConstants::Pi * sum);
1233  }
1234 
1235  /// Return the interpolated_solid_pressure
1237  {
1238  // Find number of nodes
1239  unsigned n_solid_pres = nsolid_pres();
1240  // Local shape function
1241  Shape psisp(n_solid_pres);
1242  // Find values of shape function
1243  solid_pshape(s, psisp);
1244 
1245  // Initialise value of solid_p
1246  double interpolated_solid_p = 0.0;
1247  // Loop over the local nodes and sum
1248  for (unsigned l = 0; l < n_solid_pres; l++)
1249  {
1250  interpolated_solid_p += solid_p(l) * psisp[l];
1251  }
1252 
1253  return (interpolated_solid_p);
1254  }
1255 
1256  /// Overload the output function
1257  void output(std::ostream& outfile)
1258  {
1259  FiniteElement::output(outfile);
1260  }
1261 
1262  /// Output function
1263  void output(std::ostream& outfile, const unsigned& n_plot)
1264  {
1265  // Set the output Vector
1266  Vector<double> s(2);
1267 
1268  // Tecplot header info
1269  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
1270 
1271  // Loop over element nodes
1272  for (unsigned l2 = 0; l2 < n_plot; l2++)
1273  {
1274  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
1275  for (unsigned l1 = 0; l1 < n_plot; l1++)
1276  {
1277  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
1278 
1279  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
1280  double theta = interpolated_xi(s, 1);
1281  // Output the x,y,u,v
1282  // First output x and y assuming phi = 0
1283  outfile << x_r * sin(theta) + x_theta * cos(theta) << " "
1284  << x_r * cos(theta) - x_theta * sin(theta) << " ";
1285  outfile << interpolated_solid_p(s) << " ";
1286  // Now output the true variables
1287  for (unsigned i = 0; i < 2; i++)
1288  outfile << interpolated_x(s, i) << " ";
1289  for (unsigned i = 0; i < 2; i++)
1290  outfile << interpolated_xi(s, i) << " ";
1291  outfile << std::endl;
1292  }
1293  }
1294  outfile << std::endl;
1295  }
1296 
1297  /// Overload the output function
1298  void output(FILE* file_pt)
1299  {
1300  FiniteElement::output(file_pt);
1301  }
1302 
1303 
1304  /// Output function
1305  void output(FILE* file_pt, const unsigned& n_plot)
1306  {
1307  // Set the output Vector
1308  Vector<double> s(2);
1309 
1310  // Tecplot header info
1311  fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
1312 
1313  // Loop over element nodes
1314  for (unsigned l2 = 0; l2 < n_plot; l2++)
1315  {
1316  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
1317  for (unsigned l1 = 0; l1 < n_plot; l1++)
1318  {
1319  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
1320 
1321  double x_r = interpolated_x(s, 0), x_theta = interpolated_x(s, 1);
1322  double theta = interpolated_xi(s, 1);
1323  // Output the x,y,u,v
1324  // First output x and y assuming phi = 0
1325  // outfile << x_r*sin(theta) + x_theta*cos(theta) << " " <<
1326  // x_r*cos(theta) - x_theta*sin(theta) << " ";
1327  fprintf(file_pt,
1328  "%g %g ",
1329  x_r * sin(theta) + x_theta * cos(theta),
1330  x_r * cos(theta) - x_theta * sin(theta));
1331 
1332  // Now output the true variables
1333  for (unsigned i = 0; i < 2; i++)
1334  fprintf(file_pt, "%g ", interpolated_x(s, i));
1335  for (unsigned i = 0; i < 2; i++)
1336  fprintf(file_pt, "%g ", interpolated_xi(s, i));
1337  fprintf(file_pt, "\n");
1338  }
1339  }
1340  fprintf(file_pt, "\n");
1341  }
1342  };
1343 
1344  //========================================================================
1345  /// An Element that solves the Axisymmetric principle of virtual displacements
1346  /// with separately interpolated pressure, discontinuous interpolation.
1347  //=========================================================================
1349  : public SolidQElement<2, 3>,
1351  {
1352  /// Internal index that indicates at which internal data value the
1353  /// solid pressure is stored
1355 
1356  /// Overload the access function for the solid pressure equation numbers
1357  inline int solid_p_local_eqn(const unsigned& i)
1358  {
1360  }
1361 
1362  /// Return the pressure shape functions
1363  inline void solid_pshape(const Vector<double>& s, Shape& psi) const;
1364 
1365  public:
1366  /// There is internal solid data so we can't use the automatic
1367  /// assignment of consistent initial conditions for time-dependent problems.
1369  {
1370  return true;
1371  }
1372 
1373  /// Constructor, there are 3 internal data items
1376  {
1377  // Allocate and add one Internal data object that stores 3 pressure
1378  // values
1380  }
1381 
1382  /// Return the l-th pressure value
1383  double solid_p(const unsigned& l)
1384  {
1385  return this->internal_data_pt(P_solid_internal_index)->value(l);
1386  }
1387 
1388  /// Return number of pressure values
1389  unsigned nsolid_pres() const
1390  {
1391  return 3;
1392  }
1393 
1394  /// Fix the pressure dof l to the value pvalue
1395  void fix_solid_pressure(const unsigned& l, const double& pvalue)
1396  {
1397  this->internal_data_pt(P_solid_internal_index)->pin(l);
1398  this->internal_data_pt(P_solid_internal_index)->set_value(l, pvalue);
1399  }
1400 
1401  /// Overload the output function
1402  void output(std::ostream& outfile)
1403  {
1404  FiniteElement::output(outfile);
1405  }
1406 
1407  /// Output function
1408  void output(std::ostream& outfile, const unsigned& n_plot)
1409  {
1411  }
1412 
1413 
1414  /// Overload the output function
1415  void output(FILE* file_pt)
1416  {
1417  FiniteElement::output(file_pt);
1418  }
1419 
1420  /// Output function
1421  void output(FILE* file_pt, const unsigned& n_plot)
1422  {
1424  }
1425  };
1426 
1427  /// Define the pressure shape functions
1429  const Vector<double>& s, Shape& psi) const
1430  {
1431  psi[0] = 1.0;
1432  psi[1] = s[0];
1433  psi[2] = s[1];
1434  }
1435 
1436  // Explicit definition of the face geometry for the AxisymQPVDElement element
1437  template<>
1439  : public virtual SolidQElement<1, 3>
1440  {
1441  public:
1443  };
1444 
1445 } // namespace oomph
1446 
1447 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
An element that solved the AxisymmetricPVDEquations with (diagonal) Hermite interpolation for the pos...
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt)
Overload the output function.
void output(FILE *file_pt, const unsigned &n_plot)
Output function.
AxisymDiagHermitePVDElement()
Constructor, there are no internal data points.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
An Element that solves the Axisymmetric principle of virtual displacements with separately interpolat...
unsigned P_solid_internal_index
Internal index that indicates at which internal data value the solid pressure is stored.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt, const unsigned &n_plot)
Output function.
double solid_p(const unsigned &l)
Return the l-th pressure value.
AxisymQPVDElementWithPressure()
Constructor, there are 3 internal data items.
int solid_p_local_eqn(const unsigned &i)
Overload the access function for the solid pressure equation numbers.
void output(FILE *file_pt)
Overload the output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
unsigned nsolid_pres() const
Return number of pressure values.
void solid_pshape(const Vector< double > &s, Shape &psi) const
Return the pressure shape functions.
void fix_solid_pressure(const unsigned &l, const double &pvalue)
Fix the pressure dof l to the value pvalue.
bool has_internal_solid_data()
There is internal solid data so we can't use the automatic assignment of consistent initial condition...
An element that solved the AxisymmetricPVDEquations with quadratic interpolation for the positions.
AxisymQPVDElement()
Constructor, there are no internal data points.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void output(FILE *file_pt)
Overload the output function.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt, const unsigned &n_plot)
Output function.
A class for elements that solve the equations of solid mechanics, based on the principle of virtual d...
void solid_pshape_at_knot(const unsigned &ipt, Shape &psi) const
Return the stored solid shape functions at the knots.
bool Incompressible
Boolean to determine whether the solid is incompressible or not.
void fill_in_generic_residual_contribution_axisym_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Return the residuals for the equations of solid mechanics formulated in the incompressible case!
void output(FILE *file_pt, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Return the residuals and the jacobian.
virtual unsigned nsolid_pres() const =0
Return the number of solid pressure degrees of freedom.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
AxisymmetricPVDEquationsWithPressure()
Constructor, by default the element is not incompressible.
void set_incompressible()
Set the material to be incompressible.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
virtual int solid_p_local_eqn(const unsigned &i)=0
Access function that returns the local equation number for the n-th solid pressure value.
bool is_incompressible() const
Return whether the material is incompressible.
double compute_physical_size() const
Overload/implement the size function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
virtual void solid_pshape(const Vector< double > &s, Shape &psi) const =0
Return the solid pressure shape functions.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt)
Overload the output function.
double interpolated_solid_p(const Vector< double > &s)
Return the interpolated_solid_pressure.
virtual double solid_p(const unsigned &l)=0
Return the lth solid pressures.
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma, DenseMatrix< double > &Gup, double &pressure_stress, double &kappa)
Return the stress tensor, as calculated from the constitutive law in the Near-incompresible formulati...
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma, DenseMatrix< double > &Gup, double &detG)
Return the stress tensor, as calculated from the constitutive law in the "true" incompresible formula...
void set_compressible()
Set the material to be compressible.
A class for elements that solve the equations of solid mechanics, based on the principle of virtual d...
void fill_in_contribution_to_residuals_axisym_pvd(Vector< double > &residuals)
Return the residuals for the equations of solid mechanics.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void output(FILE *file_pt)
Overload the output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals by calling the generic function.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Assign the contribution to the residual using only finite differences.
double compute_physical_size() const
Overload/implement the function to calculate the volume of the element.
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Return the stress tensor, as calculated from the constitutive law.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)=0
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:293
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
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
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n.
Definition: elements.h:2349
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data.
Definition: elements.h:267
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
SolidQHermiteElements in which we assume the local and global coordinates to be aligned so that the J...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
double lagrangian_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return Generalised Lagrangian coordinate at local node n. ‘Direction’ i, ‘Type’ k.
Definition: elements.h:3912
bool Solve_for_consistent_newmark_accel_flag
Flag to indicate which system of equations to solve when assigning initial conditions for time-depend...
Definition: elements.h:4302
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:4137
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s.
Definition: elements.cc:7104
virtual double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Return the geometric shape functions and also first derivatives w.r.t. Lagrangian coordinates at ipt-...
Definition: elements.cc:6737
virtual void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use finite differences to calculate the Jacobian entries corresponding to the solid positions....
Definition: elements.cc:6985
void fill_in_jacobian_for_newmark_accel(DenseMatrix< double > &jacobian)
Fill in the contributions of the Jacobian matrix for the consistent assignment of the initial "accele...
Definition: elements.cc:7227
SolidQElement elements are quadrilateral elements whose derivatives also include those based upon the...
Definition: Qelements.h:1742
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const double Pi
50 digits from maple
//////////////////////////////////////////////////////////////////// ////////////////////////////////...