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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_POROELASTICITY_ELEMENTS_HEADER
27#define OOMPH_POROELASTICITY_ELEMENTS_HEADER
28
29// Config header generated by autoconfig
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34#include "../generic/elements.h"
35#include "../generic/shape.h"
36
37#include "elasticity_tensor.h"
38
39namespace oomph
40{
41 /// Class implementing the generic maths of the poroelasticity
42 /// equations: linear elasticity coupled with Darcy equations (using
43 /// Raviart-Thomas elements with both edge and internal degrees of freedom)
44 template<unsigned DIM>
46 {
47 public:
48 /// Source function pointer typedef
49 typedef void (*SourceFctPt)(const double& time,
50 const Vector<double>& x,
52
53 /// Mass source function pointer typedef
54 typedef void (*MassSourceFctPt)(const double& time,
55 const Vector<double>& x,
56 double& f);
57
58 /// Constructor
69 {
70 }
71
72 /// Access function for timescale ratio (nondim density)
73 const double& lambda_sq() const
74 {
75 return *Lambda_sq_pt;
76 }
77
78 /// Access function for pointer to timescale ratio (nondim density)
79 double*& lambda_sq_pt()
80 {
81 return Lambda_sq_pt;
82 }
83
84 /// Access function for the density ratio
85 const double& density_ratio() const
86 {
87 return *Density_ratio_pt;
88 }
89
90 /// Access function for pointer to the density ratio
92 {
93 return Density_ratio_pt;
94 }
95
96 /// Access function for the nondim inverse permeability
97 const double& k_inv() const
98 {
99 return *K_inv_pt;
100 }
101
102 /// Access function for pointer to the nondim inverse permeability
103 double*& k_inv_pt()
104 {
105 return K_inv_pt;
106 }
107
108 /// Access function for alpha
109 const double& alpha() const
110 {
111 return *Alpha_pt;
112 }
113
114 /// Access function for pointer to alpha
115 double*& alpha_pt()
116 {
117 return Alpha_pt;
118 }
119
120 /// Access function for the porosity
121 const double& porosity() const
122 {
123 return *Porosity_pt;
124 }
125
126 /// Access function for pointer to the porosity
127 double*& porosity_pt()
128 {
129 return Porosity_pt;
130 }
131
132 /// Access function: Pointer to solid force function
134 {
135 return Force_solid_fct_pt;
136 }
137
138 /// Access function: Pointer to solid force function (const version)
140 {
141 return Force_solid_fct_pt;
142 }
143
144 /// Access function: Pointer to fluid force function
146 {
147 return Force_fluid_fct_pt;
148 }
149
150 /// Access function: Pointer to fluid force function (const version)
152 {
153 return Force_fluid_fct_pt;
154 }
155
156 /// Access function: Pointer to mass source function
158 {
159 return Mass_source_fct_pt;
160 }
161
162 /// Access function: Pointer to mass source function (const version)
164 {
165 return Mass_source_fct_pt;
166 }
167
168 /// Indirect access to the solid force function - returns 0 if no
169 /// forcing function has been set
170 void force_solid(const double& time,
171 const Vector<double>& x,
172 Vector<double>& b) const
173 {
174 // If no function has been set, return zero vector
175 if (Force_solid_fct_pt == 0)
176 {
177 // Get spatial dimension of element
178 unsigned n = dim();
179 for (unsigned i = 0; i < n; i++)
180 {
181 b[i] = 0.0;
182 }
183 }
184 else
185 {
186 (*Force_solid_fct_pt)(time, x, b);
187 }
188 }
189
190 /// Indirect access to the fluid forcing function - returns 0 if no
191 /// forcing function has been set
192 void force_fluid(const double& time,
193 const Vector<double>& x,
194 Vector<double>& b) const
195 {
196 // If no function has been set, return zero vector
197 if (Force_fluid_fct_pt == 0)
198 {
199 // Get spatial dimension of element
200 unsigned n = dim();
201 for (unsigned i = 0; i < n; i++)
202 {
203 b[i] = 0.0;
204 }
205 }
206 else
207 {
208 (*Force_fluid_fct_pt)(time, x, b);
209 }
210 }
211
212 /// Indirect access to the mass source function - returns 0 if no
213 /// mass source function has been set
214 void mass_source(const double& time,
215 const Vector<double>& x,
216 double& b) const
217 {
218 // If no function has been set, return zero vector
219 if (Mass_source_fct_pt == 0)
220 {
221 b = 0.0;
222 }
223 else
224 {
225 (*Mass_source_fct_pt)(time, x, b);
226 }
227 }
228
229 /// Return the pointer to the elasticity_tensor
231 {
233 }
234
235 /// Access function to the entries in the elasticity tensor
236 const double E(const unsigned& i,
237 const unsigned& j,
238 const unsigned& k,
239 const unsigned& l) const
240 {
241 return (*Elasticity_tensor_pt)(i, j, k, l);
242 }
243
244 /// Return the Cauchy stress tensor, as calculated
245 /// from the elasticity tensor at specified local coordinate
246 void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma) const;
247
248 /// Return the strain tensor
249 void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const;
250
251 /// Number of values required at node n
252 virtual unsigned required_nvalue(const unsigned& n) const = 0;
253
254 /// Return the nodal index of the n-th solid displacement unknown
255 virtual unsigned u_index(const unsigned& n) const = 0;
256
257 /// Return the equation number of the n-th edge (flux) degree of freedom
258 virtual int q_edge_local_eqn(const unsigned& n) const = 0;
259
260 /// Return the equation number of the n-th internal (moment) degree of
261 /// freedom
262 virtual int q_internal_local_eqn(const unsigned& n) const = 0;
263
264 /// Return the nodal index at which the nth edge unknown is stored
265 virtual unsigned q_edge_index(const unsigned& n) const = 0;
266
267 /// Return the index of the internal data where the q_internal
268 /// degrees of freedom are stored
269 virtual unsigned q_internal_index() const = 0;
270
271 /// Return the number of the node where the nth edge unknown is stored
272 virtual unsigned q_edge_node_number(const unsigned& n) const = 0;
273
274 /// Return the values of the edge (flux) degrees of freedom
275 virtual double q_edge(const unsigned& n) const = 0;
276
277 /// Return the values of the edge (flux) degrees of freedom at time
278 /// history level t
279 virtual double q_edge(const unsigned& t, const unsigned& n) const = 0;
280
281 /// Return the values of the internal (moment) degrees of freedom
282 virtual double q_internal(const unsigned& n) const = 0;
283
284 /// Return the values of the internal (moment) degrees of freedom at
285 /// time history level t
286 virtual double q_internal(const unsigned& t, const unsigned& n) const = 0;
287
288 /// Return the total number of computational basis functions for q
289 virtual unsigned nq_basis() const = 0;
290
291 /// Return the number of edge basis functions for q
292 virtual unsigned nq_basis_edge() const = 0;
293
294 /// Returns the local form of the q basis at local coordinate s
295 virtual void get_q_basis_local(const Vector<double>& s,
296 Shape& q_basis) const = 0;
297
298 /// Returns the local form of the q basis and dbasis/ds at local coordinate
299 /// s
301 Shape& div_q_basis_ds) const = 0;
302
303 /// Returns the transformed basis at local coordinate s
304 void get_q_basis(const Vector<double>& s, Shape& q_basis) const
305 {
306 const unsigned n_node = this->nnode();
307 Shape psi(n_node, DIM);
308 const unsigned n_q_basis = this->nq_basis();
309 Shape q_basis_local(n_q_basis, DIM);
310 this->get_q_basis_local(s, q_basis_local);
311 (void)this->transform_basis(s, q_basis_local, psi, q_basis);
312 }
313
314 /// Returns the number of gauss points along each edge of the element
315 virtual unsigned nedge_gauss_point() const = 0;
316
317 /// Returns the nth gauss point along an edge
318 virtual double edge_gauss_point(const unsigned& edge,
319 const unsigned& n) const = 0;
320
321 /// Returns the global coordinates of the nth gauss point along an edge
322 virtual void edge_gauss_point_global(const unsigned& edge,
323 const unsigned& n,
324 Vector<double>& x) const = 0;
325
326 /// Pin the nth internal q value
327 virtual void pin_q_internal_value(const unsigned& n) = 0;
328
329 /// Return the equation number of the n-th pressure degree of freedom
330 virtual int p_local_eqn(const unsigned& n) const = 0;
331
332 /// Return the nth pressure value
333 virtual double p_value(unsigned& n) const = 0;
334
335 /// Return the total number of pressure basis functions
336 virtual unsigned np_basis() const = 0;
337
338 /// Return the pressure basis
339 virtual void get_p_basis(const Vector<double>& s, Shape& p_basis) const = 0;
340
341 /// Pin the nth pressure value
342 virtual void pin_p_value(const unsigned& n, const double& p) = 0;
343
344 /// Scale the edge basis to allow arbitrary edge mappings
345 virtual void scale_basis(Shape& basis) const = 0;
346
347 /// Performs a div-conserving transformation of the vector basis
348 /// functions from the reference element to the actual element
349 double transform_basis(const Vector<double>& s,
350 const Shape& q_basis_local,
351 Shape& psi,
352 DShape& dpsi,
353 Shape& q_basis) const;
354
355 /// Performs a div-conserving transformation of the vector basis
356 /// functions from the reference element to the actual element
358 const Shape& q_basis_local,
359 Shape& psi,
360 Shape& q_basis) const
361 {
362 const unsigned n_node = this->nnode();
363 DShape dpsi(n_node, DIM);
364 return transform_basis(s, q_basis_local, psi, dpsi, q_basis);
365 }
366
367 /// Fill in contribution to residuals for the Darcy equations
369 {
372 }
373
374 /// Fill in the Jacobian matrix for the Newton method
376 DenseMatrix<double>& jacobian)
377 {
378 this->fill_in_generic_residual_contribution(residuals, jacobian, 1);
379 }
380
381 /// Calculate the FE representation of u
383 {
384 // Find number of nodes
385 unsigned n_node = nnode();
386
387 // Local shape function
388 Shape psi(n_node);
389
390 // Find values of shape function
391 shape(s, psi);
392
393 for (unsigned i = 0; i < DIM; i++)
394 {
395 // Index at which the nodal value is stored
396 unsigned u_nodal_index = u_index(i);
397
398 // Initialise value of u
399 disp[i] = 0.0;
400
401 // Loop over the local nodes and sum
402 for (unsigned l = 0; l < n_node; l++)
403 {
404 disp[i] += nodal_value(l, u_nodal_index) * psi[l];
405 }
406 }
407 }
408
409 /// Calculate the FE representation of the i-th component of u
410 double interpolated_u(const Vector<double>& s, const unsigned& i) const
411 {
412 // Find number of nodes
413 unsigned n_node = nnode();
414
415 // Local shape function
416 Shape psi(n_node);
417
418 // Find values of shape function
419 shape(s, psi);
420
421 // Get nodal index at which i-th velocity is stored
422 unsigned u_nodal_index = u_index(i);
423
424 // Initialise value of u
425 double interpolated_u = 0.0;
426
427 // Loop over the local nodes and sum
428 for (unsigned l = 0; l < n_node; l++)
429 {
430 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
431 }
432
433 return (interpolated_u);
434 }
435
436 /// Calculate the FE representation of q
438 {
439 unsigned n_q_basis = nq_basis();
440 unsigned n_q_basis_edge = nq_basis_edge();
441
442 Shape q_basis(n_q_basis, DIM);
443
444 get_q_basis(s, q_basis);
445 for (unsigned i = 0; i < DIM; i++)
446 {
447 u[i] = 0.0;
448 for (unsigned l = 0; l < n_q_basis_edge; l++)
449 {
450 u[i] += q_edge(l) * q_basis(l, i);
451 }
452 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
453 {
454 u[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
455 }
456 }
457 }
458
459 /// Calculate the FE representation of the i-th component of q
460 double interpolated_q(const Vector<double>& s, const unsigned i) const
461 {
462 unsigned n_q_basis = nq_basis();
463 unsigned n_q_basis_edge = nq_basis_edge();
464
465 Shape q_basis(n_q_basis, DIM);
466
467 get_q_basis(s, q_basis);
468 double q_i = 0.0;
469 for (unsigned l = 0; l < n_q_basis_edge; l++)
470 {
471 q_i += q_edge(l) * q_basis(l, i);
472 }
473 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
474 {
475 q_i += q_internal(l - n_q_basis_edge) * q_basis(l, i);
476 }
477
478 return q_i;
479 }
480
481 /// Calculate the FE representation of div u
482 void interpolated_div_q(const Vector<double>& s, double& div_q) const
483 {
484 // Zero the divergence
485 div_q = 0;
486
487 // Get the number of nodes, q basis function, and q edge basis functions
488 unsigned n_node = nnode();
489 const unsigned n_q_basis = nq_basis();
490 const unsigned n_q_basis_edge = nq_basis_edge();
491
492 // Storage for the divergence basis
493 Shape div_q_basis_ds(n_q_basis);
494
495 // Storage for the geometric basis and it's derivatives
496 Shape psi(n_node);
497 DShape dpsi(n_node, DIM);
498
499 // Call the geometric shape functions and their derivatives
500 this->dshape_local(s, psi, dpsi);
501
502 // Storage for the inverse of the geometric jacobian (just so we can call
503 // the local to eulerian mapping)
504 DenseMatrix<double> inverse_jacobian(DIM);
505
506 // Get the determinant of the geometric mapping
507 double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
508
509 // Get the divergence basis (wrt local coords) at local coords s
510 get_div_q_basis_local(s, div_q_basis_ds);
511
512 // Add the contribution to the divergence from the edge basis functions
513 for (unsigned l = 0; l < n_q_basis_edge; l++)
514 {
515 div_q += 1.0 / det * div_q_basis_ds(l) * q_edge(l);
516 }
517
518 // Add the contribution to the divergence from the internal basis
519 // functions
520 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
521 {
522 div_q += 1.0 / det * div_q_basis_ds(l) * q_internal(l - n_q_basis_edge);
523 }
524 }
525
526 /// Calculate the FE representation of div q and return it
528 {
529 // Temporary storage for div u
530 double div_q = 0;
531
532 // Get the intepolated divergence
533 interpolated_div_q(s, div_q);
534
535 // Return it
536 return div_q;
537 }
538
539 /// Calculate the FE representation of p
540 void interpolated_p(const Vector<double>& s, double& p) const
541 {
542 // Get the number of p basis functions
543 unsigned n_p_basis = np_basis();
544
545 // Storage for the p basis
546 Shape p_basis(n_p_basis);
547
548 // Call the p basis
549 get_p_basis(s, p_basis);
550
551 // Zero the pressure
552 p = 0;
553
554 // Add the contribution to the pressure from each basis function
555 for (unsigned l = 0; l < n_p_basis; l++)
556 {
557 p += p_value(l) * p_basis(l);
558 }
559 }
560
561 /// Calculate the FE representation of p and return it
562 double interpolated_p(const Vector<double>& s) const
563 {
564 // Temporary storage for p
565 double p = 0;
566
567 // Get the interpolated pressure
568 interpolated_p(s, p);
569
570 // Return it
571 return p;
572 }
573
574 /// du/dt at local node n
575 double du_dt(const unsigned& n, const unsigned& i) const
576 {
577 // Get the timestepper
579
580 // Storage for the derivative - initialise to 0
581 double du_dt = 0.0;
582
583 // If we are doing an unsteady solve then calculate the derivative
585 {
586 // Get the nodal index
587 const unsigned u_nodal_index = u_index(i);
588
589 // Get the number of values required to represent history
590 const unsigned n_time = time_stepper_pt->ntstorage();
591
592 // Loop over history values
593 for (unsigned t = 0; t < n_time; t++)
594 {
595 // Add the contribution to the derivative
596 du_dt +=
597 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
598 }
599 }
600
601 return du_dt;
602 }
603
604 /// d^2u/dt^2 at local node n
605 double d2u_dt2(const unsigned& n, const unsigned& i) const
606 {
607 // Get the timestepper
609
610 // Storage for the derivative - initialise to 0
611 double d2u_dt2 = 0.0;
612
613 // If we are doing an unsteady solve then calculate the derivative
615 {
616 // Get the nodal index
617 const unsigned u_nodal_index = u_index(i);
618
619 // Get the number of values required to represent history
620 const unsigned n_time = time_stepper_pt->ntstorage();
621
622 // Loop over history values
623 for (unsigned t = 0; t < n_time; t++)
624 {
625 // Add the contribution to the derivative
626 d2u_dt2 +=
627 time_stepper_pt->weight(2, t) * nodal_value(t, n, u_nodal_index);
628 }
629 }
630
631 return d2u_dt2;
632 }
633
634 /// dq_edge/dt for the n-th edge degree of freedom
635 double dq_edge_dt(const unsigned& n) const
636 {
637 unsigned node_num = q_edge_node_number(n);
638
639 // get the timestepper
641
642 // storage for the derivative - initialise to 0
643 double dq_dt = 0.0;
644
645 // if we are doing an unsteady solve then calculate the derivative
647 {
648 // get the number of values required to represent history
649 const unsigned n_time = time_stepper_pt->ntstorage();
650
651 // loop over history values
652 for (unsigned t = 0; t < n_time; t++)
653 {
654 // add the contribution to the derivative
655 dq_dt += time_stepper_pt->weight(1, t) * q_edge(t, n);
656 }
657 }
658
659 return dq_dt;
660 }
661
662 /// dq_internal/dt for the n-th internal degree of freedom
663 double dq_internal_dt(const unsigned& n) const
664 {
665 // get the internal data index for q
666 unsigned internal_index = q_internal_index();
667
668 // get the timestepper
670 internal_data_pt(internal_index)->time_stepper_pt();
671
672 // storage for the derivative - initialise to 0
673 double dq_dt = 0.0;
674
675 // if we are doing an unsteady solve then calculate the derivative
677 {
678 // get the number of values required to represent history
679 const unsigned n_time = time_stepper_pt->ntstorage();
680
681 // loop over history values
682 for (unsigned t = 0; t < n_time; t++)
683 {
684 // add the contribution to the derivative
685 dq_dt += time_stepper_pt->weight(1, t) * q_internal(t, n);
686 }
687 }
688
689 return dq_dt;
690 }
691
692 /// Set the timestepper of the q internal data object
694 {
695 unsigned q_index = q_internal_index();
696
697 this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt, false);
698 }
699
700 unsigned self_test()
701 {
702 return 0;
703 }
704
705 /// Output with default number of plot points
706 void output(std::ostream& outfile)
707 {
708 unsigned nplot = 5;
709 output(outfile, nplot);
710 }
711
712 /// Output FE representation of soln: x,y,u1,u2,div_q,p at
713 /// Nplot^DIM plot points
714 void output(std::ostream& outfile, const unsigned& nplot);
715
716 /// Output FE representation of exact soln: x,y,u1,u2,div_q,p at
717 /// Nplot^DIM plot points
718 void output_fct(std::ostream& outfile,
719 const unsigned& nplot,
721
722 /// Output FE representation of exact soln: x,y,u1,u2,div_q,p at
723 /// Nplot^DIM plot points. Unsteady version
724 void output_fct(std::ostream& outfile,
725 const unsigned& nplot,
726 const double& time,
728
729 /// Compute the error between the FE solution and the exact solution
730 /// using the H(div) norm for q and L^2 norm for p
731 void compute_error(std::ostream& outfile,
733 Vector<double>& error,
734 Vector<double>& norm);
735
736 /// Compute the error between the FE solution and the exact solution
737 /// using the H(div) norm for q and L^2 norm for p. Unsteady version
738 void compute_error(std::ostream& outfile,
740 const double& time,
741 Vector<double>& error,
742 Vector<double>& norm);
743
744 protected:
745 /// Returns the geometric basis, and the q, p and divergence basis
746 /// functions and test functions at local coordinate s
748 Shape& psi,
749 DShape& dpsi,
750 Shape& u_basis,
751 Shape& u_test,
752 DShape& du_basis_dx,
753 DShape& du_test_dx,
754 Shape& q_basis,
755 Shape& q_test,
756 Shape& p_basis,
757 Shape& p_test,
758 Shape& div_q_basis_ds,
759 Shape& div_q_test_ds) const = 0;
760
761 /// Returns the geometric basis, and the q, p and divergence basis
762 /// functions and test functions at integration point ipt
764 const unsigned& ipt,
765 Shape& psi,
766 DShape& dpsi,
767 Shape& u_basis,
768 Shape& u_test,
769 DShape& du_basis_dx,
770 DShape& du_test_dx,
771 Shape& q_basis,
772 Shape& q_test,
773 Shape& p_basis,
774 Shape& p_test,
775 Shape& div_q_basis_ds,
776 Shape& div_q_test_ds) const = 0;
777
778 // fill in residuals and, if flag==true, jacobian
780 Vector<double>& residuals, DenseMatrix<double>& jacobian, bool flag);
781
782 /// Pointer to the elasticity tensor
784
785 private:
786 /// Pointer to solid source function
788
789 /// Pointer to fluid source function
791
792 /// Pointer to the mass source function
794
795 /// Timescale ratio (non-dim. density)
797
798 /// Density ratio
800
801 /// 1/k
802 double* K_inv_pt;
803
804 /// Alpha
805 double* Alpha_pt;
806
807 /// Porosity
808 double* Porosity_pt;
809
810 /// Static default value for timescale ratio (1.0 -- for natural scaling)
812
813 /// Static default value for the density ratio
815
816 /// Static default value for 1/k
817 static double Default_k_inv_value;
818
819 /// Static default value for alpha
820 static double Default_alpha_value;
821
822 /// Static default value for the porosity
824 };
825
826} // namespace oomph
827
828#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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:406
A base class that represents the fourth-rank elasticity tensor defined such that.
A general Finite Element class.
Definition: elements.h:1313
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:2593
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...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1508
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
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:227
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
Class implementing the generic maths of the poroelasticity equations: linear elasticity coupled with ...
double dq_internal_dt(const unsigned &n) const
dq_internal/dt for the n-th internal degree of freedom
virtual int q_edge_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th edge (flux) degree of freedom.
SourceFctPt Force_fluid_fct_pt
Pointer to fluid source function.
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 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.
const double & density_ratio() const
Access function for the density ratio.
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(* SourceFctPt)(const double &time, const Vector< double > &x, Vector< double > &f)
Source function pointer typedef.
SourceFctPt force_fluid_fct_pt() const
Access function: Pointer to fluid force function (const version)
double interpolated_u(const Vector< double > &s, const unsigned &i) const
Calculate the FE representation of the i-th component of u.
virtual void pin_q_internal_value(const unsigned &n)=0
Pin the nth internal q value.
const double & k_inv() const
Access function for the nondim inverse permeability.
static double Default_porosity_value
Static default value for the porosity.
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.
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.
SourceFctPt Force_solid_fct_pt
Pointer to solid source function.
virtual unsigned q_edge_node_number(const unsigned &n) const =0
Return the number of the node where the nth edge unknown is stored.
const double & porosity() const
Access function for the porosity.
ElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
MassSourceFctPt & mass_source_fct_pt()
Access function: Pointer to mass source function.
virtual unsigned u_index(const unsigned &n) const =0
Return the nodal index of the n-th solid displacement unknown.
virtual int q_internal_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th internal (moment) degree of freedom.
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, DShape &dpsi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
static double Default_density_ratio_value
Static default value for the density ratio.
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...
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
virtual unsigned q_internal_index() const =0
Return the index of the internal data where the q_internal degrees of freedom are stored.
double dq_edge_dt(const unsigned &n) const
dq_edge/dt for the n-th edge degree of freedom
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals for the Darcy equations.
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.
double interpolated_div_q(const Vector< double > &s)
Calculate the FE representation of div q and return it.
virtual double q_edge(const unsigned &n) const =0
Return the values of the edge (flux) degrees of freedom.
void interpolated_q(const Vector< double > &s, Vector< double > &u) const
Calculate the FE representation of q.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the Jacobian matrix for the Newton method.
SourceFctPt & force_solid_fct_pt()
Access function: Pointer to solid force function.
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
const double & alpha() const
Access function for alpha.
void(* MassSourceFctPt)(const double &time, const Vector< double > &x, double &f)
Mass source function pointer typedef.
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
virtual double q_internal(const unsigned &n) const =0
Return the values of the internal (moment) degrees of freedom.
static double Default_alpha_value
Static default value for alpha.
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.
virtual double q_internal(const unsigned &t, const unsigned &n) const =0
Return the values of the internal (moment) degrees of freedom at time history level t.
virtual int p_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th pressure degree of freedom.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 ...
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
void set_q_internal_timestepper(TimeStepper *const time_stepper_pt)
Set the timestepper of the q internal data object.
virtual unsigned nedge_gauss_point() const =0
Returns the number of gauss points along each edge of the element.
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.
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Returns the transformed basis 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.
virtual unsigned required_nvalue(const unsigned &n) const =0
Number of values required at node n.
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
MassSourceFctPt mass_source_fct_pt() const
Access function: Pointer to mass source function (const version)
virtual double p_value(unsigned &n) const =0
Return the nth pressure value.
virtual unsigned q_edge_index(const unsigned &n) const =0
Return the nodal index at which the nth edge unknown is stored.
double interpolated_p(const Vector< double > &s) const
Calculate the FE representation of p and return it.
static double Default_k_inv_value
Static default value for 1/k.
double *& porosity_pt()
Access function for pointer to the porosity.
double interpolated_q(const Vector< double > &s, const unsigned i) const
Calculate the FE representation of the i-th component of q.
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div u.
void output(std::ostream &outfile)
Output with default number of plot points.
unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
double du_dt(const unsigned &n, const unsigned &i) const
du/dt at local node n
double *& k_inv_pt()
Access function for pointer to the nondim inverse permeability.
double *& density_ratio_pt()
Access function for pointer to the density ratio.
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.
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.
void interpolated_u(const Vector< double > &s, Vector< double > &disp) const
Calculate the FE representation of u.
SourceFctPt force_solid_fct_pt() const
Access function: Pointer to solid force function (const version)
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
virtual unsigned nq_basis() const =0
Return the total number of computational basis functions for q.
virtual double q_edge(const unsigned &t, const unsigned &n) const =0
Return the values of the edge (flux) degrees of freedom at time history level t.
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 ...
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
double *& alpha_pt()
Access function for pointer to alpha.
virtual double edge_gauss_point(const unsigned &edge, const unsigned &n) const =0
Returns the nth gauss point along an edge.
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Return the pressure basis.
SourceFctPt & force_fluid_fct_pt()
Access function: Pointer to fluid force function.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
//////////////////////////////////////////////////////////////////// ////////////////////////////////...