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