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