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_POROELASTICITY_ELEMENTS_HEADER
31 #define OOMPH_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 
41 #include "elasticity_tensor.h"
42 
43 namespace oomph
44 {
45 
46  /// \short Class implementing the generic maths of the poroelasticity
47  /// equations: linear elasticity coupled with Darcy equations (using
48  /// Raviart-Thomas elements with both edge and internal degrees of freedom)
49  template <unsigned DIM>
50  class PoroelasticityEquations : public virtual FiniteElement
51  {
52  public:
53 
54  /// Source function pointer typedef
55  typedef void (*SourceFctPt)(const double &time,
56  const Vector<double>& x,
57  Vector<double>& f);
58 
59  /// Mass source function pointer typedef
60  typedef void (*MassSourceFctPt)(const double &time,
61  const Vector<double>& x,
62  double& f);
63 
64  /// Constructor
74  {
75  }
76 
77  /// Access function for timescale ratio (nondim density)
78  const double& lambda_sq() const {return *Lambda_sq_pt;}
79 
80  /// Access function for pointer to timescale ratio (nondim density)
81  double* &lambda_sq_pt() {return Lambda_sq_pt;}
82 
83  /// Access function for the density ratio
84  const double& density_ratio() const {return *Density_ratio_pt;}
85 
86  /// Access function for pointer to the density ratio
87  double* &density_ratio_pt() {return Density_ratio_pt;}
88 
89  /// Access function for the nondim inverse permeability
90  const double& k_inv() const {return *K_inv_pt;}
91 
92  /// Access function for pointer to the nondim inverse permeability
93  double* &k_inv_pt() {return K_inv_pt;}
94 
95  /// Access function for alpha
96  const double& alpha() const {return *Alpha_pt;}
97 
98  /// Access function for pointer to alpha
99  double* &alpha_pt() {return Alpha_pt;}
100 
101  /// Access function for the porosity
102  const double& porosity() const {return *Porosity_pt;}
103 
104  /// Access function for pointer to the porosity
105  double* &porosity_pt() {return Porosity_pt;}
106 
107  /// Access function: Pointer to solid force function
109 
110  /// Access function: Pointer to solid force function (const version)
112 
113  /// Access function: Pointer to fluid force function
115 
116  /// Access function: Pointer to fluid force function (const version)
118 
119  /// Access function: Pointer to mass source function
121 
122  /// Access function: Pointer to mass source function (const version)
124 
125  /// \short Indirect access to the solid force function - returns 0 if no
126  /// forcing function has been set
127  void force_solid(const double& time,
128  const Vector<double>& x,
129  Vector<double>& b) const
130  {
131  // If no function has been set, return zero vector
132  if(Force_solid_fct_pt==0)
133  {
134  // Get spatial dimension of element
135  unsigned n=dim();
136  for(unsigned i=0;i<n;i++)
137  {
138  b[i] = 0.0;
139  }
140  }
141  else
142  {
143  (*Force_solid_fct_pt)(time,x,b);
144  }
145  }
146 
147  /// \short Indirect access to the fluid forcing function - returns 0 if no
148  /// forcing function has been set
149  void force_fluid(const double& time,
150  const Vector<double>& x,
151  Vector<double>& b) const
152  {
153  // If no function has been set, return zero vector
154  if(Force_fluid_fct_pt==0)
155  {
156  // Get spatial dimension of element
157  unsigned n=dim();
158  for(unsigned i=0;i<n;i++)
159  {
160  b[i] = 0.0;
161  }
162  }
163  else
164  {
165  (*Force_fluid_fct_pt)(time,x,b);
166  }
167  }
168 
169  /// \short Indirect access to the mass source function - returns 0 if no
170  /// mass source function has been set
171  void mass_source(const double& time,
172  const Vector<double>& x,
173  double& b) const
174  {
175  // If no function has been set, return zero vector
176  if(Mass_source_fct_pt==0)
177  {
178  b = 0.0;
179  }
180  else
181  {
182  (*Mass_source_fct_pt)(time,x,b);
183  }
184  }
185 
186  /// Return the pointer to the elasticity_tensor
188 
189  /// Access function to the entries in the elasticity tensor
190  const double E(const unsigned &i, const unsigned &j,
191  const unsigned &k, const unsigned &l) const
192  {
193  return (*Elasticity_tensor_pt)(i,j,k,l);
194  }
195 
196  /// \short Return the Cauchy stress tensor, as calculated
197  /// from the elasticity tensor at specified local coordinate
198  void get_stress(const Vector<double> &s,
199  DenseMatrix<double> &sigma) const;
200 
201  /// \short Return the strain tensor
202  void get_strain(const Vector<double> &s, DenseMatrix<double> &strain) const;
203 
204  /// Number of values required at node n
205  virtual unsigned required_nvalue(const unsigned &n) const = 0;
206 
207  /// Return the nodal index of the n-th solid displacement unknown
208  virtual unsigned u_index(const unsigned &n) const = 0;
209 
210  /// Return the equation number of the n-th edge (flux) degree of freedom
211  virtual int q_edge_local_eqn(const unsigned &n) const = 0;
212 
213  /// Return the equation number of the n-th internal (moment) degree of freedom
214  virtual int q_internal_local_eqn(const unsigned &n) const = 0;
215 
216  /// Return the nodal index at which the nth edge unknown is stored
217  virtual unsigned q_edge_index(const unsigned &n) const = 0;
218 
219  /// \short Return the index of the internal data where the q_internal degrees
220  /// of freedom are stored
221  virtual unsigned q_internal_index() const = 0;
222 
223  /// Return the number of the node where the nth edge unknown is stored
224  virtual unsigned q_edge_node_number(const unsigned &n) const = 0;
225 
226  /// Return the values of the edge (flux) degrees of freedom
227  virtual double q_edge(const unsigned &n) const = 0;
228 
229  /// \short Return the values of the edge (flux) degrees of freedom at time
230  /// history level t
231  virtual double q_edge(const unsigned &t,const unsigned &n) const = 0;
232 
233  /// Return the values of the internal (moment) degrees of freedom
234  virtual double q_internal(const unsigned &n) const = 0;
235 
236  /// \short Return the values of the internal (moment) degrees of freedom at
237  /// time history level t
238  virtual double q_internal(const unsigned &t,const unsigned &n) const = 0;
239 
240  /// Return the total number of computational basis functions for q
241  virtual unsigned nq_basis() const = 0;
242 
243  /// Return the number of edge basis functions for q
244  virtual unsigned nq_basis_edge() const = 0;
245 
246  /// Returns the local form of the q basis at local coordinate s
247  virtual void get_q_basis_local(const Vector<double> &s,
248  Shape &q_basis) const = 0;
249 
250  /// Returns the local form of the q basis and dbasis/ds at local coordinate s
251  virtual void get_div_q_basis_local(const Vector<double> &s,
252  Shape &div_q_basis_ds) const = 0;
253 
254  /// Returns the transformed basis at local coordinate s
256  Shape &q_basis) const
257  {
258  const unsigned n_node = this->nnode();
259  Shape psi(n_node,DIM);
260  const unsigned n_q_basis = this->nq_basis();
261  Shape q_basis_local(n_q_basis,DIM);
262  this->get_q_basis_local(s,q_basis_local);
263  (void)this->transform_basis(s,q_basis_local,psi,q_basis);
264  }
265 
266  /// Returns the number of gauss points along each edge of the element
267  virtual unsigned nedge_gauss_point() const = 0;
268 
269  /// Returns the nth gauss point along an edge
270  virtual double edge_gauss_point(const unsigned &edge,
271  const unsigned &n) const = 0;
272 
273  /// Returns the global coordinates of the nth gauss point along an edge
274  virtual void edge_gauss_point_global(const unsigned &edge,
275  const unsigned &n,
276  Vector<double> &x) const = 0;
277 
278  /// Pin the nth internal q value
279  virtual void pin_q_internal_value(const unsigned &n) = 0;
280 
281  /// Return the equation number of the n-th pressure degree of freedom
282  virtual int p_local_eqn(const unsigned &n) const = 0;
283 
284  /// Return the nth pressure value
285  virtual double p_value(unsigned &n) const = 0;
286 
287  /// Return the total number of pressure basis functions
288  virtual unsigned np_basis() const = 0;
289 
290  /// Return the pressure basis
291  virtual void get_p_basis(const Vector<double> &s,
292  Shape &p_basis) const = 0;
293 
294  /// Pin the nth pressure value
295  virtual void pin_p_value(const unsigned &n, const double &p) = 0;
296 
297  /// Scale the edge basis to allow arbitrary edge mappings
298  virtual void scale_basis(Shape& basis) const = 0;
299 
300  /// \short Performs a div-conserving transformation of the vector basis
301  /// functions from the reference element to the actual element
302  double transform_basis(const Vector<double> &s,
303  const Shape &q_basis_local,
304  Shape &psi,
305  DShape &dpsi,
306  Shape &q_basis) const;
307 
308  /// \short Performs a div-conserving transformation of the vector basis
309  /// functions from the reference element to the actual element
311  const Shape &q_basis_local,
312  Shape &psi,
313  Shape &q_basis) const
314  {
315  const unsigned n_node = this->nnode();
316  DShape dpsi(n_node,DIM);
317  return transform_basis(s,q_basis_local,psi,dpsi,q_basis);
318  }
319 
320  /// Fill in contribution to residuals for the Darcy equations
322  {
325  }
326 
327  /// Fill in the Jacobian matrix for the Newton method
329  DenseMatrix<double> &jacobian)
330  {
331  this->fill_in_generic_residual_contribution(residuals,jacobian,1);
332  }
333 
334  /// Calculate the FE representation of u
336  Vector<double>& disp) const
337  {
338  //Find number of nodes
339  unsigned n_node = nnode();
340 
341  //Local shape function
342  Shape psi(n_node);
343 
344  //Find values of shape function
345  shape(s,psi);
346 
347  for(unsigned i=0;i<DIM;i++)
348  {
349  //Index at which the nodal value is stored
350  unsigned u_nodal_index = u_index(i);
351 
352  //Initialise value of u
353  disp[i] = 0.0;
354 
355  //Loop over the local nodes and sum
356  for(unsigned l=0;l<n_node;l++)
357  {
358  disp[i] += nodal_value(l,u_nodal_index)*psi[l];
359  }
360  }
361  }
362 
363  /// Calculate the FE representation of the i-th component of u
365  const unsigned &i) const
366  {
367  //Find number of nodes
368  unsigned n_node = nnode();
369 
370  //Local shape function
371  Shape psi(n_node);
372 
373  //Find values of shape function
374  shape(s,psi);
375 
376  //Get nodal index at which i-th velocity is stored
377  unsigned u_nodal_index = u_index(i);
378 
379  //Initialise value of u
380  double interpolated_u = 0.0;
381 
382  //Loop over the local nodes and sum
383  for(unsigned l=0;l<n_node;l++)
384  {
385  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
386  }
387 
388  return(interpolated_u);
389  }
390 
391  /// Calculate the FE representation of q
393  Vector<double> &u) const
394  {
395  unsigned n_q_basis = nq_basis();
396  unsigned n_q_basis_edge = nq_basis_edge();
397 
398  Shape q_basis(n_q_basis,DIM);
399 
400  get_q_basis(s,q_basis);
401  for(unsigned i=0;i<DIM;i++)
402  {
403  u[i]=0.0;
404  for(unsigned l=0;l<n_q_basis_edge;l++)
405  {
406  u[i] += q_edge(l)*q_basis(l,i);
407  }
408  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
409  {
410  u[i] += q_internal(l-n_q_basis_edge)*q_basis(l,i);
411  }
412  }
413  }
414 
415  /// Calculate the FE representation of the i-th component of q
417  const unsigned i) const
418  {
419  unsigned n_q_basis = nq_basis();
420  unsigned n_q_basis_edge = nq_basis_edge();
421 
422  Shape q_basis(n_q_basis,DIM);
423 
424  get_q_basis(s,q_basis);
425  double q_i=0.0;
426  for(unsigned l=0;l<n_q_basis_edge;l++)
427  {
428  q_i += q_edge(l)*q_basis(l,i);
429  }
430  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
431  {
432  q_i += q_internal(l-n_q_basis_edge)*q_basis(l,i);
433  }
434 
435  return q_i;
436  }
437 
438  /// Calculate the FE representation of div u
439  void interpolated_div_q(const Vector<double> &s, double &div_q) const
440  {
441  // Zero the divergence
442  div_q=0;
443 
444  // Get the number of nodes, q basis function, and q edge basis functions
445  unsigned n_node=nnode();
446  const unsigned n_q_basis = nq_basis();
447  const unsigned n_q_basis_edge = nq_basis_edge();
448 
449  // Storage for the divergence basis
450  Shape div_q_basis_ds(n_q_basis);
451 
452  // Storage for the geometric basis and it's derivatives
453  Shape psi(n_node);
454  DShape dpsi(n_node,DIM);
455 
456  // Call the geometric shape functions and their derivatives
457  this->dshape_local(s,psi,dpsi);
458 
459  // Storage for the inverse of the geometric jacobian (just so we can call
460  // the local to eulerian mapping)
461  DenseMatrix<double> inverse_jacobian(DIM);
462 
463  // Get the determinant of the geometric mapping
464  double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
465 
466  // Get the divergence basis (wrt local coords) at local coords s
467  get_div_q_basis_local(s,div_q_basis_ds);
468 
469  // Add the contribution to the divergence from the edge basis functions
470  for(unsigned l=0;l<n_q_basis_edge;l++)
471  {
472  div_q+=1.0/det*div_q_basis_ds(l)*q_edge(l);
473  }
474 
475  // Add the contribution to the divergence from the internal basis functions
476  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
477  {
478  div_q+=1.0/det*div_q_basis_ds(l)*q_internal(l-n_q_basis_edge);
479  }
480  }
481 
482  /// Calculate the FE representation of div q and return it
484  {
485  // Temporary storage for div u
486  double div_q=0;
487 
488  // Get the intepolated divergence
489  interpolated_div_q(s,div_q);
490 
491  // Return it
492  return div_q;
493  }
494 
495  /// Calculate the FE representation of p
497  double &p) const
498  {
499  // Get the number of p basis functions
500  unsigned n_p_basis = np_basis();
501 
502  // Storage for the p basis
503  Shape p_basis(n_p_basis);
504 
505  // Call the p basis
506  get_p_basis(s,p_basis);
507 
508  // Zero the pressure
509  p=0;
510 
511  // Add the contribution to the pressure from each basis function
512  for(unsigned l=0;l<n_p_basis;l++)
513  {
514  p+=p_value(l)*p_basis(l);
515  }
516  }
517 
518  /// Calculate the FE representation of p and return it
519  double interpolated_p(const Vector<double> &s) const
520  {
521  // Temporary storage for p
522  double p=0;
523 
524  // Get the interpolated pressure
525  interpolated_p(s,p);
526 
527  // Return it
528  return p;
529  }
530 
531  /// du/dt at local node n
532  double du_dt(const unsigned &n,
533  const unsigned &i) const
534  {
535  // Get the timestepper
537 
538  // Storage for the derivative - initialise to 0
539  double du_dt=0.0;
540 
541  // If we are doing an unsteady solve then calculate the derivative
542  if(!time_stepper_pt->is_steady())
543  {
544  // Get the nodal index
545  const unsigned u_nodal_index=u_index(i);
546 
547  // Get the number of values required to represent history
548  const unsigned n_time=time_stepper_pt->ntstorage();
549 
550  // Loop over history values
551  for(unsigned t=0;t<n_time;t++)
552  {
553  //Add the contribution to the derivative
554  du_dt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
555  }
556  }
557 
558  return du_dt;
559  }
560 
561  /// d^2u/dt^2 at local node n
562  double d2u_dt2(const unsigned &n,
563  const unsigned &i) const
564  {
565  // Get the timestepper
567 
568  // Storage for the derivative - initialise to 0
569  double d2u_dt2=0.0;
570 
571  // If we are doing an unsteady solve then calculate the derivative
572  if(!time_stepper_pt->is_steady())
573  {
574  // Get the nodal index
575  const unsigned u_nodal_index=u_index(i);
576 
577  // Get the number of values required to represent history
578  const unsigned n_time=time_stepper_pt->ntstorage();
579 
580  // Loop over history values
581  for(unsigned t=0;t<n_time;t++)
582  {
583  //Add the contribution to the derivative
584  d2u_dt2+=time_stepper_pt->weight(2,t)*nodal_value(t,n,u_nodal_index);
585  }
586  }
587 
588  return d2u_dt2;
589  }
590 
591  /// dq_edge/dt for the n-th edge degree of freedom
592  double dq_edge_dt(const unsigned &n) const
593  {
594  unsigned node_num=q_edge_node_number(n);
595 
596  // get the timestepper
598 
599  // storage for the derivative - initialise to 0
600  double dq_dt=0.0;
601 
602  // if we are doing an unsteady solve then calculate the derivative
603  if(!time_stepper_pt->is_steady())
604  {
605  // get the number of values required to represent history
606  const unsigned n_time=time_stepper_pt->ntstorage();
607 
608  // loop over history values
609  for(unsigned t=0;t<n_time;t++)
610  {
611  // add the contribution to the derivative
612  dq_dt+=
613  time_stepper_pt->weight(1,t)*q_edge(t,n);
614  }
615  }
616 
617  return dq_dt;
618  }
619 
620  /// dq_internal/dt for the n-th internal degree of freedom
621  double dq_internal_dt(const unsigned &n) const
622  {
623  // get the internal data index for q
624  unsigned internal_index=q_internal_index();
625 
626  // get the timestepper
628  internal_data_pt(internal_index)->time_stepper_pt();
629 
630  // storage for the derivative - initialise to 0
631  double dq_dt=0.0;
632 
633  // if we are doing an unsteady solve then calculate the derivative
634  if(!time_stepper_pt->is_steady())
635  {
636  // get the number of values required to represent history
637  const unsigned n_time=time_stepper_pt->ntstorage();
638 
639  // loop over history values
640  for(unsigned t=0;t<n_time;t++)
641  {
642  // add the contribution to the derivative
643  dq_dt+=time_stepper_pt->weight(1,t)*q_internal(t,n);
644  }
645  }
646 
647  return dq_dt;
648  }
649 
650  /// Set the timestepper of the q internal data object
652  {
653  unsigned q_index=q_internal_index();
654 
655  this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt,false);
656  }
657 
658  unsigned self_test()
659  {
660  return 0;
661  }
662 
663  /// Output with default number of plot points
664  void output(std::ostream &outfile)
665  {
666  unsigned nplot=5;
667  output(outfile,nplot);
668  }
669 
670  /// \short Output FE representation of soln: x,y,u1,u2,div_q,p at
671  /// Nplot^DIM plot points
672  void output(std::ostream &outfile, const unsigned &nplot);
673 
674  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
675  /// Nplot^DIM plot points
676  void output_fct(std::ostream &outfile,
677  const unsigned &nplot,
679 
680  /// \short Output FE representation of exact soln: x,y,u1,u2,div_q,p at
681  /// Nplot^DIM plot points. Unsteady version
682  void output_fct(std::ostream &outfile,
683  const unsigned &nplot,
684  const double &time,
686 
687  /// \short Compute the error between the FE solution and the exact solution
688  /// using the H(div) norm for q and L^2 norm for p
689  void compute_error(std::ostream &outfile,
691  Vector<double>& error,
692  Vector<double>& norm);
693 
694  /// \short Compute the error between the FE solution and the exact solution
695  /// using the H(div) norm for q and L^2 norm for p. Unsteady version
696  void compute_error(std::ostream &outfile,
698  const double &time,
699  Vector<double>& error,
700  Vector<double>& norm);
701 
702  protected:
703 
704  /// \short Returns the geometric basis, and the q, p and divergence basis
705  /// functions and test functions at local coordinate s
706  virtual double shape_basis_test_local(const Vector<double> &s,
707  Shape &psi,
708  DShape &dpsi,
709  Shape &u_basis,
710  Shape &u_test,
711  DShape &du_basis_dx,
712  DShape &du_test_dx,
713  Shape &q_basis,
714  Shape &q_test,
715  Shape &p_basis,
716  Shape &p_test,
717  Shape &div_q_basis_ds,
718  Shape &div_q_test_ds) const = 0;
719 
720  /// \short Returns the geometric basis, and the q, p and divergence basis
721  /// functions and test functions at integration point ipt
722  virtual double shape_basis_test_local_at_knot(const unsigned &ipt,
723  Shape &psi,
724  DShape &dpsi,
725  Shape &u_basis,
726  Shape &u_test,
727  DShape &du_basis_dx,
728  DShape &du_test_dx,
729  Shape &q_basis,
730  Shape &q_test,
731  Shape &p_basis,
732  Shape &p_test,
733  Shape &div_q_basis_ds,
734  Shape &div_q_test_ds) const = 0;
735 
736  // fill in residuals and, if flag==true, jacobian
738  Vector<double> &residuals,
739  DenseMatrix<double> &jacobian,
740  bool flag);
741 
742  /// Pointer to the elasticity tensor
744 
745  private:
746 
747  /// Pointer to solid source function
749 
750  /// Pointer to fluid source function
752 
753  /// Pointer to the mass source function
755 
756  /// Timescale ratio (non-dim. density)
757  double* Lambda_sq_pt;
758 
759  /// Density ratio
761 
762  /// 1/k
763  double* K_inv_pt;
764 
765  /// Alpha
766  double* Alpha_pt;
767 
768  /// Porosity
769  double* Porosity_pt;
770 
771  /// Static default value for timescale ratio (1.0 -- for natural scaling)
772  static double Default_lambda_sq_value;
773 
774  /// Static default value for the density ratio
776 
777  /// Static default value for 1/k
778  static double Default_k_inv_value;
779 
780  /// Static default value for alpha
781  static double Default_alpha_value;
782 
783  /// Static default value for the porosity
784  static double Default_porosity_value;
785 
786  };
787 
788 } // End of oomph namespace
789 
790 #endif
791 
SourceFctPt Force_fluid_fct_pt
Pointer to fluid source function.
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div u.
virtual int p_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th pressure degree of freedom.
void force_solid(const double &time, const Vector< double > &x, Vector< double > &b) const
Indirect access to the solid force function - returns 0 if no forcing function has been set...
static double Default_k_inv_value
Static default value for 1/k.
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, DShape &dpsi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
virtual double q_edge(const unsigned &n) const =0
Return the values of the edge (flux) degrees of freedom.
void(* MassSourceFctPt)(const double &time, const Vector< double > &x, double &f)
Mass source function pointer typedef.
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
Returns the geometric basis, and the q, p and divergence basis functions and test functions at integr...
void output(std::ostream &outfile)
Output with default number of plot points.
cstr elem_len * i
Definition: cfortran.h:607
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Returns the transformed basis at local coordinate s.
virtual double p_value(unsigned &n) const =0
Return the nth pressure value.
double dq_edge_dt(const unsigned &n) const
dq_edge/dt for the n-th edge degree of freedom
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the Jacobian matrix for the Newton method.
virtual unsigned required_nvalue(const unsigned &n) const =0
Number of values required at node n.
SourceFctPt force_fluid_fct_pt() const
Access function: Pointer to fluid force function (const version)
const double E(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Access function to the entries in the elasticity tensor.
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
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 & k_inv() const
Access function for the nondim inverse permeability.
A general Finite Element class.
Definition: elements.h:1271
virtual unsigned nedge_gauss_point() const =0
Returns the number of gauss points along each edge of the element.
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
static double Default_porosity_value
Static default value for the porosity.
MassSourceFctPt mass_source_fct_pt() const
Access function: Pointer to mass source function (const version)
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^DIM plot points.
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
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
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
void(* SourceFctPt)(const double &time, const Vector< double > &x, Vector< double > &f)
Source function pointer typedef.
void interpolated_q(const Vector< double > &s, Vector< double > &u) const
Calculate the FE representation of q.
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.
const double & density_ratio() const
Access function for the density ratio.
virtual int q_internal_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th internal (moment) degree of freedom.
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
unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
double interpolated_q(const Vector< double > &s, const unsigned i) const
Calculate the FE representation of the i-th component of q.
SourceFctPt & force_solid_fct_pt()
Access function: Pointer to solid force function.
SourceFctPt & force_fluid_fct_pt()
Access function: Pointer to fluid force function.
virtual unsigned u_index(const unsigned &n) const =0
Return the nodal index of the n-th solid displacement unknown.
virtual int q_edge_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th edge (flux) degree of freedom.
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
Returns the geometric basis, and the q, p and divergence basis functions and test functions at local ...
MassSourceFctPt Mass_source_fct_pt
Pointer to the mass source function.
static double Default_alpha_value
Static default value for alpha.
const double & alpha() const
Access function for alpha.
MassSourceFctPt & mass_source_fct_pt()
Access function: Pointer to mass source function.
void set_q_internal_timestepper(TimeStepper *const time_stepper_pt)
Set the timestepper of the q internal data object.
double d2u_dt2(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
virtual void pin_p_value(const unsigned &n, const double &p)=0
Pin the nth pressure value.
static char t char * s
Definition: cfortran.h:572
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...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Return the pressure basis.
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
virtual void pin_q_internal_value(const unsigned &n)=0
Pin the nth internal q value.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
ElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
virtual unsigned q_internal_index() const =0
Return the index of the internal data where the q_internal degrees of freedom are stored...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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 ...
Class implementing the generic maths of the poroelasticity equations: linear elasticity coupled with ...
double interpolated_u(const Vector< double > &s, const unsigned &i) const
Calculate the FE representation of the i-th component of u.
virtual double edge_gauss_point(const unsigned &edge, const unsigned &n) const =0
Returns the nth gauss point along an edge.
virtual void edge_gauss_point_global(const unsigned &edge, const unsigned &n, Vector< double > &x) const =0
Returns the global coordinates of the nth gauss point along an edge.
double du_dt(const unsigned &n, const unsigned &i) const
du/dt at local node n
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
double dq_internal_dt(const unsigned &n) const
dq_internal/dt for the n-th internal degree of freedom
SourceFctPt Force_solid_fct_pt
Pointer to solid source function.
virtual double q_internal(const unsigned &n) const =0
Return the values of the internal (moment) degrees of freedom.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
double interpolated_div_q(const Vector< double > &s)
Calculate the FE representation of div q and return it.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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 *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
double * Density_ratio_pt
Density ratio.
void force_fluid(const double &time, const Vector< double > &x, Vector< double > &b) const
Indirect access to the fluid forcing function - returns 0 if no forcing function has been set...
double *& density_ratio_pt()
Access function for pointer to the density ratio.
virtual unsigned q_edge_node_number(const unsigned &n) const =0
Return the number of the node where the nth edge unknown is stored.
double *& alpha_pt()
Access function for pointer to alpha.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals for the Darcy equations.
void interpolated_u(const Vector< double > &s, Vector< double > &disp) const
Calculate the FE representation of u.
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
virtual void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const =0
Returns the local form of the q basis and dbasis/ds at local coordinate s.
virtual void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const =0
Returns the local form of the q basis at local coordinate s.
const double & porosity() const
Access function for the porosity.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
virtual unsigned nq_basis() const =0
Return the total number of computational basis functions for q.
double interpolated_p(const Vector< double > &s) const
Calculate the FE representation of p and return it.
SourceFctPt force_solid_fct_pt() const
Access function: Pointer to solid force function (const version)
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...
double *& k_inv_pt()
Access function for pointer to the nondim inverse permeability.
static double Default_density_ratio_value
Static default value for the density ratio.
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
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
virtual unsigned q_edge_index(const unsigned &n) const =0
Return the nodal index at which the nth edge unknown is stored.
double *& porosity_pt()
Access function for pointer to the porosity.
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 ...