RefineableQDPVDElement.h
Go to the documentation of this file.
1 //
2 // Created by Thomas Weinhart on 29/04/2022.
3 //
4 
5 #ifndef MERCURYDPM_REFINEABLEQDPVDELEMENT_H
6 #define MERCURYDPM_REFINEABLEQDPVDELEMENT_H
7 #include "refineable_solid_elements.h"
8 
9 namespace oomph
10 {
11 
17 template<unsigned DIM, unsigned NNODE_1D>
18 class RefineableQDPVDElement : public RefineableQPVDElement<DIM,NNODE_1D>
19 {
20 public:
21 
23 
25  Vector<double>& residuals,
26  DenseMatrix<double>& jacobian,
27  const unsigned& flag )
28 
29  {
30 #ifdef PARANOID
31  // Check if the constitutive equation requires the explicit imposition of an
32  // incompressibility constraint
33  if (this->Constitutive_law_pt->requires_incompressibility_constraint())
34  {
35  throw OomphLibError("RefineablePVDEquations cannot be used with "
36  "incompressible constitutive laws.",
37  OOMPH_CURRENT_FUNCTION,
38  OOMPH_EXCEPTION_LOCATION);
39  }
40 #endif
41 
42  // Simply set up initial condition?
43  if (this->Solid_ic_pt != 0)
44  {
45  this->get_residuals_for_solid_ic(residuals);
46  return;
47  }
48 
49  // Find out how many nodes there are
50  const unsigned n_node = this->nnode();
51 
52  // Find out how many positional dofs there are
53  const unsigned n_position_type = this->nnodal_position_type();
54 
55  // Integers to store local equation numbers
56  int local_eqn = 0;
57 
58  // Timescale ratio (non-dim density)
59  double lambda_sq = this->lambda_sq();
60 
61  // Time factor
62  double time_factor = 0.0;
63  double time_factor1 = 0.0;
64  if(lambda_sq>0)
65  {
66  time_factor = this->node_pt(0)->position_time_stepper_pt()->weight(2, 0);
67  time_factor1 = this->node_pt(0)->position_time_stepper_pt()->weight(1, 0);
68  }
69 
70  // Set up memory for the shape functions
71  Shape psi(n_node, n_position_type);
72  DShape dpsidxi(n_node, n_position_type, DIM);
73 
74  // Set the value of n_intpt -- the number of integration points
75  const unsigned n_intpt = this->integral_pt()->nweight();
76 
77  // Set the vector to hold the local coordinates in the element
78  Vector<double> s(DIM);
79 
80  // Loop over the integration points
81  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
82  {
83  // Assign the values of s
84  for (unsigned i = 0; i < DIM; ++i)
85  {
86  s[i] = this->integral_pt()->knot(ipt, i);
87  }
88 
89  // Get the integral weight
90  double w = this->integral_pt()->weight(ipt);
91 
92  // Call the derivatives of the shape functions (and get Jacobian)
93  double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
94 
95  // Calculate interpolated values of the derivative of global position
96  // wrt lagrangian coordinates
97  DenseMatrix<double> interpolated_G(DIM, DIM, 0.0);
98 
99  // Setup memory for accelerations
100  Vector<double> accel(DIM, 0.0);
101 
102  // Storage for Lagrangian coordinates (initialised to zero)
103  Vector<double> interpolated_xi(DIM, 0.0);
104 
105  // Calculate displacements and derivatives and lagrangian coordinates
106  for (unsigned l = 0; l < n_node; l++)
107  {
108  // Loop over positional dofs
109  for (unsigned k = 0; k < n_position_type; k++)
110  {
111  double psi_ = psi(l, k);
112  // Loop over displacement components (deformed position)
113  for (unsigned i = 0; i < DIM; i++)
114  {
115  // Calculate the Lagrangian coordinates and the accelerations
116  interpolated_xi[i] += this->lagrangian_position_gen(l, k, i) * psi_;
117 
118  // Only compute accelerations if inertia is switched on
119  // otherwise the timestepper might not be able to
120  // work out dx_gen_dt(2,...)
121  if ((lambda_sq > 0.0) && (this->Unsteady))
122  {
123  accel[i] += this->dnodal_position_gen_dt(2, l, k, i) * psi_;
124  accel[i] += this->dissipation_ *
125  this->dnodal_position_gen_dt(1, l, k, i) * psi_; //TW
126  }
127 
128  // Loop over derivative directions
129  for (unsigned j = 0; j < DIM; j++)
130  {
131  interpolated_G(j, i) +=
132  this->nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
133  }
134  }
135  }
136  }
137 
138  // Get isotropic growth factor
139  double gamma = 1.0;
140  this->get_isotropic_growth(ipt, s, interpolated_xi, gamma);
141 
142  // Get body force at current time
143  Vector<double> b(DIM);
144  this->body_force(interpolated_xi, b);
145 
146  // We use Cartesian coordinates as the reference coordinate
147  // system. In this case the undeformed metric tensor is always
148  // the identity matrix -- stretched by the isotropic growth
149  double diag_entry = pow(gamma, 2.0 / double(DIM));
150  DenseMatrix<double> g(DIM);
151  for (unsigned i = 0; i < DIM; i++)
152  {
153  for (unsigned j = 0; j < DIM; j++)
154  {
155  if (i == j)
156  {
157  g(i, j) = diag_entry;
158  }
159  else
160  {
161  g(i, j) = 0.0;
162  }
163  }
164  }
165 
166  // Premultiply the undeformed volume ratio (from the isotropic
167  // growth), the weights and the Jacobian
168  double W = gamma * w * J;
169 
170  // Declare and calculate the deformed metric tensor
171  DenseMatrix<double> G(DIM);
172 
173  // Assign values of G
174  for (unsigned i = 0; i < DIM; i++)
175  {
176  // Do upper half of matrix
177  for (unsigned j = i; j < DIM; j++)
178  {
179  // Initialise G(i,j) to zero
180  G(i, j) = 0.0;
181  // Now calculate the dot product
182  for (unsigned k = 0; k < DIM; k++)
183  {
184  G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
185  }
186  }
187  // Matrix is symmetric so just copy lower half
188  for (unsigned j = 0; j < i; j++)
189  {
190  G(i, j) = G(j, i);
191  }
192  }
193 
194  // Now calculate the stress tensor from the constitutive law
195  DenseMatrix<double> sigma(DIM);
196  this->get_stress(g, G, sigma);
197 
198  // Get stress derivative by FD only needed for Jacobian
199  //-----------------------------------------------------
200 
201  // Stress derivative
202  RankFourTensor<double> d_stress_dG(DIM, DIM, DIM, DIM, 0.0);
203  // Derivative of metric tensor w.r.t. to nodal coords
204  RankFiveTensor<double> d_G_dX(
205  n_node, n_position_type, DIM, DIM, DIM, 0.0);
206 
207  // Get Jacobian too?
208  if (flag == 1)
209  {
210  // Derivative of metric tensor w.r.t. to discrete positional dofs
211  // NOTE: Since G is symmetric we only compute the upper triangle
212  // and DO NOT copy the entries across. Subsequent computations
213  // must (and, in fact, do) therefore only operate with upper
214  // triangular entries
215  for (unsigned ll = 0; ll < n_node; ll++)
216  {
217  for (unsigned kk = 0; kk < n_position_type; kk++)
218  {
219  for (unsigned ii = 0; ii < DIM; ii++)
220  {
221  for (unsigned aa = 0; aa < DIM; aa++)
222  {
223  for (unsigned bb = aa; bb < DIM; bb++)
224  {
225  d_G_dX(ll, kk, ii, aa, bb) =
226  interpolated_G(aa, ii) * dpsidxi(ll, kk, bb) +
227  interpolated_G(bb, ii) * dpsidxi(ll, kk, aa);
228  }
229  }
230  }
231  }
232  }
233 
234  // Get the "upper triangular"
235  // entries of the derivatives of the stress tensor with
236  // respect to G
237  this->get_d_stress_dG_upper(g, G, sigma, d_stress_dG);
238  }
239 
240 
241  // Add pre-stress
242  for (unsigned i = 0; i < DIM; i++)
243  {
244  for (unsigned j = 0; j < DIM; j++)
245  {
246  sigma(i, j) += this->prestress(i, j, interpolated_xi);
247  }
248  }
249 
250  //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
251  // DISPLACEMENTS========
252 
253 
254  // Default setting for non-hanging node
255  unsigned n_master = 1;
256  double hang_weight = 1.0;
257 
258  // Loop over the test functions, nodes of the element
259  for (unsigned l = 0; l < n_node; l++)
260  {
261  // Get pointer to local node l
262  Node* local_node_pt = this->node_pt(l);
263 
264  // Cache hang status
265  bool is_hanging = local_node_pt->is_hanging();
266 
267  // If the node is a hanging node
268  if (is_hanging)
269  {
270  n_master = local_node_pt->hanging_pt()->nmaster();
271  }
272  // Otherwise the node is its own master
273  else
274  {
275  n_master = 1;
276  }
277 
278 
279  // Storage for local equation numbers at node indexed by
280  // type and direction
281  DenseMatrix<int> position_local_eqn_at_node(n_position_type, DIM);
282 
283  // Loop over the master nodes
284  for (unsigned m = 0; m < n_master; m++)
285  {
286  if (is_hanging)
287  {
288  // Find the equation numbers
289  position_local_eqn_at_node = this->local_position_hang_eqn(
290  local_node_pt->hanging_pt()->master_node_pt(m));
291 
292  // Find the hanging node weight
293  hang_weight = local_node_pt->hanging_pt()->master_weight(m);
294  }
295  else
296  {
297  // Loop of types of dofs
298  for (unsigned k = 0; k < n_position_type; k++)
299  {
300  // Loop over the displacement components
301  for (unsigned i = 0; i < DIM; i++)
302  {
303  position_local_eqn_at_node(k, i) = this->position_local_eqn(l, k, i);
304  }
305  }
306 
307  // Hang weight is one
308  hang_weight = 1.0;
309  }
310 
311  // Loop of types of dofs
312  for (unsigned k = 0; k < n_position_type; k++)
313  {
314  // Offset for faster access
315  const unsigned offset5 = dpsidxi.offset(l, k);
316 
317  // Loop over the displacement components
318  for (unsigned i = 0; i < DIM; i++)
319  {
320  local_eqn = position_local_eqn_at_node(k, i);
321 
322  /*IF it's not a boundary condition*/
323  if (local_eqn >= 0)
324  {
325  // Initialise the contribution
326  double sum = 0.0;
327 
328  // Acceleration and body force
329  sum += (lambda_sq * accel[i] - b[i]) * psi(l, k);
330 
331  // Stress term
332  for (unsigned a = 0; a < DIM; a++)
333  {
334  unsigned count = offset5;
335  for (unsigned b = 0; b < DIM; b++)
336  {
337  // Add the stress terms to the residuals
338  sum += sigma(a, b) * interpolated_G(a, i) *
339  dpsidxi.raw_direct_access(count);
340  ++count;
341  }
342  }
343  residuals[local_eqn] += W * sum * hang_weight;
344 
345 
346  // Get Jacobian too?
347  if (flag == 1)
348  {
349  // Offset for faster access in general stress loop
350  const unsigned offset1 = d_G_dX.offset(l, k, i);
351 
352  // Default setting for non-hanging node
353  unsigned nn_master = 1;
354  double hhang_weight = 1.0;
355 
356  // Loop over the nodes of the element again
357  for (unsigned ll = 0; ll < n_node; ll++)
358  {
359  // Get pointer to local node ll
360  Node* llocal_node_pt = this->node_pt(ll);
361 
362  // Cache hang status
363  bool iis_hanging = llocal_node_pt->is_hanging();
364 
365  // If the node is a hanging node
366  if (iis_hanging)
367  {
368  nn_master = llocal_node_pt->hanging_pt()->nmaster();
369  }
370  // Otherwise the node is its own master
371  else
372  {
373  nn_master = 1;
374  }
375 
376 
377  // Storage for local unknown numbers at node indexed by
378  // type and direction
379  DenseMatrix<int> position_local_unk_at_node(n_position_type,
380  DIM);
381 
382  // Loop over the master nodes
383  for (unsigned mm = 0; mm < nn_master; mm++)
384  {
385  if (iis_hanging)
386  {
387  // Find the unknown numbers
388  position_local_unk_at_node = this->local_position_hang_eqn(
389  llocal_node_pt->hanging_pt()->master_node_pt(mm));
390 
391  // Find the hanging node weight
392  hhang_weight =
393  llocal_node_pt->hanging_pt()->master_weight(mm);
394  }
395  else
396  {
397  // Loop of types of dofs
398  for (unsigned kk = 0; kk < n_position_type; kk++)
399  {
400  // Loop over the displacement components
401  for (unsigned ii = 0; ii < DIM; ii++)
402  {
403  position_local_unk_at_node(kk, ii) =
404  this->position_local_eqn(ll, kk, ii);
405  }
406  }
407 
408  // Hang weight is one
409  hhang_weight = 1.0;
410  }
411 
412 
413  // Loop of types of dofs again
414  for (unsigned kk = 0; kk < n_position_type; kk++)
415  {
416  // Loop over the displacement components again
417  for (unsigned ii = 0; ii < DIM; ii++)
418  {
419  // Get the number of the unknown
420  int local_unknown =
421  position_local_unk_at_node(kk, ii);
422 
423 
424  /*IF it's not a boundary condition*/
425  if (local_unknown >= 0)
426  {
427  // Offset for faster access in general stress loop
428  const unsigned offset2 = d_G_dX.offset(ll, kk, ii);
429  const unsigned offset4 = dpsidxi.offset(ll, kk);
430 
431 
432  // General stress term
433  //--------------------
434  double sum = 0.0;
435  unsigned count1 = offset1;
436  for (unsigned a = 0; a < DIM; a++)
437  {
438  // Bump up direct access because we're only
439  // accessing upper triangle
440  count1 += a;
441  for (unsigned b = a; b < DIM; b++)
442  {
443  double factor =
444  d_G_dX.raw_direct_access(count1);
445  if (a == b) factor *= 0.5;
446 
447  // Offset for faster access
448  unsigned offset3 = d_stress_dG.offset(a, b);
449  unsigned count2 = offset2;
450  unsigned count3 = offset3;
451 
452  for (unsigned aa = 0; aa < DIM; aa++)
453  {
454  // Bump up direct access because we're only
455  // accessing upper triangle
456  count2 += aa;
457  count3 += aa;
458 
459  // Only upper half of derivatives w.r.t.
460  // symm tensor
461  for (unsigned bb = aa; bb < DIM; bb++)
462  {
463  sum +=
464  factor *
465  d_stress_dG.raw_direct_access(count3) *
466  d_G_dX.raw_direct_access(count2);
467  ++count2;
468  ++count3;
469  }
470  }
471  ++count1;
472  }
473  }
474 
475  // Multiply by weight and add contribution
476  // (Add directly because this bit is nonsymmetric)
477  jacobian(local_eqn, local_unknown) +=
478  sum * W * hang_weight * hhang_weight;
479 
480  // Only upper triangle (no separate test for bc as
481  // local_eqn is already nonnegative)
482  if ((i == ii) && (local_unknown >= local_eqn))
483  {
484  // Initialise the contribution
485  double sum = 0.0;
486 
487  // Inertia term
488  sum += lambda_sq * time_factor * psi(ll, kk) *
489  psi(l, k);
490  sum += lambda_sq * this->dissipation_
491  * time_factor1 * psi(ll,kk) * psi(l,k); //TW
492 
493  // Stress term
494  unsigned count4 = offset4;
495  for (unsigned a = 0; a < DIM; a++)
496  {
497  // Cache term
498  const double factor =
499  dpsidxi.raw_direct_access(count4); // ll ,kk
500  ++count4;
501 
502  unsigned count5 = offset5;
503  for (unsigned b = 0; b < DIM; b++)
504  {
505  sum +=
506  sigma(a, b) * factor *
507  dpsidxi.raw_direct_access(count5); // l ,k
508  ++count5;
509  }
510  }
511 
512  // Multiply by weights to form contribution
513  double sym_entry =
514  sum * W * hang_weight * hhang_weight;
515  // Add contribution to jacobian
516  jacobian(local_eqn, local_unknown) += sym_entry;
517  // Add to lower triangular entries
518  if (local_eqn != local_unknown)
519  {
520  jacobian(local_unknown, local_eqn) += sym_entry;
521  }
522  }
523  } // End of if not boundary condition
524  }
525  }
526  }
527  }
528  }
529 
530  } // End of if not boundary condition
531 
532  } // End of loop over coordinate directions
533  } // End of loop over type of dof
534  } // End of loop over master nodes
535  } // End of loop over nodes
536  } // End of loop over integration points
537  }
538 
539  double dissipation_ = 0.0; //TW
540 
541 public:
542 
543  void setDissipation(double dissipation) {
544  dissipation_ = dissipation;
545  }
546 }; // end class
547 
548 //==============================================================
550 //==============================================================
551 template<unsigned NNODE_1D>
553  : public virtual SolidQElement<1, NNODE_1D>
554 {
555 public:
556  // Make sure that we call the constructor of the SolidQElement
557  // Only the Intel compiler seems to need this!
558  FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
559 };
560 
561 //==============================================================
563 //==============================================================
564 template<unsigned NNODE_1D>
566  : public virtual PointElement
567 {
568 public:
569  // Make sure that we call the constructor of the SolidQElement
570  // Only the Intel compiler seems to need this!
572 };
573 
574 
575 //==============================================================
577 //==============================================================
578 template<unsigned NNODE_1D>
580  : public virtual SolidQElement<2, NNODE_1D>
581 {
582 public:
583  // Make sure that we call the constructor of the SolidQElement
584  // Only the Intel compiler seems to need this!
585  FaceGeometry() : SolidQElement<2, NNODE_1D>() {}
586 };
587 
588 //==============================================================
590 //==============================================================
591 template<unsigned NNODE_1D>
593  : public virtual SolidQElement<1, NNODE_1D>
594 {
595 public:
596  // Make sure that we call the constructor of the SolidQElement
597  // Only the Intel compiler seems to need this!
598  FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
599 };
600 
601 } // end namespace oomph
602 
603 #endif//MERCURYDPM_REFINEABLEQDPVDELEMENT_H
FaceGeometry()
Definition: RefineableQDPVDElement.h:558
FaceGeometry()
Definition: RefineableQDPVDElement.h:585
Definition: RefineableQDPVDElement.h:19
void setDissipation(double dissipation)
Definition: RefineableQDPVDElement.h:543
RefineableQDPVDElement()
Definition: RefineableQDPVDElement.h:22
double dissipation_
Definition: RefineableQDPVDElement.h:539
void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Definition: RefineableQDPVDElement.h:24
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
Definition: TwenteMeshGluing.cpp:61
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:137
Definition: AnisotropicHookean.h:31