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 /// Static default value for timescale ratio (1.0 -- for natural scaling)
31 template<unsigned DIM>
33
34 /// Static default value for the density ratio
35 template<unsigned DIM>
37
38 /// Static default value for 1/k
39 template<unsigned DIM>
41
42 /// Static default value for alpha
43 template<unsigned DIM>
45
46 /// Static default value for the porosity
47 template<unsigned DIM>
49
50 //======================================================================
51 /// Compute the strain tensor at local coordinate s
52 //======================================================================
53 template<unsigned DIM>
55 const Vector<double>& s, DenseMatrix<double>& strain) const
56 {
57#ifdef PARANOID
58 if ((strain.ncol() != DIM) || (strain.nrow() != DIM))
59 {
60 std::ostringstream error_message;
61 error_message << "Strain matrix is " << strain.ncol() << " x "
62 << strain.nrow() << ", but dimension of the equations is "
63 << DIM << std::endl;
64 throw OomphLibError(
65 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
66 }
67
68 /*
69 //Find out how many position types there are
70 unsigned n_position_type = this->nnodal_position_type();
71
72 if(n_position_type != 1)
73 {
74 throw OomphLibError(
75 "PoroElasticity is not yet implemented for more than one position type",
76 OOMPH_CURRENT_FUNCTION,
77 OOMPH_EXCEPTION_LOCATION);
78 }
79 */
80#endif
81
82
83 // Find out how many nodes there are in the element
84 unsigned n_node = nnode();
85
86 // Find the indices at which the local velocities are stored
87 unsigned u_nodal_index[DIM];
88 for (unsigned i = 0; i < DIM; i++)
89 {
90 u_nodal_index[i] = u_index(i);
91 }
92
93 // Set up memory for the shape and derivative functions
94 Shape psi(n_node);
95 DShape dpsidx(n_node, DIM);
96
97 // Call the derivatives of the shape functions
98 (void)dshape_eulerian(s, psi, dpsidx);
99
100 // Calculate interpolated values of the derivative of global position
101 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
102
103 // Loop over nodes
104 for (unsigned l = 0; l < n_node; l++)
105 {
106 // Loop over velocity components
107 for (unsigned i = 0; i < DIM; i++)
108 {
109 // Get the nodal value
110 const double u_value = this->nodal_value(l, u_nodal_index[i]);
111
112 // Loop over derivative directions
113 for (unsigned j = 0; j < DIM; j++)
114 {
115 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
116 }
117 }
118 }
119
120 /// Now fill in the entries of the strain tensor
121 for (unsigned i = 0; i < DIM; i++)
122 {
123 // Do upper half of matrix
124 // Note that j must be signed here for the comparison test to work
125 // Also i must be cast to an int
126 for (int j = (DIM - 1); j >= static_cast<int>(i); j--)
127 {
128 // Off diagonal terms
129 if (static_cast<int>(i) != j)
130 {
131 strain(i, j) =
132 0.5 * (interpolated_dudx(i, j) + interpolated_dudx(j, i));
133 }
134 // Diagonal terms will including growth factor when it comes back in
135 else
136 {
137 strain(i, i) = interpolated_dudx(i, i);
138 }
139 }
140 // Matrix is symmetric so just copy lower half
141 for (int j = (i - 1); j >= 0; j--)
142 {
143 strain(i, j) = strain(j, i);
144 }
145 }
146 }
147
148
149 //======================================================================
150 /// Compute the Cauchy stress tensor at local coordinate s for
151 /// displacement formulation.
152 //======================================================================
153 template<unsigned DIM>
155 const Vector<double>& s, DenseMatrix<double>& stress) const
156 {
157#ifdef PARANOID
158 if ((stress.ncol() != DIM) || (stress.nrow() != DIM))
159 {
160 std::ostringstream error_message;
161 error_message << "Stress matrix is " << stress.ncol() << " x "
162 << stress.nrow() << ", but dimension of the equations is "
163 << DIM << std::endl;
164 throw OomphLibError(
165 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
166 }
167#endif
168
169 // Get strain
170 DenseMatrix<double> strain(DIM, DIM);
171 this->get_strain(s, strain);
172
173 // Now fill in the entries of the stress tensor without exploiting
174 // symmetry -- sorry too lazy. This fct is only used for
175 // postprocessing anyway...
176 for (unsigned i = 0; i < DIM; i++)
177 {
178 for (unsigned j = 0; j < DIM; j++)
179 {
180 stress(i, j) = 0.0;
181 for (unsigned k = 0; k < DIM; k++)
182 {
183 for (unsigned l = 0; k < DIM; k++)
184 {
185 stress(i, j) += this->E(i, j, k, l) * strain(k, l);
186 }
187 }
188 }
189 }
190 }
191
192
193 /// Performs a div-conserving transformation of the vector basis
194 /// functions from the reference element to the actual element
195 template<unsigned DIM>
197 const Vector<double>& s,
198 const Shape& q_basis_local,
199 Shape& psi,
200 DShape& dpsi,
201 Shape& q_basis) const
202 {
203 // Call the (geometric) shape functions and their derivatives
204 //(void)this->dshape_eulerian(s,psi,dpsi);
205 this->dshape_local(s, psi, dpsi);
206
207 // Storage for the (geometric) jacobian and its inverse
208 DenseMatrix<double> jacobian(DIM), inverse_jacobian(DIM);
209
210 // Get the jacobian of the geometric mapping and its determinant
211 double det = local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
212
213 // Transform the derivative of the geometric basis so that it's w.r.t.
214 // global coordinates
215 this->transform_derivatives(inverse_jacobian, dpsi);
216
217 // Get the number of computational basis vectors
218 const unsigned n_q_basis = this->nq_basis();
219
220 // Loop over the basis vectors
221 for (unsigned l = 0; l < n_q_basis; l++)
222 {
223 // Loop over the spatial components
224 for (unsigned i = 0; i < DIM; i++)
225 {
226 // Zero the basis
227 q_basis(l, i) = 0.0;
228 }
229 }
230
231 // Loop over the spatial components
232 for (unsigned i = 0; i < DIM; i++)
233 {
234 // And again
235 for (unsigned j = 0; j < DIM; j++)
236 {
237 // Get the element of the jacobian (must transpose it due to different
238 // conventions) and divide by the determinant
239 double jac_trans = jacobian(j, i) / det;
240
241 // Loop over the computational basis vectors
242 for (unsigned l = 0; l < n_q_basis; l++)
243 {
244 // Apply the appropriate div-conserving mapping
245 q_basis(l, i) += jac_trans * q_basis_local(l, j);
247 }
248 }
250 // Scale the basis by the ratio of the length of the edge to the length of
251 // the corresponding edge on the reference element
252 scale_basis(q_basis);
253
254 return det;
255 }
256
257 /// Output FE representation of soln: x,y,u1,u2,q1,q2,div_q,p at
258 /// Nplot^DIM plot points
259 template<unsigned DIM>
260 void PoroelasticityEquations<DIM>::output(std::ostream& outfile,
261 const unsigned& nplot)
262 {
263 // Vector of local coordinates
264 Vector<double> s(DIM);
265
266 // Tecplot header info
267 outfile << tecplot_zone_string(nplot);
268
269 // Loop over plot points
270 unsigned num_plot_points = nplot_points(nplot);
271
272 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
273 {
274 // Get local coordinates of plot point
275 get_s_plot(iplot, nplot, s);
276
277 // Output the components of the position
278 for (unsigned i = 0; i < DIM; i++)
279 {
280 outfile << interpolated_x(s, i) << " ";
281 }
282
283 // Output the components of the FE representation of u at s
284 for (unsigned i = 0; i < DIM; i++)
285 {
286 outfile << interpolated_u(s, i) << " ";
287 }
288
289 // Output the components of the FE representation of q at s
290 for (unsigned i = 0; i < DIM; i++)
291 {
292 outfile << interpolated_q(s, i) << " ";
293 }
294
295 // Output FE representation of div u at s
296 outfile << interpolated_div_q(s) << " ";
297
298 // Output FE representation of p at s
299 outfile << interpolated_p(s) << " ";
300
301 const unsigned n_node = this->nnode();
302
303 Shape psi(n_node);
304 shape(s, psi);
305
306 Vector<double> interpolated_du_dt(DIM, 0.0);
307
308 for (unsigned i = 0; i < DIM; i++)
309 {
310 for (unsigned l = 0; l < n_node; l++)
311 {
312 interpolated_du_dt[i] += du_dt(l, i) * psi(l);
313 }
314 outfile << interpolated_du_dt[i] << " ";
315 }
316
317 Vector<double> interpolated_d2u_dt2(DIM, 0.0);
318
319 for (unsigned i = 0; i < DIM; i++)
320 {
321 for (unsigned l = 0; l < n_node; l++)
322 {
323 interpolated_d2u_dt2[i] += d2u_dt2(l, i) * psi(l);
324 }
325 outfile << interpolated_d2u_dt2[i] << " ";
326 }
327
328 const unsigned n_q_basis = this->nq_basis();
329 const unsigned n_q_basis_edge = this->nq_basis_edge();
330
331 Shape q_basis(n_q_basis, DIM), q_basis_local(n_q_basis, DIM);
332 this->get_q_basis_local(s, q_basis_local);
333 (void)this->transform_basis(s, q_basis_local, psi, q_basis);
334
335 Vector<double> interpolated_dq_dt(DIM, 0.0);
336
337 for (unsigned i = 0; i < DIM; i++)
338 {
339 for (unsigned l = 0; l < n_q_basis_edge; l++)
340 {
341 interpolated_dq_dt[i] += dq_edge_dt(l) * q_basis(l, i);
342 }
343 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
344 {
345 interpolated_dq_dt[i] +=
346 dq_internal_dt(l - n_q_basis_edge) * q_basis(l, i);
347 }
348 outfile << interpolated_dq_dt[i] << " ";
350
351 outfile << std::endl;
352 }
353
354 // Write tecplot footer (e.g. FE connectivity lists)
355 this->write_tecplot_zone_footer(outfile, nplot);
356 }
357
358 /// Output FE representation of exact soln: x,y,u1,u2,q1,q2,div_q,p at
359 /// Nplot^DIM plot points
360 template<unsigned DIM>
362 std::ostream& outfile,
363 const unsigned& nplot,
365 {
366 // Vector of local coordinates
367 Vector<double> s(DIM);
368
369 // Vector for coordintes
370 Vector<double> x(DIM);
371
372 // Tecplot header info
373 outfile << this->tecplot_zone_string(nplot);
374
375 // Exact solution Vector
376 Vector<double> exact_soln(2 * DIM + 2);
377
378 // Loop over plot points
379 unsigned num_plot_points = this->nplot_points(nplot);
380
381 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
382 {
383 // Get local coordinates of plot point
384 this->get_s_plot(iplot, nplot, s);
385
386 // Get x position as Vector
387 this->interpolated_x(s, x);
388
389 // Get exact solution at this point
390 (*exact_soln_pt)(x, exact_soln);
391
392 // Output x,y,q_exact,p_exact,div_q_exact
393 for (unsigned i = 0; i < DIM; i++)
394 {
395 outfile << x[i] << " ";
396 }
397 for (unsigned i = 0; i < 2 * DIM + 2; i++)
398 {
399 outfile << exact_soln[i] << " ";
400 }
401 outfile << std::endl;
402 }
403
404 // Write tecplot footer (e.g. FE connectivity lists)
405 this->write_tecplot_zone_footer(outfile, nplot);
406 }
407
408 /// Output FE representation of exact soln: x,y,u1,u2,div_q,p at
409 /// Nplot^DIM plot points. Unsteady version
410 template<unsigned DIM>
412 std::ostream& outfile,
413 const unsigned& nplot,
414 const double& time,
416 {
417 // Vector of local coordinates
418 Vector<double> s(DIM);
419
420 // Vector for coordintes
421 Vector<double> x(DIM);
422
423 // Tecplot header info
424 outfile << this->tecplot_zone_string(nplot);
425
426 // Exact solution Vector
427 Vector<double> exact_soln(2 * DIM + 2);
428
429 // Loop over plot points
430 unsigned num_plot_points = this->nplot_points(nplot);
431
432 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
433 {
434 // Get local coordinates of plot point
435 this->get_s_plot(iplot, nplot, s);
436
437 // Get x position as Vector
438 this->interpolated_x(s, x);
439
440 // Get exact solution at this point
441 (*exact_soln_pt)(time, x, exact_soln);
442
443 // Output x,y,q_exact,p_exact,div_q_exact
444 for (unsigned i = 0; i < DIM; i++)
445 {
446 outfile << x[i] << " ";
447 }
448 for (unsigned i = 0; i < 2 * DIM + 2; i++)
449 {
450 outfile << exact_soln[i] << " ";
451 }
452 outfile << std::endl;
453 }
454
455 // Write tecplot footer (e.g. FE connectivity lists)
456 this->write_tecplot_zone_footer(outfile, nplot);
457 }
458
459 /// Compute the error between the FE solution and the exact solution
460 /// using the H(div) norm for q and L^2 norm for u and p
461 template<unsigned DIM>
463 std::ostream& outfile,
465 Vector<double>& error,
466 Vector<double>& norm)
467 {
468 for (unsigned i = 0; i < 3; i++)
469 {
470 error[i] = 0.0;
471 norm[i] = 0.0;
472 }
473
474 // Vector of local coordinates
475 Vector<double> s(DIM);
476
477 // Vector for coordinates
478 Vector<double> x(DIM);
479
480 // Set the value of n_intpt
481 unsigned n_intpt = this->integral_pt()->nweight();
482
483 outfile << "ZONE" << std::endl;
484
485 // Exact solution Vector (u,v,[w])
486 Vector<double> exact_soln(2 * DIM + 2);
487
488 // Loop over the integration points
489 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
490 {
491 // Assign values of s
492 for (unsigned i = 0; i < DIM; i++)
493 {
494 s[i] = this->integral_pt()->knot(ipt, i);
495 }
496
497 // Get the integral weight
498 double w = this->integral_pt()->weight(ipt);
499
500 // Get jacobian of mapping
501 double J = this->J_eulerian(s);
502
503 // Premultiply the weights and the Jacobian
504 double W = w * J;
505
506 // Get x position as Vector
507 this->interpolated_x(s, x);
508
509 // Get exact solution at this point
510 (*exact_soln_pt)(x, exact_soln);
511
512 // Displacement error
513 for (unsigned i = 0; i < DIM; i++)
514 {
515 norm[0] += exact_soln[i] * exact_soln[i] * W;
516 // Error due to q_i
517 error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
518 (exact_soln[i] - this->interpolated_u(s, i)) * W;
519 }
520
521 // Flux error
522 for (unsigned i = 0; i < DIM; i++)
523 {
524 norm[1] += exact_soln[DIM + i] * exact_soln[DIM + i] * W;
525 // Error due to q_i
526 error[1] += (exact_soln[DIM + i] - this->interpolated_q(s, i)) *
527 (exact_soln[DIM + i] - this->interpolated_q(s, i)) * W;
528 }
529
530 // Flux divergence error
531 norm[1] += exact_soln[2 * DIM] * exact_soln[2 * DIM] * W;
532 error[1] += (exact_soln[2 * DIM] - interpolated_div_q(s)) *
533 (exact_soln[2 * DIM] - interpolated_div_q(s)) * W;
534
535 // Pressure error
536 norm[2] += exact_soln[2 * DIM + 1] * exact_soln[2 * DIM + 1] * W;
537 error[2] += (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) *
538 (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) * W;
539
540 // Output x,y,[z]
541 for (unsigned i = 0; i < DIM; i++)
542 {
543 outfile << x[i] << " ";
544 }
545
546 // Output u_1_error,u_2_error,...
547 for (unsigned i = 0; i < DIM; i++)
548 {
549 outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
550 }
551
552 // Output q_1_error,q_2_error,...
553 for (unsigned i = 0; i < DIM; i++)
554 {
555 outfile << exact_soln[DIM + i] - this->interpolated_q(s, i) << " ";
556 }
557
558 // Output p_error
559 outfile << exact_soln[2 * DIM + 1] - this->interpolated_p(s) << " ";
560
561 outfile << std::endl;
562 }
563 }
564
565 /// Compute the error between the FE solution and the exact solution
566 /// using the H(div) norm for u and L^2 norm for p. Unsteady version
567 template<unsigned DIM>
569 std::ostream& outfile,
571 const double& time,
572 Vector<double>& error,
573 Vector<double>& norm)
574 {
575 for (unsigned i = 0; i < 3; i++)
576 {
577 error[i] = 0.0;
578 norm[i] = 0.0;
579 }
580
581 // Vector of local coordinates
582 Vector<double> s(DIM);
583
584 // Vector for coordinates
585 Vector<double> x(DIM);
586
587 // Set the value of n_intpt
588 unsigned n_intpt = this->integral_pt()->nweight();
589
590 outfile << "ZONE" << std::endl;
591
592 // Exact solution Vector (u,v,[w])
593 Vector<double> exact_soln(2 * DIM + 2);
594
595 // Loop over the integration points
596 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
597 {
598 // Assign values of s
599 for (unsigned i = 0; i < DIM; i++)
600 {
601 s[i] = this->integral_pt()->knot(ipt, i);
602 }
603
604 // Get the integral weight
605 double w = this->integral_pt()->weight(ipt);
606
607 // Get jacobian of mapping
608 double J = this->J_eulerian(s);
609
610 // Premultiply the weights and the Jacobian
611 double W = w * J;
612
613 // Get x position as Vector
614 this->interpolated_x(s, x);
615
616 // Get exact solution at this point
617 (*exact_soln_pt)(time, x, exact_soln);
618
619 // Displacement error
620 for (unsigned i = 0; i < DIM; i++)
621 {
622 norm[0] += exact_soln[i] * exact_soln[i] * W;
623 // Error due to q_i
624 error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
625 (exact_soln[i] - this->interpolated_u(s, i)) * W;
626 }
627
628 // Flux error
629 for (unsigned i = 0; i < DIM; i++)
630 {
631 norm[1] += exact_soln[DIM + i] * exact_soln[DIM + i] * W;
632 // Error due to q_i
633 error[1] += (exact_soln[DIM + i] - this->interpolated_q(s, i)) *
634 (exact_soln[DIM + i] - this->interpolated_q(s, i)) * W;
635 }
636
637 // Flux divergence error
638 norm[1] += exact_soln[2 * DIM] * exact_soln[2 * DIM] * W;
639 error[1] += (exact_soln[2 * DIM] - interpolated_div_q(s)) *
640 (exact_soln[2 * DIM] - interpolated_div_q(s)) * W;
641
642 // Pressure error
643 norm[2] += exact_soln[2 * DIM + 1] * exact_soln[2 * DIM + 1] * W;
644 error[2] += (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) *
645 (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) * W;
646
647 // Output x,y,[z]
648 for (unsigned i = 0; i < DIM; i++)
649 {
650 outfile << x[i] << " ";
651 }
652
653 // Output u_1_error,u_2_error,...
654 for (unsigned i = 0; i < DIM; i++)
655 {
656 outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
657 }
658
659 // Output q_1_error,q_2_error,...
660 for (unsigned i = 0; i < DIM; i++)
661 {
662 outfile << exact_soln[DIM + i] - this->interpolated_q(s, i) << " ";
663 }
664
665 // Output p_error
666 outfile << exact_soln[2 * DIM + 1] - this->interpolated_p(s) << " ";
667
668 outfile << std::endl;
669 }
670 }
671
672 /// Fill in residuals and, if flag==true, jacobian
673 template<unsigned DIM>
675 Vector<double>& residuals, DenseMatrix<double>& jacobian, bool flag)
676 {
677 // Get the number of geometric nodes, total number of basis functions,
678 // and number of edges basis functions
679 const unsigned n_node = nnode();
680 const unsigned n_q_basis = nq_basis();
681 const unsigned n_q_basis_edge = nq_basis_edge();
682 const unsigned n_p_basis = np_basis();
683
684 // Storage for the geometric and computational bases and the test functions
685 Shape psi(n_node), u_basis(n_node), u_test(n_node), q_basis(n_q_basis, DIM),
686 q_test(n_q_basis, DIM), p_basis(n_p_basis), p_test(n_p_basis),
687 div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
688
689 DShape dpsidx(n_node, DIM), du_basis_dx(n_node, DIM),
690 du_test_dx(n_node, DIM);
691
692 // Get the number of integration points
693 unsigned n_intpt = integral_pt()->nweight();
694
695 // Storage for the local coordinates
696 Vector<double> s(DIM);
697
698 // Storage for the elasticity source function
699 Vector<double> f_solid(DIM);
700
701 // Storage for the source function
702 Vector<double> f_fluid(DIM);
703
704 // Storage for the mass source function
705 double mass_source_local;
706
707 // Storage for Lambda_sq
708 double lambda_sq = this->lambda_sq();
709
710 // Get the value of 1/k
711 double k_inv = this->k_inv();
712
713 // Get the value of alpha
714 double alpha = this->alpha();
715
716 // Get the value of the porosity
717 double porosity = this->porosity();
719 // Get the density ratio
720 double density_ratio = this->density_ratio();
721
722 // Precompute the ratio of fluid density to combined density
723 double rho_f_over_rho =
724 density_ratio / (1 + porosity * (density_ratio - 1));
725
726 // Get continuous time from timestepper of first node
727 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
728
729 // Local equation and unknown numbers
730 int local_eqn = 0, local_unknown = 0;
732 // Loop over the integration points
733 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
734 {
735 // Find the local coordinates at the integration point
736 for (unsigned i = 0; i < DIM; i++)
737 {
738 s[i] = integral_pt()->knot(ipt, i);
739 }
740
741 // Get the weight of the intetgration point
742 double w = integral_pt()->weight(ipt);
743
744 // Call the basis functions and test functions and get the
745 // (geometric) jacobian of the current element
746 double J = shape_basis_test_local_at_knot(ipt,
747 psi,
748 dpsidx,
749 u_basis,
750 u_test,
751 du_basis_dx,
752 du_test_dx,
753 q_basis,
754 q_test,
755 p_basis,
756 p_test,
757 div_q_basis_ds,
758 div_q_test_ds);
759
760 // Storage for interpolated values
761 Vector<double> interpolated_x(DIM, 0.0);
762 Vector<double> interpolated_u(DIM, 0.0);
763 DenseMatrix<double> interpolated_du_dx(DIM, DIM, 0.0);
764 double interpolated_div_du_dt_dx = 0.0;
765 Vector<double> interpolated_d2u_dt2(DIM, 0.0);
766 Vector<double> interpolated_q(DIM, 0.0);
767 double interpolated_div_q_ds = 0.0;
768 Vector<double> interpolated_dq_dt(DIM, 0.0);
769 double interpolated_p = 0.0;
770
771 // loop over the vector basis functions to find interpolated x
772 for (unsigned l = 0; l < n_node; l++)
773 {
774 // Loop over the geometric basis functions
775 for (unsigned i = 0; i < DIM; i++)
776 {
777 interpolated_x[i] += nodal_position(l, i) * psi(l);
778
779 interpolated_d2u_dt2[i] += this->d2u_dt2(l, i) * u_basis(l);
780
781 // Get the nodal displacements
782 const double u_value = this->raw_nodal_value(l, u_index(i));
783 interpolated_u[i] += u_value * u_basis(l);
784
785 // Loop over derivative directions
786 for (unsigned j = 0; j < DIM; j++)
787 {
788 interpolated_du_dx(i, j) += u_value * du_basis_dx(l, j);
789 }
790
791 // divergence of the time derivative of the solid displacement
792 interpolated_div_du_dt_dx += this->du_dt(l, i) * du_basis_dx(l, i);
793 }
794 }
795
796 // loop over the nodes and use the geometric basis functions to find the
797 // interpolated flux and its time derivative
798 for (unsigned i = 0; i < DIM; i++)
799 {
800 // Loop over the edge basis vectors
801 for (unsigned l = 0; l < n_q_basis_edge; l++)
802 {
803 interpolated_q[i] += q_edge(l) * q_basis(l, i);
804 interpolated_dq_dt[i] += dq_edge_dt(l) * q_basis(l, i);
805 }
806 // Loop over the internal basis vectors
807 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
808 {
809 interpolated_q[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
810 interpolated_dq_dt[i] +=
811 dq_internal_dt(l - n_q_basis_edge) * q_basis(l, i);
812 }
813 }
814
815 // loop over the pressure basis and find the interpolated pressure
816 for (unsigned l = 0; l < n_p_basis; l++)
817 {
818 interpolated_p += p_value(l) * p_basis(l);
819 }
820
821 // loop over the u edge divergence basis and the u internal divergence
822 // basis to find interpolated div u
823 for (unsigned l = 0; l < n_q_basis_edge; l++)
824 {
825 interpolated_div_q_ds += q_edge(l) * div_q_basis_ds(l);
826 }
827 for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
828 {
829 interpolated_div_q_ds +=
830 q_internal(l - n_q_basis_edge) * div_q_basis_ds(l);
831 }
832
833 // Get the solid forcing
834 this->force_solid(time, interpolated_x, f_solid);
835
836 // Get the fluid forcing
837 this->force_fluid(time, interpolated_x, f_fluid);
838
839 // Get the mass source function
840 this->mass_source(time, interpolated_x, mass_source_local);
841
842 // Premultiply the weights and the Jacobian
843 double W = w * J;
844
845 // Linear elasticity:
846
847 for (unsigned l = 0; l < n_node; l++)
848 {
849 for (unsigned a = 0; a < DIM; a++)
850 {
851 local_eqn = this->nodal_local_eqn(l, u_index(a));
852
853 if (local_eqn >= 0)
854 {
855 residuals[local_eqn] +=
856 (lambda_sq * (interpolated_d2u_dt2[a] +
857 rho_f_over_rho * interpolated_dq_dt[a]) -
858 f_solid[a]) *
859 u_test(l) * W;
860
861 // Stress term
862 for (unsigned b = 0; b < DIM; b++)
863 {
864 if (a == b)
865 {
866 residuals[local_eqn] -=
867 alpha * interpolated_p * du_test_dx(l, b) * W;
868 }
869
870 for (unsigned c = 0; c < DIM; c++)
871 {
872 for (unsigned d = 0; d < DIM; d++)
873 {
874 // Add the stress terms to the residuals
875 residuals[local_eqn] += this->E(a, b, c, d) *
876 interpolated_du_dx(c, d) *
877 du_test_dx(l, b) * W;
878 }
879 }
880 }
881 // Jacobian entries
882 if (flag)
883 {
884 // d(u_eqn_l,a)/d(U_l2,c)
885 for (unsigned l2 = 0; l2 < n_node; l2++)
886 {
887 for (unsigned c = 0; c < DIM; c++)
888 {
889 local_unknown = this->nodal_local_eqn(l2, u_index(c));
890 if (local_unknown >= 0)
891 {
892 if (a == c)
893 {
894 jacobian(local_eqn, local_unknown) +=
895 lambda_sq *
896 this->node_pt(l2)->time_stepper_pt()->weight(2, 0) *
897 u_basis(l2) * u_test(l) * W;
898 }
899
900 for (unsigned b = 0; b < DIM; b++)
901 {
902 for (unsigned d = 0; d < DIM; d++)
903 {
904 jacobian(local_eqn, local_unknown) +=
905 this->E(a, b, c, d) * du_basis_dx(l2, d) *
906 du_test_dx(l, b) * W;
907 }
908 }
909 }
910 }
911 }
912
913 // d(u_eqn_l,a)/d(Q_l2)
914 for (unsigned l2 = 0; l2 < n_q_basis; l2++)
915 {
916 TimeStepper* timestepper_pt = 0;
917
918 if (l2 < n_q_basis_edge)
919 {
920 local_unknown = q_edge_local_eqn(l2);
921 timestepper_pt =
922 this->node_pt(q_edge_node_number(l2))->time_stepper_pt();
923 }
924 else // n_q_basis_edge <= l < n_basis
925 {
926 local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
927 timestepper_pt = this->internal_data_pt(q_internal_index())
928 ->time_stepper_pt();
929 }
930
931 if (local_unknown >= 0)
932 {
933 jacobian(local_eqn, local_unknown) +=
934 lambda_sq * rho_f_over_rho * timestepper_pt->weight(1, 0) *
935 q_basis(l2, a) * u_test(l) * W;
936 }
937 }
938
939 // d(u_eqn_l,a)/d(P_l2)
940 for (unsigned l2 = 0; l2 < n_p_basis; l2++)
941 {
942 local_unknown = p_local_eqn(l2);
943 if (local_unknown >= 0)
944 {
945 jacobian(local_eqn, local_unknown) -=
946 alpha * p_basis(l2) * du_test_dx(l, a) * W;
947 }
948 }
949 } // End of Jacobian entries
950 } // End of if not boundary condition
951 } // End of loop over dimensions
952 } // End of loop over u test functions
953
954
955 // Darcy:
956 // Loop over the test functions
957 for (unsigned l = 0; l < n_q_basis; l++)
958 {
959 if (l < n_q_basis_edge)
960 {
961 local_eqn = q_edge_local_eqn(l);
962 }
963 else // n_q_basis_edge <= l < n_basis
964 {
965 local_eqn = q_internal_local_eqn(l - n_q_basis_edge);
966 }
967
968 // If it's not a boundary condition
969 if (local_eqn >= 0)
970 {
971 for (unsigned i = 0; i < DIM; i++)
972 {
973 residuals[local_eqn] +=
974 (k_inv * interpolated_q[i] - rho_f_over_rho * f_fluid[i] +
975 rho_f_over_rho * lambda_sq *
976 (interpolated_d2u_dt2[i] +
977 (1.0 / porosity) * interpolated_dq_dt[i])) *
978 q_test(l, i) * w * J;
979 }
980
981 // deliberately no jacobian factor in this integral
982 residuals[local_eqn] -= (interpolated_p * div_q_test_ds(l)) * w;
983
984 // Jacobian entries
985 if (flag)
986 {
987 // d(q_eqn_l)/d(U_l2,c)
988 for (unsigned l2 = 0; l2 < n_node; l2++)
989 {
990 for (unsigned c = 0; c < DIM; c++)
991 {
992 local_unknown = this->nodal_local_eqn(l2, u_index(c));
993 if (local_unknown >= 0)
994 {
995 jacobian(local_eqn, local_unknown) +=
996 rho_f_over_rho * lambda_sq *
997 this->node_pt(l2)->time_stepper_pt()->weight(2, 0) *
998 u_basis(l2) * q_test(l, c) * W;
999 }
1000 }
1001 }
1002
1003 // d(q_eqn_l)/d(Q_l2)
1004 for (unsigned l2 = 0; l2 < n_q_basis; l2++)
1005 {
1006 TimeStepper* timestepper_pt = 0;
1007
1008 if (l2 < n_q_basis_edge)
1009 {
1010 local_unknown = q_edge_local_eqn(l2);
1011 timestepper_pt =
1012 this->node_pt(q_edge_node_number(l2))->time_stepper_pt();
1013 }
1014 else // n_q_basis_edge <= l < n_basis
1015 {
1016 local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
1017 timestepper_pt =
1018 this->internal_data_pt(q_internal_index())->time_stepper_pt();
1019 }
1020
1021 if (local_unknown >= 0)
1022 {
1023 for (unsigned c = 0; c < DIM; c++)
1024 {
1025 jacobian(local_eqn, local_unknown) +=
1026 q_basis(l2, c) * q_test(l, c) *
1027 (k_inv + rho_f_over_rho * lambda_sq *
1028 timestepper_pt->weight(1, 0) / porosity) *
1029 W;
1030 }
1031 }
1032 }
1033
1034 // d(q_eqn_l)/d(P_l2)
1035 for (unsigned l2 = 0; l2 < n_p_basis; l2++)
1036 {
1037 local_unknown = p_local_eqn(l2);
1038
1039 // deliberately no jacobian factor in this term
1040 jacobian(local_eqn, local_unknown) -=
1041 p_basis(l2) * div_q_test_ds(l) * w;
1042 }
1043 } // End of Jacobian entries
1044 } // End of if not boundary condition
1045 } // End of loop over q test functions
1046
1047 // loop over pressure test functions
1048 for (unsigned l = 0; l < n_p_basis; l++)
1049 {
1050 // get the local equation number
1051 local_eqn = p_local_eqn(l);
1052
1053 // If it's not a boundary condition
1054 if (local_eqn >= 0)
1055 {
1056 // solid divergence term
1057 residuals[local_eqn] +=
1058 alpha * interpolated_div_du_dt_dx * p_test(l) * w * J;
1059 // fluid divergence term - deliberately no jacobian factor in this
1060 // term
1061 residuals[local_eqn] += interpolated_div_q_ds * p_test(l) * w;
1062 // mass source term
1063 residuals[local_eqn] -= mass_source_local * p_test(l) * w * J;
1064
1065 // Jacobian entries
1066 if (flag)
1067 {
1068 // d(p_eqn_l)/d(U_l2,c)
1069 for (unsigned l2 = 0; l2 < n_node; l2++)
1070 {
1071 for (unsigned c = 0; c < DIM; c++)
1072 {
1073 local_unknown = this->nodal_local_eqn(l2, u_index(c));
1074
1075 if (local_unknown >= 0)
1076 {
1077 jacobian(local_eqn, local_unknown) +=
1078 alpha * this->node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1079 du_basis_dx(l2, c) * p_test(l) * W;
1080 }
1081 }
1082 }
1083
1084 // d(p_eqn_l)/d(Q_l2)
1085 for (unsigned l2 = 0; l2 < n_q_basis; l2++)
1086 {
1087 if (l2 < n_q_basis_edge)
1088 {
1089 local_unknown = q_edge_local_eqn(l2);
1090 }
1091 else // n_q_basis_edge <= l < n_basis
1092 {
1093 local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
1094 }
1095
1096 if (local_unknown >= 0)
1097 {
1098 jacobian(local_eqn, local_unknown) +=
1099 p_test(l) * div_q_basis_ds(l2) * w;
1100 }
1101 }
1102 } // End of Jacobian entries
1103 } // End of if not boundary condition
1104 } // End of loop over p test functions
1105 } // End of loop over integration points
1106 }
1107
1108 // Force building of templates
1109 template class PoroelasticityEquations<2>;
1110
1111} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
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
An OomphLibError object which should be thrown when an run-time error is encountered....
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output FE representation of exact soln: x,y,u1,u2,div_q,p at Nplot^DIM plot points.
static double Default_porosity_value
Static default value for the porosity.
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 ...
static double Default_density_ratio_value
Static default value for the density ratio.
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
static double Default_alpha_value
Static default value for alpha.
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 ...
static double Default_k_inv_value
Static default value for 1/k.
void output(std::ostream &outfile)
Output with default number of plot points.
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:564
//////////////////////////////////////////////////////////////////// ////////////////////////////////...