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