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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for elements that solve the principle of virtual
27// equations of solid mechanics
28
29#include "solid_elements.h"
30
31
32namespace 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
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
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>;
2491 template class PVDEquationsBase<2>;
2492 template class PVDEquations<2>;
2493 template class PVDEquationsWithPressure<2>;
2494
2495 template class QPVDElementWithPressure<3>;
2497 template class PVDEquationsBase<3>;
2498 template class PVDEquations<3>;
2499 template class PVDEquationsWithPressure<3>;
2500
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...