axisym_poroelasticity_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_AXISYM_POROELASTICITY_ELEMENTS_HEADER
31 #define OOMPH_AXISYM_POROELASTICITY_ELEMENTS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 #include "../generic/elements.h"
39 #include "../generic/shape.h"
40 #include "../generic/error_estimator.h"
41 #include "../generic/projection.h"
42 
43 
44 namespace oomph
45 {
46 
47 //=========================================================================
48  /// \short Class implementing the generic maths of the axisym poroelasticity
49  /// equations: axisym linear elasticity coupled with axisym Darcy
50  /// equations (using Raviart-Thomas elements with both edge and internal
51  /// degrees of freedom) including inertia in both.
52 //=========================================================================
54  public virtual ElementWithZ2ErrorEstimator
55  {
56  public:
57 
58  /// Source function pointer typedef
59  typedef void (*SourceFctPt)(const double &time,
60  const Vector<double>& x,
61  Vector<double>& f);
62 
63  /// Mass source function pointer typedef
64  typedef void (*MassSourceFctPt)(const double &time,
65  const Vector<double>& x,
66  double& f);
67 
68  /// Constructor
74  Nu_pt(0),
82  {
83  }
84 
85  /// \short Access function to non-dim Young's modulus (ratio of actual
86  /// Young's modulus to reference stress used to non-dim the equations.
87  /// (const version)
88  const double& youngs_modulus() const {return *Youngs_modulus_pt;}
89 
90  /// \short Pointer to non-dim Young's modulus (ratio of actual
91  /// Young's modulus to reference stress used to non-dim the equations.
92  double* &youngs_modulus_pt() {return Youngs_modulus_pt;}
93 
94  ///Access function for Poisson's ratio
95  const double& nu() const
96  {
97 #ifdef PARANOID
98  if (Nu_pt==0)
99  {
100  std::ostringstream error_message;
101  error_message << "No pointer to Poisson's ratio set. Please set one!\n";
102  throw OomphLibError(
103  error_message.str(),
104  OOMPH_CURRENT_FUNCTION,
105  OOMPH_EXCEPTION_LOCATION);
106  }
107 #endif
108  return *Nu_pt;
109  }
110 
111  /// Access function for pointer to Poisson's ratio
112  double* &nu_pt() {return Nu_pt;}
113 
114  /// Access function for timescale ratio (nondim density)
115  const double& lambda_sq() const {return *Lambda_sq_pt;}
116 
117  /// Access function for pointer to timescale ratio (nondim density)
118  double* &lambda_sq_pt() {return Lambda_sq_pt;}
119 
120  /// Access function for the density ratio (fluid to solid)
121  const double& density_ratio() const {return *Density_ratio_pt;}
122 
123  /// Access function for pointer to the density ratio (fluid to solid)
124  double* &density_ratio_pt() {return Density_ratio_pt;}
125 
126  /// Access function for the nondim permeability
127  const double& permeability() const {return *Permeability_pt;}
128 
129  /// Access function for pointer to the nondim permeability
130  double* &permeability_pt() {return Permeability_pt;}
131 
132 
133  /// \short Access function for the ratio of the material's actual permeability
134  /// to the permeability used in the non-dimensionalisation of the equations
135  const double& permeability_ratio() const {return *Permeability_ratio_pt;}
136 
137  /// Access function for pointer to ratio of the material's actual permeability
138  /// to the permeability used in the non-dimensionalisation of the equations
140 
141  /// Access function for alpha, the Biot parameter
142  const double& alpha() const {return *Alpha_pt;}
143 
144  /// Access function for pointer to alpha, the Biot parameter
145  double* &alpha_pt() {return Alpha_pt;}
146 
147  /// Access function for the porosity
148  const double& porosity() const {return *Porosity_pt;}
149 
150  /// Access function for pointer to the porosity
151  double* &porosity_pt() {return Porosity_pt;}
152 
153  /// Access function: Pointer to solid body force function
155 
156  /// Access function: Pointer to solid body force function (const version)
158 
159  /// Access function: Pointer to fluid force function
161 
162  /// Access function: Pointer to fluid force function (const version)
164 
165  /// Access function: Pointer to mass source function
167 
168  /// Access function: Pointer to mass source function (const version)
170 
171  /// \short Indirect access to the solid body force function - returns 0 if no
172  /// forcing function has been set
173  void solid_body_force(const double& time,
174  const Vector<double>& x,
175  Vector<double>& b) const
176  {
177  // If no function has been set, return zero vector
179  {
180  // Get spatial dimension of element
181  unsigned n=dim();
182  for(unsigned i=0;i<n;i++)
183  {
184  b[i] = 0.0;
185  }
186  }
187  else
188  {
189  (*Solid_body_force_fct_pt)(time,x,b);
190  }
191  }
192 
193  /// \short Indirect access to the fluid body force function - returns 0 if no
194  /// forcing function has been set
195  void fluid_body_force(const double& time,
196  const Vector<double>& x,
197  Vector<double>& b) const
198  {
199  // If no function has been set, return zero vector
201  {
202  // Get spatial dimension of element
203  unsigned n=dim();
204  for(unsigned i=0;i<n;i++)
205  {
206  b[i] = 0.0;
207  }
208  }
209  else
210  {
211  (*Fluid_body_force_fct_pt)(time,x,b);
212  }
213  }
214 
215  /// \short Indirect access to the mass source function - returns 0 if no
216  /// mass source function has been set
217  void mass_source(const double& time,
218  const Vector<double>& x,
219  double& b) const
220  {
221  // If no function has been set, return zero vector
222  if(Mass_source_fct_pt==0)
223  {
224  b = 0.0;
225  }
226  else
227  {
228  (*Mass_source_fct_pt)(time,x,b);
229  }
230  }
231 
232  /// Number of values required at node n
233  virtual unsigned required_nvalue(const unsigned &n) const = 0;
234 
235  /// Return the nodal index of the j-th solid displacement unknown
236  virtual unsigned u_index_axisym_poroelasticity(const unsigned &j) const = 0;
237 
238  /// Return the equation number of the j-th edge (flux) degree of freedom
239  virtual int q_edge_local_eqn(const unsigned &j) const = 0;
240 
241  /// Return the equation number of the j-th internal degree of freedom
242  virtual int q_internal_local_eqn(const unsigned &j) const = 0;
243 
244  /// \short Return vector of pointers to the Data objects that store the
245  /// edge flux values
246  virtual Vector<Data*> q_edge_data_pt() const=0;
247 
248  /// Return pointer to the Data object that stores the internal flux values
249  virtual Data* q_internal_data_pt() const=0;
250 
251  /// Return the nodal index at which the jth edge unknown is stored
252  virtual unsigned q_edge_index(const unsigned &j) const = 0;
253 
254  /// \short Return the index of the internal data where the q internal degrees
255  /// of freedom are stored
256  virtual unsigned q_internal_index() const = 0;
257 
258  /// Return the number of the node where the jth edge unknown is stored
259  virtual unsigned q_edge_node_number(const unsigned &j) const = 0;
260 
261  /// Return the values of the j-th edge (flux) degree of freedom
262  virtual double q_edge(const unsigned &j) const = 0;
263 
264  /// \short Return the values of the j-th edge (flux) degree of freedom at time
265  /// history level t
266  virtual double q_edge(const unsigned &t,const unsigned &j) const = 0;
267 
268  /// Return the face index associated with j-th edge flux degree of freedom
269  virtual unsigned face_index_of_q_edge_basis_fct(const unsigned& j) const=0;
270 
271  /// Return the face index associated with specified edge
272  virtual unsigned face_index_of_edge(const unsigned& j) const=0;
273 
274  /// \short Compute the face element coordinates of the nth flux interpolation
275  /// point along an edge
277  const unsigned &edge,
278  const unsigned &n,
279  Vector<double> &s) const=0;
280 
281  /// Return the values of the j-th internal degree of freedom
282  virtual double q_internal(const unsigned &j) const = 0;
283 
284  /// \short Return the values of the j-th internal degree of freedom at
285  /// time history level t
286  virtual double q_internal(const unsigned &t,const unsigned &j) const = 0;
287 
288  /// Set the values of the j-th edge (flux) degree of freedom
289  virtual void set_q_edge(const unsigned &j, const double& value)=0;
290 
291  /// Set the values of the j-th internal degree of freedom
292  virtual void set_q_internal(const unsigned &j, const double& value)=0;
293 
294  /// \short Set the values of the j-th edge (flux) degree of freedom at
295  /// time history level t
296  virtual void set_q_edge(const unsigned &t, const unsigned &j,
297  const double& value)=0;
298 
299  /// \short Set the values of the j-th internal degree of freedom at
300  /// time history level t
301  virtual void set_q_internal(const unsigned &t, const unsigned &j,
302  const double& value)=0;
303 
304  /// Return the total number of computational basis functions for q
305  virtual unsigned nq_basis() const
306  {
307  return nq_basis_edge()+nq_basis_internal();
308  }
309 
310  /// Return the number of edge basis functions for q
311  virtual unsigned nq_basis_edge() const = 0;
312 
313  /// Return the number of internal basis functions for q
314  virtual unsigned nq_basis_internal() const = 0;
315 
316  /// Comute the local form of the q basis at local coordinate s
317  virtual void get_q_basis_local(const Vector<double> &s,
318  Shape &q_basis) const = 0;
319 
320  /// Compute the local form of the q basis and dbasis/ds at local coordinate s
321  virtual void get_div_q_basis_local(const Vector<double> &s,
322  Shape &div_q_basis_ds) const = 0;
323 
324  /// Compute the transformed basis at local coordinate s
326  Shape &q_basis) const
327  {
328  const unsigned n_node = this->nnode();
329  Shape psi(n_node,2);
330  const unsigned n_q_basis = this->nq_basis();
331  Shape q_basis_local(n_q_basis,2);
332  this->get_q_basis_local(s,q_basis_local);
333  (void)this->transform_basis(s,q_basis_local,psi,q_basis);
334  }
335 
336  /// \short Returns the number of flux_interpolation points along each edge of
337  /// the element
338  virtual unsigned nedge_flux_interpolation_point() const = 0;
339 
340  /// \short Returns the local coordinate of the jth flux_interpolation point
341  /// along the specified edge
342  virtual Vector<double> edge_flux_interpolation_point(const unsigned &edge,
343  const unsigned &j)
344  const = 0;
345 
346  /// \short Compute the global coordinates of the jth flux_interpolation
347  /// point along an edge
348  virtual void edge_flux_interpolation_point_global(const unsigned &edge,
349  const unsigned &j,
350  Vector<double> &x)
351  const = 0;
352 
353  /// Pin the jth internal q value and set it to specified value
354  virtual void pin_q_internal_value(const unsigned &j,
355  const double& value) = 0;
356 
357  /// Pin the j-th edge (flux) degree of freedom and set it to specified value
358  virtual void pin_q_edge_value(const unsigned &j, const double& value) = 0;
359 
360  /// Return the equation number of the j-th pressure degree of freedom
361  virtual int p_local_eqn(const unsigned &j) const = 0;
362 
363  /// Return the jth pressure value
364  virtual double p_value(const unsigned &j) const = 0;
365 
366  /// Return the total number of pressure basis functions
367  virtual unsigned np_basis() const = 0;
368 
369  /// Compute the pressure basis
370  virtual void get_p_basis(const Vector<double> &s,
371  Shape &p_basis) const = 0;
372 
373  /// Pin the jth pressure value and set it to p
374  virtual void pin_p_value(const unsigned &j, const double &p) = 0;
375 
376  /// Set the jth pressure value
377  virtual void set_p_value(const unsigned &j, const double& value)=0;
378 
379  /// Return pointer to the Data object that stores the pressure values
380  virtual Data* p_data_pt() const=0;
381 
382  /// Scale the edge basis to allow arbitrary edge mappings
383  virtual void scale_basis(Shape& basis) const = 0;
384 
385  /// \short Performs a div-conserving transformation of the vector basis
386  /// functions from the reference element to the actual element
387  double transform_basis(const Vector<double> &s,
388  const Shape &q_basis_local,
389  Shape &psi,
390  DShape &dpsi,
391  Shape &q_basis) const;
392 
393  /// \short Performs a div-conserving transformation of the vector basis
394  /// functions from the reference element to the actual element
396  const Shape &q_basis_local,
397  Shape &psi,
398  Shape &q_basis) const
399  {
400  const unsigned n_node = this->nnode();
401  DShape dpsi(n_node,2);
402  return transform_basis(s,q_basis_local,psi,dpsi,q_basis);
403  }
404 
405  /// Fill in contribution to residuals for the Darcy equations
407  {
410  }
411 
412  /// Fill in the Jacobian matrix for the Newton method
414  DenseMatrix<double> &jacobian)
415  {
416  this->fill_in_generic_residual_contribution(residuals,jacobian,1);
417  }
418 
419 
420  /// Calculate the FE representation of the divergence of the
421  /// skeleton velocity, div(du/dt), and its
422  /// components: 1/r diff(r*du_r/dt,r) and diff(du_z/dt,z).
424  Vector<double>& div_dudt_components) const
425  {
426  //Find number of nodes
427  unsigned n_node = nnode();
428 
429  //Local shape function
430  Shape psi(n_node);
431  DShape dpsidx(n_node,2);
432 
433  //Find values of shape function
434  dshape_eulerian(s,psi,dpsidx);
435 
436  // Local coordinates
437  double r=interpolated_x(s,0);
438 
439  // Assemble the "cartesian-like" contributions
440  for(unsigned i=0;i<2;i++)
441  {
442  // Initialise
443  div_dudt_components[i] = 0.0;
444 
445  // Loop over the local nodes and sum the "cartesian-like"
446  // contributions
447  for(unsigned l=0;l<n_node;l++)
448  {
449  div_dudt_components[i] += du_dt(l,i)*dpsidx(l,i);
450  }
451  }
452 
453  // Radial skeleton veloc
454  double dur_dt=0.0;
455  for(unsigned l=0;l<n_node;l++)
456  {
457  dur_dt+=du_dt(l,0)*psi(l);
458  }
459 
460  // Add geometric component to radial contribution
461  div_dudt_components[0]+=dur_dt/r;
462 
463  // Return sum
464  return div_dudt_components[0]+div_dudt_components[1];
465  }
466 
467 
468  /// Calculate the FE representation of the divergence of the
469  /// skeleton displ, div(u), and its
470  /// components: 1/r diff(r*u_r,r) and diff(u_z,z).
472  Vector<double>& div_u_components) const
473  {
474  //Find number of nodes
475  unsigned n_node = nnode();
476 
477  //Local shape function
478  Shape psi(n_node);
479  DShape dpsidx(n_node,2);
480 
481  //Find values of shape function
482  dshape_eulerian(s,psi,dpsidx);
483 
484  // Local coordinates
485  double r=interpolated_x(s,0);
486 
487  // Assemble the "cartesian-like" contributions
488  for(unsigned i=0;i<2;i++)
489  {
490  //Get nodal index at which i-th velocity is stored
491  unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
492 
493  // Initialise
494  div_u_components[i] = 0.0;
495 
496  // Loop over the local nodes and sum the "cartesian-like"
497  // contributions
498  for(unsigned l=0;l<n_node;l++)
499  {
500  div_u_components[i] += nodal_value(l,u_nodal_index)*dpsidx(l,i);
501  }
502  }
503 
504  // Radial skeleton displ
505  double ur=0.0;
506  for(unsigned l=0;l<n_node;l++)
507  {
509  }
510 
511  // Add geometric component to radial contribution
512  div_u_components[0]+=ur/r;
513 
514  // Return sum
515  return div_u_components[0]+div_u_components[1];
516  }
517 
518 
519 
520  /// Calculate the FE representation of u
522  Vector<double>& disp) const
523  {
524  //Find number of nodes
525  unsigned n_node = nnode();
526 
527  //Local shape function
528  Shape psi(n_node);
529 
530  //Find values of shape function
531  shape(s,psi);
532 
533  for(unsigned i=0;i<2;i++)
534  {
535  //Index at which the nodal value is stored
536  unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
537 
538  //Initialise value of u
539  disp[i] = 0.0;
540 
541  //Loop over the local nodes and sum
542  for(unsigned l=0;l<n_node;l++)
543  {
544  disp[i] += nodal_value(l,u_nodal_index)*psi[l];
545  }
546  }
547  }
548 
549  /// Calculate the FE representation of the i-th component of u
551  const unsigned &i) const
552  {
553  //Find number of nodes
554  unsigned n_node = nnode();
555 
556  //Local shape function
557  Shape psi(n_node);
558 
559  //Find values of shape function
560  shape(s,psi);
561 
562  //Get nodal index at which i-th velocity is stored
563  unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
564 
565  //Initialise value of u
566  double interpolated_u = 0.0;
567 
568  //Loop over the local nodes and sum
569  for(unsigned l=0;l<n_node;l++)
570  {
571  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
572  }
573 
574  return(interpolated_u);
575  }
576 
577 
578  /// \short Calculate the FE representation of the i-th component of u
579  /// at time level t (t=0: current)
580  double interpolated_u(const unsigned& t,
581  const Vector<double> &s,
582  const unsigned &i) const
583  {
584  //Find number of nodes
585  unsigned n_node = nnode();
586 
587  //Local shape function
588  Shape psi(n_node);
589 
590  //Find values of shape function
591  shape(s,psi);
592 
593  //Get nodal index at which i-th velocity is stored
594  unsigned u_nodal_index = u_index_axisym_poroelasticity(i);
595 
596  //Initialise value of u
597  double interpolated_u = 0.0;
598 
599  //Loop over the local nodes and sum
600  for(unsigned l=0;l<n_node;l++)
601  {
602  interpolated_u += nodal_value(t,l,u_nodal_index)*psi[l];
603  }
604 
605  return(interpolated_u);
606  }
607 
608 
609  /// Calculate the FE representation of du_dt
611  Vector<double>& du_dt) const
612  {
613  //Find number of nodes
614  unsigned n_node = nnode();
615 
616  //Local shape function
617  Shape psi(n_node);
618 
619  //Find values of shape function
620  shape(s,psi);
621 
622  for(unsigned i=0;i<2;i++)
623  {
624  //Initialise value of u
625  du_dt[i] = 0.0;
626 
627  //Loop over the local nodes and sum
628  for(unsigned l=0;l<n_node;l++)
629  {
630  du_dt[i] += this->du_dt(l,i)*psi[l];
631  }
632  }
633  }
634 
635  /// Calculate the FE representation of q
637  Vector<double> &q) const
638  {
639  unsigned n_q_basis = nq_basis();
640  unsigned n_q_basis_edge = nq_basis_edge();
641  Shape q_basis(n_q_basis,2);
642  get_q_basis(s,q_basis);
643  for(unsigned i=0;i<2;i++)
644  {
645  q[i]=0.0;
646  for(unsigned l=0;l<n_q_basis_edge;l++)
647  {
648  q[i] += q_edge(l)*q_basis(l,i);
649  }
650  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
651  {
652  q[i] += q_internal(l-n_q_basis_edge)*q_basis(l,i);
653  }
654  }
655  }
656 
657  /// \short Calculate the FE representation of q
658  /// at time level t (t=0: current)
659  void interpolated_q(const unsigned& t,
660  const Vector<double> &s,
661  Vector<double> &q) const
662  {
663  unsigned n_q_basis = nq_basis();
664  unsigned n_q_basis_edge = nq_basis_edge();
665  Shape q_basis(n_q_basis,2);
666  get_q_basis(s,q_basis);
667  for(unsigned i=0;i<2;i++)
668  {
669  q[i]=0.0;
670  for(unsigned l=0;l<n_q_basis_edge;l++)
671  {
672  q[i] += q_edge(t,l)*q_basis(l,i);
673  }
674  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
675  {
676  q[i] += q_internal(t,l-n_q_basis_edge)*q_basis(l,i);
677  }
678  }
679  }
680 
681  /// Calculate the FE representation of the i-th component of q
683  const unsigned i) const
684  {
685  unsigned n_q_basis = nq_basis();
686  unsigned n_q_basis_edge = nq_basis_edge();
687 
688  Shape q_basis(n_q_basis,2);
689 
690  get_q_basis(s,q_basis);
691  double q_i=0.0;
692  for(unsigned l=0;l<n_q_basis_edge;l++)
693  {
694  q_i += q_edge(l)*q_basis(l,i);
695  }
696  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
697  {
698  q_i += q_internal(l-n_q_basis_edge)*q_basis(l,i);
699  }
700 
701  return q_i;
702  }
703 
704  /// Calculate the FE representation of the i-th component of q
705  /// at time level t (t=0: current)
706  double interpolated_q(const unsigned& t,
707  const Vector<double> &s,
708  const unsigned i) const
709  {
710  unsigned n_q_basis = nq_basis();
711  unsigned n_q_basis_edge = nq_basis_edge();
712 
713  Shape q_basis(n_q_basis,2);
714 
715  get_q_basis(s,q_basis);
716  double q_i=0.0;
717  for(unsigned l=0;l<n_q_basis_edge;l++)
718  {
719  q_i += q_edge(t,l)*q_basis(l,i);
720  }
721  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
722  {
723  q_i += q_internal(t,l-n_q_basis_edge)*q_basis(l,i);
724  }
725 
726  return q_i;
727  }
728 
729 
730  /// Calculate the FE representation of div u
731  void interpolated_div_q(const Vector<double> &s, double &div_q) const
732  {
733  // Zero the divergence
734  div_q=0;
735 
736  // Get the number of nodes, q basis function, and q edge basis functions
737  unsigned n_node=nnode();
738  const unsigned n_q_basis = nq_basis();
739  const unsigned n_q_basis_edge = nq_basis_edge();
740 
741  // Storage for the divergence basis
742  Shape div_q_basis_ds(n_q_basis);
743 
744  // Storage for the geometric basis and it's derivatives
745  Shape psi(n_node);
746  DShape dpsi(n_node,2);
747 
748  // Call the geometric shape functions and their derivatives
749  this->dshape_local(s,psi,dpsi);
750 
751  // Storage for the inverse of the geometric jacobian (just so we can call
752  // the local to eulerian mapping)
753  DenseMatrix<double> inverse_jacobian(2);
754 
755  // Get the determinant of the geometric mapping
756  double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
757 
758  // Get the divergence basis (wrt local coords) at local coords s
759  get_div_q_basis_local(s,div_q_basis_ds);
760 
761  // Add the contribution to the divergence from the edge basis functions
762  for(unsigned l=0;l<n_q_basis_edge;l++)
763  {
764  div_q+=1.0/det*div_q_basis_ds(l)*q_edge(l);
765  }
766 
767  // Add the contribution to the divergence from the internal basis functions
768  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
769  {
770  div_q+=1.0/det*div_q_basis_ds(l)*q_internal(l-n_q_basis_edge);
771  }
772 
773  // Extra term due to cylindrical coords
774  if(std::abs(interpolated_x(s,0)) > 1e-10)
775  {
776  div_q+=interpolated_q(s,0)/interpolated_x(s,0);
777  }
778  }
779 
780 
781 
782 
783  /// Calculate the FE representation of div q and return it
784  double interpolated_div_q(const Vector<double> &s) const
785  {
786  // Temporary storage for div q
787  double div_q=0;
788 
789  // Get the intepolated divergence
790  interpolated_div_q(s,div_q);
791 
792  // Return it
793  return div_q;
794  }
795 
796  /// Calculate the FE representation of p
798  double &p) const
799  {
800  // Get the number of p basis functions
801  unsigned n_p_basis = np_basis();
802 
803  // Storage for the p basis
804  Shape p_basis(n_p_basis);
805 
806  // Call the p basis
807  get_p_basis(s,p_basis);
808 
809  // Zero the pressure
810  p=0;
811 
812  // Add the contribution to the pressure from each basis function
813  for(unsigned l=0;l<n_p_basis;l++)
814  {
815  p+=p_value(l)*p_basis(l);
816  }
817  }
818 
819  /// Calculate the FE representation of p and return it
820  double interpolated_p(const Vector<double> &s) const
821  {
822  // Temporary storage for p
823  double p=0;
824 
825  // Get the interpolated pressure
826  interpolated_p(s,p);
827 
828  // Return it
829  return p;
830  }
831 
832  /// du/dt at local node n
833  double du_dt(const unsigned &n,
834  const unsigned &i) const
835  {
836  // Get the timestepper
838 
839  // Storage for the derivative - initialise to 0
840  double du_dt=0.0;
841 
842  // If we are doing an unsteady solve then calculate the derivative
843  if(!time_stepper_pt->is_steady())
844  {
845  // Get the nodal index
846  const unsigned u_nodal_index=u_index_axisym_poroelasticity(i);
847 
848  // Get the number of values required to represent history
849  const unsigned n_time=time_stepper_pt->ntstorage();
850 
851  // Loop over history values
852  for(unsigned t=0;t<n_time;t++)
853  {
854  //Add the contribution to the derivative
855  du_dt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
856  }
857  }
858 
859  return du_dt;
860  }
861 
862  /// d^2u/dt^2 at local node n
863  double d2u_dt2(const unsigned &n,
864  const unsigned &i) const
865  {
866  // Get the timestepper
868 
869  // Storage for the derivative - initialise to 0
870  double d2u_dt2=0.0;
871 
872  // If we are doing an unsteady solve then calculate the derivative
873  if(!time_stepper_pt->is_steady())
874  {
875  // Get the nodal index
876  const unsigned u_nodal_index=u_index_axisym_poroelasticity(i);
877 
878  // Get the number of values required to represent history
879  const unsigned n_time=time_stepper_pt->ntstorage();
880 
881  // Loop over history values
882  for(unsigned t=0;t<n_time;t++)
883  {
884  //Add the contribution to the derivative
885  d2u_dt2+=time_stepper_pt->weight(2,t)*nodal_value(t,n,u_nodal_index);
886  }
887  }
888 
889  return d2u_dt2;
890  }
891 
892  /// dq_edge/dt for the n-th edge degree of freedom
893  double dq_edge_dt(const unsigned &n) const
894  {
895  unsigned node_num=q_edge_node_number(n);
896 
897  // get the timestepper
899 
900  // storage for the derivative - initialise to 0
901  double dq_dt=0.0;
902 
903  // if we are doing an unsteady solve then calculate the derivative
904  if(!time_stepper_pt->is_steady())
905  {
906  // get the number of values required to represent history
907  const unsigned n_time=time_stepper_pt->ntstorage();
908 
909  // loop over history values
910  for(unsigned t=0;t<n_time;t++)
911  {
912  // add the contribution to the derivative
913  dq_dt+=
914  time_stepper_pt->weight(1,t)*q_edge(t,n);
915  }
916  }
917 
918  return dq_dt;
919  }
920 
921  /// dq_internal/dt for the n-th internal degree of freedom
922  double dq_internal_dt(const unsigned &n) const
923  {
924  // get the internal data index for q
925  unsigned internal_index=q_internal_index();
926 
927  // get the timestepper
929  internal_data_pt(internal_index)->time_stepper_pt();
930 
931  // storage for the derivative - initialise to 0
932  double dq_dt=0.0;
933 
934  // if we are doing an unsteady solve then calculate the derivative
935  if(!time_stepper_pt->is_steady())
936  {
937  // get the number of values required to represent history
938  const unsigned n_time=time_stepper_pt->ntstorage();
939 
940  // loop over history values
941  for(unsigned t=0;t<n_time;t++)
942  {
943  // add the contribution to the derivative
944  dq_dt+=time_stepper_pt->weight(1,t)*q_internal(t,n);
945  }
946  }
947 
948  return dq_dt;
949  }
950 
951  /// Set the timestepper of the q internal data object
953  {
954  unsigned q_index=q_internal_index();
955  this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt,false);
956  }
957 
958 
959  /// Is Darcy flow switched off?
961  {
962  return Darcy_is_switched_off;
963  }
964 
965 
966  /// Switch off Darcy flow
968  {
970 
971  // Pin pressures and set them to zero
972  double p=0.0;
973  unsigned np=np_basis();
974  for (unsigned j=0;j<np;j++)
975  {
976  pin_p_value(j,p);
977  }
978 
979  // Pin internal flux data and set it to zero
980  double q=0.0;
981  unsigned nq=nq_basis_internal();
982  for (unsigned j=0;j<nq;j++)
983  {
985  }
986 
987  // Pin edge flux data and set it to zero
988  nq=nq_basis_edge();
989  for (unsigned j=0;j<nq;j++)
990  {
991  pin_q_edge_value(j,q);
992  }
993  }
994 
995  /// Self test
996  unsigned self_test()
997  {
998  return 0;
999  }
1000 
1001 
1002  /// \short Number of scalars/fields output by this element. Reimplements
1003  /// broken virtual function in base class.
1004  unsigned nscalar_paraview() const
1005  {
1006  return 8;
1007  }
1008 
1009  /// \short Write values of the i-th scalar field at the plot points. Needs
1010  /// to be implemented for each new specific element type.
1011  void scalar_value_paraview(std::ofstream& file_out,
1012  const unsigned& i,
1013  const unsigned& nplot) const
1014  {
1015  // Vector of local coordinates
1016  Vector<double> s(2);
1017 
1018  // Loop over plot points
1019  unsigned num_plot_points=nplot_points_paraview(nplot);
1020  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1021  {
1022 
1023  // Get local coordinates of plot point
1024  get_s_plot(iplot,nplot,s);
1025 
1026  // Skeleton velocity
1027  Vector<double> du_dt(2);
1028  interpolated_du_dt(s,du_dt);
1029 
1030  // Displacements
1031  if(i<2)
1032  {
1033  file_out << interpolated_u(s,i) << std::endl;
1034  }
1035  // Flux
1036  else if(i<4)
1037  {
1038  file_out << interpolated_q(s,i-2) << std::endl;
1039  }
1040  // Divergence of flux
1041  else if (i==4)
1042  {
1043  file_out << interpolated_div_q(s) << std::endl;
1044  }
1045  else if (i==5)
1046  {
1047  file_out << interpolated_p(s) << std::endl;
1048  }
1049  else if (i==6)
1050  {
1051  file_out << du_dt[0] << std::endl;
1052  }
1053  else if (i==7)
1054  {
1055  file_out << du_dt[1] << std::endl;
1056  }
1057  // Never get here
1058  else
1059  {
1060  std::stringstream error_stream;
1061  error_stream
1062  << "Axisymmetric poroelasticity elements only store 6 fields "
1063  << std::endl;
1064  throw OomphLibError(
1065  error_stream.str(),
1066  OOMPH_CURRENT_FUNCTION,
1067  OOMPH_EXCEPTION_LOCATION);
1068  }
1069  }
1070  }
1071 
1072  /// \short Name of the i-th scalar field. Default implementation
1073  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
1074  /// overloaded with more meaningful names in specific elements.
1075  std::string scalar_name_paraview(const unsigned& i) const
1076  {
1077  switch (i)
1078  {
1079  case 0:
1080  return "radial displacement";
1081  break;
1082 
1083  case 1:
1084  return "axial displacement";
1085  break;
1086 
1087  case 2:
1088  return "radial flux";
1089  break;
1090 
1091  case 3:
1092  return "axial flux";
1093  break;
1094 
1095  case 4:
1096  return "divergence of flux";
1097  break;
1098 
1099  case 5:
1100  return "pore pressure";
1101  break;
1102 
1103  case 6:
1104  return "radial skeleton velocity";
1105  break;
1106 
1107  case 7:
1108  return "axial skeleton velocity";
1109  break;
1110 
1111  default:
1112 
1113  std::stringstream error_stream;
1114  error_stream
1115  << "AxisymmetricPoroelasticityEquations only store 8 fields "
1116  << std::endl;
1117  throw OomphLibError(
1118  error_stream.str(),
1119  OOMPH_CURRENT_FUNCTION,
1120  OOMPH_EXCEPTION_LOCATION);
1121 
1122  // Dummy return
1123  return " ";
1124  }
1125  }
1126 
1127  /// \short Output solution in data vector at local cordinates s:
1128  /// r,z,u_r,u_z,q_r,q_z,div_q,p,durdt,duzdt
1130  {
1131  // Output the components of the position
1132  for(unsigned i=0;i<2;i++)
1133  {
1134  data.push_back(interpolated_x(s,i));
1135  }
1136 
1137  // Output the components of the FE representation of u at s
1138  for(unsigned i=0;i<2;i++)
1139  {
1140  data.push_back(interpolated_u(s,i));
1141  }
1142 
1143  // Output the components of the FE representation of q at s
1144  for(unsigned i=0;i<2;i++)
1145  {
1146  data.push_back(interpolated_q(s,i));
1147  }
1148 
1149  // Output FE representation of div u at s
1150  data.push_back(interpolated_div_q(s));
1151 
1152  // Output FE representation of p at s
1153  data.push_back(interpolated_p(s));
1154 
1155  // Skeleton velocity
1156  Vector<double> du_dt(2);
1157  interpolated_du_dt(s,du_dt);
1158  data.push_back(du_dt[0]);
1159  data.push_back(du_dt[1]);
1160 
1161  // Get divergence of skeleton velocity and its components
1162  Vector<double> div_dudt_components(2);
1163  data.push_back(interpolated_div_du_dt(s,div_dudt_components));
1164  data.push_back(div_dudt_components[0]);
1165  data.push_back(div_dudt_components[1]);
1166 
1167  // Get divergence of skeleton displacement and its components
1168  Vector<double> div_u_components(2);
1169  data.push_back(interpolated_div_u(s,div_u_components));
1170  data.push_back(div_u_components[0]);
1171  data.push_back(div_u_components[1]);
1172 
1173  }
1174 
1175 
1176  /// Output with default number of plot points
1177  void output(std::ostream &outfile)
1178  {
1179  unsigned nplot=5;
1180  output(outfile,nplot);
1181  }
1182 
1183  /// \short Output FE representation of soln: x,y,u1,u2,div_q,p at
1184  /// Nplot^2 plot points
1185  void output(std::ostream &outfile, const unsigned &nplot);
1186 
1187  /// \short Output incl. projection of fluxes into direction of
1188  /// the specified unit vector
1189  void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot,
1190  const Vector<double> unit_normal);
1191 
1192  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
1193  /// Nplot^2 plot points
1194  void output_fct(std::ostream &outfile,
1195  const unsigned &nplot,
1197 
1198  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
1199  /// Nplot^2 plot points. Unsteady version
1200  void output_fct(std::ostream &outfile,
1201  const unsigned &nplot,
1202  const double &time,
1204 
1205  /// \short Compute the error between the FE solution and the exact solution
1206  /// using the H(div) norm for q and L^2 norm for p
1207  void compute_error(std::ostream &outfile,
1209  Vector<double>& error,
1210  Vector<double>& norm);
1211 
1212  /// \short Compute the error between the FE solution and the exact solution
1213  /// using the H(div) norm for q and L^2 norm for p. Unsteady version
1214  void compute_error(std::ostream &outfile,
1216  const double &time,
1217  Vector<double>& error,
1218  Vector<double>& norm);
1219 
1220 
1221 
1222  // Z2 stuff -- currently based on Darcy flux
1223 
1224  /// Number off flux terms for Z2 error estimator (use Darcy flux)
1226  {
1227  return 2;
1228  }
1229 
1230  /// Z2 flux (use Darcy flux)
1232  {
1233  interpolated_q(s,flux);
1234  }
1235 
1236  protected:
1237 
1238  /// \short Compute the geometric basis, and the q, p and divergence basis
1239  /// functions and test functions at local coordinate s
1240  virtual double shape_basis_test_local(const Vector<double> &s,
1241  Shape &psi,
1242  DShape &dpsi,
1243  Shape &u_basis,
1244  Shape &u_test,
1245  DShape &du_basis_dx,
1246  DShape &du_test_dx,
1247  Shape &q_basis,
1248  Shape &q_test,
1249  Shape &p_basis,
1250  Shape &p_test,
1251  Shape &div_q_basis_ds,
1252  Shape &div_q_test_ds) const = 0;
1253 
1254  /// \short Compute the geometric basis, and the q, p and divergence basis
1255  /// functions and test functions at integration point ipt
1256  virtual double shape_basis_test_local_at_knot(const unsigned &ipt,
1257  Shape &psi,
1258  DShape &dpsi,
1259  Shape &u_basis,
1260  Shape &u_test,
1261  DShape &du_basis_dx,
1262  DShape &du_test_dx,
1263  Shape &q_basis,
1264  Shape &q_test,
1265  Shape &p_basis,
1266  Shape &p_test,
1267  Shape &div_q_basis_ds,
1268  Shape &div_q_test_ds) const = 0;
1269 
1270  // fill in residuals and, if flag==true, jacobian
1272  Vector<double> &residuals,
1273  DenseMatrix<double> &jacobian,
1274  bool flag);
1275 
1276  private:
1277 
1278  /// Pointer to solid body force function
1280 
1281  /// Pointer to fluid source function
1283 
1284  /// Pointer to the mass source function
1286 
1287  /// Pointer to the nondim Young's modulus
1289 
1290  /// Pointer to Poisson's ratio
1291  double* Nu_pt;
1292 
1293  /// Timescale ratio (non-dim. density)
1294  double* Lambda_sq_pt;
1295 
1296  /// Density ratio
1298 
1299  /// permeability
1301 
1302  /// \short Ratio of the material's actual permeability to the permeability
1303  /// used in the non-dimensionalisation of the equations
1305 
1306  /// Alpha -- the biot parameter
1307  double* Alpha_pt;
1308 
1309  /// Porosity
1310  double* Porosity_pt;
1311 
1312  /// Boolean to record that darcy has been switched off
1314 
1315  /// Static default value for Young's modulus (1.0 -- for natural
1316  /// scaling, i.e. all stresses have been non-dimensionalised by
1317  /// the same reference Young's modulus. Setting the "non-dimensional"
1318  /// Young's modulus (obtained by de-referencing Youngs_modulus_pt)
1319  /// to a number larger than one means that the material is stiffer
1320  /// than assumed in the non-dimensionalisation.
1322 
1323  /// Static default value for timescale ratio
1325 
1326  /// Static default value for the density ratio
1328 
1329  /// \short Static default value for the permeability (1.0 for natural scaling;
1330  /// i.e. timescale is given by L^2/(k^* E)
1332 
1333  /// \short Static default value for the ratio of the material's actual
1334  /// permeability to the permeability used to non-dimensionalise the
1335  /// equations
1337 
1338  /// Static default value for alpha, the biot parameter
1339  static double Default_alpha_value;
1340 
1341  /// Static default value for the porosity
1342  static double Default_porosity_value;
1343 
1344  };
1345 
1346 
1347 ////////////////////////////////////////////////////////////////////////
1348 ////////////////////////////////////////////////////////////////////////
1349 ////////////////////////////////////////////////////////////////////////
1350 
1351 
1352 //==========================================================
1353 /// Axisymmetric poro elasticity upgraded to become projectable
1354 //==========================================================
1355  template<class AXISYMMETRIC_POROELASTICITY_ELEMENT>
1357  public virtual ProjectableElement<AXISYMMETRIC_POROELASTICITY_ELEMENT>,
1358  public virtual ProjectableElementBase
1359  {
1360 
1361  public:
1362 
1363  /// \short Constructor [this was only required explicitly
1364  /// from gcc 4.5.2 onwards...]
1366 
1367  /// \short Specify the values associated with field fld.
1368  /// The information is returned in a vector of pairs which comprise
1369  /// the Data object and the value within it, that correspond to field fld.
1371  {
1372 
1373 #ifdef PARANOID
1374  if (fld>3)
1375  {
1376  std::stringstream error_stream;
1377  error_stream
1378  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1379  << fld << " is illegal \n";
1380  throw OomphLibError(
1381  error_stream.str(),
1382  OOMPH_CURRENT_FUNCTION,
1383  OOMPH_EXCEPTION_LOCATION);
1384  }
1385 #endif
1386 
1387  // Create the vector
1388  Vector<std::pair<Data*,unsigned> > data_values;
1389 
1390 
1391  // Pressure
1392  //---------
1393  if (fld==2)
1394  {
1395  Data* dat_pt=this->p_data_pt();
1396  unsigned nvalue=dat_pt->nvalue();
1397  for (unsigned i=0;i<nvalue;i++)
1398  {
1399  data_values.push_back(std::make_pair(dat_pt,i));
1400  }
1401  }
1402  // Darcy flux
1403  //-----------
1404  else if (fld==3)
1405  {
1406  Vector<Data*> edge_dat_pt=this->q_edge_data_pt();
1407  unsigned n=edge_dat_pt.size();
1408  for (unsigned j=0;j<n;j++)
1409  {
1410  unsigned nvalue=this->nedge_flux_interpolation_point();
1411  for (unsigned i=0;i<nvalue;i++)
1412  {
1413  data_values.push_back(std::make_pair(edge_dat_pt[j],
1414  this->q_edge_index(i)));
1415  }
1416  }
1417  if (this->nq_basis_internal()>0)
1418  {
1419  Data* int_dat_pt=this->q_internal_data_pt();
1420  unsigned nvalue=int_dat_pt->nvalue();
1421  for (unsigned i=0;i<nvalue;i++)
1422  {
1423  data_values.push_back(std::make_pair(int_dat_pt,i));
1424  }
1425  }
1426  }
1427  // Solid displacements
1428  else
1429  {
1430  // Loop over all nodes
1431  unsigned nnod=this->nnode();
1432  for (unsigned j=0;j<nnod;j++)
1433  {
1434  // Add the data value associated with the displacement components
1435  // (stored first)
1436  data_values.push_back(std::make_pair(
1437  this->node_pt(j),
1438  this->u_index_axisym_poroelasticity(fld)));
1439  }
1440  }
1441 
1442  // Return the vector
1443  return data_values;
1444  }
1445 
1446  /// \short Number of fields to be projected: 4 (two displacement components,
1447  /// pressure, Darcy flux)
1449  {
1450  return 4;
1451  }
1452 
1453  /// \short Number of history values to be stored for fld-th field.
1454  /// (Note: count includes current value!)
1455  unsigned nhistory_values_for_projection(const unsigned &fld)
1456  {
1457 #ifdef PARANOID
1458  if (fld>3)
1459  {
1460  std::stringstream error_stream;
1461  error_stream
1462  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1463  << fld << " is illegal\n";
1464  throw OomphLibError(
1465  error_stream.str(),
1466  OOMPH_CURRENT_FUNCTION,
1467  OOMPH_EXCEPTION_LOCATION);
1468  }
1469 #endif
1470 
1471  // Displacements -- infer from first (vertex) node
1472  unsigned return_value=this->node_pt(0)->ntstorage();
1473 
1474  // Pressure: No history values (just present one!)
1475  if (fld==2) return_value=1;
1476 
1477  // Flux: infer from first midside node
1478  if (fld==3)
1479  {
1480  return_value=this->node_pt(3)->ntstorage();
1481  }
1482  return return_value;
1483  }
1484 
1485  ///\short Number of positional history values
1486  /// (Note: count includes current value!)
1488  {
1489  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1490  }
1491 
1492  /// \short Return Jacobian of mapping and shape functions of field fld
1493  /// at local coordinate s
1494  double jacobian_and_shape_of_field(const unsigned &fld,
1495  const Vector<double> &s,
1496  Shape &psi)
1497  {
1498 #ifdef PARANOID
1499  if (fld>3)
1500  {
1501  std::stringstream error_stream;
1502  error_stream
1503  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1504  << fld << " is illegal.\n";
1505  throw OomphLibError(
1506  error_stream.str(),
1507  OOMPH_CURRENT_FUNCTION,
1508  OOMPH_EXCEPTION_LOCATION);
1509  }
1510 #endif
1511 
1512 
1513  // Get the number of geometric nodes, total number of basis functions,
1514  // and number of edges basis functions
1515  const unsigned n_dim=this->dim();
1516  const unsigned n_node = this->nnode();
1517  const unsigned n_q_basis = this->nq_basis();
1518  const unsigned n_p_basis = this->np_basis();
1519 
1520  // Storage for the geometric and computational bases and the test functions
1521  Shape psi_geom(n_node), q_basis(n_q_basis,n_dim), q_test(n_q_basis,n_dim),
1522  p_basis(n_p_basis), p_test(n_p_basis),
1523  div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
1524  DShape dpsidx_geom(n_node,n_dim);
1525  Shape u_basis(n_node);
1526  Shape u_test(n_node);
1527  DShape du_basis_dx(n_node,n_dim);
1528  DShape du_test_dx(n_node,n_dim);
1529  double J= this->shape_basis_test_local(s,
1530  psi_geom,
1531  dpsidx_geom,
1532  u_basis,
1533  u_test,
1534  du_basis_dx,
1535  du_test_dx,
1536  q_basis,
1537  q_test,
1538  p_basis,
1539  p_test,
1540  div_q_basis_ds,
1541  div_q_test_ds);
1542  // Pressure basis functions
1543  if (fld==2)
1544  {
1545  unsigned n=p_basis.nindex1();
1546  for (unsigned i=0;i<n;i++)
1547  {
1548  psi[i]=p_basis[i];
1549  }
1550  }
1551  // Flux basis functions
1552  else if (fld==3)
1553  {
1554  unsigned n=q_basis.nindex1();
1555  unsigned m=q_basis.nindex2();
1556  for (unsigned i=0;i<n;i++)
1557  {
1558  for (unsigned j=0;j<m;j++)
1559  {
1560  psi(i,j)=q_basis(i,j);
1561  }
1562  }
1563  }
1564  // Displacement components
1565  else
1566  {
1567  for (unsigned j=0;j<n_node;j++)
1568  {
1569  psi[j]=u_basis[j];
1570  }
1571  }
1572 
1573  return J;
1574  }
1575 
1576 
1577 
1578  /// \short Return interpolated field fld at local coordinate s, at time level
1579  /// t (t=0: present; t>0: history values)
1580  double get_field(const unsigned &t,
1581  const unsigned &fld,
1582  const Vector<double>& s)
1583  {
1584 #ifdef PARANOID
1585  if (fld>3)
1586  {
1587  std::stringstream error_stream;
1588  error_stream
1589  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1590  << fld << " is illegal\n";
1591  throw OomphLibError(
1592  error_stream.str(),
1593  OOMPH_CURRENT_FUNCTION,
1594  OOMPH_EXCEPTION_LOCATION);
1595  }
1596 #endif
1597 
1598  double return_value=0.0;
1599 
1600  // Pressure
1601  if (fld==2)
1602  {
1603  // No time-dependence in here
1604 #ifdef PARANOID
1605  if (t!=0)
1606  {
1607  throw OomphLibError(
1608  "Pressure in AxisymmetricPoroelasticity has no time-dep.!",
1609  OOMPH_CURRENT_FUNCTION,
1610  OOMPH_EXCEPTION_LOCATION);
1611  }
1612 #endif
1613  this->interpolated_p(s,return_value);
1614  }
1615  // Darcy flux -- doesn't really work as it's a vector field
1616  else if (fld==3)
1617  {
1618  throw OomphLibError(
1619  "Don't call this function for AxisymmetricPoroelasticity!",
1620  OOMPH_CURRENT_FUNCTION,
1621  OOMPH_EXCEPTION_LOCATION);
1622  }
1623  // Displacement components
1624  else
1625  {
1626  return_value=this->interpolated_u(t,s,fld);
1627  }
1628  return return_value;
1629  }
1630 
1631 
1632  ///Return number of values in field fld
1633  unsigned nvalue_of_field(const unsigned &fld)
1634  {
1635 #ifdef PARANOID
1636  if (fld>3)
1637  {
1638  std::stringstream error_stream;
1639  error_stream
1640  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1641  << fld << " is illegal\n";
1642  throw OomphLibError(
1643  error_stream.str(),
1644  OOMPH_CURRENT_FUNCTION,
1645  OOMPH_EXCEPTION_LOCATION);
1646  }
1647 #endif
1648 
1649  unsigned return_value=0;
1650  // Pressure
1651  if (fld==2)
1652  {
1653  return_value=this->np_basis();
1654  }
1655  // Darcy flux
1656  else if (fld==3)
1657  {
1658  return_value=this->nq_basis();
1659  }
1660  // Displacements
1661  else
1662  {
1663  return_value=this->nnode();
1664  }
1665 
1666  return return_value;
1667  }
1668 
1669 
1670  ///Return local equation number of value j in field fld.
1671  int local_equation(const unsigned &fld,
1672  const unsigned &j)
1673  {
1674 #ifdef PARANOID
1675  if (fld>3)
1676  {
1677  std::stringstream error_stream;
1678  error_stream
1679  << "AxisymmetricPoroelasticity elements only store 4 fields so fld = "
1680  << fld << " is illegal\n";
1681  throw OomphLibError(
1682  error_stream.str(),
1683  OOMPH_CURRENT_FUNCTION,
1684  OOMPH_EXCEPTION_LOCATION);
1685  }
1686 #endif
1687 
1688  int return_value=0;
1689 
1690  // Pressure
1691  if (fld==2)
1692  {
1693  return_value=this->p_local_eqn(j);
1694  }
1695  // Darcy flux
1696  else if (fld==3)
1697  {
1698  unsigned nedge=this->nq_basis_edge();
1699  if (j<nedge)
1700  {
1701  return_value=this->q_edge_local_eqn(j);
1702  }
1703  else
1704  {
1705  return_value=this->q_internal_local_eqn(j-nedge);
1706  }
1707  }
1708  // Displacement
1709  else
1710  {
1711  return_value=
1712  this->nodal_local_eqn(j,this->u_index_axisym_poroelasticity(fld));
1713  }
1714 
1715  return return_value;
1716  }
1717 
1718 
1719  /// \short Output FE representation of soln as in underlying element
1720  void output(std::ostream &outfile, const unsigned &nplot)
1721  {
1723  }
1724 
1725 
1726  /// \short Residual for the projection step. Flag indicates if we
1727  /// want the Jacobian (1) or not (0). Virtual so it can be
1728  /// overloaded if necessary
1730  DenseMatrix<double> &jacobian,
1731  const unsigned& flag)
1732  {
1733  //Am I a solid element
1734  SolidFiniteElement* solid_el_pt =
1735  dynamic_cast<SolidFiniteElement*>(this);
1736 
1737  unsigned n_dim=this->dim();
1738 
1739  //Allocate storage for local coordinates
1740  Vector<double> s(n_dim);
1741 
1742  //Current field
1743  unsigned fld=this->Projected_field;
1744 
1745  //Number of nodes
1746  const unsigned n_node = this->nnode();
1747 
1748  //Number of positional dofs
1749  const unsigned n_position_type = this->nnodal_position_type();
1750 
1751  //Number of dof for current field
1752  const unsigned n_value=nvalue_of_field(fld);
1753 
1754  //Set the value of n_intpt
1755  const unsigned n_intpt = this->integral_pt()->nweight();
1756 
1757  //Loop over the integration points
1758  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1759  {
1760  // Get the local coordinates of Gauss point
1761  for(unsigned i=0;i<n_dim;i++) s[i] = this->integral_pt()->knot(ipt,i);
1762 
1763  //Get the integral weight
1764  double w = this->integral_pt()->weight(ipt);
1765 
1766  // Find same point in base mesh using external storage
1767  FiniteElement* other_el_pt=0;
1768  other_el_pt=this->external_element_pt(0,ipt);
1769  Vector<double> other_s(n_dim);
1770  other_s=this->external_element_local_coord(0,ipt);
1773  other_el_pt);
1774 
1775  //Storage for the local equation and local unknown
1776  int local_eqn=0, local_unknown=0;
1777 
1778  //Now set up the three different projection problems
1779  switch(Projection_type)
1780  {
1781  case Lagrangian:
1782  {
1783  //If we don't have a solid element there is a problem
1784  if(solid_el_pt==0)
1785  {
1786  throw OomphLibError(
1787  "Trying to project Lagrangian coordinates in non-SolidElement\n",
1788  OOMPH_CURRENT_FUNCTION,
1789  OOMPH_EXCEPTION_LOCATION);
1790  }
1791 
1792  //Find the position shape functions
1793  Shape psi(n_node,n_position_type);
1794  //Get the position shape functions
1795  this->shape(s,psi);
1796  //Get the jacobian
1797  double J = this->J_eulerian(s);
1798 
1799  //Premultiply the weights and the Jacobian
1800  double W = w*J;
1801 
1802  //Get the value of the current position of the 0th coordinate
1803  //in the current element
1804  //at the current time level (which is the unkown)
1805  double interpolated_xi_proj = this->interpolated_x(s,0);
1806 
1807  //Get the Lagrangian position in the other element
1808  double interpolated_xi_bar=
1809  dynamic_cast<SolidFiniteElement*>(cast_el_pt)
1810  ->interpolated_xi(other_s,Projected_lagrangian);
1811 
1812  //Now loop over the nodes and position dofs
1813  for(unsigned l=0;l<n_node;++l)
1814  {
1815  //Loop over position unknowns
1816  for(unsigned k=0;k<n_position_type;++k)
1817  {
1818  //The local equation is going to be the positional local equation
1819  local_eqn = solid_el_pt->position_local_eqn(l,k,0);
1820 
1821  //Now assemble residuals
1822  if(local_eqn >= 0)
1823  {
1824  //calculate residuals
1825  residuals[local_eqn] +=
1826  (interpolated_xi_proj - interpolated_xi_bar)*psi(l,k)*W;
1827 
1828  //Calculate the jacobian
1829  if(flag==1)
1830  {
1831  for(unsigned l2=0;l2<n_node;l2++)
1832  {
1833  //Loop over position dofs
1834  for(unsigned k2=0;k2<n_position_type;k2++)
1835  {
1836  local_unknown =
1837  solid_el_pt->position_local_eqn(l2,k2,0);
1838  if(local_unknown >= 0)
1839  {
1840  jacobian(local_eqn,local_unknown)
1841  += psi(l2,k2)*psi(l,k)*W;
1842  }
1843  }
1844  }
1845  } //end of jacobian
1846  }
1847  }
1848  }
1849  } //End of Lagrangian coordinate case
1850 
1851  break;
1852 
1853  //Now the coordinate history case
1854  case Coordinate:
1855  {
1856  //Find the position shape functions
1857  Shape psi(n_node,n_position_type);
1858  //Get the position shape functions
1859  this->shape(s,psi);
1860  //Get the jacobian
1861  double J = this->J_eulerian(s);
1862 
1863  //Premultiply the weights and the Jacobian
1864  double W = w*J;
1865 
1866  //Get the value of the current position in the current element
1867  //at the current time level (which is the unkown)
1868  double interpolated_x_proj = 0.0;
1869  //If we are a solid element read it out directly from the data
1870  if(solid_el_pt!=0)
1871  {
1872  interpolated_x_proj = this->interpolated_x(s,Projected_coordinate);
1873  }
1874  //Otherwise it's the field value at the current time
1875  else
1876  {
1877  interpolated_x_proj = this->get_field(0,fld,s);
1878  }
1879 
1880  //Get the position in the other element
1881  double interpolated_x_bar=
1883  other_s,Projected_coordinate);
1884 
1885  //Now loop over the nodes and position dofs
1886  for(unsigned l=0;l<n_node;++l)
1887  {
1888  //Loop over position unknowns
1889  for(unsigned k=0;k<n_position_type;++k)
1890  {
1891  //If I'm a solid element
1892  if(solid_el_pt!=0)
1893  {
1894  //The local equation is going to be the positional local equation
1895  local_eqn =
1896  solid_el_pt->position_local_eqn(l,k,Projected_coordinate);
1897  }
1898  //Otherwise just pick the local unknown of a field to
1899  //project into
1900  else
1901  {
1902  //Complain if using generalised position types
1903  //but this is not a solid element, because the storage
1904  //is then not clear
1905  if(n_position_type!=1)
1906  {
1907  throw OomphLibError(
1908  "Trying to project generalised positions not in SolidElement\n",
1909  OOMPH_CURRENT_FUNCTION,
1910  OOMPH_EXCEPTION_LOCATION);
1911  }
1912  local_eqn = local_equation(fld,l);
1913  }
1914 
1915  //Now assemble residuals
1916  if(local_eqn >= 0)
1917  {
1918  //calculate residuals
1919  residuals[local_eqn] +=
1920  (interpolated_x_proj - interpolated_x_bar)*psi(l,k)*W;
1921 
1922  //Calculate the jacobian
1923  if(flag==1)
1924  {
1925  for(unsigned l2=0;l2<n_node;l2++)
1926  {
1927  //Loop over position dofs
1928  for(unsigned k2=0;k2<n_position_type;k2++)
1929  {
1930  //If I'm a solid element
1931  if(solid_el_pt!=0)
1932  {
1933  local_unknown = solid_el_pt
1935  }
1936  else
1937  {
1938  local_unknown = local_equation(fld,l2);
1939  }
1940 
1941  if(local_unknown >= 0)
1942  {
1943  jacobian(local_eqn,local_unknown)
1944  += psi(l2,k2)*psi(l,k)*W;
1945  }
1946  }
1947  }
1948  } //end of jacobian
1949  }
1950  }
1951  }
1952  } //End of coordinate case
1953  break;
1954 
1955  //Now the value case
1956  case Value:
1957  {
1958 
1959  // Pressure or displacements -- "normal" procedure
1960  if (fld<=2)
1961  {
1962  //Field shape function
1963  Shape psi(n_value);
1964 
1965  //Calculate jacobian and shape functions for that field
1966  double J=jacobian_and_shape_of_field(fld,s,psi);
1967 
1968  //Premultiply the weights and the Jacobian
1969  double W = w*J;
1970 
1971  //Value of field in current element at current time level
1972  //(the unknown)
1973  unsigned now=0;
1974  double interpolated_value_proj = this->get_field(now,fld,s);
1975 
1976  //Value of the interpolation of element located in base mesh
1977  double interpolated_value_bar =
1978  cast_el_pt->get_field(Time_level_for_projection,fld,other_s);
1979 
1980  //Loop over dofs of field fld
1981  for(unsigned l=0;l<n_value;l++)
1982  {
1983  local_eqn = local_equation(fld,l);
1984  if(local_eqn >= 0)
1985  {
1986  //calculate residuals
1987  residuals[local_eqn] +=
1988  (interpolated_value_proj - interpolated_value_bar)*psi[l]*W;
1989 
1990  //Calculate the jacobian
1991  if(flag==1)
1992  {
1993  for(unsigned l2=0;l2<n_value;l2++)
1994  {
1995  local_unknown = local_equation(fld,l2);
1996  if(local_unknown >= 0)
1997  {
1998  jacobian(local_eqn,local_unknown)
1999  += psi[l2]*psi[l]*W;
2000  }
2001  }
2002  } //end of jacobian
2003  }
2004  }
2005  }
2006  // Flux -- need inner product
2007  else if (fld==3)
2008  {
2009 
2010  //Field shape function
2011  Shape psi(n_value,n_dim);
2012 
2013  //Calculate jacobian and shape functions for that field
2014  double J=jacobian_and_shape_of_field(fld,s,psi);
2015 
2016  //Premultiply the weights and the Jacobian
2017  double W = w*J;
2018 
2019  //Value of flux in current element at current time level
2020  //(the unknown)
2021  unsigned now=0;
2022  Vector<double> q_proj(n_dim);
2023  this->interpolated_q(now,s,q_proj);
2024 
2025  //Value of the interpolation of element located in base mesh
2026  Vector<double> q_bar(n_dim);
2027  cast_el_pt->interpolated_q(Time_level_for_projection,other_s,q_bar);
2028 
2029  //Loop over dofs of field fld
2030  for(unsigned l=0;l<n_value;l++)
2031  {
2032  local_eqn = local_equation(fld,l);
2033  if(local_eqn >= 0)
2034  {
2035  // Loop over spatial dimension for dot product
2036  for (unsigned i=0;i<n_dim;i++)
2037  {
2038  // Add to residuals
2039  residuals[local_eqn] +=
2040  (q_proj[i] - q_bar[i])*psi(l,i)*W;
2041 
2042  //Calculate the jacobian
2043  if(flag==1)
2044  {
2045  for(unsigned l2=0;l2<n_value;l2++)
2046  {
2047  local_unknown = local_equation(fld,l2);
2048  if(local_unknown >= 0)
2049  {
2050  jacobian(local_eqn,local_unknown)
2051  += psi(l2,i)*psi(l,i)*W;
2052  }
2053  }
2054  } //end of jacobian
2055  }
2056  }
2057  }
2058  }
2059  else
2060  {
2061  throw OomphLibError(
2062  "Wrong field specified. This should never happen\n",
2063  OOMPH_CURRENT_FUNCTION,
2064  OOMPH_EXCEPTION_LOCATION);
2065  }
2066 
2067 
2068  break;
2069 
2070  default:
2071  throw OomphLibError(
2072  "Wrong type specificied in Projection_type. This should never happen\n",
2073  OOMPH_CURRENT_FUNCTION,
2074  OOMPH_EXCEPTION_LOCATION);
2075  }
2076  } //End of the switch statement
2077 
2078  }//End of loop over ipt
2079 
2080  }//End of residual_for_projection function
2081 
2082  };
2083 
2084 
2085 //=======================================================================
2086 /// Face geometry for element is the same as that for the underlying
2087 /// wrapped element
2088 //=======================================================================
2089  template<class ELEMENT>
2091  : public virtual FaceGeometry<ELEMENT>
2092  {
2093  public:
2094  FaceGeometry() : FaceGeometry<ELEMENT>() {}
2095  };
2096 
2097 
2098 
2099 ////////////////////////////////////////////////////////////////////////
2100 ////////////////////////////////////////////////////////////////////////
2101 ////////////////////////////////////////////////////////////////////////
2102 
2103 
2104 
2105 
2106 
2107 
2108 
2109 } // End of oomph namespace
2110 
2111 #endif
2112 
virtual void pin_q_internal_value(const unsigned &j, const double &value)=0
Pin the jth internal q value and set it to specified value.
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...
Projection_Type Projection_type
Variable to indicate if we're projecting the history values of the nodal coordinates (Coordinate) the...
Definition: projection.h:85
virtual void pin_q_edge_value(const unsigned &j, const double &value)=0
Pin the j-th edge (flux) degree of freedom and set it to specified value.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
virtual unsigned q_edge_index(const unsigned &j) const =0
Return the nodal index at which the jth edge unknown is stored.
unsigned num_Z2_flux_terms()
Number off flux terms for Z2 error estimator (use Darcy flux)
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.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
virtual Vector< double > edge_flux_interpolation_point(const unsigned &edge, const unsigned &j) const =0
Returns the local coordinate of the jth flux_interpolation point along the specified edge...
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.
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
double interpolated_div_u(const Vector< double > &s, Vector< double > &div_u_components) const
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 ...
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
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
SourceFctPt & solid_body_force_fct_pt()
Access function: Pointer to solid body force function.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
double interpolated_div_du_dt(const Vector< double > &s, Vector< double > &div_dudt_components) const
virtual unsigned nq_basis_internal() const =0
Return the number of internal basis functions for q.
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
Template-free Base class for projectable elements.
Definition: projection.h:59
unsigned nfields_for_projection()
Number of fields to be projected: 4 (two displacement components, pressure, Darcy flux) ...
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
double interpolated_div_q(const Vector< double > &s) const
Calculate the FE representation of div q and return it.
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
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
char t
Definition: cfortran.h:572
bool darcy_is_switched_off()
Is Darcy flow switched off?
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double *& alpha_pt()
Access function for pointer to alpha, the Biot parameter.
Axisymmetric poro elasticity upgraded to become projectable.
double interpolated_q(const unsigned &t, const Vector< double > &s, const unsigned i) const
unsigned Time_level_for_projection
Time level we are projecting (0=current values; >0: history values)
Definition: projection.h:72
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
double interpolated_u(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Calculate the FE representation of the i-th component of u at time level t (t=0: current) ...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
const double & porosity() const
Access function for the porosity.
double *& density_ratio_pt()
Access function for pointer to the density ratio (fluid to solid)
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting...
Definition: projection.h:76
e
Definition: cfortran.h:575
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
virtual int p_local_eqn(const unsigned &j) const =0
Return the equation number of the j-th pressure degree of freedom.
virtual void edge_flux_interpolation_point_global(const unsigned &edge, const unsigned &j, Vector< double > &x) const =0
Compute the global coordinates of the jth flux_interpolation point along an edge. ...
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
virtual Vector< Data * > q_edge_data_pt() const =0
Return vector of pointers to the Data objects that store the edge flux values.
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 void set_q_internal(const unsigned &j, const double &value)=0
Set the values of the j-th internal degree of freedom.
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Compute the transformed basis at local coordinate s.
SourceFctPt Fluid_body_force_fct_pt
Pointer to fluid source function.
void interpolated_q(const Vector< double > &s, Vector< double > &q) const
Calculate the FE representation of q.
void set_q_internal_timestepper(TimeStepper *const time_stepper_pt)
Set the timestepper of the q internal data object.
virtual void pin_p_value(const unsigned &j, const double &p)=0
Pin the jth pressure value and set it to p.
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.
double interpolated_q(const Vector< double > &s, const unsigned i) const
Calculate the FE representation of the i-th component of q.
void(* SourceFctPt)(const double &time, const Vector< double > &x, Vector< double > &f)
Source function pointer typedef.
ProjectableAxisymmetricPoroelasticityElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
virtual unsigned face_index_of_q_edge_basis_fct(const unsigned &j) const =0
Return the face index associated with j-th edge flux degree of freedom.
double interpolated_p(const Vector< double > &s) const
Calculate the FE representation of p and return it.
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 scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
SourceFctPt fluid_body_force_fct_pt() const
Access function: Pointer to fluid force function (const version)
void(* MassSourceFctPt)(const double &time, const Vector< double > &x, double &f)
Mass source function pointer typedef.
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
virtual unsigned face_index_of_edge(const unsigned &j) const =0
Return the face index associated with specified edge.
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.
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
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Compute the pressure basis.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void set_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set a new timestepper by resizing the appropriate storage. If already assigned the equation numbering...
Definition: nodes.cc:392
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Z2 flux (use Darcy flux)
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
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.
virtual Data * q_internal_data_pt() const =0
Return pointer to the Data object that stores the internal flux values.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals for the Darcy equations.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
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
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
SourceFctPt Solid_body_force_fct_pt
Pointer to solid body force function.
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
double *& youngs_modulus_pt()
Pointer to non-dim Young's modulus (ratio of actual Young's modulus to reference stress used to non-d...
virtual void set_p_value(const unsigned &j, const double &value)=0
Set the jth pressure value.
static double Default_density_ratio_value
Static default value for the density ratio.
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
void output(std::ostream &outfile)
Output with default number of plot points.
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:262
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
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln as in underlying element.
void interpolated_q(const unsigned &t, const Vector< double > &s, Vector< double > &q) const
Calculate the FE representation of q at time level t (t=0: current)
void residual_for_projection(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Residual for the projection step. Flag indicates if we want the Jacobian (1) or not (0)...
virtual int q_internal_local_eqn(const unsigned &j) const =0
Return the equation number of the j-th internal degree of freedom.
virtual unsigned u_index_axisym_poroelasticity(const unsigned &j) const =0
Return the nodal index of the j-th solid displacement unknown.
virtual void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const =0
Compute the local form of the q basis and dbasis/ds at local coordinate s.
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.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
double interpolated_u(const Vector< double > &s, const unsigned &i) const
Calculate the FE representation of the i-th component of u.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
MassSourceFctPt Mass_source_fct_pt
Pointer to the mass source function.
SourceFctPt & fluid_body_force_fct_pt()
Access function: Pointer to fluid force function.
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...
double * Youngs_modulus_pt
Pointer to the nondim Young's modulus.
virtual void face_local_coordinate_of_flux_interpolation_point(const unsigned &edge, const unsigned &n, Vector< double > &s) const =0
Compute the face element coordinates of the nth flux interpolation point along an edge...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
double *& nu_pt()
Access function for pointer to Poisson's ratio.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
void interpolated_du_dt(const Vector< double > &s, Vector< double > &du_dt) const
Calculate the FE representation of du_dt.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
virtual unsigned nedge_flux_interpolation_point() const =0
Returns the number of flux_interpolation points along each edge of the element.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
virtual Data * p_data_pt() const =0
Return pointer to the Data object that stores the pressure values.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
virtual double get_field(const unsigned &time, const unsigned &fld, const Vector< double > &s)=0
Return the fld-th field at local coordinates s at time-level time (time=0: current value; time>0: his...
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.
MassSourceFctPt mass_source_fct_pt() const
Access function: Pointer to mass source function (const version)
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
double *& permeability_pt()
Access function for pointer to the nondim permeability.
static double Default_permeability_ratio_value
Static default value for the ratio of the material's actual permeability to the permeability used to ...
unsigned Projected_field
Field that is currently being projected.
Definition: projection.h:69
virtual void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const =0
Comute the local form of the q basis at local coordinate s.
Class implementing the generic maths of the axisym poroelasticity equations: axisym linear elasticity...
void point_output_data(const Vector< double > &s, Vector< double > &data)
Output solution in data vector at local cordinates s: r,z,u_r,u_z,q_r,q_z,div_q,p,durdt,duzdt.
void interpolated_u(const Vector< double > &s, Vector< double > &disp) const
Calculate the FE representation of u.
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
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 unsigned required_nvalue(const unsigned &n) const =0
Number of values required at node n.
double *& porosity_pt()
Access function for pointer to 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.
bool Darcy_is_switched_off
Boolean to record that darcy has been switched off.
void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot, const Vector< double > unit_normal)
Output incl. projection of fluxes into direction of the specified unit vector.
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:3881
virtual void set_q_edge(const unsigned &j, const double &value)=0
Set the values 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.
MassSourceFctPt & mass_source_fct_pt()
Access function: Pointer to mass source function.
SolidFiniteElement class.
Definition: elements.h:3320
virtual double q_internal(const unsigned &j) const =0
Return the values of the j-th internal degree of freedom.
double * Permeability_ratio_pt
Ratio of the material's actual permeability to the permeability used in the non-dimensionalisation of...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the Jacobian matrix for the Newton method.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
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
SourceFctPt solid_body_force_fct_pt() const
Access function: Pointer to solid body force function (const version)
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.
Definition: projection.h:80
virtual double shape_basis_test_local(const Vector< double > &s, 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 local ...