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