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