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