solid_elements.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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 //Non-inline functions for elements that solve the principle of virtual
31 //equations of solid mechanics
32 
33 #include "solid_elements.h"
34 
35 
36 namespace oomph
37 {
38 
39 
40 /// Static default value for timescale ratio (1.0 -- for natural scaling)
41 template <unsigned DIM>
42 double PVDEquationsBase<DIM>::Default_lambda_sq_value=1.0;
43 
44 
45 
46 //////////////////////////////////////////////////////////////////
47 //////////////////////////////////////////////////////////////////
48 //////////////////////////////////////////////////////////////////
49 
50 //======================================================================
51 /// Compute the strain tensor at local coordinate s
52 //======================================================================
53 template<unsigned DIM>
55  DenseMatrix<double> &strain) const
56 {
57 #ifdef PARANOID
58  if ((strain.ncol()!=DIM)||(strain.nrow()!=DIM))
59  {
60  std::ostringstream error_message;
61  error_message << "Strain matrix is " << strain.ncol() << " x "
62  << strain.nrow() << ", but dimension of the equations is "
63  << DIM << std::endl;
64  throw OomphLibError(error_message.str(),
65  OOMPH_CURRENT_FUNCTION,
66  OOMPH_EXCEPTION_LOCATION);
67  }
68 #endif
69 
70  //Find out how many nodes there are in the element
71  const unsigned n_node = nnode();
72 
73  //Find out how many position types there are
74  const unsigned n_position_type = this->nnodal_position_type();
75 
76  //Set up memory for the shape and test functions
77  Shape psi(n_node,n_position_type);
78  DShape dpsidxi(n_node,n_position_type,DIM);
79 
80  //Call the derivatives of the shape functions
81  (void) dshape_lagrangian(s,psi,dpsidxi);
82 
83  //Calculate interpolated values of the derivative of global position
84  DenseMatrix<double> interpolated_G(DIM);
85 
86  //Initialise to zero
87  for(unsigned i=0;i<DIM;i++)
88  {
89  for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}
90  }
91 
92  //Storage for Lagrangian coordinates (initialised to zero)
93  Vector<double> interpolated_xi(DIM,0.0);
94 
95  //Loop over nodes
96  for(unsigned l=0;l<n_node;l++)
97  {
98  //Loop over the positional dofs
99  for(unsigned k=0;k<n_position_type;k++)
100  {
101  //Loop over velocity components
102  for(unsigned i=0;i<DIM;i++)
103  {
104  //Calculate the Lagrangian coordinates
105  interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);
106 
107  //Loop over derivative directions
108  for(unsigned j=0;j<DIM;j++)
109  {
110  interpolated_G(j,i) +=
111  this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
112  }
113  }
114  }
115  }
116 
117  //Get isotropic growth factor
118  double gamma=1.0;
119  //Dummy integration point
120  unsigned ipt=0;
121  get_isotropic_growth(ipt,s,interpolated_xi,gamma);
122 
123  // We use Cartesian coordinates as the reference coordinate
124  // system. In this case the undeformed metric tensor is always
125  // the identity matrix -- stretched by the isotropic growth
126  double diag_entry=pow(gamma,2.0/double(DIM));
127  DenseMatrix<double> g(DIM);
128  for(unsigned i=0;i<DIM;i++)
129  {
130  for(unsigned j=0;j<DIM;j++)
131  {
132  if(i==j) {g(i,j) = diag_entry;} else {g(i,j) = 0.0;}
133  }
134  }
135 
136 
137 
138  //Declare and calculate the deformed metric tensor
139  DenseMatrix<double> G(DIM);
140 
141  //Assign values of G
142  for(unsigned i=0;i<DIM;i++)
143  {
144  //Do upper half of matrix
145  for(unsigned j=i;j<DIM;j++)
146  {
147  //Initialise G(i,j) to zero
148  G(i,j) = 0.0;
149  //Now calculate the dot product
150  for(unsigned k=0;k<DIM;k++)
151  {
152  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
153  }
154  }
155  //Matrix is symmetric so just copy lower half
156  for(unsigned j=0;j<i;j++)
157  {
158  G(i,j) = G(j,i);
159  }
160  }
161 
162  //Fill in the strain tensor
163  for(unsigned i=0;i<DIM;i++)
164  {
165  for(unsigned j=0;j<DIM;j++) {strain(i,j)= 0.5*(G(i,j) - g(i,j));}
166  }
167 
168 }
169 
170 //=======================================================================
171 /// Compute the residuals for the discretised principle of
172 /// virtual displacements.
173 //=======================================================================
174 template <unsigned DIM>
177  DenseMatrix<double> &jacobian,
178  const unsigned& flag)
179 {
180 
181 #ifdef PARANOID
182  // Check if the constitutive equation requires the explicit imposition of an
183  // incompressibility constraint
185  {
186  throw OomphLibError(
187  "PVDEquations cannot be used with incompressible constitutive laws.",
188  OOMPH_CURRENT_FUNCTION,
189  OOMPH_EXCEPTION_LOCATION);
190  }
191 #endif
192 
193  // Simply set up initial condition?
194  if (this->Solid_ic_pt!=0)
195  {
196  this->fill_in_residuals_for_solid_ic(residuals);
197  return;
198  }
199 
200  //Find out how many nodes there are
201  const unsigned n_node = this->nnode();
202 
203  //Find out how many positional dofs there are
204  const unsigned n_position_type = this->nnodal_position_type();
205 
206  //Set up memory for the shape functions
207  Shape psi(n_node,n_position_type);
208  DShape dpsidxi(n_node,n_position_type,DIM);
209 
210  //Set the value of Nintpt -- the number of integration points
211  const unsigned n_intpt = this->integral_pt()->nweight();
212 
213  //Set the vector to hold the local coordinates in the element
214  Vector<double> s(DIM);
215 
216  // Timescale ratio (non-dim density)
217  double lambda_sq = this->lambda_sq();
218 
219  // Time factor
220  double time_factor=0.0;
221  if (lambda_sq>0)
222  {
223  time_factor=this->node_pt(0)->position_time_stepper_pt()->weight(2,0);
224  }
225 
226  //Integer to store the local equation number
227  int local_eqn=0;
228 
229  //Loop over the integration points
230  for(unsigned ipt=0;ipt<n_intpt;ipt++)
231  {
232  //Assign the values of s
233  for(unsigned i=0;i<DIM;++i) {s[i] = this->integral_pt()->knot(ipt,i);}
234 
235  //Get the integral weight
236  double w = this->integral_pt()->weight(ipt);
237 
238  //Call the derivatives of the shape functions (and get Jacobian)
239  double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
240 
241  //Calculate interpolated values of the derivative of global position
242  //wrt lagrangian coordinates
243  DenseMatrix<double> interpolated_G(DIM);
244 
245  // Setup memory for accelerations
246  Vector<double> accel(DIM);
247 
248  //Initialise to zero
249  for(unsigned i=0;i<DIM;i++)
250  {
251  // Initialise acclerations
252  accel[i]=0.0;
253  for(unsigned j=0;j<DIM;j++)
254  {
255  interpolated_G(i,j) = 0.0;
256  }
257  }
258 
259  //Storage for Lagrangian coordinates (initialised to zero)
260  Vector<double> interpolated_xi(DIM,0.0);
261 
262  //Calculate displacements and derivatives and lagrangian coordinates
263  for(unsigned l=0;l<n_node;l++)
264  {
265  //Loop over positional dofs
266  for(unsigned k=0;k<n_position_type;k++)
267  {
268  double psi_ = psi(l,k);
269  //Loop over displacement components (deformed position)
270  for(unsigned i=0;i<DIM;i++)
271  {
272  //Calculate the Lagrangian coordinates and the accelerations
273  interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi_;
274 
275  // Only compute accelerations if inertia is switched on
276  if ((lambda_sq>0.0)&&(this->Unsteady))
277  {
278  accel[i] += this->dnodal_position_gen_dt(2,l,k,i)*psi_;
279  }
280 
281  //Loop over derivative directions
282  for(unsigned j=0;j<DIM;j++)
283  {
284  interpolated_G(j,i) +=
285  this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
286  }
287  }
288  }
289  }
290 
291  //Get isotropic growth factor
292  double gamma=1.0;
293  this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);
294 
295 
296  //Get body force at current time
297  Vector<double> b(DIM);
298  this->body_force(interpolated_xi,b);
299 
300  // We use Cartesian coordinates as the reference coordinate
301  // system. In this case the undeformed metric tensor is always
302  // the identity matrix -- stretched by the isotropic growth
303  double diag_entry=pow(gamma,2.0/double(DIM));
304  DenseMatrix<double> g(DIM);
305  for(unsigned i=0;i<DIM;i++)
306  {
307  for(unsigned j=0;j<DIM;j++)
308  {
309  if(i==j) {g(i,j) = diag_entry;}
310  else {g(i,j) = 0.0;}
311  }
312  }
313 
314  //Premultiply the undeformed volume ratio (from the isotropic
315  // growth), the weights and the Jacobian
316  double W = gamma*w*J;
317 
318  //Declare and calculate the deformed metric tensor
319  DenseMatrix<double> G(DIM);
320 
321  //Assign values of G
322  for(unsigned i=0;i<DIM;i++)
323  {
324  //Do upper half of matrix
325  for(unsigned j=i;j<DIM;j++)
326  {
327  //Initialise G(i,j) to zero
328  G(i,j) = 0.0;
329  //Now calculate the dot product
330  for(unsigned k=0;k<DIM;k++)
331  {
332  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
333  }
334  }
335  //Matrix is symmetric so just copy lower half
336  for(unsigned j=0;j<i;j++)
337  {
338  G(i,j) = G(j,i);
339  }
340  }
341 
342  //Now calculate the stress tensor from the constitutive law
343  DenseMatrix<double> sigma(DIM);
344  get_stress(g,G,sigma);
345 
346  // Add pre-stress
347  for(unsigned i=0;i<DIM;i++)
348  {
349  for(unsigned j=0;j<DIM;j++)
350  {
351  sigma(i,j) += this->prestress(i,j,interpolated_xi);
352  }
353  }
354 
355  // Get stress derivative by FD only needed for Jacobian
356  //-----------------------------------------------------
357 
358  // Stress derivative
359  RankFourTensor<double> d_stress_dG(DIM,DIM,DIM,DIM,0.0);
360  // Derivative of metric tensor w.r.t. to nodal coords
361  RankFiveTensor<double> d_G_dX(n_node,n_position_type,DIM,DIM,DIM,0.0);
362 
363  // Get Jacobian too?
364  if (flag==1)
365  {
366  // Derivative of metric tensor w.r.t. to discrete positional dofs
367  // NOTE: Since G is symmetric we only compute the upper triangle
368  // and DO NOT copy the entries across. Subsequent computations
369  // must (and, in fact, do) therefore only operate with upper
370  // triangular entries
371  for (unsigned ll=0;ll<n_node;ll++)
372  {
373  for (unsigned kk=0;kk<n_position_type;kk++)
374  {
375  for (unsigned ii=0;ii<DIM;ii++)
376  {
377  for (unsigned aa=0;aa<DIM;aa++)
378  {
379  for (unsigned bb=aa;bb<DIM;bb++)
380  {
381  d_G_dX(ll,kk,ii,aa,bb)=
382  interpolated_G(aa,ii)*dpsidxi(ll,kk,bb)+
383  interpolated_G(bb,ii)*dpsidxi(ll,kk,aa);
384  }
385  }
386  }
387  }
388  }
389 
390  //Get the "upper triangular" entries of the derivatives of the stress
391  //tensor with respect to G
392  this->get_d_stress_dG_upper(g,G,sigma,d_stress_dG);
393  }
394 
395 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
396 
397  //Loop over the test functions, nodes of the element
398  for(unsigned l=0;l<n_node;l++)
399  {
400  //Loop of types of dofs
401  for(unsigned k=0;k<n_position_type;k++)
402  {
403  // Offset for faster access
404  const unsigned offset5=dpsidxi.offset(l ,k);
405 
406  //Loop over the displacement components
407  for(unsigned i=0;i<DIM;i++)
408  {
409  //Get the equation number
410  local_eqn = this->position_local_eqn(l,k,i);
411 
412  /*IF it's not a boundary condition*/
413  if(local_eqn >= 0)
414  {
415  //Initialise contribution to sum
416  double sum=0.0;
417 
418  // Acceleration and body force
419  sum+=(lambda_sq*accel[i]-b[i])*psi(l,k);
420 
421  // Stress term
422  for(unsigned a=0;a<DIM;a++)
423  {
424  unsigned count=offset5;
425  for(unsigned b=0;b<DIM;b++)
426  {
427  //Add the stress terms to the residuals
428  sum+=sigma(a,b)*interpolated_G(a,i)*
429  dpsidxi.raw_direct_access(count);
430  ++count;
431  }
432  }
433  residuals[local_eqn] += W*sum;
434 
435  // Get Jacobian too?
436  if (flag==1)
437  {
438  // Offset for faster access in general stress loop
439  const unsigned offset1=d_G_dX.offset( l, k, i);
440 
441  //Loop over the nodes of the element again
442  for(unsigned ll=0;ll<n_node;ll++)
443  {
444  //Loop of types of dofs again
445  for(unsigned kk=0;kk<n_position_type;kk++)
446  {
447  //Loop over the displacement components again
448  for(unsigned ii=0;ii<DIM;ii++)
449  {
450  //Get the number of the unknown
451  int local_unknown = this->position_local_eqn(ll,kk,ii);
452 
453  /*IF it's not a boundary condition*/
454  if(local_unknown >= 0)
455  {
456  // Offset for faster access in general stress loop
457  const unsigned offset2=d_G_dX.offset( ll, kk, ii);
458  const unsigned offset4=dpsidxi.offset(ll, kk);
459 
460  // General stress term
461  //--------------------
462  double sum=0.0;
463  unsigned count1=offset1;
464  for(unsigned a=0;a<DIM;a++)
465  {
466  // Bump up direct access because we're only
467  // accessing upper triangle
468  count1+=a;
469  for(unsigned b=a;b<DIM;b++)
470  {
471  double factor=d_G_dX.raw_direct_access(count1);
472  if (a==b) factor*=0.5;
473 
474  // Offset for faster access
475  unsigned offset3=d_stress_dG.offset(a,b);
476  unsigned count2=offset2;
477  unsigned count3=offset3;
478 
479  for(unsigned aa=0;aa<DIM;aa++)
480  {
481  // Bump up direct access because we're only
482  // accessing upper triangle
483  count2+=aa;
484  count3+=aa;
485 
486  // Only upper half of derivatives w.r.t. symm tensor
487  for(unsigned bb=aa;bb<DIM;bb++)
488  {
489  sum+=factor*
490  d_stress_dG.raw_direct_access(count3)*
491  d_G_dX.raw_direct_access(count2);
492  ++count2;
493  ++count3;
494  }
495  }
496  ++count1;
497  }
498 
499  }
500 
501  // Multiply by weight and add contribution
502  // (Add directly because this bit is nonsymmetric)
503  jacobian(local_eqn,local_unknown)+=sum*W;
504 
505  // Only upper triangle (no separate test for bc as
506  // local_eqn is already nonnegative)
507  if((i==ii) && (local_unknown >= local_eqn))
508  {
509  //Initialise contribution
510  double sum=0.0;
511 
512  // Inertia term
513  sum+=lambda_sq*time_factor*psi(ll,kk)*psi(l,k);
514 
515  // Stress term
516  unsigned count4=offset4;
517  for(unsigned a=0;a<DIM;a++)
518  {
519  //Cache term
520  const double factor=
521  dpsidxi.raw_direct_access(count4);// ll ,kk
522  ++count4;
523 
524  unsigned count5=offset5;
525  for(unsigned b=0;b<DIM;b++)
526  {
527  sum+=sigma(a,b)*factor*
528  dpsidxi.raw_direct_access(count5); // l ,k
529  ++count5;
530  }
531  }
532  //Add contribution to jacobian
533  jacobian(local_eqn,local_unknown) += sum*W;
534  //Add to lower triangular section
535  if(local_eqn != local_unknown)
536  {
537  jacobian(local_unknown,local_eqn) += sum*W;
538  }
539  }
540 
541  } //End of if not boundary condition
542  }
543  }
544  }
545  }
546 
547  } //End of if not boundary condition
548 
549  } //End of loop over coordinate directions
550  } //End of loop over type of dof
551  } //End of loop over shape functions
552  } //End of loop over integration points
553 }
554 
555 
556 //=======================================================================
557 /// Output: x,y,[z],xi0,xi1,[xi2],gamma
558 //=======================================================================
559 template <unsigned DIM>
560 void PVDEquations<DIM>::output(std::ostream &outfile, const unsigned &n_plot)
561 {
562 
563  Vector<double> x(DIM);
564  Vector<double> xi(DIM);
565  Vector<double> s(DIM);
566 
567  // Tecplot header info
568  outfile << this->tecplot_zone_string(n_plot);
569 
570  // Loop over plot points
571  unsigned num_plot_points=this->nplot_points(n_plot);
572  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
573  {
574  // Get local coordinates of plot point
575  this->get_s_plot(iplot,n_plot,s);
576 
577  // Get Eulerian and Lagrangian coordinates
578  this->interpolated_x(s,x);
579  this->interpolated_xi(s,xi);
580 
581  // Get isotropic growth
582  double gamma;
583  // Dummy integration point
584  unsigned ipt=0;
585  this->get_isotropic_growth(ipt,s,xi,gamma);
586 
587  //Output the x,y,..
588  for(unsigned i=0;i<DIM;i++)
589  {outfile << x[i] << " ";}
590 
591  // Output xi0,xi1,..
592  for(unsigned i=0;i<DIM;i++)
593  {outfile << xi[i] << " ";}
594 
595  // Output growth
596  outfile << gamma;
597  outfile << std::endl;
598  }
599 
600 
601  // Write tecplot footer (e.g. FE connectivity lists)
602  this->write_tecplot_zone_footer(outfile,n_plot);
603  outfile << std::endl;
604 }
605 
606 
607 
608 
609 //=======================================================================
610 /// C-style output: x,y,[z],xi0,xi1,[xi2],gamma
611 //=======================================================================
612 template <unsigned DIM>
613 void PVDEquations<DIM>::output(FILE* file_pt, const unsigned &n_plot)
614  {
615  //Set output Vector
616  Vector<double> s(DIM);
617  Vector<double> x(DIM);
618  Vector<double> xi(DIM);
619 
620  switch(DIM)
621  {
622 
623  case 2:
624 
625  //Tecplot header info
626  //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
627  fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
628 
629  //Loop over element nodes
630  for(unsigned l2=0;l2<n_plot;l2++)
631  {
632  s[1] = -1.0 + l2*2.0/(n_plot-1);
633  for(unsigned l1=0;l1<n_plot;l1++)
634  {
635  s[0] = -1.0 + l1*2.0/(n_plot-1);
636 
637  // Get Eulerian and Lagrangian coordinates
638  this->interpolated_x(s,x);
639  this->interpolated_xi(s,xi);
640 
641  // Get isotropic growth
642  double gamma;
643  // Dummy integration point
644  unsigned ipt=0;
645  this->get_isotropic_growth(ipt,s,xi,gamma);
646 
647  //Output the x,y,..
648  for(unsigned i=0;i<DIM;i++)
649  {
650  //outfile << x[i] << " ";
651  fprintf(file_pt,"%g ",x[i]);
652  }
653  // Output xi0,xi1,..
654  for(unsigned i=0;i<DIM;i++)
655  {
656  //outfile << xi[i] << " ";
657  fprintf(file_pt,"%g ",xi[i]);
658  }
659  // Output growth
660  //outfile << gamma << " ";
661  //outfile << std::endl;
662  fprintf(file_pt,"%g \n",gamma);
663  }
664  }
665  //outfile << std::endl;
666  fprintf(file_pt,"\n");
667 
668  break;
669 
670  case 3:
671 
672  //Tecplot header info
673  //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
674  fprintf(file_pt,"ZONE I=%i, J=%i, K=%i \n",n_plot,n_plot,n_plot);
675 
676  //Loop over element nodes
677  for(unsigned l3=0;l3<n_plot;l3++)
678  {
679  s[2] = -1.0 + l3*2.0/(n_plot-1);
680  for(unsigned l2=0;l2<n_plot;l2++)
681  {
682  s[1] = -1.0 + l2*2.0/(n_plot-1);
683  for(unsigned l1=0;l1<n_plot;l1++)
684  {
685  s[0] = -1.0 + l1*2.0/(n_plot-1);
686 
687  // Get Eulerian and Lagrangian coordinates
688  this->interpolated_x(s,x);
689  this->interpolated_xi(s,xi);
690 
691  // Get isotropic growth
692  double gamma;
693  // Dummy integration point
694  unsigned ipt=0;
695  this->get_isotropic_growth(ipt,s,xi,gamma);
696 
697  //Output the x,y,z
698  for(unsigned i=0;i<DIM;i++)
699  {
700  //outfile << x[i] << " ";
701  fprintf(file_pt,"%g ",x[i]);
702  }
703  // Output xi0,xi1,xi2
704  for(unsigned i=0;i<DIM;i++)
705  {
706  //outfile << xi[i] << " ";
707  fprintf(file_pt,"%g ",xi[i]);
708  }
709  // Output growth
710  //outfile << gamma << " ";
711  //outfile << std::endl;
712  fprintf(file_pt,"%g \n",gamma);
713  }
714  }
715  }
716  //outfile << std::endl;
717  fprintf(file_pt,"\n");
718 
719  break;
720 
721  default:
722  std::ostringstream error_message;
723  error_message << "No output routine for PVDEquations<" <<
724  DIM << "> elements -- write it yourself!" << std::endl;
725  throw OomphLibError(error_message.str(),
726  OOMPH_CURRENT_FUNCTION,
727  OOMPH_EXCEPTION_LOCATION);
728  }
729  }
730 
731 
732 //=======================================================================
733 /// Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components
734 //=======================================================================
735 template <unsigned DIM>
736 void PVDEquations<DIM>::extended_output(std::ostream &outfile,
737  const unsigned &n_plot)
738 {
739 
740  Vector<double> x(DIM);
741  Vector<double> xi(DIM);
742  Vector<double> s(DIM);
743  DenseMatrix<double> stress_or_strain(DIM,DIM);
744 
745  // Tecplot header info
746  outfile << this->tecplot_zone_string(n_plot);
747 
748  // Loop over plot points
749  unsigned num_plot_points=this->nplot_points(n_plot);
750  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
751  {
752  // Get local coordinates of plot point
753  this->get_s_plot(iplot,n_plot,s);
754 
755  // Get Eulerian and Lagrangian coordinates
756  this->interpolated_x(s,x);
757  this->interpolated_xi(s,xi);
758 
759  // Get isotropic growth
760  double gamma;
761  // Dummy integration point
762  unsigned ipt=0;
763  this->get_isotropic_growth(ipt,s,xi,gamma);
764 
765  //Output the x,y,..
766  for(unsigned i=0;i<DIM;i++)
767  {outfile << x[i] << " ";}
768 
769  // Output xi0,xi1,..
770  for(unsigned i=0;i<DIM;i++)
771  {outfile << xi[i] << " ";}
772 
773  // Output growth
774  outfile << gamma << " ";
775 
776  //get the strain
777  this->get_strain(s,stress_or_strain);
778  for(unsigned i=0;i<DIM;i++)
779  {
780  for(unsigned j=0;j<=i;j++)
781  {
782  outfile << stress_or_strain(j,i) << " " ;
783  }
784  }
785 
786  //get the stress
787  this->get_stress(s,stress_or_strain);
788  for(unsigned i=0;i<DIM;i++)
789  {
790  for(unsigned j=0;j<=i;j++)
791  {
792  outfile << stress_or_strain(j,i) << " " ;
793  }
794  }
795 
796 
797  outfile << std::endl;
798  }
799 
800 
801  // Write tecplot footer (e.g. FE connectivity lists)
802  this->write_tecplot_zone_footer(outfile,n_plot);
803  outfile << std::endl;
804 }
805 
806 
807 //=======================================================================
808 /// Get potential (strain) and kinetic energy
809 //=======================================================================
810 template <unsigned DIM>
811 void PVDEquationsBase<DIM>::get_energy(double &pot_en, double &kin_en)
812 {
813  // Initialise
814  pot_en=0;
815  kin_en=0;
816 
817  //Set the value of n_intpt
818  unsigned n_intpt = this->integral_pt()->nweight();
819 
820  //Set the Vector to hold local coordinates
821  Vector<double> s(DIM);
822 
823  //Find out how many nodes there are
824  const unsigned n_node = this->nnode();
825 
826  //Find out how many positional dofs there are
827  const unsigned n_position_type = this->nnodal_position_type();
828 
829  //Set up memory for the shape functions
830  Shape psi(n_node,n_position_type);
831  DShape dpsidxi(n_node,n_position_type,DIM);
832 
833  // Timescale ratio (non-dim density)
834  double lambda_sq = this->lambda_sq();
835 
836  //Loop over the integration points
837  for(unsigned ipt=0;ipt<n_intpt;ipt++)
838  {
839  //Assign values of s
840  for(unsigned i=0;i<DIM;i++) { s[i] = this->integral_pt()->knot(ipt,i); }
841 
842  //Get the integral weight
843  double w = this->integral_pt()->weight(ipt);
844 
845  // Call the derivatives of the shape functions and get Jacobian
846  double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
847 
848  //Storage for Lagrangian coordinates and velocity (initialised to zero)
849  Vector<double> interpolated_xi(DIM,0.0);
850  Vector<double> veloc(DIM,0.0);
851 
852  //Calculate lagrangian coordinates
853  for(unsigned l=0;l<n_node;l++)
854  {
855  //Loop over positional dofs
856  for(unsigned k=0;k<n_position_type;k++)
857  {
858  //Loop over displacement components (deformed position)
859  for(unsigned i=0;i<DIM;i++)
860  {
861  //Calculate the Lagrangian coordinates
862  interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);
863 
864  //Calculate the velocity components (if unsteady solve)
865  if (this->Unsteady)
866  {
867  veloc[i] += this->dnodal_position_gen_dt(l,k,i)*psi(l,k);
868  }
869  }
870  }
871  }
872 
873  //Get isotropic growth factor
874  double gamma=1.0;
875  this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);
876 
877  //Premultiply the undeformed volume ratio (from the isotropic
878  // growth), the weights and the Jacobian
879  double W = gamma*w*J;
880 
881  DenseMatrix<double> sigma(DIM,DIM);
882  DenseMatrix<double> strain(DIM,DIM);
883 
884  //Now calculate the stress tensor from the constitutive law
885  this->get_stress(s,sigma);
886 
887  // Add pre-stress
888  for(unsigned i=0;i<DIM;i++)
889  {
890  for(unsigned j=0;j<DIM;j++)
891  {
892  sigma(i,j) += prestress(i,j,interpolated_xi);
893  }
894  }
895 
896  //get the strain
897  this->get_strain(s,strain);
898 
899  // Initialise
900  double local_pot_en=0;
901  double veloc_sq=0;
902 
903  // Compute integrals
904  for(unsigned i=0;i<DIM;i++)
905  {
906  for(unsigned j=0;j<DIM;j++)
907  {
908  local_pot_en+=sigma(i,j)*strain(i,j);
909  }
910  veloc_sq += veloc[i]*veloc[i];
911  }
912 
913  pot_en+=0.5*local_pot_en*W;
914  kin_en+=lambda_sq*0.5*veloc_sq*W;
915  }
916 }
917 
918 
919 
920 //=======================================================================
921 /// Compute the contravariant second Piola Kirchoff stress at a given local
922 /// coordinate. Note: this replicates a lot of code that is already
923 /// coontained in get_residuals() but without sacrificing efficiency
924 /// (re-computing the shape functions several times) or creating
925 /// helper functions with horrendous interfaces (to pass all the
926 /// functions which shouldn't be recomputed) about this is
927 /// unavoidable.
928 //=======================================================================
929 template <unsigned DIM>
931  DenseMatrix<double> &sigma)
932 {
933  //Find out how many nodes there are
934  unsigned n_node = this->nnode();
935 
936  //Find out how many positional dofs there are
937  unsigned n_position_type = this->nnodal_position_type();
938 
939  //Set up memory for the shape functions
940  Shape psi(n_node,n_position_type);
941  DShape dpsidxi(n_node,n_position_type,DIM);
942 
943  //Call the derivatives of the shape functions (ignore Jacobian)
944  (void) this->dshape_lagrangian(s,psi,dpsidxi);
945 
946  // Lagrangian coordinates
947  Vector<double> xi(DIM);
948  this->interpolated_xi(s,xi);
949 
950  //Get isotropic growth factor
951  double gamma;
952  //Dummy integration point
953  unsigned ipt=0;
954  this->get_isotropic_growth(ipt,s,xi,gamma);
955 
956  // We use Cartesian coordinates as the reference coordinate
957  // system. In this case the undeformed metric tensor is always
958  // the identity matrix -- stretched by the isotropic growth
959  double diag_entry=pow(gamma,2.0/double(DIM));
960  DenseMatrix<double> g(DIM);
961  for(unsigned i=0;i<DIM;i++)
962  {
963  for(unsigned j=0;j<DIM;j++)
964  {
965  if(i==j) {g(i,j) = diag_entry;}
966  else {g(i,j) = 0.0;}
967  }
968  }
969 
970 
971  //Calculate interpolated values of the derivative of global position
972  //wrt lagrangian coordinates
973  DenseMatrix<double> interpolated_G(DIM);
974 
975  //Initialise to zero
976  for(unsigned i=0;i<DIM;i++)
977  {for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}}
978 
979  //Calculate displacements and derivatives
980  for(unsigned l=0;l<n_node;l++)
981  {
982  //Loop over positional dofs
983  for(unsigned k=0;k<n_position_type;k++)
984  {
985  //Loop over displacement components (deformed position)
986  for(unsigned i=0;i<DIM;i++)
987  {
988  //Loop over derivative directions
989  for(unsigned j=0;j<DIM;j++)
990  {
991  interpolated_G(j,i) +=
992  this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
993  }
994  }
995  }
996  }
997 
998  //Declare and calculate the deformed metric tensor
999  DenseMatrix<double> G(DIM);
1000  //Assign values of G
1001  for(unsigned i=0;i<DIM;i++)
1002  {
1003  //Do upper half of matrix
1004  //Note that j must be signed here for the comparison test to work
1005  //Also i must be cast to an int
1006  for(int j=(DIM-1);j>=static_cast<int>(i);j--)
1007  {
1008  //Initialise G(i,j) to zero
1009  G(i,j) = 0.0;
1010  //Now calculate the dot product
1011  for(unsigned k=0;k<DIM;k++)
1012  {
1013  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
1014  }
1015  }
1016  //Matrix is symmetric so just copy lower half
1017  for(int j=(i-1);j>=0;j--)
1018  {
1019  G(i,j) = G(j,i);
1020  }
1021  }
1022 
1023  //Now calculate the stress tensor from the constitutive law
1024  get_stress(g,G,sigma);
1025 
1026 }
1027 
1028 
1029 
1030 //////////////////////////////////////////////////////////////////
1031 //////////////////////////////////////////////////////////////////
1032 //////////////////////////////////////////////////////////////////
1033 
1034 
1035 
1036 //=======================================================================
1037 /// Compute principal stress vectors and (scalar) principal stresses
1038 /// at specified local coordinate: \c principal_stress_vector(i,j)
1039 /// is the j-th component of the i-th principal stress vector.
1040 //=======================================================================
1041 template <unsigned DIM>
1043  const Vector<double> &s, DenseMatrix<double>& principal_stress_vector,
1044  Vector<double>& principal_stress)
1045 {
1046  // Compute contravariant ("upper") 2nd Piola Kirchhoff stress
1047  DenseDoubleMatrix sigma(DIM,DIM);
1048  get_stress(s,sigma);
1049 
1050  // Lagrangian coordinates
1051  Vector<double> xi(DIM);
1052  this->interpolated_xi(s,xi);
1053 
1054  // Add pre-stress
1055  for(unsigned i=0;i<DIM;i++)
1056  {
1057  for(unsigned j=0;j<DIM;j++)
1058  {
1059  sigma(i,j) += this->prestress(i,j,xi);
1060  }
1061  }
1062 
1063  // Get covariant base vectors in deformed configuration
1064  DenseMatrix<double> lower_deformed_basis(DIM);
1065  get_deformed_covariant_basis_vectors(s,lower_deformed_basis);
1066 
1067  // Work out covariant ("lower") metric tensor
1068  DenseDoubleMatrix lower_metric(DIM);
1069  for (unsigned i=0;i<DIM;i++)
1070  {
1071  for (unsigned j=0;j<DIM;j++)
1072  {
1073  lower_metric(i,j)=0.0;
1074  for (unsigned k=0;k<DIM;k++)
1075  {
1076  lower_metric(i,j)+=
1077  lower_deformed_basis(i,k)*lower_deformed_basis(j,k);
1078  }
1079  }
1080  }
1081 
1082  // Work out cartesian components of contravariant ("upper") basis vectors
1083  DenseMatrix<double> upper_deformed_basis(DIM);
1084 
1085  // Loop over RHSs
1086  Vector<double> rhs(DIM);
1087  Vector<double> aux(DIM);
1088  for (unsigned k=0;k<DIM;k++)
1089  {
1090 
1091  for (unsigned l=0;l<DIM;l++)
1092  {
1093  rhs[l] = lower_deformed_basis(l,k);
1094  }
1095 
1096  lower_metric.solve(rhs,aux);
1097 
1098  for (unsigned l=0;l<DIM;l++)
1099  {
1100  upper_deformed_basis(l,k) = aux[l];
1101  }
1102  }
1103 
1104  // Eigenvalues (=principal stresses) and eigenvectors
1105  DenseMatrix<double> ev(DIM);
1106 
1107  // Get eigenvectors of contravariant 2nd Piola Kirchoff stress
1108  sigma.eigenvalues_by_jacobi(principal_stress,ev);
1109 
1110  // ev(j,i) is the i-th component of the j-th eigenvector
1111  // relative to the deformed "lower variance" basis!
1112  // Work out cartesian components of eigenvectors by multiplying
1113  // the "lower variance components" by these "upper variance" basis
1114  // vectors
1115 
1116  // Loop over cartesian compnents
1117  for(unsigned i=0;i<DIM;i++)
1118  {
1119  // Initialise the row
1120  for(unsigned j=0;j<DIM;j++) {principal_stress_vector(j,i)=0.0;}
1121 
1122  // Loop over basis vectors
1123  for(unsigned j=0;j<DIM;j++)
1124  {
1125  for(unsigned k=0;k<DIM;k++)
1126  {
1127  principal_stress_vector(j,i) += upper_deformed_basis(k,i)*ev(j,k);
1128  }
1129  }
1130  }
1131 
1132  // Scaling factor to turn these vectors into unit vectors
1133  Vector<double> norm(DIM);
1134  for (unsigned i=0;i<DIM;i++)
1135  {
1136  norm[i]=0.0;
1137  for (unsigned j=0;j<DIM;j++)
1138  {
1139  norm[i] += pow(principal_stress_vector(i,j),2);
1140  }
1141  norm[i] = sqrt(norm[i]);
1142  }
1143 
1144 
1145  // Scaling and then multiplying by eigenvalue gives the principal stress
1146  // vectors
1147  for(unsigned i=0;i<DIM;i++)
1148  {
1149  for (unsigned j=0;j<DIM;j++)
1150  {
1151  principal_stress_vector(j,i) = ev(j,i)/norm[j]*principal_stress[j];
1152  }
1153  }
1154 
1155 }
1156 
1157 
1158 
1159 //=======================================================================
1160 /// Return the deformed covariant basis vectors
1161 /// at specified local coordinate: \c def_covariant_basis(i,j)
1162 /// is the j-th component of the i-th basis vector.
1163 //=======================================================================
1164 template <unsigned DIM>
1166  const Vector<double> &s, DenseMatrix<double>& def_covariant_basis)
1167 {
1168 
1169  //Find out how many nodes there are
1170  unsigned n_node = nnode();
1171 
1172  //Find out how many positional dofs there are
1173  unsigned n_position_type = this->nnodal_position_type();
1174 
1175  //Set up memory for the shape functions
1176  Shape psi(n_node,n_position_type);
1177  DShape dpsidxi(n_node,n_position_type,DIM);
1178 
1179 
1180  //Call the derivatives of the shape functions (ignore Jacobian)
1181  (void) dshape_lagrangian(s,psi,dpsidxi);
1182 
1183 
1184  //Initialise to zero
1185  for(unsigned i=0;i<DIM;i++)
1186  {for(unsigned j=0;j<DIM;j++) {def_covariant_basis(i,j) = 0.0;}}
1187 
1188  //Calculate displacements and derivatives
1189  for(unsigned l=0;l<n_node;l++)
1190  {
1191  //Loop over positional dofs
1192  for(unsigned k=0;k<n_position_type;k++)
1193  {
1194  //Loop over displacement components (deformed position)
1195  for(unsigned i=0;i<DIM;i++)
1196  {
1197  //Loop over derivative directions (i.e. base vectors)
1198  for(unsigned j=0;j<DIM;j++)
1199  {
1200  def_covariant_basis(j,i) +=
1201  nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
1202  }
1203  }
1204  }
1205  }
1206 }
1207 
1208 
1209 //////////////////////////////////////////////////////////////////////
1210 //////////////////////////////////////////////////////////////////////
1211 //////////////////////////////////////////////////////////////////////
1212 
1213 //=====================================================================
1214 /// "Magic" number that indicates that the solid pressure is not stored
1215 /// at a node. It is a negative number that cannot be -1 because that is
1216 /// used to represent the positional hanging scheme in Hanging_pt objects
1217 //======================================================================
1218 template<unsigned DIM>
1220 
1221 
1222 
1223 
1224 
1225 
1226 //=======================================================================
1227 /// Fill in element's contribution to the elemental
1228 /// residual vector and/or Jacobian matrix.
1229 /// flag=0: compute only residual vector
1230 /// flag=1: compute both, fully analytically
1231 /// flag=2: compute both, using FD for the derivatives w.r.t. to the
1232 /// discrete displacment dofs.
1233 /// flag=3: compute residuals, jacobian (full analytic) and mass matrix
1234 /// flag=4: compute residuals, jacobian (FD for derivatives w.r.t.
1235 /// displacements) and mass matrix
1236 //=======================================================================
1237 template <unsigned DIM>
1240  Vector<double> &residuals,DenseMatrix<double> &jacobian,
1241  DenseMatrix<double> &mass_matrix, const unsigned& flag)
1242 {
1243 
1244 #ifdef PARANOID
1245  // Check if the constitutive equation requires the explicit imposition of an
1246  // incompressibility constraint
1248  (!Incompressible))
1249  {
1250  throw OomphLibError(
1251  "The constitutive law requires the use of the incompressible formulation by setting the element's member function set_incompressible()",
1252  OOMPH_CURRENT_FUNCTION,
1253  OOMPH_EXCEPTION_LOCATION);
1254  }
1255 #endif
1256 
1257 
1258  // Simply set up initial condition?
1259  if (this->Solid_ic_pt!=0)
1260  {
1261  this->get_residuals_for_solid_ic(residuals);
1262  return;
1263  }
1264 
1265  //Find out how many nodes there are
1266  const unsigned n_node = this->nnode();
1267 
1268  //Find out how many position types of dof there are
1269  const unsigned n_position_type = this->nnodal_position_type();
1270 
1271  //Find out how many pressure dofs there are
1272  const unsigned n_solid_pres = this->npres_solid();
1273 
1274  //Set up memory for the shape functions
1275  Shape psi(n_node,n_position_type);
1276  DShape dpsidxi(n_node,n_position_type,DIM);
1277 
1278  //Set up memory for the pressure shape functions
1279  Shape psisp(n_solid_pres);
1280 
1281  //Set the value of n_intpt
1282  const unsigned n_intpt = this->integral_pt()->nweight();
1283 
1284  //Set the vector to hold the local coordinates in the element
1285  Vector<double> s(DIM);
1286 
1287  // Timescale ratio (non-dim density)
1288  double lambda_sq = this->lambda_sq();
1289 
1290  // Time factor
1291  double time_factor=0.0;
1292  if (lambda_sq>0)
1293  {
1294  time_factor=this->node_pt(0)->position_time_stepper_pt()->weight(2,0);
1295  }
1296 
1297  //Integers to hold the local equation and unknown numbers
1298  int local_eqn=0, local_unknown=0;
1299 
1300  //Loop over the integration points
1301  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1302  {
1303  //Assign the values of s
1304  for(unsigned i=0;i<DIM;++i)
1305  {
1306  s[i] = this->integral_pt()->knot(ipt,i);
1307  }
1308 
1309  //Get the integral weight
1310  double w = this->integral_pt()->weight(ipt);
1311 
1312  //Call the derivatives of the shape functions
1313  double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
1314 
1315  //Call the pressure shape functions
1316  solid_pshape_at_knot(ipt,psisp);
1317 
1318  //Storage for Lagrangian coordinates (initialised to zero)
1319  Vector<double> interpolated_xi(DIM,0.0);
1320 
1321  // Deformed tangent vectors
1322  DenseMatrix<double> interpolated_G(DIM);
1323 
1324  // Setup memory for accelerations
1325  Vector<double> accel(DIM);
1326 
1327  //Initialise to zero
1328  for(unsigned i=0;i<DIM;i++)
1329  {
1330  // Initialise acclerations
1331  accel[i]=0.0;
1332  for(unsigned j=0;j<DIM;j++)
1333  {
1334  interpolated_G(i,j) = 0.0;
1335  }
1336  }
1337 
1338  //Calculate displacements and derivatives and lagrangian coordinates
1339  for(unsigned l=0;l<n_node;l++)
1340  {
1341  //Loop over positional dofs
1342  for(unsigned k=0;k<n_position_type;k++)
1343  {
1344  //Loop over displacement components (deformed position)
1345  for(unsigned i=0;i<DIM;i++)
1346  {
1347  //Calculate the lagrangian coordinates and the accelerations
1348  interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi(l,k);
1349 
1350  // Only compute accelerations if inertia is switched on
1351  // otherwise the timestepper might not be able to
1352  // work out dx_gen_dt(2,...)
1353  if ((lambda_sq>0.0)&&(this->Unsteady))
1354  {
1355  accel[i] += this->dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
1356  }
1357 
1358  //Loop over derivative directions
1359  for(unsigned j=0;j<DIM;j++)
1360  {
1361  interpolated_G(j,i) +=
1362  this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
1363  }
1364  }
1365  }
1366  }
1367 
1368  //Get isotropic growth factor
1369  double gamma=1.0;
1370  this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);
1371 
1372  //Get body force at current time
1373  Vector<double> b(DIM);
1374  this->body_force(interpolated_xi,b);
1375 
1376  // We use Cartesian coordinates as the reference coordinate
1377  // system. In this case the undeformed metric tensor is always
1378  // the identity matrix -- stretched by the isotropic growth
1379  double diag_entry=pow(gamma,2.0/double(DIM));
1380  DenseMatrix<double> g(DIM);
1381  for(unsigned i=0;i<DIM;i++)
1382  {
1383  for(unsigned j=0;j<DIM;j++)
1384  {
1385  if(i==j) {g(i,j) = diag_entry;}
1386  else {g(i,j) = 0.0;}
1387  }
1388  }
1389 
1390  //Premultiply the undeformed volume ratio (from the isotropic
1391  // growth), the weights and the Jacobian
1392  double W = gamma*w*J;
1393 
1394  //Calculate the interpolated solid pressure
1395  double interpolated_solid_p=0.0;
1396  for(unsigned l=0;l<n_solid_pres;l++)
1397  {
1398  interpolated_solid_p += solid_p(l)*psisp[l];
1399  }
1400 
1401 
1402  //Declare and calculate the deformed metric tensor
1403  DenseMatrix<double> G(DIM);
1404 
1405  //Assign values of G
1406  for(unsigned i=0;i<DIM;i++)
1407  {
1408  //Do upper half of matrix
1409  for(unsigned j=i;j<DIM;j++)
1410  {
1411  //Initialise G(i,j) to zero
1412  G(i,j) = 0.0;
1413  //Now calculate the dot product
1414  for(unsigned k=0;k<DIM;k++)
1415  {
1416  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
1417  }
1418  }
1419  //Matrix is symmetric so just copy lower half
1420  for(unsigned j=0;j<i;j++)
1421  {
1422  G(i,j) = G(j,i);
1423  }
1424  }
1425 
1426  //Now calculate the deviatoric stress and all pressure-related
1427  //quantitites
1428  DenseMatrix<double> sigma(DIM,DIM), sigma_dev(DIM,DIM), Gup(DIM,DIM);
1429  double detG = 0.0;
1430  double gen_dil=0.0;
1431  double inv_kappa=0.0;
1432 
1433  // Get stress derivative by FD only needed for Jacobian
1434 
1435  // Stress etc derivatives
1436  RankFourTensor<double> d_stress_dG(DIM,DIM,DIM,DIM,0.0);
1437  DenseMatrix<double> d_detG_dG(DIM,DIM,0.0);
1438  DenseMatrix<double> d_gen_dil_dG(DIM,DIM,0.0);
1439 
1440  // Derivative of metric tensor w.r.t. to nodal coords
1441  RankFiveTensor<double> d_G_dX(n_node,n_position_type,DIM,DIM,DIM,0.0);
1442 
1443  // Get Jacobian too?
1444  if ((flag==1) || (flag==3))
1445  {
1446  // Derivative of metric tensor w.r.t. to discrete positional dofs
1447  // NOTE: Since G is symmetric we only compute the upper triangle
1448  // and DO NOT copy the entries across. Subsequent computations
1449  // must (and, in fact, do) therefore only operate with upper
1450  // triangular entries
1451  for (unsigned ll=0;ll<n_node;ll++)
1452  {
1453  for (unsigned kk=0;kk<n_position_type;kk++)
1454  {
1455  for (unsigned ii=0;ii<DIM;ii++)
1456  {
1457  for (unsigned aa=0;aa<DIM;aa++)
1458  {
1459  for (unsigned bb=aa;bb<DIM;bb++)
1460  {
1461  d_G_dX(ll,kk,ii,aa,bb)=
1462  interpolated_G(aa,ii)*dpsidxi(ll,kk,bb)+
1463  interpolated_G(bb,ii)*dpsidxi(ll,kk,aa);
1464  }
1465  }
1466  }
1467  }
1468  }
1469  }
1470 
1471 
1472  // Incompressible: Compute the deviatoric part of the stress tensor, the
1473  // contravariant deformed metric tensor and the determinant
1474  // of the deformed covariant metric tensor.
1475  if(Incompressible)
1476  {
1477  get_stress(g,G,sigma_dev,Gup,detG);
1478 
1479  // Get full stress
1480  for (unsigned a=0;a<DIM;a++)
1481  {
1482  for (unsigned b=0;b<DIM;b++)
1483  {
1484  sigma(a,b)=sigma_dev(a,b)-interpolated_solid_p*Gup(a,b);
1485  }
1486  }
1487 
1488  // Get Jacobian too?
1489  if ((flag==1) || (flag==3))
1490  {
1491  //Get the "upper triangular" entries of the derivatives of the stress
1492  //tensor with respect to G
1493  this->
1494  get_d_stress_dG_upper(g,G,sigma,detG,interpolated_solid_p,
1495  d_stress_dG,d_detG_dG);
1496  }
1497  }
1498  // Nearly incompressible: Compute the deviatoric part of the
1499  // stress tensor, the contravariant deformed metric tensor,
1500  // the generalised dilatation and the inverse bulk modulus.
1501  else
1502  {
1503  get_stress(g,G,sigma_dev,Gup,gen_dil,inv_kappa);
1504 
1505  // Get full stress
1506  for (unsigned a=0;a<DIM;a++)
1507  {
1508  for (unsigned b=0;b<DIM;b++)
1509  {
1510  sigma(a,b)=sigma_dev(a,b)-interpolated_solid_p*Gup(a,b);
1511  }
1512  }
1513 
1514  // Get Jacobian too?
1515  if ((flag==1) || (flag==3))
1516  {
1517  //Get the "upper triangular" entries of the derivatives of the stress
1518  //tensor with respect to G
1519  this->get_d_stress_dG_upper(g,G,sigma,gen_dil,inv_kappa,
1520  interpolated_solid_p,
1521  d_stress_dG,d_gen_dil_dG);
1522  }
1523  }
1524 
1525  // Add pre-stress
1526  for(unsigned i=0;i<DIM;i++)
1527  {
1528  for(unsigned j=0;j<DIM;j++)
1529  {
1530  sigma(i,j) += this->prestress(i,j,interpolated_xi);
1531  }
1532  }
1533 
1534 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
1535 
1536  //Loop over the test functions, nodes of the element
1537  for(unsigned l=0;l<n_node;l++)
1538  {
1539  //Loop over the types of dof
1540  for(unsigned k=0;k<n_position_type;k++)
1541  {
1542  // Offset for faster access
1543  const unsigned offset5=dpsidxi.offset(l ,k);
1544 
1545  //Loop over the displacement components
1546  for(unsigned i=0;i<DIM;i++)
1547  {
1548  //Get the equation number
1549  local_eqn = this->position_local_eqn(l,k,i);
1550 
1551  /*IF it's not a boundary condition*/
1552  if(local_eqn >= 0)
1553  {
1554  //Initialise the contribution
1555  double sum=0.0;
1556 
1557  // Acceleration and body force
1558  sum+=(lambda_sq*accel[i]-b[i])*psi(l,k);
1559 
1560  // Stress term
1561  for(unsigned a=0;a<DIM;a++)
1562  {
1563  unsigned count=offset5;
1564  for(unsigned b=0;b<DIM;b++)
1565  {
1566  //Add the stress terms to the residuals
1567  sum+=sigma(a,b)*interpolated_G(a,i)*
1568  dpsidxi.raw_direct_access(count);
1569  ++count;
1570  }
1571  }
1572  residuals[local_eqn] += W*sum;
1573 
1574  //Add in the mass matrix terms
1575  if(flag > 2)
1576  {
1577  //Loop over the nodes of the element again
1578  for(unsigned ll=0;ll<n_node;ll++)
1579  {
1580  //Loop of types of dofs again
1581  for(unsigned kk=0;kk<n_position_type;kk++)
1582  {
1583  //Get the number of the unknown
1584  int local_unknown = this->position_local_eqn(ll,kk,i);
1585 
1586  /*IF it's not a boundary condition*/
1587  if(local_unknown >= 0)
1588  {
1589  mass_matrix(local_eqn,local_unknown) +=
1590  lambda_sq*psi(l,k)*psi(ll,kk)*W;
1591  }
1592  }
1593  }
1594  }
1595 
1596  //Add in the jacobian terms
1597  if((flag==1) || (flag==3))
1598  {
1599  // Offset for faster access in general stress loop
1600  const unsigned offset1=d_G_dX.offset( l, k, i);
1601 
1602  //Loop over the nodes of the element again
1603  for(unsigned ll=0;ll<n_node;ll++)
1604  {
1605  //Loop of types of dofs again
1606  for(unsigned kk=0;kk<n_position_type;kk++)
1607  {
1608  //Loop over the displacement components again
1609  for(unsigned ii=0;ii<DIM;ii++)
1610  {
1611  //Get the number of the unknown
1612  int local_unknown = this->position_local_eqn(ll,kk,ii);
1613 
1614  /*IF it's not a boundary condition*/
1615  if(local_unknown >= 0)
1616  {
1617 
1618  // Offset for faster access in general stress loop
1619  const unsigned offset2=d_G_dX.offset( ll, kk, ii);
1620  const unsigned offset4=dpsidxi.offset(ll, kk);
1621 
1622  // General stress term
1623  //--------------------
1624  double sum=0.0;
1625  unsigned count1=offset1;
1626  for(unsigned a=0;a<DIM;a++)
1627  {
1628  // Bump up direct access because we're only
1629  // accessing upper triangle
1630  count1+=a;
1631  for(unsigned b=a;b<DIM;b++)
1632  {
1633  double factor=d_G_dX.raw_direct_access(count1);
1634  if (a==b) factor*=0.5;
1635 
1636  // Offset for faster access
1637  unsigned offset3=d_stress_dG.offset(a,b);
1638  unsigned count2=offset2;
1639  unsigned count3=offset3;
1640 
1641  for(unsigned aa=0;aa<DIM;aa++)
1642  {
1643  // Bump up direct access because we're only
1644  // accessing upper triangle
1645  count2+=aa;
1646  count3+=aa;
1647 
1648  // Only upper half of derivatives w.r.t. symm tensor
1649  for(unsigned bb=aa;bb<DIM;bb++)
1650  {
1651  sum+=factor*
1652  d_stress_dG.raw_direct_access(count3)*
1653  d_G_dX.raw_direct_access(count2);
1654  ++count2;
1655  ++count3;
1656  }
1657  }
1658  ++count1;
1659  }
1660 
1661  }
1662 
1663  // Multiply by weight and add contribution
1664  // (Add directly because this bit is nonsymmetric)
1665  jacobian(local_eqn,local_unknown)+=sum*W;
1666 
1667  // Only upper triangle (no separate test for bc as
1668  // local_eqn is already nonnegative)
1669  if((i==ii) && (local_unknown >= local_eqn))
1670  {
1671  //Initialise the contribution
1672  double sum=0.0;
1673 
1674  // Inertia term
1675  sum+=lambda_sq*time_factor*psi(ll,kk)*psi(l,k);
1676 
1677  // Stress term
1678  unsigned count4=offset4;
1679  for(unsigned a=0;a<DIM;a++)
1680  {
1681  //Cache term
1682  const double factor=
1683  dpsidxi.raw_direct_access(count4);// ll ,kk
1684  ++count4;
1685 
1686  unsigned count5=offset5;
1687  for(unsigned b=0;b<DIM;b++)
1688  {
1689  sum+=sigma(a,b)*factor*
1690  dpsidxi.raw_direct_access(count5); // l ,k
1691  ++count5;
1692  }
1693  }
1694 
1695  //Add to jacobian
1696  jacobian(local_eqn,local_unknown) += sum*W;
1697  //Add to lower triangular parts
1698  if(local_eqn!=local_unknown)
1699  {
1700  jacobian(local_unknown,local_eqn) += sum*W;
1701  }
1702  }
1703  } //End of if not boundary condition
1704  }
1705  }
1706  }
1707  }
1708 
1709  // Derivatives w.r.t. pressure dofs
1710  if (flag>0)
1711  {
1712  //Loop over the pressure dofs for unknowns
1713  for(unsigned l2=0;l2<n_solid_pres;l2++)
1714  {
1715  local_unknown = this->solid_p_local_eqn(l2);
1716 
1717  //If it's not a boundary condition
1718  if(local_unknown >= 0)
1719  {
1720 
1721  //Add the pressure terms to the jacobian
1722  for(unsigned a=0;a<DIM;a++)
1723  {
1724  for(unsigned b=0;b<DIM;b++)
1725  {
1726  jacobian(local_eqn,local_unknown) -=
1727  psisp[l2]*Gup(a,b)*
1728  interpolated_G(a,i)*dpsidxi(l,k,b)*W;
1729  }
1730  }
1731  }
1732  }
1733  } //End of Jacobian terms
1734 
1735  } //End of if not boundary condition
1736  } //End of loop over coordinate directions
1737  } //End of loop over types of dof
1738  } //End of loop over shape functions
1739 
1740  //==============CONSTRAINT EQUATIONS FOR PRESSURE=====================
1741 
1742  //Now loop over the pressure degrees of freedom
1743  for(unsigned l=0;l<n_solid_pres;l++)
1744  {
1745  local_eqn = this->solid_p_local_eqn(l);
1746 
1747  // Pinned (unlikely, actually) or real dof?
1748  if(local_eqn >= 0)
1749  {
1750  //For true incompressibility we need to conserve volume
1751  //so the determinant of the deformed metric tensor
1752  //needs to be equal to that of the undeformed one, which
1753  //is equal to the volumetric growth factor
1754  if(Incompressible)
1755  {
1756  residuals[local_eqn] += (detG - gamma)*psisp[l]*W;
1757 
1758 
1759  // Get Jacobian too?
1760  if ((flag==1) || (flag==3))
1761  {
1762  //Loop over the nodes of the element again
1763  for(unsigned ll=0;ll<n_node;ll++)
1764  {
1765  //Loop of types of dofs again
1766  for(unsigned kk=0;kk<n_position_type;kk++)
1767  {
1768  //Loop over the displacement components again
1769  for(unsigned ii=0;ii<DIM;ii++)
1770  {
1771  //Get the number of the unknown
1772  int local_unknown = this->position_local_eqn(ll,kk,ii);
1773 
1774  /*IF it's not a boundary condition*/
1775  if(local_unknown >= 0)
1776  {
1777 
1778  // Offset for faster access
1779  const unsigned offset=d_G_dX.offset( ll, kk, ii);
1780 
1781  // General stress term
1782  double sum=0.0;
1783  unsigned count=offset;
1784  for(unsigned aa=0;aa<DIM;aa++)
1785  {
1786  // Bump up direct access because we're only
1787  // accessing upper triangle
1788  count+=aa;
1789 
1790  // Only upper half
1791  for(unsigned bb=aa;bb<DIM;bb++)
1792  {
1793  sum+=d_detG_dG(aa,bb)*
1794  d_G_dX.raw_direct_access(count)*psisp(l);
1795  ++count;
1796  }
1797  }
1798  jacobian(local_eqn,local_unknown)+=sum*W;
1799  }
1800  }
1801  }
1802  }
1803 
1804  //No Jacobian terms due to pressure since it does not feature
1805  //in the incompressibility constraint
1806  }
1807  }
1808  //Nearly incompressible: (Neg.) pressure given by product of
1809  //bulk modulus and generalised dilatation
1810  else
1811  {
1812  residuals[local_eqn] +=
1813  (inv_kappa*interpolated_solid_p + gen_dil)*psisp[l]*W;
1814 
1815  //Add in the jacobian terms
1816  if ((flag==1) || (flag==3))
1817  {
1818  //Loop over the nodes of the element again
1819  for(unsigned ll=0;ll<n_node;ll++)
1820  {
1821  //Loop of types of dofs again
1822  for(unsigned kk=0;kk<n_position_type;kk++)
1823  {
1824  //Loop over the displacement components again
1825  for(unsigned ii=0;ii<DIM;ii++)
1826  {
1827  //Get the number of the unknown
1828  int local_unknown = this->position_local_eqn(ll,kk,ii);
1829 
1830  /*IF it's not a boundary condition*/
1831  if(local_unknown >= 0)
1832  {
1833  // Offset for faster access
1834  const unsigned offset=d_G_dX.offset( ll, kk, ii);
1835 
1836  // General stress term
1837  double sum=0.0;
1838  unsigned count=offset;
1839  for(unsigned aa=0;aa<DIM;aa++)
1840  {
1841  // Bump up direct access because we're only
1842  // accessing upper triangle
1843  count+=aa;
1844 
1845  // Only upper half
1846  for(unsigned bb=aa;bb<DIM;bb++)
1847  {
1848  sum+=d_gen_dil_dG(aa,bb)*
1849  d_G_dX.raw_direct_access(count)*psisp(l);
1850  ++count;
1851  }
1852  }
1853  jacobian(local_eqn,local_unknown)+=sum*W;
1854  }
1855  }
1856  }
1857  }
1858  }
1859 
1860  // Derivatives w.r.t. pressure dofs
1861  if(flag>0)
1862  {
1863  //Loop over the pressure nodes again
1864  for(unsigned l2=0;l2<n_solid_pres;l2++)
1865  {
1866  local_unknown = this->solid_p_local_eqn(l2);
1867  //If not pinnned
1868  if(local_unknown >= 0)
1869  {
1870  jacobian(local_eqn,local_unknown)
1871  += inv_kappa*psisp[l2]*psisp[l]*W;
1872  }
1873  }
1874  } //End of jacobian terms
1875 
1876  } //End of else
1877 
1878  } //End of if not boundary condition
1879  }
1880 
1881  } //End of loop over integration points
1882 }
1883 
1884 
1885 
1886 //=======================================================================
1887 /// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
1888 //=======================================================================
1889 template <unsigned DIM>
1890 void PVDEquationsWithPressure<DIM>::output(std::ostream &outfile,
1891  const unsigned &n_plot)
1892 {
1893  //Set output Vector
1894  Vector<double> s(DIM);
1895  Vector<double> x(DIM);
1896  Vector<double> xi(DIM);
1897 
1898  // Tecplot header info
1899  outfile << this->tecplot_zone_string(n_plot);
1900 
1901  // Loop over plot points
1902  unsigned num_plot_points=this->nplot_points(n_plot);
1903  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1904  {
1905  // Get local coordinates of plot point
1906  this->get_s_plot(iplot,n_plot,s);
1907 
1908  // Get Eulerian and Lagrangian coordinates
1909  this->interpolated_x(s,x);
1910  this->interpolated_xi(s,xi);
1911 
1912  // Get isotropic growth
1913  double gamma;
1914  // Dummy integration point
1915  unsigned ipt=0;
1916  this->get_isotropic_growth(ipt,s,xi,gamma);
1917 
1918  //Output the x,y,..
1919  for(unsigned i=0;i<DIM;i++)
1920  {outfile << x[i] << " ";}
1921 
1922  // Output xi0,xi1,..
1923  for(unsigned i=0;i<DIM;i++)
1924  {outfile << xi[i] << " ";}
1925 
1926  // Output growth
1927  outfile << gamma << " ";
1928  // Output pressure
1929  outfile << interpolated_solid_p(s) << " ";
1930  outfile << "\n";
1931  }
1932 
1933  // Write tecplot footer (e.g. FE connectivity lists)
1934  this->write_tecplot_zone_footer(outfile,n_plot);
1935 }
1936 
1937 
1938 
1939 
1940 //=======================================================================
1941 /// C-stsyle output: x,y,[z],xi0,xi1,[xi2],p,gamma
1942 //=======================================================================
1943 template <unsigned DIM>
1945  const unsigned &n_plot)
1946 {
1947  //Set output Vector
1948  Vector<double> s(DIM);
1949  Vector<double> x(DIM);
1950  Vector<double> xi(DIM);
1951 
1952  switch(DIM)
1953  {
1954  case 2:
1955  //Tecplot header info
1956  //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
1957  fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
1958 
1959  //Loop over element nodes
1960  for(unsigned l2=0;l2<n_plot;l2++)
1961  {
1962  s[1] = -1.0 + l2*2.0/(n_plot-1);
1963  for(unsigned l1=0;l1<n_plot;l1++)
1964  {
1965  s[0] = -1.0 + l1*2.0/(n_plot-1);
1966 
1967  // Get Eulerian and Lagrangian coordinates
1968  this->interpolated_x(s,x);
1969  this->interpolated_xi(s,xi);
1970 
1971  // Get isotropic growth
1972  double gamma;
1973  // Dummy integration point
1974  unsigned ipt=0;
1975  this->get_isotropic_growth(ipt,s,xi,gamma);
1976 
1977  //Output the x,y,..
1978  for(unsigned i=0;i<DIM;i++)
1979  {
1980  //outfile << x[i] << " ";
1981  fprintf(file_pt,"%g ",x[i]);
1982  }
1983  // Output xi0,xi1,..
1984  for(unsigned i=0;i<DIM;i++)
1985  {
1986  //outfile << xi[i] << " ";
1987  fprintf(file_pt,"%g ",xi[i]);
1988  }
1989  // Output growth
1990  //outfile << gamma << " ";
1991  fprintf(file_pt,"%g ",gamma);
1992 
1993  // Output pressure
1994  //outfile << interpolated_solid_p(s) << " ";
1995  //outfile << std::endl;
1996  fprintf(file_pt,"%g \n",interpolated_solid_p(s));
1997 
1998  }
1999  }
2000 
2001  break;
2002 
2003  case 3:
2004  //Tecplot header info
2005  //outfile << "ZONE I=" << n_plot
2006  // << ", J=" << n_plot
2007  // << ", K=" << n_plot << std::endl;
2008  fprintf(file_pt,"ZONE I=%i, J=%i, K=%i \n",n_plot,n_plot,n_plot);
2009 
2010  //Loop over element nodes
2011  for(unsigned l3=0;l3<n_plot;l3++)
2012  {
2013  s[2] = -1.0 + l3*2.0/(n_plot-1);
2014  for(unsigned l2=0;l2<n_plot;l2++)
2015  {
2016  s[1] = -1.0 + l2*2.0/(n_plot-1);
2017  for(unsigned l1=0;l1<n_plot;l1++)
2018  {
2019  s[0] = -1.0 + l1*2.0/(n_plot-1);
2020 
2021  // Get Eulerian and Lagrangian coordinates
2022  this->interpolated_x(s,x);
2023  this->interpolated_xi(s,xi);
2024 
2025  // Get isotropic growth
2026  double gamma;
2027  // Dummy integration point
2028  unsigned ipt=0;
2029  this->get_isotropic_growth(ipt,s,xi,gamma);
2030 
2031  //Output the x,y,..
2032  for(unsigned i=0;i<DIM;i++)
2033  {
2034  //outfile << x[i] << " ";
2035  fprintf(file_pt,"%g ",x[i]);
2036  }
2037  // Output xi0,xi1,..
2038  for(unsigned i=0;i<DIM;i++)
2039  {
2040  //outfile << xi[i] << " ";
2041  fprintf(file_pt,"%g ",xi[i]);
2042  }
2043  // Output growth
2044  //outfile << gamma << " ";
2045  fprintf(file_pt,"%g ",gamma);
2046 
2047  // Output pressure
2048  //outfile << interpolated_solid_p(s) << " ";
2049  //outfile << std::endl;
2050  fprintf(file_pt,"%g \n",interpolated_solid_p(s));
2051  }
2052  }
2053  }
2054  break;
2055 
2056  default:
2057  std::ostringstream error_message;
2058  error_message << "No output routine for PVDEquationsWithPressure<" <<
2059  DIM << "> elements. Write it yourself!" << std::endl;
2060  throw OomphLibError(error_message.str(),
2061  OOMPH_CURRENT_FUNCTION,
2062  OOMPH_EXCEPTION_LOCATION);
2063  }
2064 }
2065 
2066 
2067 //=======================================================================
2068 /// Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components
2069 //=======================================================================
2070 template <unsigned DIM>
2072  const unsigned &n_plot)
2073 {
2074 
2075  Vector<double> x(DIM);
2076  Vector<double> xi(DIM);
2077  Vector<double> s(DIM);
2078  DenseMatrix<double> stress_or_strain(DIM,DIM);
2079 
2080  // Tecplot header info
2081  outfile << this->tecplot_zone_string(n_plot);
2082 
2083  // Loop over plot points
2084  unsigned num_plot_points=this->nplot_points(n_plot);
2085  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
2086  {
2087  // Get local coordinates of plot point
2088  this->get_s_plot(iplot,n_plot,s);
2089 
2090  // Get Eulerian and Lagrangian coordinates
2091  this->interpolated_x(s,x);
2092  this->interpolated_xi(s,xi);
2093 
2094  // Get isotropic growth
2095  double gamma;
2096  // Dummy integration point
2097  unsigned ipt=0;
2098  this->get_isotropic_growth(ipt,s,xi,gamma);
2099 
2100  //Output the x,y,..
2101  for(unsigned i=0;i<DIM;i++)
2102  {outfile << x[i] << " ";}
2103 
2104  // Output xi0,xi1,..
2105  for(unsigned i=0;i<DIM;i++)
2106  {outfile << xi[i] << " ";}
2107 
2108  // Output growth
2109  outfile << gamma << " ";
2110 
2111  //Output pressure
2112  outfile << interpolated_solid_p(s) << " ";
2113 
2114  //get the strain
2115  this->get_strain(s,stress_or_strain);
2116  for(unsigned i=0;i<DIM;i++)
2117  {
2118  for(unsigned j=0;j<=i;j++)
2119  {
2120  outfile << stress_or_strain(j,i) << " " ;
2121  }
2122  }
2123 
2124  //get the stress
2125  this->get_stress(s,stress_or_strain);
2126  for(unsigned i=0;i<DIM;i++)
2127  {
2128  for(unsigned j=0;j<=i;j++)
2129  {
2130  outfile << stress_or_strain(j,i) << " " ;
2131  }
2132  }
2133 
2134 
2135  outfile << std::endl;
2136  }
2137 
2138 
2139  // Write tecplot footer (e.g. FE connectivity lists)
2140  this->write_tecplot_zone_footer(outfile,n_plot);
2141  outfile << std::endl;
2142 }
2143 
2144 
2145 //=======================================================================
2146 /// Compute the diagonal of the velocity mass matrix for LSC
2147 /// preconditioner.
2148 //=======================================================================
2149 template <unsigned DIM>
2151  Vector<double> &mass_diag)
2152 {
2153  // Resize and initialise
2154  mass_diag.assign(this->ndof(),0.0);
2155 
2156  // find out how many nodes there are
2157  unsigned n_node = this->nnode();
2158 
2159  //Find out how many position types of dof there are
2160  const unsigned n_position_type = this->nnodal_position_type();
2161 
2162  //Set up memory for the shape functions
2163  Shape psi(n_node,n_position_type);
2164  DShape dpsidxi(n_node,n_position_type,DIM);
2165 
2166  //Number of integration points
2167  unsigned n_intpt = this->integral_pt()->nweight();
2168 
2169  //Integer to store the local equations no
2170  int local_eqn=0;
2171 
2172  //Loop over the integration points
2173  for(unsigned ipt=0; ipt<n_intpt; ipt++)
2174  {
2175  //Get the integral weight
2176  double w = this->integral_pt()->weight(ipt);
2177 
2178  //Call the derivatives of the shape functions
2179  double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
2180 
2181  //Premultiply weights and Jacobian
2182  double W = w*J;
2183 
2184  //Loop over the nodes
2185  for(unsigned l=0; l<n_node; l++)
2186  {
2187  //Loop over the types of dof
2188  for(unsigned k=0;k<n_position_type;k++)
2189  {
2190  //Loop over the directions
2191  for(unsigned i=0; i<DIM; i++)
2192  {
2193  //Get the equation number
2194  local_eqn = this->position_local_eqn(l,k,i);
2195 
2196  //If not a boundary condition
2197  if(local_eqn >= 0)
2198  {
2199  //Add the contribution
2200  mass_diag[local_eqn] += pow(psi(l,k),2) * W;
2201  } //End of if not boundary condition statement
2202  }//End of loop over dimension
2203  } // End of dof type
2204  } //End of loop over basis functions
2205 
2206  }
2207 
2208 }
2209 
2210 
2211 
2212 //=======================================================================
2213 /// Compute the contravariant second Piola Kirchoff stress at a given local
2214 /// coordinate. Note: this replicates a lot of code that is already
2215 /// coontained in get_residuals() but without sacrificing efficiency
2216 /// (re-computing the shape functions several times) or creating
2217 /// helper functions with horrendous interfaces (to pass all the
2218 /// functions which shouldn't be recomputed) about this is
2219 /// unavoidable.
2220 //=======================================================================
2221 template <unsigned DIM>
2223  DenseMatrix<double> &sigma)
2224 {
2225  //Find out how many nodes there are
2226  unsigned n_node = this->nnode();
2227 
2228  //Find out how many positional dofs there are
2229  unsigned n_position_type = this->nnodal_position_type();
2230 
2231  //Find out how many pressure dofs there are
2232  unsigned n_solid_pres = this->npres_solid();
2233 
2234  //Set up memory for the shape functions
2235  Shape psi(n_node,n_position_type);
2236  DShape dpsidxi(n_node,n_position_type,DIM);
2237 
2238  //Set up memory for the pressure shape functions
2239  Shape psisp(n_solid_pres);
2240 
2241  //Find values of shape function
2242  solid_pshape(s,psisp);
2243 
2244  //Call the derivatives of the shape functions (ignore Jacobian)
2245  (void) this->dshape_lagrangian(s,psi,dpsidxi);
2246 
2247  // Lagrangian coordinates
2248  Vector<double> xi(DIM);
2249  this->interpolated_xi(s,xi);
2250 
2251  //Get isotropic growth factor
2252  double gamma;
2253  //Dummy integration point
2254  unsigned ipt=0;
2255  this->get_isotropic_growth(ipt,s,xi,gamma);
2256 
2257  // We use Cartesian coordinates as the reference coordinate
2258  // system. In this case the undeformed metric tensor is always
2259  // the identity matrix -- stretched by the isotropic growth
2260  double diag_entry=pow(gamma,2.0/double(DIM));
2261  DenseMatrix<double> g(DIM);
2262  for(unsigned i=0;i<DIM;i++)
2263  {
2264  for(unsigned j=0;j<DIM;j++)
2265  {
2266  if(i==j) {g(i,j) = diag_entry;}
2267  else {g(i,j) = 0.0;}
2268  }
2269  }
2270 
2271 
2272  //Calculate interpolated values of the derivative of global position
2273  //wrt lagrangian coordinates
2274  DenseMatrix<double> interpolated_G(DIM);
2275 
2276  //Initialise to zero
2277  for(unsigned i=0;i<DIM;i++)
2278  {for(unsigned j=0;j<DIM;j++) {interpolated_G(i,j) = 0.0;}}
2279 
2280  //Calculate displacements and derivatives
2281  for(unsigned l=0;l<n_node;l++)
2282  {
2283  //Loop over positional dofs
2284  for(unsigned k=0;k<n_position_type;k++)
2285  {
2286  //Loop over displacement components (deformed position)
2287  for(unsigned i=0;i<DIM;i++)
2288  {
2289  //Loop over derivative directions
2290  for(unsigned j=0;j<DIM;j++)
2291  {
2292  interpolated_G(j,i) += this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
2293  }
2294  }
2295  }
2296  }
2297 
2298  //Declare and calculate the deformed metric tensor
2299  DenseMatrix<double> G(DIM);
2300 
2301  //Assign values of G
2302  for(unsigned i=0;i<DIM;i++)
2303  {
2304  //Do upper half of matrix
2305  //Note that j must be signed here for the comparison test to work
2306  //Also i must be cast to an int
2307  for(int j=(DIM-1);j>=static_cast<int>(i);j--)
2308  {
2309  //Initialise G(i,j) to zero
2310  G(i,j) = 0.0;
2311  //Now calculate the dot product
2312  for(unsigned k=0;k<DIM;k++)
2313  {
2314  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
2315  }
2316  }
2317  //Matrix is symmetric so just copy lower half
2318  for(int j=(i-1);j>=0;j--)
2319  {
2320  G(i,j) = G(j,i);
2321  }
2322  }
2323 
2324 
2325  //Calculate the interpolated solid pressure
2326  double interpolated_solid_p=0.0;
2327  for(unsigned l=0;l<n_solid_pres;l++)
2328  {
2329  interpolated_solid_p += solid_p(l)*psisp[l];
2330  }
2331 
2332  //Now calculate the deviatoric stress and all pressure-related
2333  //quantitites
2334  DenseMatrix<double> sigma_dev(DIM), Gup(DIM);
2335  double detG = 0.0;
2336  double gen_dil=0.0;
2337  double inv_kappa=0.0;
2338 
2339  // Incompressible: Compute the deviatoric part of the stress tensor, the
2340  // contravariant deformed metric tensor and the determinant
2341  // of the deformed covariant metric tensor.
2342 
2343  if(Incompressible)
2344  {
2345  get_stress(g,G,sigma_dev,Gup,detG);
2346  }
2347  // Nearly incompressible: Compute the deviatoric part of the
2348  // stress tensor, the contravariant deformed metric tensor,
2349  // the generalised dilatation and the inverse bulk modulus.
2350  else
2351  {
2352  get_stress(g,G,sigma_dev,Gup,gen_dil,inv_kappa);
2353  }
2354 
2355  // Get complete stress
2356  for (unsigned i=0;i<DIM;i++)
2357  {
2358  for (unsigned j=0;j<DIM;j++)
2359  {
2360  sigma(i,j)=-interpolated_solid_p*Gup(i,j)+sigma_dev(i,j);
2361  }
2362  }
2363 
2364 }
2365 
2366 
2367 //////////////////////////////////////////////////////////////////////
2368 //////////////////////////////////////////////////////////////////////
2369 //////////////////////////////////////////////////////////////////////
2370 
2371 
2372 
2373 //====================================================================
2374 /// Data for the number of Variables at each node
2375 //====================================================================
2376 template<>
2378 {1,0,1,0,0,0,1,0,1};
2379 
2380 //==========================================================================
2381 /// Conversion from pressure dof to Node number at which pressure is stored
2382 //==========================================================================
2383 template<>
2385 {0,2,6,8};
2386 
2387 //====================================================================
2388 /// Data for the number of Variables at each node
2389 //====================================================================
2390 template<>
2392 {1,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,1};
2393 
2394 //==========================================================================
2395 /// Conversion from pressure dof to Node number at which pressure is stored
2396 //==========================================================================
2397 template<>
2399 {0,2,6,8,18,20,24,26};
2400 
2401 
2402 ///////////////////////////////////////////////////////////////////////////
2403 ///////////////////////////////////////////////////////////////////////////
2404 ///////////////////////////////////////////////////////////////////////////
2405 
2406 //=======================================================================
2407 /// Data for the number of variables at each node
2408 //=======================================================================
2409 template<>
2410 const unsigned TPVDElementWithContinuousPressure<2>::
2411 Initial_Nvalue[6]={1,1,1,0,0,0};
2412 
2413 //=======================================================================
2414 /// Data for the pressure conversion array
2415 //=======================================================================
2416 template<>
2418 
2419 //=======================================================================
2420 /// Data for the number of variables at each node
2421 //=======================================================================
2422 template<>
2424 ={1,1,1,1,0,0,0,0,0,0};
2425 
2426 //=======================================================================
2427 /// Data for the pressure conversion array
2428 //=======================================================================
2429 template<>
2431 
2432 
2433 
2434 //Instantiate the required elements
2435 template class QPVDElementWithPressure<2>;
2437 template class PVDEquationsBase<2>;
2438 template class PVDEquations<2>;
2439 template class PVDEquationsWithPressure<2>;
2440 
2441 template class QPVDElementWithPressure<3>;
2443 template class PVDEquationsBase<3>;
2444 template class PVDEquations<3>;
2445 template class PVDEquationsWithPressure<3>;
2446 
2449 
2450 }
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void get_energy(double &pot_en, double &kin_en)
Get potential (strain) and kinetic energy.
virtual void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
cstr elem_len * i
Definition: cfortran.h:607
double & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format. WARNING: Only for experienced users. Only use this if raw speed is of the essence, as in the solid mechanics problems.
Definition: shape.h:468
unsigned offset(const unsigned long &i, const unsigned long &j) const
Caculate the offset in flat-packed C-style, column-major format, required for a given i...
Definition: shape.h:482
A Rank 4 Tensor class.
Definition: matrices.h:1625
virtual void fill_in_generic_residual_contribution_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Returns the residuals for the discretised principle of virtual displacements, formulated in the incom...
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)...
Definition: matrices.cc:55
unsigned offset(const unsigned long &i, const unsigned long &j) const
Caculate the offset in flat-packed C-style, column-major format, required for a given i...
Definition: matrices.h:1960
T & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format. WARNING: Only for experienced users. Only use this if raw speed is of the essence, as in the solid mechanics problems.
Definition: matrices.h:2339
T & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format. WARNING: Only for experienced users. Only use this if raw speed is of the essence, as in the solid mechanics problems.
Definition: matrices.h:1946
unsigned offset(const unsigned long &i, const unsigned long &j, const unsigned long &k) const
Caculate the offset in flat-packed Cy-style, column-major format, required for a given i...
Definition: matrices.h:2354
A Rank 5 Tensor class.
Definition: matrices.h:1976
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
static char t char * s
Definition: cfortran.h:572
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)
Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified loc...
void extended_output(std::ostream &outfile, const unsigned &n_plot)
Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress components.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
void get_deformed_covariant_basis_vectors(const Vector< double > &s, DenseMatrix< double > &def_covariant_basis)
Return the deformed covariant basis vectors at specified local coordinate: def_covariant_basis(i,j) is the j-th component of the i-th basis vector.
void extended_output(std::ostream &outfile, const unsigned &n_plot)
Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress.
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)
Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified loc...
void get_mass_matrix_diagonal(Vector< double > &mass_diag)
Compute the diagonal of the displacement mass matrix for LSC preconditioner.
void output(std::ostream &outfile)
Output: x,y,[z],xi0,xi1,[xi2],gamma.
void eigenvalues_by_jacobi(Vector< double > &eigen_val, DenseMatrix< double > &eigen_vect) const
Determine eigenvalues and eigenvectors, using Jacobi rotations. Only for symmetric matrices...
Definition: matrices.cc:231
void output(std::ostream &outfile)
Output: x,y,[z],xi0,xi1,[xi2],p,gamma.
ConstitutiveLaw * Constitutive_law_pt
Create constitutive law (for smoothing by nonlinear elasticity)
Definition: mesh_smooth.h:63
virtual bool requires_incompressibility_constraint()=0
Pure virtual function in which the user must declare if the constitutive equation requires an incompr...
void get_principal_stress(const Vector< double > &s, DenseMatrix< double > &principal_stress_vector, Vector< double > &principal_stress)
Compute principal stress vectors and (scalar) principal stresses at specified local coordinate...