axisym_poroelasticity_elements.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
27 
28 namespace oomph
29 {
30  //===================================================================
31  /// Static default value for Young's modulus (1.0 -- for natural
32  /// scaling, i.e. all stresses have been non-dimensionalised by
33  /// the same reference Young's modulus. Setting the "non-dimensional"
34  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
35  /// to a number larger than one means that the material is stiffer
36  /// than assumed in the non-dimensionalisation.
37  //===================================================================
39  1.0;
40 
41  //========================================================================
42  /// Static default value for timescale ratio (1.0)
43  //===================================================================
45 
46  //===================================================================
47  /// Static default value for the density ratio (fluid to solid)
48  //===================================================================
50 
51  //===================================================================
52  /// Static default value for permeability (1.0 for natural scaling
53  /// i.e. timescale is given by L^2/(k^* E)
54  //===================================================================
56 
57 
58  //===================================================================
59  /// Static default value for ratio of the material's actual
60  /// permeability to the permeability used in the non-dimensionalisastion
61  /// of the equations
62  //===================================================================
64  1.0;
65 
66  //===================================================================
67  /// Static default value for alpha, the Biot parameter
68  //===================================================================
70 
71  //===================================================================
72  /// Static default value for the porosity
73  //===================================================================
75 
76  //======================================================================
77  /// Performs a div-conserving transformation of the vector basis
78  /// functions from the reference element to the actual element
79  //======================================================================
81  const Vector<double>& s,
82  const Shape& q_basis_local,
83  Shape& psi,
84  DShape& dpsi,
85  Shape& q_basis) const
86  {
87  // Call the (geometric) shape functions and their derivatives
88  this->dshape_local(s, psi, dpsi);
89 
90  // Storage for the (geometric) jacobian and its inverse
91  DenseMatrix<double> jacobian(2), inverse_jacobian(2);
92 
93  // Get the jacobian of the geometric mapping and its determinant
94  double det = local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
95 
96  // Transform the derivative of the geometric basis so that it's w.r.t.
97  // global coordinates
98  this->transform_derivatives(inverse_jacobian, dpsi);
99 
100  // Get the number of computational basis vectors
101  const unsigned n_q_basis = this->nq_basis();
102 
103  // Loop over the basis vectors
104  for (unsigned l = 0; l < n_q_basis; l++)
105  {
106  // Loop over the spatial components
107  for (unsigned i = 0; i < 2; i++)
108  {
109  // Zero the basis
110  q_basis(l, i) = 0.0;
111  }
112  }
113 
114  // Loop over the spatial components
115  for (unsigned i = 0; i < 2; i++)
116  {
117  // And again
118  for (unsigned j = 0; j < 2; j++)
119  {
120  // Get the element of the jacobian (must transpose it due to different
121  // conventions) and divide by the determinant
122  double jac_trans = jacobian(j, i) / det;
123 
124  // Loop over the computational basis vectors
125  for (unsigned l = 0; l < n_q_basis; l++)
126  {
127  // Apply the appropriate div-conserving mapping
128  q_basis(l, i) += jac_trans * q_basis_local(l, j);
129  }
130  }
131  }
132 
133  // Scale the basis by the ratio of the length of the edge to the length of
134  // the corresponding edge on the reference element
135  scale_basis(q_basis);
136 
137  return det;
138  }
139 
140 
141  //========================================================================
142  /// Output FE representation of soln:
143  /// x,y,u1,u2,q1,q2,div_q,p at Nplot^2 plot points
144  //========================================================================
146  const unsigned& nplot)
147  {
148  // Vector of local coordinates
149  Vector<double> s(2);
150 
151  // Skeleton velocity
153 
154  // Tecplot header info
155  outfile << tecplot_zone_string(nplot);
156 
157  // Loop over plot points
158  unsigned num_plot_points = nplot_points(nplot);
159 
160  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
161  {
162  // Get local coordinates of plot point
163  get_s_plot(iplot, nplot, s);
164 
165  // Output the components of the position
166  for (unsigned i = 0; i < 2; i++)
167  {
168  outfile << interpolated_x(s, i) << " ";
169  }
170 
171  // Output the components of the FE representation of skeleton displ. at s
172  for (unsigned i = 0; i < 2; i++)
173  {
174  outfile << interpolated_u(s, i) << " "; // soln 0 and 1
175  }
176 
177  // Output the components of the FE representation of q at s
178  for (unsigned i = 0; i < 2; i++)
179  {
180  outfile << interpolated_q(s, i) << " "; // soln 2 and 3
181  }
182 
183  // Output FE representation of div u at s
184  outfile << interpolated_div_q(s) << " "; // soln 4
185 
186  // Output FE representation of p at s
187  outfile << interpolated_p(s) << " "; // soln 5
188 
189  // Skeleton velocity
191  outfile << du_dt[0] << " "; // soln 6
192  outfile << du_dt[1] << " "; // soln 7
193 
194  outfile << std::endl;
195  }
196 
197  // Write tecplot footer (e.g. FE connectivity lists)
198  this->write_tecplot_zone_footer(outfile, nplot);
199  }
200 
201  //============================================================================
202  /// Output FE representation of exact soln at
203  /// Nplot^2 plot points
204  //============================================================================
206  std::ostream& outfile,
207  const unsigned& nplot,
209  {
210  // Vector of local coordinates
211  Vector<double> s(2);
212 
213  // Vector for coordintes
214  Vector<double> x(2);
215 
216  // Tecplot header info
217  outfile << this->tecplot_zone_string(nplot);
218 
219  // Exact solution Vector
220  Vector<double> exact_soln(13);
221 
222  // Loop over plot points
223  unsigned num_plot_points = this->nplot_points(nplot);
224 
225  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
226  {
227  // Get local coordinates of plot point
228  this->get_s_plot(iplot, nplot, s);
229 
230  // Get x position as Vector
231  this->interpolated_x(s, x);
232 
233  // Get exact solution at this point
234  (*exact_soln_pt)(x, exact_soln);
235 
236  // Output
237  for (unsigned i = 0; i < 2; i++)
238  {
239  outfile << x[i] << " ";
240  }
241  for (unsigned i = 0; i < 13; i++)
242  {
243  outfile << exact_soln[i] << " ";
244  }
245  outfile << std::endl;
246  }
247 
248  // Write tecplot footer (e.g. FE connectivity lists)
249  this->write_tecplot_zone_footer(outfile, nplot);
250  }
251 
252  //========================================================================
253  /// Output FE representation of exact soln at
254  /// Nplot^2 plot points. Unsteady version
255  //========================================================================
257  std::ostream& outfile,
258  const unsigned& nplot,
259  const double& time,
261  {
262  // Vector of local coordinates
263  Vector<double> s(2);
264 
265  // Vector for coordintes
266  Vector<double> x(2);
267 
268  // Tecplot header info
269  outfile << this->tecplot_zone_string(nplot);
270 
271  // Exact solution Vector
272  Vector<double> exact_soln(13);
273 
274  // Loop over plot points
275  unsigned num_plot_points = this->nplot_points(nplot);
276 
277  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
278  {
279  // Get local coordinates of plot point
280  this->get_s_plot(iplot, nplot, s);
281 
282  // Get x position as Vector
283  this->interpolated_x(s, x);
284 
285  // Get exact solution at this point
286  (*exact_soln_pt)(time, x, exact_soln);
287 
288  // Output
289  for (unsigned i = 0; i < 2; i++)
290  {
291  outfile << x[i] << " ";
292  }
293  for (unsigned i = 0; i < 13; i++)
294  {
295  outfile << exact_soln[i] << " ";
296  }
297  outfile << std::endl;
298  }
299 
300  // Write tecplot footer (e.g. FE connectivity lists)
301  this->write_tecplot_zone_footer(outfile, nplot);
302  }
303 
304  //========================================================================
305  /// Compute the error between the FE solution and the exact solution
306  /// using the H(div) norm for q and L^2 norm for u and p
307  //========================================================================
309  std::ostream& outfile,
311  Vector<double>& error,
312  Vector<double>& norm)
313  {
314  for (unsigned i = 0; i < 3; i++)
315  {
316  error[i] = 0.0;
317  norm[i] = 0.0;
318  }
319 
320  // Vector of local coordinates
321  Vector<double> s(2);
322 
323  // Vector for coordinates
324  Vector<double> x(2);
325 
326  // Set the value of n_intpt
327  unsigned n_intpt = this->integral_pt()->nweight();
328 
329  outfile << "ZONE" << std::endl;
330 
331  // Exact solution Vector
332  Vector<double> exact_soln(13);
333 
334  // Loop over the integration points
335  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
336  {
337  // Assign values of s
338  for (unsigned i = 0; i < 2; i++)
339  {
340  s[i] = this->integral_pt()->knot(ipt, i);
341  }
342 
343  // Get the integral weight
344  double w = this->integral_pt()->weight(ipt);
345 
346  // Get jacobian of mapping
347  double J = this->J_eulerian(s);
348 
349  // Premultiply the weights and the Jacobian
350  double W = w * J;
351 
352  // Get x position as Vector
353  this->interpolated_x(s, x);
354 
355  // Get exact solution at this point
356  (*exact_soln_pt)(x, exact_soln);
357 
358  // Displacement error
359  for (unsigned i = 0; i < 2; i++)
360  {
361  norm[0] += exact_soln[i] * exact_soln[i] * W;
362  // Error due to q_i
363  error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
364  (exact_soln[i] - this->interpolated_u(s, i)) * W;
365  }
366 
367  // Flux error
368  for (unsigned i = 0; i < 2; i++)
369  {
370  norm[1] += exact_soln[2 + i] * exact_soln[2 + i] * W;
371  // Error due to q_i
372  error[1] += (exact_soln[2 + i] - this->interpolated_q(s, i)) *
373  (exact_soln[2 + i] - this->interpolated_q(s, i)) * W;
374  }
375 
376  // Flux divergence error
377  norm[1] += exact_soln[2 * 2] * exact_soln[2 * 2] * W;
378  error[1] += (exact_soln[2 * 2] - interpolated_div_q(s)) *
379  (exact_soln[2 * 2] - interpolated_div_q(s)) * W;
380 
381  // Pressure error
382  norm[2] += exact_soln[2 * 2 + 1] * exact_soln[2 * 2 + 1] * W;
383  error[2] += (exact_soln[2 * 2 + 1] - this->interpolated_p(s)) *
384  (exact_soln[2 * 2 + 1] - this->interpolated_p(s)) * W;
385 
386  // Output x,y,[z]
387  for (unsigned i = 0; i < 2; i++)
388  {
389  outfile << x[i] << " ";
390  }
391 
392  // Output u_1_error,u_2_error,...
393  for (unsigned i = 0; i < 2; i++)
394  {
395  outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
396  }
397 
398  // Output q_1_error,q_2_error,...
399  for (unsigned i = 0; i < 2; i++)
400  {
401  outfile << exact_soln[2 + i] - this->interpolated_q(s, i) << " ";
402  }
403 
404  // Output p_error
405  outfile << exact_soln[2 * 2 + 1] - this->interpolated_p(s) << " ";
406 
407  outfile << std::endl;
408  }
409  }
410 
411  //========================================================================
412  /// Compute the error between the FE solution and the exact solution
413  /// using the H(div) norm for u and L^2 norm for p. Unsteady version
414  //========================================================================
416  std::ostream& outfile,
418  const double& time,
419  Vector<double>& error,
420  Vector<double>& norm)
421  {
422  for (unsigned i = 0; i < 3; i++)
423  {
424  error[i] = 0.0;
425  norm[i] = 0.0;
426  }
427 
428  // Vector of local coordinates
429  Vector<double> s(2);
430 
431  // Vector for coordinates
432  Vector<double> x(2);
433 
434  // Set the value of n_intpt
435  unsigned n_intpt = this->integral_pt()->nweight();
436 
437  outfile << "ZONE" << std::endl;
438 
439  // Exact solution Vector
440  Vector<double> exact_soln(13);
441 
442  // Loop over the integration points
443  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
444  {
445  // Assign values of s
446  for (unsigned i = 0; i < 2; i++)
447  {
448  s[i] = this->integral_pt()->knot(ipt, i);
449  }
450 
451  // Get the integral weight
452  double w = this->integral_pt()->weight(ipt);
453 
454  // Get jacobian of mapping
455  double J = this->J_eulerian(s);
456 
457  // Premultiply the weights and the Jacobian
458  double W = w * J;
459 
460  // Get x position as Vector
461  this->interpolated_x(s, x);
462 
463  // Get exact solution at this point
464  (*exact_soln_pt)(time, x, exact_soln);
465 
466  // Displacement error
467  for (unsigned i = 0; i < 2; i++)
468  {
469  norm[0] += exact_soln[i] * exact_soln[i] * W;
470  // Error due to q_i
471  error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
472  (exact_soln[i] - this->interpolated_u(s, i)) * W;
473  }
474 
475  // Flux error
476  for (unsigned i = 0; i < 2; i++)
477  {
478  norm[1] += exact_soln[2 + i] * exact_soln[2 + i] * W;
479  // Error due to q_i
480  error[1] += (exact_soln[2 + i] - this->interpolated_q(s, i)) *
481  (exact_soln[2 + i] - this->interpolated_q(s, i)) * W;
482  }
483 
484  // Flux divergence error
485  norm[1] += exact_soln[2 * 2] * exact_soln[2 * 2] * W;
486  error[1] += (exact_soln[2 * 2] - interpolated_div_q(s)) *
487  (exact_soln[2 * 2] - interpolated_div_q(s)) * W;
488 
489  // Pressure error
490  norm[2] += exact_soln[2 * 2 + 1] * exact_soln[2 * 2 + 1] * W;
491  error[2] += (exact_soln[2 * 2 + 1] - this->interpolated_p(s)) *
492  (exact_soln[2 * 2 + 1] - this->interpolated_p(s)) * W;
493 
494  // Output x,y,[z]
495  for (unsigned i = 0; i < 2; i++)
496  {
497  outfile << x[i] << " ";
498  }
499 
500  // Output u_1_error,u_2_error,...
501  for (unsigned i = 0; i < 2; i++)
502  {
503  outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
504  }
505 
506  // Output q_1_error,q_2_error,...
507  for (unsigned i = 0; i < 2; i++)
508  {
509  outfile << exact_soln[2 + i] - this->interpolated_q(s, i) << " ";
510  }
511 
512  // Output p_error
513  outfile << exact_soln[2 * 2 + 1] - this->interpolated_p(s) << " ";
514 
515  outfile << std::endl;
516  }
517  }
518 
519  //========================================================================
520  /// Fill in residuals and, if flag==true, jacobian
521  //========================================================================
524  DenseMatrix<double>& jacobian,
525  bool flag)
526  {
527  // Get the number of geometric nodes, total number of basis functions,
528  // and number of edges basis functions
529  const unsigned n_node = nnode();
530  const unsigned n_q_basis = nq_basis();
531  const unsigned n_q_basis_edge = nq_basis_edge();
532  const unsigned n_p_basis = np_basis();
533 
534  // Storage for the geometric and computational bases and the test functions
535  Shape psi(n_node), u_basis(n_node), u_test(n_node), q_basis(n_q_basis, 2),
536  q_test(n_q_basis, 2), p_basis(n_p_basis), p_test(n_p_basis),
537  div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
538 
539  DShape dpsidx(n_node, 2), du_basis_dx(n_node, 2), du_test_dx(n_node, 2);
540 
541  // Get the number of integration points
542  unsigned n_intpt = integral_pt()->nweight();
543 
544  // Storage for the local coordinates
545  Vector<double> s(2);
546 
547  // Storage for the elasticity source function
548  Vector<double> f_solid(2);
549 
550  // Storage for the source function
551  Vector<double> f_fluid(2);
552 
553  // Storage for the mass source function
554  double mass_source_local = 0.0;
555 
556  // Get elastic parameters
557  double nu_local = this->nu();
558  double youngs_modulus_local = this->youngs_modulus();
559 
560  // Obtain Lame parameters from Young's modulus and Poisson's ratio
561  double lambda = youngs_modulus_local * nu_local / (1.0 + nu_local) /
562  (1.0 - 2.0 * nu_local);
563 
564  double mu = youngs_modulus_local / 2.0 / (1.0 + nu_local);
565 
566  // Storage for Lambda_sq
567  double lambda_sq = this->lambda_sq();
568 
569  // Get the value of permeability
570  double local_permeability = this->permeability();
571 
572  // Ratio of the material's permeability to the permeability used
573  // to non-dimensionalise the equations
574  double local_permeability_ratio = this->permeability_ratio();
575 
576  // Get the value of alpha
577  double alpha = this->alpha();
578 
579  // Get the value of the porosity
580  double porosity = this->porosity();
581 
582  // Get the density ratio
583  double density_ratio = this->density_ratio();
584 
585  // Precompute the ratio of fluid density to combined density
586  double rho_f_over_rho =
587  density_ratio / (1.0 + porosity * (density_ratio - 1.0));
588 
589  // Get continuous time from timestepper of first node
590  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
591 
592  // Local equation and unknown numbers
593  int local_eqn = 0, local_unknown = 0;
594 
595  // Loop over the integration points
596  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
597  {
598  // Find the local coordinates at the integration point
599  for (unsigned i = 0; i < 2; i++)
600  {
601  s[i] = integral_pt()->knot(ipt, i);
602  }
603 
604  // Get the weight of the intetgration point
605  double w = integral_pt()->weight(ipt);
606 
607  // Call the basis functions and test functions and get the
608  // (geometric) jacobian of the current element
609  double J = shape_basis_test_local_at_knot(ipt,
610  psi,
611  dpsidx,
612  u_basis,
613  u_test,
614  du_basis_dx,
615  du_test_dx,
616  q_basis,
617  q_test,
618  p_basis,
619  p_test,
620  div_q_basis_ds,
621  div_q_test_ds);
622 
623  // Storage for interpolated values
626  DenseMatrix<double> interpolated_du_dx(2, 2, 0.0);
627  double interpolated_div_du_dt_dx = 0.0;
628  double interpolated_du_r_dt = 0.0;
629  Vector<double> interpolated_d2u_dt2(2, 0.0);
631  double interpolated_div_q_ds = 0.0;
632  Vector<double> interpolated_dq_dt(2, 0.0);
633  double interpolated_p = 0.0;
634 
635  // loop over geometric basis functions to find interpolated x
636  for (unsigned l = 0; l < n_node; l++)
637  {
638  // Loop over the geometric basis functions
639  for (unsigned i = 0; i < 2; i++)
640  {
641  interpolated_x[i] += nodal_position(l, i) * psi(l);
642  interpolated_d2u_dt2[i] += this->d2u_dt2(l, i) * u_basis(l);
643 
644  // Get the nodal displacements
645  const double u_value =
647  interpolated_u[i] += u_value * u_basis(l);
648 
649  // Loop over derivative directions
650  for (unsigned j = 0; j < 2; j++)
651  {
652  interpolated_du_dx(i, j) += u_value * du_basis_dx(l, j);
653  }
654 
655  // divergence of the time derivative of the solid displacement
656  interpolated_div_du_dt_dx += this->du_dt(l, i) * du_basis_dx(l, i);
657  }
658 
659  // r-component of the solid velocity
660  interpolated_du_r_dt += du_dt(l, 0) * u_basis(l);
661  }
662 
663  // loop over the nodes and use the vector basis functions to find the
664  // interpolated flux and its time derivative
665  for (unsigned i = 0; i < 2; i++)
666  {
667  // Loop over the edge basis vectors
668  for (unsigned l = 0; l < n_q_basis_edge; l++)
669  {
670  interpolated_q[i] += q_edge(l) * q_basis(l, i);
671  interpolated_dq_dt[i] += dq_edge_dt(l) * q_basis(l, i);
672  }
673  // Loop over the internal basis vectors
674  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
675  {
676  interpolated_q[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
677  interpolated_dq_dt[i] +=
678  dq_internal_dt(l - n_q_basis_edge) * q_basis(l, i);
679  }
680  }
681 
682  // loop over the pressure basis and find the interpolated pressure
683  for (unsigned l = 0; l < n_p_basis; l++)
684  {
685  interpolated_p += p_value(l) * p_basis(l);
686  }
687 
688  // loop over the q edge divergence basis and the q internal divergence
689  // basis to find interpolated div q
690  for (unsigned l = 0; l < n_q_basis_edge; l++)
691  {
692  interpolated_div_q_ds += q_edge(l) * div_q_basis_ds(l);
693  }
694  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
695  {
696  interpolated_div_q_ds +=
697  q_internal(l - n_q_basis_edge) * div_q_basis_ds(l);
698  }
699 
700  // Get the solid body force
701  this->solid_body_force(time, interpolated_x, f_solid);
702 
703  // Get the fluid nody force
704  this->fluid_body_force(time, interpolated_x, f_fluid);
705 
706  // Get the mass source function
707  this->mass_source(time, interpolated_x, mass_source_local);
708 
709  double r = interpolated_x[0];
710 
711  // Linear elasticity:
712  //-------------------
713 
714  double u_r = interpolated_u[0];
715  double du_r_dr = interpolated_du_dx(0, 0);
716  double du_r_dz = interpolated_du_dx(0, 1);
717  double du_z_dr = interpolated_du_dx(1, 0);
718  double du_z_dz = interpolated_du_dx(1, 1);
719 
720  // Storage for terms of Jacobian
721  double G_r = 0, G_z = 0;
722 
723  for (unsigned l = 0; l < n_node; l++)
724  {
725  for (unsigned a = 0; a < 2; a++)
726  {
727  local_eqn =
729 
730  if (local_eqn >= 0)
731  {
732  residuals[local_eqn] +=
733  (lambda_sq *
734  (interpolated_d2u_dt2[a] +
735  rho_f_over_rho * local_permeability * interpolated_dq_dt[a]) -
736  f_solid[a]) *
737  u_test(l) * r * w * J;
738 
739  // r-equation
740  if (a == 0)
741  {
742  residuals[local_eqn] +=
743  (mu * (2.0 * du_r_dr * du_test_dx(l, 0) +
744  du_test_dx(l, 1) * (du_r_dz + du_z_dr) +
745  2.0 * u_test(l) / pow(r, 2) * (u_r)) +
746  (lambda * (du_r_dr + u_r / r + du_z_dz) -
747  alpha * interpolated_p) *
748  (du_test_dx(l, 0) + u_test(l) / r)) *
749  r * w * J;
750  }
751  else if (a == 1)
752  {
753  residuals[local_eqn] +=
754  (mu * (du_test_dx(l, 0) * (du_r_dz + du_z_dr) +
755  2.0 * du_z_dz * du_test_dx(l, 1)) +
756  (lambda * (du_r_dr + u_r / r + du_z_dz) -
757  alpha * interpolated_p) *
758  du_test_dx(l, 1)) *
759  r * w * J;
760  }
761  // error: a should be 0 or 1
762  else
763  {
764  throw OomphLibError("a should equal 0 or 1",
765  OOMPH_CURRENT_FUNCTION,
766  OOMPH_EXCEPTION_LOCATION);
767  }
768 
769  // Jacobian entries
770  if (flag)
771  {
772  // d(u_eqn_l,a)/d(U_l2,c)
773  for (unsigned l2 = 0; l2 < n_node; l2++)
774  {
775  if (a == 0)
776  {
777  G_r =
778  (mu * (2.0 * du_basis_dx(l2, 0) * du_test_dx(l, 0) +
779  du_test_dx(l, 1) * du_basis_dx(l2, 1) +
780  2 * u_test(l) / pow(r, 2) * u_basis(l2)) +
781  (lambda * (du_basis_dx(l2, 0) + u_basis(l2) / r)) *
782  (du_test_dx(l, 0) + u_test(l) / r) +
783  lambda_sq * node_pt(l2)->time_stepper_pt()->weight(2, 0) *
784  u_basis(l2) * u_test(l)) *
785  r * w * J;
786 
787  G_z = (mu * du_test_dx(l, 1) * du_basis_dx(l2, 0) +
788  lambda * du_basis_dx(l2, 1) *
789  (du_test_dx(l, 0) + u_test(l) / r)) *
790  r * w * J;
791  }
792  else if (a == 1)
793  {
794  G_r = (mu * du_test_dx(l, 0) * du_basis_dx(l2, 1) +
795  lambda * (du_basis_dx(l2, 0) + u_basis(l2) / r) *
796  du_test_dx(l, 1)) *
797  r * w * J;
798  G_z =
799  (mu * (du_test_dx(l, 0) * du_basis_dx(l2, 0) +
800  2.0 * du_basis_dx(l2, 1) * du_test_dx(l, 1)) +
801  lambda * du_basis_dx(l2, 1) * du_test_dx(l, 1) +
802  lambda_sq * node_pt(l2)->time_stepper_pt()->weight(2, 0) *
803  u_basis(l2) * u_test(l)) *
804  r * w * J;
805  }
806 
807  for (unsigned c = 0; c < 2; c++)
808  {
809  local_unknown =
811  if (local_unknown >= 0)
812  {
813  if (c == 0)
814  {
815  jacobian(local_eqn, local_unknown) += G_r;
816  }
817  else if (c == 1)
818  {
819  jacobian(local_eqn, local_unknown) += G_z;
820  }
821  }
822  }
823  }
824 
825  // d(u_eqn_l,a)/d(Q_l2)
826  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
827  {
828  TimeStepper* timestepper_pt = 0;
829 
830  if (l2 < n_q_basis_edge)
831  {
832  local_unknown = q_edge_local_eqn(l2);
833  timestepper_pt =
835  }
836  else // n_q_basis_edge <= l < n_basis
837  {
838  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
839  timestepper_pt = this->internal_data_pt(q_internal_index())
840  ->time_stepper_pt();
841  }
842 
843  if (local_unknown >= 0)
844  {
845  jacobian(local_eqn, local_unknown) +=
846  lambda_sq * rho_f_over_rho * local_permeability *
847  timestepper_pt->weight(1, 0) * q_basis(l2, a) * u_test(l) *
848  r * w * J;
849  }
850  }
851 
852  // d(u_eqn_l,a)/d(P_l2)
853  for (unsigned l2 = 0; l2 < n_p_basis; l2++)
854  {
855  local_unknown = p_local_eqn(l2);
856  if (local_unknown >= 0)
857  {
858  if (a == 0)
859  {
860  jacobian(local_eqn, local_unknown) -=
861  alpha * p_basis(l2) * (du_test_dx(l, 0) + u_test(l) / r) *
862  r * w * J;
863  }
864  else if (a == 1)
865  {
866  jacobian(local_eqn, local_unknown) -=
867  alpha * p_basis(l2) * du_test_dx(l, 1) * r * w * J;
868  }
869  }
870  }
871  } // End of Jacobian entries
872  } // End of if not boundary condition
873  } // End of loop over dimensions
874  } // End of loop over u test functions
875 
876 
877  // Darcy:
878  //-------
879 
880  // Loop over the test functions
881  for (unsigned l = 0; l < n_q_basis; l++)
882  {
883  if (l < n_q_basis_edge)
884  {
885  local_eqn = q_edge_local_eqn(l);
886  }
887  else // n_q_basis_edge <= l < n_basis
888  {
889  local_eqn = q_internal_local_eqn(l - n_q_basis_edge);
890  }
891 
892  // If it's not a boundary condition
893  if (local_eqn >= 0)
894  {
895  for (unsigned i = 0; i < 2; i++)
896  {
897  residuals[local_eqn] +=
898  (rho_f_over_rho * lambda_sq *
899  (interpolated_d2u_dt2[i] +
900  (local_permeability / porosity) * interpolated_dq_dt[i]) +
901  interpolated_q[i] / local_permeability_ratio -
902  rho_f_over_rho * f_fluid[i]) *
903  q_test(l, i) * r * w * J;
904  }
905 
906  // deliberately no jacobian factor in this integral
907  residuals[local_eqn] -= interpolated_p * div_q_test_ds(l) * r * w;
908 
909  // deliberately no r factor in this integral
910  residuals[local_eqn] -= interpolated_p * q_test(l, 0) * w * J;
911 
912  // Jacobian entries
913  if (flag)
914  {
915  // d(q_eqn_l)/d(U_l2,c)
916  for (unsigned l2 = 0; l2 < n_node; l2++)
917  {
918  for (unsigned c = 0; c < 2; c++)
919  {
920  local_unknown =
922  if (local_unknown >= 0)
923  {
924  jacobian(local_eqn, local_unknown) +=
925  rho_f_over_rho * lambda_sq *
926  this->node_pt(l2)->time_stepper_pt()->weight(2, 0) *
927  u_basis(l2) * q_test(l, c) * r * w * J;
928  }
929  }
930  }
931 
932  // d(q_eqn_l)/d(Q_l2)
933  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
934  {
935  TimeStepper* timestepper_pt = 0;
936 
937  if (l2 < n_q_basis_edge)
938  {
939  local_unknown = q_edge_local_eqn(l2);
940  timestepper_pt =
942  }
943  else // n_q_basis_edge <= l < n_basis
944  {
945  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
946  timestepper_pt =
948  }
949 
950  if (local_unknown >= 0)
951  {
952  for (unsigned c = 0; c < 2; c++)
953  {
954  jacobian(local_eqn, local_unknown) +=
955  q_basis(l2, c) * q_test(l, c) *
956  (1.0 / local_permeability_ratio +
957  rho_f_over_rho * lambda_sq * local_permeability *
958  timestepper_pt->weight(1, 0) / porosity) *
959  r * w * J;
960  }
961  }
962  }
963 
964  // d(q_eqn_l)/d(P_l2)
965  for (unsigned l2 = 0; l2 < n_p_basis; l2++)
966  {
967  local_unknown = p_local_eqn(l2);
968 
969  if (local_unknown >= 0)
970  {
971  jacobian(local_eqn, local_unknown) -=
972  p_basis(l2) * (div_q_test_ds(l) * r + q_test(l, 0) * J) * w;
973  }
974  }
975  } // End of Jacobian entries
976  } // End of if not boundary condition
977  } // End of loop over q test functions
978 
979  // loop over pressure test functions
980  for (unsigned l = 0; l < n_p_basis; l++)
981  {
982  // get the local equation number
983  local_eqn = p_local_eqn(l);
984 
985  // If it's not a boundary condition
986  if (local_eqn >= 0)
987  {
988  // solid divergence term
989  residuals[local_eqn] +=
990  alpha * interpolated_div_du_dt_dx * p_test(l) * r * w * J;
991  residuals[local_eqn] +=
992  alpha * interpolated_du_r_dt * p_test(l) * w * J;
993  // deliberately no jacobian factor in this integral
994  residuals[local_eqn] +=
995  local_permeability * interpolated_div_q_ds * p_test(l) * r * w;
996  // deliberately no r factor in this integral
997  residuals[local_eqn] +=
998  local_permeability * interpolated_q[0] * p_test(l) * w * J;
999  residuals[local_eqn] -= mass_source_local * p_test(l) * r * w * J;
1000 
1001  // Jacobian entries
1002  if (flag)
1003  {
1004  // d(p_eqn_l)/d(U_l2,c)
1005  for (unsigned l2 = 0; l2 < n_node; l2++)
1006  {
1007  for (unsigned c = 0; c < 2; c++)
1008  {
1009  local_unknown =
1011 
1012  if (local_unknown >= 0)
1013  {
1014  jacobian(local_eqn, local_unknown) +=
1015  alpha * this->node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1016  du_basis_dx(l2, c) * p_test(l) * r * w * J;
1017 
1018  // Extra term due to cylindrical coordinate system
1019  if (c == 0)
1020  {
1021  jacobian(local_eqn, local_unknown) +=
1022  alpha *
1023  this->node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1024  u_basis(l2) * p_test(l) * w * J;
1025  }
1026  }
1027  }
1028  }
1029 
1030  // d(p_eqn_l)/d(Q_l2)
1031  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
1032  {
1033  if (l2 < n_q_basis_edge)
1034  {
1035  local_unknown = q_edge_local_eqn(l2);
1036  }
1037  else // n_q_basis_edge <= l < n_basis
1038  {
1039  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
1040  }
1041 
1042  if (local_unknown >= 0)
1043  {
1044  jacobian(local_eqn, local_unknown) +=
1045  (div_q_basis_ds(l2) * r + q_basis(l2, 0) * J) *
1046  local_permeability * p_test(l) * w;
1047  }
1048  }
1049  } // End of Jacobian entries
1050  } // End of if not boundary condition
1051  } // End of loop over p test functions
1052  } // End of loop over integration points
1053  }
1054 
1055 
1056 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
static double Default_lambda_sq_value
Static default value for timescale ratio.
void mass_source(const double &time, const Vector< double > &x, double &b) const
Indirect access to the mass source function - returns 0 if no mass source function has been set.
static double Default_youngs_modulus_value
Static default value for Young's modulus (1.0 – for natural scaling, i.e. all stresses have been non-...
virtual unsigned q_internal_index() const =0
Return the index of the internal data where the q internal degrees of freedom are stored.
double dq_internal_dt(const unsigned &n) const
dq_internal/dt for the n-th internal degree of freedom
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 ...
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.
virtual unsigned u_index_axisym_poroelasticity(const unsigned &j) const =0
Return the nodal index of the j-th solid displacement unknown.
const double & porosity() const
Access function for the porosity.
virtual int q_edge_local_eqn(const unsigned &j) const =0
Return the equation number of the j-th edge (flux) degree of freedom.
void interpolated_u(const Vector< double > &s, Vector< double > &disp) const
Calculate the FE representation of u.
const double & youngs_modulus() const
Access function to non-dim Young's modulus (ratio of actual Young's modulus to reference stress used ...
virtual unsigned q_edge_node_number(const unsigned &j) const =0
Return the number of the node where the jth edge unknown is stored.
virtual int q_internal_local_eqn(const unsigned &j) const =0
Return the equation number of the j-th internal degree of freedom.
virtual double shape_basis_test_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsi, Shape &u_basis, Shape &u_test, DShape &du_basis_dx, DShape &du_test_dx, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
Compute the geometric basis, and the q, p and divergence basis functions and test functions at integr...
void interpolated_q(const Vector< double > &s, Vector< double > &q) const
Calculate the FE representation of q.
const double & permeability_ratio() const
Access function for the ratio of the material's actual permeability to the permeability used in the n...
static double Default_porosity_value
Static default value for the porosity.
double du_dt(const unsigned &n, const unsigned &i) const
du/dt at local node n
static double Default_permeability_value
Static default value for the permeability (1.0 for natural scaling; i.e. timescale is given by L^2/(k...
static double Default_permeability_ratio_value
Static default value for the ratio of the material's actual permeability to the permeability used to ...
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div u.
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, DShape &dpsi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
virtual unsigned nq_basis() const
Return the total number of computational basis functions for q.
const double & density_ratio() const
Access function for the density ratio (fluid to solid)
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
virtual double q_internal(const unsigned &j) const =0
Return the values of the j-th internal degree of freedom.
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
void interpolated_du_dt(const Vector< double > &s, Vector< double > &du_dt) const
Calculate the FE representation of du_dt.
virtual double p_value(const unsigned &j) const =0
Return the jth pressure value.
const double & permeability() const
Access function for the nondim permeability.
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
double d2u_dt2(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
static double Default_density_ratio_value
Static default value for the density ratio.
const double & alpha() const
Access function for alpha, the Biot parameter.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output FE representation of exact soln: x,y,u1,u2,div_q,p at Nplot^2 plot points.
void fluid_body_force(const double &time, const Vector< double > &x, Vector< double > &b) const
Indirect access to the fluid body force function - returns 0 if no forcing function has been set.
static double Default_alpha_value
Static default value for alpha, the biot parameter.
void solid_body_force(const double &time, const Vector< double > &x, Vector< double > &b) const
Indirect access to the solid body force function - returns 0 if no forcing function has been set.
virtual double q_edge(const unsigned &j) const =0
Return the values of the j-th edge (flux) degree of freedom.
const double & nu() const
Access function for Poisson's ratio.
double dq_edge_dt(const unsigned &n) const
dq_edge/dt for the n-th edge degree of freedom
void output(std::ostream &outfile)
Output with default number of plot points.
virtual int p_local_eqn(const unsigned &j) const =0
Return the equation number of the j-th pressure degree of freedom.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4103
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t.local coordinates to derivatives w.r.t the coordinates used to assemble the ...
Definition: elements.cc:2833
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1508
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2576
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
//////////////////////////////////////////////////////////////////// ////////////////////////////////...