refineable_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 member functions and static member data for refineable solid
27// mechanics elements
28
30
31namespace oomph
32{
33 //====================================================================
34 /// Residuals for Refineable QPVDElements
35 //====================================================================
36 template<unsigned DIM>
39 DenseMatrix<double>& jacobian,
40 const unsigned& flag)
41 {
42#ifdef PARANOID
43 // Check if the constitutive equation requires the explicit imposition of an
44 // incompressibility constraint
46 {
47 throw OomphLibError("RefineablePVDEquations cannot be used with "
48 "incompressible constitutive laws.",
49 OOMPH_CURRENT_FUNCTION,
50 OOMPH_EXCEPTION_LOCATION);
51 }
52#endif
53
54 // Simply set up initial condition?
55 if (Solid_ic_pt != 0)
56 {
57 get_residuals_for_solid_ic(residuals);
58 return;
59 }
60
61 // Find out how many nodes there are
62 const unsigned n_node = nnode();
63
64 // Find out how many positional dofs there are
65 const unsigned n_position_type = this->nnodal_position_type();
66
67 // Integers to store local equation numbers
68 int local_eqn = 0;
69
70 // Timescale ratio (non-dim density)
71 double lambda_sq = this->lambda_sq();
72
73 // Time factor
74 double time_factor = 0.0;
75 if (lambda_sq > 0)
76 {
77 time_factor = this->node_pt(0)->position_time_stepper_pt()->weight(2, 0);
78 }
79
80 // Set up memory for the shape functions
81 Shape psi(n_node, n_position_type);
82 DShape dpsidxi(n_node, n_position_type, DIM);
83
84 // Set the value of n_intpt -- the number of integration points
85 const unsigned n_intpt = integral_pt()->nweight();
86
87 // Set the vector to hold the local coordinates in the element
88 Vector<double> s(DIM);
89
90 // Loop over the integration points
91 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
92 {
93 // Assign the values of s
94 for (unsigned i = 0; i < DIM; ++i)
95 {
96 s[i] = integral_pt()->knot(ipt, i);
97 }
98
99 // Get the integral weight
100 double w = integral_pt()->weight(ipt);
101
102 // Call the derivatives of the shape functions (and get Jacobian)
103 double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
104
105 // Calculate interpolated values of the derivative of global position
106 // wrt lagrangian coordinates
107 DenseMatrix<double> interpolated_G(DIM, DIM, 0.0);
108
109 // Setup memory for accelerations
110 Vector<double> accel(DIM, 0.0);
111
112 // Storage for Lagrangian coordinates (initialised to zero)
113 Vector<double> interpolated_xi(DIM, 0.0);
114
115 // Calculate displacements and derivatives and lagrangian coordinates
116 for (unsigned l = 0; l < n_node; l++)
117 {
118 // Loop over positional dofs
119 for (unsigned k = 0; k < n_position_type; k++)
120 {
121 double psi_ = psi(l, k);
122 // Loop over displacement components (deformed position)
123 for (unsigned i = 0; i < DIM; i++)
124 {
125 // Calculate the Lagrangian coordinates and the accelerations
126 interpolated_xi[i] += lagrangian_position_gen(l, k, i) * psi_;
127
128 // Only compute accelerations if inertia is switched on
129 // otherwise the timestepper might not be able to
130 // work out dx_gen_dt(2,...)
131 if ((lambda_sq > 0.0) && (this->Unsteady))
132 {
133 accel[i] += dnodal_position_gen_dt(2, l, k, i) * psi_;
134 }
135
136 // Loop over derivative directions
137 for (unsigned j = 0; j < DIM; j++)
138 {
139 interpolated_G(j, i) +=
140 this->nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
141 }
142 }
143 }
144 }
145
146 // Get isotropic growth factor
147 double gamma = 1.0;
148 this->get_isotropic_growth(ipt, s, interpolated_xi, gamma);
149
150 // Get body force at current time
151 Vector<double> b(DIM);
152 this->body_force(interpolated_xi, b);
153
154 // We use Cartesian coordinates as the reference coordinate
155 // system. In this case the undeformed metric tensor is always
156 // the identity matrix -- stretched by the isotropic growth
157 double diag_entry = pow(gamma, 2.0 / double(DIM));
158 DenseMatrix<double> g(DIM);
159 for (unsigned i = 0; i < DIM; i++)
160 {
161 for (unsigned j = 0; j < DIM; j++)
162 {
163 if (i == j)
164 {
165 g(i, j) = diag_entry;
166 }
167 else
168 {
169 g(i, j) = 0.0;
170 }
171 }
172 }
173
174 // Premultiply the undeformed volume ratio (from the isotropic
175 // growth), the weights and the Jacobian
176 double W = gamma * w * J;
177
178 // Declare and calculate the deformed metric tensor
179 DenseMatrix<double> G(DIM);
180
181 // Assign values of G
182 for (unsigned i = 0; i < DIM; i++)
183 {
184 // Do upper half of matrix
185 for (unsigned j = i; j < DIM; j++)
186 {
187 // Initialise G(i,j) to zero
188 G(i, j) = 0.0;
189 // Now calculate the dot product
190 for (unsigned k = 0; k < DIM; k++)
191 {
192 G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
193 }
194 }
195 // Matrix is symmetric so just copy lower half
196 for (unsigned j = 0; j < i; j++)
197 {
198 G(i, j) = G(j, i);
199 }
200 }
201
202 // Now calculate the stress tensor from the constitutive law
203 DenseMatrix<double> sigma(DIM);
204 this->get_stress(g, G, sigma);
205
206 // Get stress derivative by FD only needed for Jacobian
207 //-----------------------------------------------------
208
209 // Stress derivative
210 RankFourTensor<double> d_stress_dG(DIM, DIM, DIM, DIM, 0.0);
211 // Derivative of metric tensor w.r.t. to nodal coords
213 n_node, n_position_type, DIM, DIM, DIM, 0.0);
214
215 // Get Jacobian too?
216 if (flag == 1)
217 {
218 // Derivative of metric tensor w.r.t. to discrete positional dofs
219 // NOTE: Since G is symmetric we only compute the upper triangle
220 // and DO NOT copy the entries across. Subsequent computations
221 // must (and, in fact, do) therefore only operate with upper
222 // triangular entries
223 for (unsigned ll = 0; ll < n_node; ll++)
224 {
225 for (unsigned kk = 0; kk < n_position_type; kk++)
226 {
227 for (unsigned ii = 0; ii < DIM; ii++)
228 {
229 for (unsigned aa = 0; aa < DIM; aa++)
230 {
231 for (unsigned bb = aa; bb < DIM; bb++)
232 {
233 d_G_dX(ll, kk, ii, aa, bb) =
234 interpolated_G(aa, ii) * dpsidxi(ll, kk, bb) +
235 interpolated_G(bb, ii) * dpsidxi(ll, kk, aa);
236 }
237 }
238 }
239 }
240 }
241
242 // Get the "upper triangular"
243 // entries of the derivatives of the stress tensor with
244 // respect to G
245 this->get_d_stress_dG_upper(g, G, sigma, d_stress_dG);
246 }
247
248
249 // Add pre-stress
250 for (unsigned i = 0; i < DIM; i++)
251 {
252 for (unsigned j = 0; j < DIM; j++)
253 {
254 sigma(i, j) += this->prestress(i, j, interpolated_xi);
255 }
256 }
257
258 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
259 // DISPLACEMENTS========
260
261
262 // Default setting for non-hanging node
263 unsigned n_master = 1;
264 double hang_weight = 1.0;
265
266 // Loop over the test functions, nodes of the element
267 for (unsigned l = 0; l < n_node; l++)
268 {
269 // Get pointer to local node l
270 Node* local_node_pt = node_pt(l);
271
272 // Cache hang status
273 bool is_hanging = local_node_pt->is_hanging();
274
275 // If the node is a hanging node
276 if (is_hanging)
277 {
278 n_master = local_node_pt->hanging_pt()->nmaster();
279 }
280 // Otherwise the node is its own master
281 else
282 {
283 n_master = 1;
284 }
285
286
287 // Storage for local equation numbers at node indexed by
288 // type and direction
289 DenseMatrix<int> position_local_eqn_at_node(n_position_type, DIM);
290
291 // Loop over the master nodes
292 for (unsigned m = 0; m < n_master; m++)
293 {
294 if (is_hanging)
295 {
296 // Find the equation numbers
297 position_local_eqn_at_node = local_position_hang_eqn(
298 local_node_pt->hanging_pt()->master_node_pt(m));
299
300 // Find the hanging node weight
301 hang_weight = local_node_pt->hanging_pt()->master_weight(m);
302 }
303 else
304 {
305 // Loop of types of dofs
306 for (unsigned k = 0; k < n_position_type; k++)
307 {
308 // Loop over the displacement components
309 for (unsigned i = 0; i < DIM; i++)
310 {
311 position_local_eqn_at_node(k, i) = position_local_eqn(l, k, i);
312 }
313 }
314
315 // Hang weight is one
316 hang_weight = 1.0;
317 }
318
319 // Loop of types of dofs
320 for (unsigned k = 0; k < n_position_type; k++)
321 {
322 // Offset for faster access
323 const unsigned offset5 = dpsidxi.offset(l, k);
324
325 // Loop over the displacement components
326 for (unsigned i = 0; i < DIM; i++)
327 {
328 local_eqn = position_local_eqn_at_node(k, i);
329
330 /*IF it's not a boundary condition*/
331 if (local_eqn >= 0)
332 {
333 // Initialise the contribution
334 double sum = 0.0;
335
336 // Acceleration and body force
337 sum += (lambda_sq * accel[i] - b[i]) * psi(l, k);
338
339 // Stress term
340 for (unsigned a = 0; a < DIM; a++)
341 {
342 unsigned count = offset5;
343 for (unsigned b = 0; b < DIM; b++)
344 {
345 // Add the stress terms to the residuals
346 sum += sigma(a, b) * interpolated_G(a, i) *
347 dpsidxi.raw_direct_access(count);
348 ++count;
349 }
350 }
351 residuals[local_eqn] += W * sum * hang_weight;
352
353
354 // Get Jacobian too?
355 if (flag == 1)
356 {
357 // Offset for faster access in general stress loop
358 const unsigned offset1 = d_G_dX.offset(l, k, i);
359
360 // Default setting for non-hanging node
361 unsigned nn_master = 1;
362 double hhang_weight = 1.0;
363
364 // Loop over the nodes of the element again
365 for (unsigned ll = 0; ll < n_node; ll++)
366 {
367 // Get pointer to local node ll
368 Node* llocal_node_pt = node_pt(ll);
369
370 // Cache hang status
371 bool iis_hanging = llocal_node_pt->is_hanging();
372
373 // If the node is a hanging node
374 if (iis_hanging)
375 {
376 nn_master = llocal_node_pt->hanging_pt()->nmaster();
377 }
378 // Otherwise the node is its own master
379 else
380 {
381 nn_master = 1;
382 }
383
384
385 // Storage for local unknown numbers at node indexed by
386 // type and direction
387 DenseMatrix<int> position_local_unk_at_node(n_position_type,
388 DIM);
389
390 // Loop over the master nodes
391 for (unsigned mm = 0; mm < nn_master; mm++)
392 {
393 if (iis_hanging)
394 {
395 // Find the unknown numbers
396 position_local_unk_at_node = local_position_hang_eqn(
397 llocal_node_pt->hanging_pt()->master_node_pt(mm));
398
399 // Find the hanging node weight
400 hhang_weight =
401 llocal_node_pt->hanging_pt()->master_weight(mm);
402 }
403 else
404 {
405 // Loop of types of dofs
406 for (unsigned kk = 0; kk < n_position_type; kk++)
407 {
408 // Loop over the displacement components
409 for (unsigned ii = 0; ii < DIM; ii++)
410 {
411 position_local_unk_at_node(kk, ii) =
412 position_local_eqn(ll, kk, ii);
413 }
414 }
415
416 // Hang weight is one
417 hhang_weight = 1.0;
418 }
419
420
421 // Loop of types of dofs again
422 for (unsigned kk = 0; kk < n_position_type; kk++)
423 {
424 // Loop over the displacement components again
425 for (unsigned ii = 0; ii < DIM; ii++)
426 {
427 // Get the number of the unknown
428 int local_unknown =
429 position_local_unk_at_node(kk, ii);
430
431
432 /*IF it's not a boundary condition*/
433 if (local_unknown >= 0)
434 {
435 // Offset for faster access in general stress loop
436 const unsigned offset2 = d_G_dX.offset(ll, kk, ii);
437 const unsigned offset4 = dpsidxi.offset(ll, kk);
438
439
440 // General stress term
441 //--------------------
442 double sum = 0.0;
443 unsigned count1 = offset1;
444 for (unsigned a = 0; a < DIM; a++)
445 {
446 // Bump up direct access because we're only
447 // accessing upper triangle
448 count1 += a;
449 for (unsigned b = a; b < DIM; b++)
450 {
451 double factor =
452 d_G_dX.raw_direct_access(count1);
453 if (a == b) factor *= 0.5;
454
455 // Offset for faster access
456 unsigned offset3 = d_stress_dG.offset(a, b);
457 unsigned count2 = offset2;
458 unsigned count3 = offset3;
459
460 for (unsigned aa = 0; aa < DIM; aa++)
461 {
462 // Bump up direct access because we're only
463 // accessing upper triangle
464 count2 += aa;
465 count3 += aa;
466
467 // Only upper half of derivatives w.r.t.
468 // symm tensor
469 for (unsigned bb = aa; bb < DIM; bb++)
470 {
471 sum +=
472 factor *
473 d_stress_dG.raw_direct_access(count3) *
474 d_G_dX.raw_direct_access(count2);
475 ++count2;
476 ++count3;
477 }
478 }
479 ++count1;
480 }
481 }
482
483 // Multiply by weight and add contribution
484 // (Add directly because this bit is nonsymmetric)
485 jacobian(local_eqn, local_unknown) +=
486 sum * W * hang_weight * hhang_weight;
487
488 // Only upper triangle (no separate test for bc as
489 // local_eqn is already nonnegative)
490 if ((i == ii) && (local_unknown >= local_eqn))
491 {
492 // Initialise the contribution
493 double sum = 0.0;
494
495 // Inertia term
496 sum += lambda_sq * time_factor * psi(ll, kk) *
497 psi(l, k);
498
499 // Stress term
500 unsigned count4 = offset4;
501 for (unsigned a = 0; a < DIM; a++)
502 {
503 // Cache term
504 const double factor =
505 dpsidxi.raw_direct_access(count4); // ll ,kk
506 ++count4;
507
508 unsigned count5 = offset5;
509 for (unsigned b = 0; b < DIM; b++)
510 {
511 sum +=
512 sigma(a, b) * factor *
513 dpsidxi.raw_direct_access(count5); // l ,k
514 ++count5;
515 }
516 }
517
518 // Multiply by weights to form contribution
519 double sym_entry =
520 sum * W * hang_weight * hhang_weight;
521 // Add contribution to jacobian
522 jacobian(local_eqn, local_unknown) += sym_entry;
523 // Add to lower triangular entries
524 if (local_eqn != local_unknown)
525 {
526 jacobian(local_unknown, local_eqn) += sym_entry;
527 }
528 }
529 } // End of if not boundary condition
530 }
531 }
532 }
533 }
534 }
535
536 } // End of if not boundary condition
537
538 } // End of loop over coordinate directions
539 } // End of loop over type of dof
540 } // End of loop over master nodes
541 } // End of loop over nodes
542 } // End of loop over integration points
543 }
544
545
546 //=======================================================================
547 /// Compute the diagonal of the velocity mass matrix for LSC
548 /// preconditioner.
549 //=======================================================================
550 template<unsigned DIM>
552 Vector<double>& mass_diag)
553 {
554 // Resize and initialise
555 mass_diag.assign(this->ndof(), 0.0);
556
557 // find out how many nodes there are
558 unsigned n_node = this->nnode();
559
560 // Find out how many position types of dof there are
561 const unsigned n_position_type = this->nnodal_position_type();
562
563 // Set up memory for the shape functions
564 Shape psi(n_node, n_position_type);
565 DShape dpsidxi(n_node, n_position_type, DIM);
566
567 // Number of integration points
568 unsigned n_intpt = this->integral_pt()->nweight();
569
570 // Integer to store the local equations no
571 int local_eqn = 0;
572
573 // Loop over the integration points
574 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
575 {
576 // Get the integral weight
577 double w = this->integral_pt()->weight(ipt);
578
579 // Call the derivatives of the shape functions
580 double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
581
582 // Premultiply weights and Jacobian
583 double W = w * J;
584
585 unsigned n_master = 1;
586 double hang_weight = 1.0;
587
588 // Loop over the nodes
589 for (unsigned l = 0; l < n_node; l++)
590 {
591 // Get pointer to local node l
592 Node* local_node_pt = node_pt(l);
593
594 // Cache hang status
595 bool is_hanging = local_node_pt->is_hanging();
596
597 // If the node is a hanging node
598 if (is_hanging)
599 {
600 n_master = local_node_pt->hanging_pt()->nmaster();
601 }
602 // Otherwise the node is its own master
603 else
604 {
605 n_master = 1;
606 }
607
608 // Storage for local equation numbers at node indexed by
609 // type and direction
610 DenseMatrix<int> position_local_eqn_at_node(n_position_type, DIM);
611
612 // Loop over the master nodes
613 for (unsigned m = 0; m < n_master; m++)
614 {
615 if (is_hanging)
616 {
617 // Find the equation numbers
618 position_local_eqn_at_node = local_position_hang_eqn(
619 local_node_pt->hanging_pt()->master_node_pt(m));
620
621 // Find the hanging node weight
622 hang_weight = local_node_pt->hanging_pt()->master_weight(m);
623 }
624 else
625 {
626 // Loop of types of dofs
627 for (unsigned k = 0; k < n_position_type; k++)
628 {
629 // Loop over the displacement components
630 for (unsigned i = 0; i < DIM; i++)
631 {
632 position_local_eqn_at_node(k, i) = position_local_eqn(l, k, i);
633 }
634 }
635
636 // Hang weight is one
637 hang_weight = 1.0;
638 }
639
640 // Loop over the types of dof
641 for (unsigned k = 0; k < n_position_type; k++)
642 {
643 // Loop over the directions
644 for (unsigned i = 0; i < DIM; i++)
645 {
646 // Get the equation number
647 local_eqn = position_local_eqn_at_node(k, i);
648
649 // If not a boundary condition
650 if (local_eqn >= 0)
651 {
652 // Add the contribution
653 mass_diag[local_eqn] += pow(psi(l, k) * hang_weight, 2) * W;
654 } // End of if not boundary condition statement
655 } // End of loop over dimension
656 } // End of dof type
657 } // End of loop over master nodes
658 } // End of loop over basis functions (nodes)
659 } // End integration loop
660 }
661
662
663 //===========================================================================
664 /// Fill in element's contribution to the elemental
665 /// residual vector and/or Jacobian matrix.
666 /// flag=0: compute only residual vector
667 /// flag=1: compute both, fully analytically
668 /// flag=2: compute both, using FD for the derivatives w.r.t. to the
669 /// discrete displacment dofs.
670 /// flag=3: compute residuals, jacobian (full analytic) and mass matrix
671 /// flag=4: compute residuals, jacobian (FD for derivatives w.r.t.
672 /// displacements) and mass matrix
673 //==========================================================================
674 template<unsigned DIM>
677 Vector<double>& residuals,
678 DenseMatrix<double>& jacobian,
679 DenseMatrix<double>& mass_matrix,
680 const unsigned& flag)
681 {
682#ifdef PARANOID
683 // Check if the constitutive equation requires the explicit imposition of an
684 // incompressibility constraint
686 (!this->Incompressible))
687 {
688 throw OomphLibError("The constitutive law requires the use of the "
689 "incompressible formulation by setting the element's "
690 "member function set_incompressible()",
691 OOMPH_CURRENT_FUNCTION,
692 OOMPH_EXCEPTION_LOCATION);
693 }
694#endif
695
696
697 // Simply set up initial condition?
698 if (Solid_ic_pt != 0)
699 {
700 get_residuals_for_solid_ic(residuals);
701 return;
702 }
703
704 // Find out how many nodes there are
705 const unsigned n_node = nnode();
706
707 // Find out how many position types of dof there are
708 const unsigned n_position_type = this->nnodal_position_type();
709
710 // Find out how many pressure dofs there are
711 const unsigned n_solid_pres = this->npres_solid();
712
713 // Find out the index of the solid dof
714 const int solid_p_index = this->solid_p_nodal_index();
715
716 // Local array of booleans that is true if the l-th pressure value is
717 // hanging This is an optimization because it avoids repeated virtual
718 // function calls
719 bool solid_pressure_dof_is_hanging[n_solid_pres];
720
721 // If the solid pressure is stored at a node
722 if (solid_p_index >= 0)
723 {
724 // Read out whether the solid pressure is hanging
725 for (unsigned l = 0; l < n_solid_pres; ++l)
726 {
727 solid_pressure_dof_is_hanging[l] =
728 solid_pressure_node_pt(l)->is_hanging(solid_p_index);
729 }
730 }
731 // Otherwise the pressure is not stored at a node and so
732 // it cannot hang
733 else
734 {
735 for (unsigned l = 0; l < n_solid_pres; ++l)
736 {
737 solid_pressure_dof_is_hanging[l] = false;
738 }
739 }
740
741 // Integer for storage of local equation and unknown numbers
742 int local_eqn = 0, local_unknown = 0;
743
744 // Timescale ratio (non-dim density)
745 double lambda_sq = this->lambda_sq();
746
747
748 // Time factor
749 double time_factor = 0.0;
750 if (lambda_sq > 0)
751 {
752 time_factor = this->node_pt(0)->position_time_stepper_pt()->weight(2, 0);
753 }
754
755
756 // Set up memory for the shape functions
757 Shape psi(n_node, n_position_type);
758 DShape dpsidxi(n_node, n_position_type, DIM);
759
760 // Set up memory for the pressure shape functions
761 Shape psisp(n_solid_pres);
762
763 // Set the value of n_intpt
764 const unsigned n_intpt = integral_pt()->nweight();
765
766 // Set the vector to hold the local coordinates in the element
767 Vector<double> s(DIM);
768
769 // Loop over the integration points
770 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
771 {
772 // Assign the values of s
773 for (unsigned i = 0; i < DIM; ++i)
774 {
775 s[i] = integral_pt()->knot(ipt, i);
776 }
777
778 // Get the integral weight
779 double w = integral_pt()->weight(ipt);
780
781 // Call the derivatives of the shape functions
782 double J = dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
783
784 // Call the pressure shape functions
785 this->solid_pshape_at_knot(ipt, psisp);
786
787 // Storage for Lagrangian coordinates (initialised to zero)
788 Vector<double> interpolated_xi(DIM, 0.0);
789
790 // Deformed tangent vectors
791 DenseMatrix<double> interpolated_G(DIM, DIM, 0.0);
792
793 // Setup memory for accelerations
794 Vector<double> accel(DIM, 0.0);
795
796 // Calculate displacements and derivatives and lagrangian coordinates
797 for (unsigned l = 0; l < n_node; l++)
798 {
799 // Loop over positional dofs
800 for (unsigned k = 0; k < n_position_type; k++)
801 {
802 double psi_ = psi(l, k);
803 // Loop over displacement components (deformed position)
804 for (unsigned i = 0; i < DIM; i++)
805 {
806 // Calculate the lagrangian coordinates and the accelerations
807 interpolated_xi[i] += lagrangian_position_gen(l, k, i) * psi_;
808
809 // Only compute accelerations if inertia is switched on
810 // otherwise the timestepper might not be able to
811 // work out dx_gen_dt(2,...)
812 if ((lambda_sq > 0.0) && (this->Unsteady))
813 {
814 accel[i] += dnodal_position_gen_dt(2, l, k, i) * psi_;
815 }
816
817 // Loop over derivative directions
818 for (unsigned j = 0; j < DIM; j++)
819 {
820 interpolated_G(j, i) +=
821 nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
822 }
823 }
824 }
825 }
826
827 // Get isotropic growth factor
828 double gamma = 1.0;
829 this->get_isotropic_growth(ipt, s, interpolated_xi, gamma);
830
831 // Get body force at current time
832 Vector<double> b(DIM);
833 this->body_force(interpolated_xi, b);
834
835 // We use Cartesian coordinates as the reference coordinate
836 // system. In this case the undeformed metric tensor is always
837 // the identity matrix -- stretched by the isotropic growth
838 double diag_entry = pow(gamma, 2.0 / double(DIM));
839 DenseMatrix<double> g(DIM);
840 for (unsigned i = 0; i < DIM; i++)
841 {
842 for (unsigned j = 0; j < DIM; j++)
843 {
844 if (i == j)
845 {
846 g(i, j) = diag_entry;
847 }
848 else
849 {
850 g(i, j) = 0.0;
851 }
852 }
853 }
854
855 // Premultiply the undeformed volume ratio (from the isotropic
856 // growth), the weights and the Jacobian
857 double W = gamma * w * J;
858
859 // Calculate the interpolated solid pressure
860 double interpolated_solid_p = 0.0;
861 for (unsigned l = 0; l < n_solid_pres; l++)
862 {
863 interpolated_solid_p += this->solid_p(l) * psisp[l];
864 }
865
866
867 // Declare and calculate the deformed metric tensor
868 DenseMatrix<double> G(DIM);
869
870 // Assign values of G
871 for (unsigned i = 0; i < DIM; i++)
872 {
873 // Do upper half of matrix
874 for (unsigned j = i; j < DIM; j++)
875 {
876 // Initialise G(i,j) to zero
877 G(i, j) = 0.0;
878 // Now calculate the dot product
879 for (unsigned k = 0; k < DIM; k++)
880 {
881 G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
882 }
883 }
884 // Matrix is symmetric so just copy lower half
885 for (unsigned j = 0; j < i; j++)
886 {
887 G(i, j) = G(j, i);
888 }
889 }
890
891 // Now calculate the deviatoric stress and all pressure-related
892 // quantitites
893 DenseMatrix<double> sigma(DIM, DIM), sigma_dev(DIM, DIM), Gup(DIM, DIM);
894 double detG = 0.0;
895 double gen_dil = 0.0;
896 double inv_kappa = 0.0;
897
898 // Get stress derivative by FD only needed for Jacobian
899
900 // Stress etc derivatives
901 RankFourTensor<double> d_stress_dG(DIM, DIM, DIM, DIM, 0.0);
902 // RankFourTensor<double> d_Gup_dG(DIM,DIM,DIM,DIM,0.0);
903 DenseMatrix<double> d_detG_dG(DIM, DIM, 0.0);
904 DenseMatrix<double> d_gen_dil_dG(DIM, DIM, 0.0);
905
906 // Derivative of metric tensor w.r.t. to nodal coords
908 n_node, n_position_type, DIM, DIM, DIM, 0.0);
909
910 // Get Jacobian too?
911 if ((flag == 1) || (flag == 3))
912 {
913 // Derivative of metric tensor w.r.t. to discrete positional dofs
914 // NOTE: Since G is symmetric we only compute the upper triangle
915 // and DO NOT copy the entries across. Subsequent computations
916 // must (and, in fact, do) therefore only operate with upper
917 // triangular entries
918 for (unsigned ll = 0; ll < n_node; ll++)
919 {
920 for (unsigned kk = 0; kk < n_position_type; kk++)
921 {
922 for (unsigned ii = 0; ii < DIM; ii++)
923 {
924 for (unsigned aa = 0; aa < DIM; aa++)
925 {
926 for (unsigned bb = aa; bb < DIM; bb++)
927 {
928 d_G_dX(ll, kk, ii, aa, bb) =
929 interpolated_G(aa, ii) * dpsidxi(ll, kk, bb) +
930 interpolated_G(bb, ii) * dpsidxi(ll, kk, aa);
931 }
932 }
933 }
934 }
935 }
936 }
937
938
939 // Incompressible: Compute the deviatoric part of the stress tensor, the
940 // contravariant deformed metric tensor and the determinant
941 // of the deformed covariant metric tensor.
942 if (this->Incompressible)
943 {
944 this->get_stress(g, G, sigma_dev, Gup, detG);
945
946 // Get full stress
947 for (unsigned a = 0; a < DIM; a++)
948 {
949 for (unsigned b = 0; b < DIM; b++)
950 {
951 sigma(a, b) = sigma_dev(a, b) - interpolated_solid_p * Gup(a, b);
952 }
953 }
954
955 // Get Jacobian too?
956 if ((flag == 1) || (flag == 3))
957 {
958 // Get the "upper triangular" entries of the
959 // derivatives of the stress tensor with
960 // respect to G
961 this->get_d_stress_dG_upper(
962 g, G, sigma, detG, interpolated_solid_p, d_stress_dG, d_detG_dG);
963 }
964 }
965 // Nearly incompressible: Compute the deviatoric part of the
966 // stress tensor, the contravariant deformed metric tensor,
967 // the generalised dilatation and the inverse bulk modulus.
968 else
969 {
970 this->get_stress(g, G, sigma_dev, Gup, gen_dil, inv_kappa);
971
972 // Get full stress
973 for (unsigned a = 0; a < DIM; a++)
974 {
975 for (unsigned b = 0; b < DIM; b++)
976 {
977 sigma(a, b) = sigma_dev(a, b) - interpolated_solid_p * Gup(a, b);
978 }
979 }
980
981 // Get Jacobian too?
982 if ((flag == 1) || (flag == 3))
983 {
984 // Get the "upper triangular" entries of the derivatives
985 // of the stress tensor with
986 // respect to G
987 this->get_d_stress_dG_upper(g,
988 G,
989 sigma,
990 gen_dil,
991 inv_kappa,
992 interpolated_solid_p,
993 d_stress_dG,
994 d_gen_dil_dG);
995 }
996 }
997
998 // Add pre-stress
999 for (unsigned i = 0; i < DIM; i++)
1000 {
1001 for (unsigned j = 0; j < DIM; j++)
1002 {
1003 sigma(i, j) += this->prestress(i, j, interpolated_xi);
1004 }
1005 }
1006
1007 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
1008 // DISPLACEMENTS========
1009
1010 unsigned n_master = 1;
1011 double hang_weight = 1.0;
1012
1013 // Loop over the test functions, nodes of the element
1014 for (unsigned l = 0; l < n_node; l++)
1015 {
1016 // Get pointer to local node l
1017 Node* local_node_pt = node_pt(l);
1018
1019 // Cache hang status
1020 bool is_hanging = local_node_pt->is_hanging();
1021
1022 // If the node is a hanging node
1023 if (is_hanging)
1024 {
1025 n_master = local_node_pt->hanging_pt()->nmaster();
1026 }
1027 // Otherwise the node is its own master
1028 else
1029 {
1030 n_master = 1;
1031 }
1032
1033
1034 // Storage for local equation numbers at node indexed by
1035 // type and direction
1036 DenseMatrix<int> position_local_eqn_at_node(n_position_type, DIM);
1037
1038 // Loop over the master nodes
1039 for (unsigned m = 0; m < n_master; m++)
1040 {
1041 if (is_hanging)
1042 {
1043 // Find the equation numbers
1044 position_local_eqn_at_node = local_position_hang_eqn(
1045 local_node_pt->hanging_pt()->master_node_pt(m));
1046
1047 // Find the hanging node weight
1048 hang_weight = local_node_pt->hanging_pt()->master_weight(m);
1049 }
1050 else
1051 {
1052 // Loop of types of dofs
1053 for (unsigned k = 0; k < n_position_type; k++)
1054 {
1055 // Loop over the displacement components
1056 for (unsigned i = 0; i < DIM; i++)
1057 {
1058 position_local_eqn_at_node(k, i) = position_local_eqn(l, k, i);
1059 }
1060 }
1061
1062 // Hang weight is one
1063 hang_weight = 1.0;
1064 }
1065
1066
1067 // Loop of types of dofs
1068 for (unsigned k = 0; k < n_position_type; k++)
1069 {
1070 // Offset for faster access
1071 const unsigned offset5 = dpsidxi.offset(l, k);
1072
1073 // Loop over the displacement components
1074 for (unsigned i = 0; i < DIM; i++)
1075 {
1076 local_eqn = position_local_eqn_at_node(k, i);
1077
1078 /*IF it's not a boundary condition*/
1079 if (local_eqn >= 0)
1080 {
1081 // Initialise contribution to zero
1082 double sum = 0.0;
1083
1084 // Acceleration and body force
1085 sum += (lambda_sq * accel[i] - b[i]) * psi(l, k);
1086
1087 // Stress term
1088 for (unsigned a = 0; a < DIM; a++)
1089 {
1090 unsigned count = offset5;
1091 for (unsigned b = 0; b < DIM; b++)
1092 {
1093 // Add the stress terms to the residuals
1094 sum += sigma(a, b) * interpolated_G(a, i) *
1095 dpsidxi.raw_direct_access(count);
1096 ++count;
1097 }
1098 }
1099 residuals[local_eqn] += W * sum * hang_weight;
1100
1101
1102 // Get the mass matrix
1103 // This involves another loop over the points
1104 // because the jacobian may NOT be being calculated analytically
1105 // It could be made more efficient in th event that
1106 // we eventually decide not (never) to
1107 // use finite differences.
1108 if (flag > 2)
1109 {
1110 // Default setting for non-hanging node
1111 unsigned nn_master = 1;
1112 double hhang_weight = 1.0;
1113
1114 // Loop over the nodes of the element again
1115 for (unsigned ll = 0; ll < n_node; ll++)
1116 {
1117 // Get pointer to local node ll
1118 Node* llocal_node_pt = node_pt(ll);
1119
1120 // Cache hang status
1121 bool iis_hanging = llocal_node_pt->is_hanging();
1122
1123 // If the node is a hanging node
1124 if (iis_hanging)
1125 {
1126 nn_master = llocal_node_pt->hanging_pt()->nmaster();
1127 }
1128 // Otherwise the node is its own master
1129 else
1130 {
1131 nn_master = 1;
1132 }
1133
1134
1135 // Storage for local unknown numbers at node indexed by
1136 // type and direction
1137 DenseMatrix<int> position_local_unk_at_node(n_position_type,
1138 DIM);
1139
1140 // Loop over the master nodes
1141 for (unsigned mm = 0; mm < nn_master; mm++)
1142 {
1143 if (iis_hanging)
1144 {
1145 // Find the unknown numbers
1146 position_local_unk_at_node = local_position_hang_eqn(
1147 llocal_node_pt->hanging_pt()->master_node_pt(mm));
1148
1149 // Find the hanging node weight
1150 hhang_weight =
1151 llocal_node_pt->hanging_pt()->master_weight(mm);
1152 }
1153 else
1154 {
1155 // Loop of types of dofs
1156 for (unsigned kk = 0; kk < n_position_type; kk++)
1157 {
1158 // Loop over the displacement components
1159 for (unsigned ii = 0; ii < DIM; ii++)
1160 {
1161 position_local_unk_at_node(kk, ii) =
1162 position_local_eqn(ll, kk, ii);
1163 }
1164 }
1165
1166 // Hang weight is one
1167 hhang_weight = 1.0;
1168 }
1169
1170
1171 // Loop of types of dofs again
1172 for (unsigned kk = 0; kk < n_position_type; kk++)
1173 {
1174 // Get the number of the unknown
1175 int local_unknown = position_local_unk_at_node(kk, i);
1176
1177 /*IF it's not a boundary condition*/
1178 if (local_unknown >= 0)
1179 {
1180 mass_matrix(local_eqn, local_unknown) +=
1181 lambda_sq * psi(l, k) * psi(ll, kk) * hang_weight *
1182 hhang_weight * W;
1183 }
1184 }
1185 }
1186 }
1187 }
1188
1189
1190 // Get Jacobian too?
1191 if ((flag == 1) || (flag == 3))
1192 {
1193 // Offset for faster access in general stress loop
1194 const unsigned offset1 = d_G_dX.offset(l, k, i);
1195
1196 // Default setting for non-hanging node
1197 unsigned nn_master = 1;
1198 double hhang_weight = 1.0;
1199
1200 // Loop over the nodes of the element again
1201 for (unsigned ll = 0; ll < n_node; ll++)
1202 {
1203 // Get pointer to local node ll
1204 Node* llocal_node_pt = node_pt(ll);
1205
1206 // Cache hang status
1207 bool iis_hanging = llocal_node_pt->is_hanging();
1208
1209 // If the node is a hanging node
1210 if (iis_hanging)
1211 {
1212 nn_master = llocal_node_pt->hanging_pt()->nmaster();
1213 }
1214 // Otherwise the node is its own master
1215 else
1216 {
1217 nn_master = 1;
1218 }
1219
1220
1221 // Storage for local unknown numbers at node indexed by
1222 // type and direction
1223 DenseMatrix<int> position_local_unk_at_node(n_position_type,
1224 DIM);
1225
1226 // Loop over the master nodes
1227 for (unsigned mm = 0; mm < nn_master; mm++)
1228 {
1229 if (iis_hanging)
1230 {
1231 // Find the unknown numbers
1232 position_local_unk_at_node = local_position_hang_eqn(
1233 llocal_node_pt->hanging_pt()->master_node_pt(mm));
1234
1235 // Find the hanging node weight
1236 hhang_weight =
1237 llocal_node_pt->hanging_pt()->master_weight(mm);
1238 }
1239 else
1240 {
1241 // Loop of types of dofs
1242 for (unsigned kk = 0; kk < n_position_type; kk++)
1243 {
1244 // Loop over the displacement components
1245 for (unsigned ii = 0; ii < DIM; ii++)
1246 {
1247 position_local_unk_at_node(kk, ii) =
1248 position_local_eqn(ll, kk, ii);
1249 }
1250 }
1251
1252 // Hang weight is one
1253 hhang_weight = 1.0;
1254 }
1255
1256
1257 // Loop of types of dofs again
1258 for (unsigned kk = 0; kk < n_position_type; kk++)
1259 {
1260 // Loop over the displacement components again
1261 for (unsigned ii = 0; ii < DIM; ii++)
1262 {
1263 // Get the number of the unknown
1264 int local_unknown =
1265 position_local_unk_at_node(kk, ii);
1266
1267 /*IF it's not a boundary condition*/
1268 if (local_unknown >= 0)
1269 {
1270 // Offset for faster access in general stress loop
1271 const unsigned offset2 = d_G_dX.offset(ll, kk, ii);
1272 const unsigned offset4 = dpsidxi.offset(ll, kk);
1273
1274
1275 // General stress term
1276 //--------------------
1277 double sum = 0.0;
1278 unsigned count1 = offset1;
1279 for (unsigned a = 0; a < DIM; a++)
1280 {
1281 // Bump up direct access because we're only
1282 // accessing upper triangle
1283 count1 += a;
1284 for (unsigned b = a; b < DIM; b++)
1285 {
1286 double factor =
1287 d_G_dX.raw_direct_access(count1);
1288 if (a == b) factor *= 0.5;
1289
1290 // Offset for faster access
1291 unsigned offset3 = d_stress_dG.offset(a, b);
1292 unsigned count2 = offset2;
1293 unsigned count3 = offset3;
1294
1295 for (unsigned aa = 0; aa < DIM; aa++)
1296 {
1297 // Bump up direct access because we're only
1298 // accessing upper triangle
1299 count2 += aa;
1300 count3 += aa;
1301
1302 // Only upper half of derivatives w.r.t.
1303 // symm tensor
1304 for (unsigned bb = aa; bb < DIM; bb++)
1305 {
1306 sum +=
1307 factor *
1308 d_stress_dG.raw_direct_access(count3) *
1309 d_G_dX.raw_direct_access(count2);
1310 ++count2;
1311 ++count3;
1312 }
1313 }
1314 ++count1;
1315 }
1316 }
1317
1318 // Multiply by weight and add contribution
1319 // (Add directly because this bit is nonsymmetric)
1320 jacobian(local_eqn, local_unknown) +=
1321 sum * W * hang_weight * hhang_weight;
1322
1323 // Only upper triangle (no separate test for bc as
1324 // local_eqn is already nonnegative)
1325 if ((i == ii) && (local_unknown >= local_eqn))
1326 {
1327 // Initialise contribution
1328 double sum = 0.0;
1329
1330 // Inertia term
1331 sum += lambda_sq * time_factor * psi(ll, kk) *
1332 psi(l, k);
1333
1334 // Stress term
1335 unsigned count4 = offset4;
1336 for (unsigned a = 0; a < DIM; a++)
1337 {
1338 // Cache term
1339 const double factor =
1340 dpsidxi.raw_direct_access(count4); // ll ,kk
1341 ++count4;
1342
1343 unsigned count5 = offset5;
1344 for (unsigned b = 0; b < DIM; b++)
1345 {
1346 sum +=
1347 sigma(a, b) * factor *
1348 dpsidxi.raw_direct_access(count5); // l ,k
1349 ++count5;
1350 }
1351 }
1352
1353 // Multiply by weights to form contribution
1354 double sym_entry =
1355 sum * W * hang_weight * hhang_weight;
1356 // Add contribution to jacobian
1357 jacobian(local_eqn, local_unknown) += sym_entry;
1358 // Add to lower triangular entries
1359 if (local_eqn != local_unknown)
1360 {
1361 jacobian(local_unknown, local_eqn) += sym_entry;
1362 }
1363 }
1364 } // End of if not boundary condition
1365 }
1366 }
1367 }
1368 }
1369 }
1370
1371 // Can add in the pressure jacobian terms
1372 if (flag > 0)
1373 {
1374 // Loop over the pressure nodes
1375 for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1376 {
1377 unsigned n_master2 = 1;
1378 double hang_weight2 = 1.0;
1379 HangInfo* hang_info2_pt = 0;
1380
1381 bool is_hanging2 = solid_pressure_dof_is_hanging[l2];
1382 if (is_hanging2)
1383 {
1384 // Get the HangInfo object associated with the
1385 // hanging solid pressure
1386 hang_info2_pt =
1387 solid_pressure_node_pt(l2)->hanging_pt(solid_p_index);
1388
1389 n_master2 = hang_info2_pt->nmaster();
1390 }
1391 else
1392 {
1393 n_master2 = 1;
1394 }
1395
1396 // Loop over all the master nodes
1397 for (unsigned m2 = 0; m2 < n_master2; m2++)
1398 {
1399 if (is_hanging2)
1400 {
1401 // Get the equation numbers at the master node
1402 local_unknown = local_hang_eqn(
1403 hang_info2_pt->master_node_pt(m2), solid_p_index);
1404
1405 // Find the hanging node weight at the node
1406 hang_weight2 = hang_info2_pt->master_weight(m2);
1407 }
1408 else
1409 {
1410 local_unknown = this->solid_p_local_eqn(l2);
1411 hang_weight2 = 1.0;
1412 }
1413
1414 // If it's not a boundary condition
1415 if (local_unknown >= 0)
1416 {
1417 // Add the pressure terms to the jacobian
1418 for (unsigned a = 0; a < DIM; a++)
1419 {
1420 for (unsigned b = 0; b < DIM; b++)
1421 {
1422 jacobian(local_eqn, local_unknown) -=
1423 psisp[l2] * Gup(a, b) * interpolated_G(a, i) *
1424 dpsidxi(l, k, b) * W * hang_weight * hang_weight2;
1425 }
1426 }
1427 }
1428 } // End of loop over master nodes
1429 } // End of loop over pressure dofs
1430 } // End of Jacobian terms
1431
1432 } // End of if not boundary condition
1433 }
1434 }
1435 } // End of loop of over master nodes
1436
1437 } // End of loop over nodes
1438
1439 //==============CONSTRAINT EQUATIONS FOR PRESSURE=====================
1440
1441 // Now loop over the pressure degrees of freedom
1442 for (unsigned l = 0; l < n_solid_pres; l++)
1443 {
1444 bool is_hanging = solid_pressure_dof_is_hanging[l];
1445
1446 unsigned n_master = 1;
1447 double hang_weight = 1.0;
1448 HangInfo* hang_info_pt = 0;
1449
1450 // If the node is a hanging node
1451 if (is_hanging)
1452 {
1453 // Get a pointer to the HangInfo object associated with the
1454 // solid pressure (stored at solid_p_index)
1455 hang_info_pt = solid_pressure_node_pt(l)->hanging_pt(solid_p_index);
1456
1457 // Number of master nodes
1458 n_master = hang_info_pt->nmaster();
1459 }
1460 // Otherwise the node is its own master
1461 else
1462 {
1463 n_master = 1;
1464 }
1465
1466 // Loop over all the master nodes
1467 // Note that the pressure is stored at the inded solid_p_index
1468 for (unsigned m = 0; m < n_master; m++)
1469 {
1470 if (is_hanging)
1471 {
1472 // Get the equation numbers at the master node
1473 local_eqn =
1474 local_hang_eqn(hang_info_pt->master_node_pt(m), solid_p_index);
1475
1476 // Find the hanging node weight at the node
1477 hang_weight = hang_info_pt->master_weight(m);
1478 }
1479 else
1480 {
1481 local_eqn = this->solid_p_local_eqn(l);
1482 }
1483
1484 // Pinned (unlikely, actually) or real dof?
1485 if (local_eqn >= 0)
1486 {
1487 // For true incompressibility we need to conserve volume
1488 // so the determinant of the deformed metric tensor
1489 // needs to be equal to that of the undeformed one, which
1490 // is equal to the volumetric growth factor
1491 if (this->Incompressible)
1492 {
1493 residuals[local_eqn] +=
1494 (detG - gamma) * psisp[l] * W * hang_weight;
1495
1496 // Get Jacobian too?
1497 if ((flag == 1) || (flag == 3))
1498 {
1499 // Default setting for non-hanging node
1500 unsigned nn_master = 1;
1501 double hhang_weight = 1.0;
1502
1503 // Loop over the nodes of the element again
1504 for (unsigned ll = 0; ll < n_node; ll++)
1505 {
1506 // Get pointer to local node ll
1507 Node* llocal_node_pt = node_pt(ll);
1508
1509 // Cache hang status
1510 bool iis_hanging = llocal_node_pt->is_hanging();
1511
1512 // If the node is a hanging node
1513 if (iis_hanging)
1514 {
1515 nn_master = llocal_node_pt->hanging_pt()->nmaster();
1516 }
1517 // Otherwise the node is its own master
1518 else
1519 {
1520 nn_master = 1;
1521 }
1522
1523 // Storage for local unknown numbers at node indexed by
1524 // type and direction
1525 DenseMatrix<int> position_local_unk_at_node(n_position_type,
1526 DIM);
1527
1528 // Loop over the master nodes
1529 for (unsigned mm = 0; mm < nn_master; mm++)
1530 {
1531 if (iis_hanging)
1532 {
1533 // Find the unknown numbers
1534 position_local_unk_at_node = local_position_hang_eqn(
1535 llocal_node_pt->hanging_pt()->master_node_pt(mm));
1536
1537 // Find the hanging node weight
1538 hhang_weight =
1539 llocal_node_pt->hanging_pt()->master_weight(mm);
1540 }
1541 else
1542 {
1543 // Loop of types of dofs
1544 for (unsigned kk = 0; kk < n_position_type; kk++)
1545 {
1546 // Loop over the displacement components
1547 for (unsigned ii = 0; ii < DIM; ii++)
1548 {
1549 position_local_unk_at_node(kk, ii) =
1550 position_local_eqn(ll, kk, ii);
1551 }
1552 }
1553
1554 // Hang weight is one
1555 hhang_weight = 1.0;
1556 }
1557
1558
1559 // Loop of types of dofs again
1560 for (unsigned kk = 0; kk < n_position_type; kk++)
1561 {
1562 // Loop over the displacement components again
1563 for (unsigned ii = 0; ii < DIM; ii++)
1564 {
1565 // Get the number of the unknown
1566 int local_unknown = position_local_unk_at_node(kk, ii);
1567
1568 /*IF it's not a boundary condition*/
1569 if (local_unknown >= 0)
1570 {
1571 // Offset for faster access
1572 const unsigned offset = d_G_dX.offset(ll, kk, ii);
1573
1574 // General stress term
1575 double sum = 0.0;
1576 unsigned count = offset;
1577 for (unsigned aa = 0; aa < DIM; aa++)
1578 {
1579 // Bump up direct access because we're only
1580 // accessing upper triangle
1581 count += aa;
1582
1583 // Only upper half
1584 for (unsigned bb = aa; bb < DIM; bb++)
1585 {
1586 sum += d_detG_dG(aa, bb) *
1587 d_G_dX.raw_direct_access(count) * psisp(l);
1588 ++count;
1589 }
1590 }
1591 jacobian(local_eqn, local_unknown) +=
1592 sum * W * hang_weight * hhang_weight;
1593 }
1594 }
1595 }
1596 }
1597 }
1598
1599 // No Jacobian terms due to pressure since it does not feature
1600 // in the incompressibility constraint
1601 }
1602 }
1603 // Nearly incompressible: (Neg.) pressure given by product of
1604 // bulk modulus and generalised dilatation
1605 else
1606 {
1607 residuals[local_eqn] +=
1608 (inv_kappa * interpolated_solid_p + gen_dil) * psisp[l] * W *
1609 hang_weight;
1610
1611 // Add in the jacobian terms
1612 if ((flag == 1) || (flag == 3))
1613 {
1614 // Default setting for non-hanging node
1615 unsigned nn_master = 1;
1616 double hhang_weight = 1.0;
1617
1618 // Loop over the nodes of the element again
1619 for (unsigned ll = 0; ll < n_node; ll++)
1620 {
1621 // Get pointer to local node ll
1622 Node* llocal_node_pt = node_pt(ll);
1623
1624 // Cache hang status
1625 bool iis_hanging = llocal_node_pt->is_hanging();
1626
1627 // If the node is a hanging node
1628 if (iis_hanging)
1629 {
1630 nn_master = llocal_node_pt->hanging_pt()->nmaster();
1631 }
1632 // Otherwise the node is its own master
1633 else
1634 {
1635 nn_master = 1;
1636 }
1637
1638
1639 // Storage for local unknown numbers at node indexed by
1640 // type and direction
1641 DenseMatrix<int> position_local_unk_at_node(n_position_type,
1642 DIM);
1643
1644 // Loop over the master nodes
1645 for (unsigned mm = 0; mm < nn_master; mm++)
1646 {
1647 if (iis_hanging)
1648 {
1649 // Find the unknown numbers
1650 position_local_unk_at_node = local_position_hang_eqn(
1651 llocal_node_pt->hanging_pt()->master_node_pt(mm));
1652
1653 // Find the hanging node weight
1654 hhang_weight =
1655 llocal_node_pt->hanging_pt()->master_weight(mm);
1656 }
1657 else
1658 {
1659 // Loop of types of dofs
1660 for (unsigned kk = 0; kk < n_position_type; kk++)
1661 {
1662 // Loop over the displacement components
1663 for (unsigned ii = 0; ii < DIM; ii++)
1664 {
1665 position_local_unk_at_node(kk, ii) =
1666 position_local_eqn(ll, kk, ii);
1667 }
1668 }
1669
1670 // Hang weight is one
1671 hhang_weight = 1.0;
1672 }
1673
1674
1675 // Loop of types of dofs again
1676 for (unsigned kk = 0; kk < n_position_type; kk++)
1677 {
1678 // Loop over the displacement components again
1679 for (unsigned ii = 0; ii < DIM; ii++)
1680 {
1681 // Get the number of the unknown
1682 int local_unknown = position_local_unk_at_node(kk, ii);
1683
1684 /*IF it's not a boundary condition*/
1685 if (local_unknown >= 0)
1686 {
1687 // Offset for faster access
1688 const unsigned offset = d_G_dX.offset(ll, kk, ii);
1689
1690 // General stress term
1691 double sum = 0.0;
1692 unsigned count = offset;
1693 for (unsigned aa = 0; aa < DIM; aa++)
1694 {
1695 // Bump up direct access because we're only
1696 // accessing upper triangle
1697 count += aa;
1698
1699 // Only upper half
1700 for (unsigned bb = aa; bb < DIM; bb++)
1701 {
1702 sum += d_gen_dil_dG(aa, bb) *
1703 d_G_dX.raw_direct_access(count) * psisp(l);
1704 ++count;
1705 }
1706 }
1707 jacobian(local_eqn, local_unknown) +=
1708 sum * W * hang_weight * hhang_weight;
1709 }
1710 }
1711 }
1712 }
1713 }
1714 }
1715
1716
1717 // Add in the pressure jacobian terms
1718 if (flag > 0)
1719 {
1720 // Loop over the pressure nodes again
1721 for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1722 {
1723 bool is_hanging2 = solid_pressure_dof_is_hanging[l2];
1724
1725 unsigned n_master2 = 1;
1726 double hang_weight2 = 1.0;
1727 HangInfo* hang_info2_pt = 0;
1728
1729 if (is_hanging2)
1730 {
1731 // Get pointer to hang info object
1732 // Note that the pressure is stored at
1733 // the index solid_p_index
1734 hang_info2_pt =
1735 solid_pressure_node_pt(l2)->hanging_pt(solid_p_index);
1736
1737 n_master2 = hang_info2_pt->nmaster();
1738 }
1739 else
1740 {
1741 n_master2 = 1;
1742 }
1743
1744 // Loop over all the master nodes
1745 for (unsigned m2 = 0; m2 < n_master2; m2++)
1746 {
1747 if (is_hanging2)
1748 {
1749 // Get the equation numbers at the master node
1750 local_unknown = local_hang_eqn(
1751 hang_info2_pt->master_node_pt(m2), solid_p_index);
1752
1753 // Find the hanging node weight at the node
1754 hang_weight2 = hang_info2_pt->master_weight(m2);
1755 }
1756 else
1757 {
1758 local_unknown = this->solid_p_local_eqn(l2);
1759 hang_weight2 = 1.0;
1760 }
1761
1762 // If it's not a boundary condition
1763 if (local_unknown >= 0)
1764 {
1765 jacobian(local_eqn, local_unknown) +=
1766 inv_kappa * psisp[l2] * psisp[l] * W * hang_weight *
1767 hang_weight2;
1768 }
1769
1770 } // End of loop over master nodes
1771 } // End of loop over pressure dofs
1772 } // End of pressure Jacobian
1773
1774
1775 } // End of nearly incompressible case
1776 } // End of if not boundary condition
1777 } // End of loop over master nodes
1778 } // End of loop over pressure dofs
1779 } // End of loop over integration points
1780 }
1781
1782
1783 //====================================================================
1784 /// Forcing building of required templates
1785 //====================================================================
1786 template class RefineablePVDEquations<2>;
1787 template class RefineablePVDEquations<3>;
1788
1791
1792} // 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 that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
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
Class for Refineable solid mechanics elements in near-incompressible/ incompressible formulation,...
void fill_in_generic_residual_contribution_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
void get_mass_matrix_diagonal(Vector< double > &mass_diag)
Compute the diagonal of the displacement mass matrix for LSC preconditioner.
Class for Refineable PVD equations.
void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Call the residuals including hanging node cases.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
ConstitutiveLaw * Constitutive_law_pt
Create constitutive law (for smoothing by nonlinear elasticity)
Definition: mesh_smooth.h:55
//////////////////////////////////////////////////////////////////// ////////////////////////////////...