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