darcy_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: 1099 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-07 09:05:59 +0000 (Thu, 07 Jan 2016) $
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//====================================================================
30 #include "darcy_elements.h"
31 
32 namespace oomph
33 {
34 
35  //=====================================================================
36  /// \short Performs a div-conserving transformation of the vector basis
37  /// functions from the reference element to the actual element
38  //=====================================================================
39  template<unsigned DIM>
41  const Shape &q_basis_local,
42  Shape &psi,
43  Shape &q_basis) const
44  {
45  // Get the number of nodes in the element
46  const unsigned n_node = this->nnode();
47 
48  // Storage for derivatives of the (geometric) shape functions, and call
49  // the shape functions
50  DShape dpsi(n_node,DIM);
51  this->dshape_local(s,psi,dpsi);
52 
53  // Storage for the (geometric) jacobian and its inverse
54  DenseMatrix<double> jacobian(DIM), inverse_jacobian(DIM);
55 
56  // Get the jacobian of the geometric mapping and its determinant
57  double det = local_to_eulerian_mapping(dpsi,jacobian,inverse_jacobian);
58 
59  // Get the number of computational basis vectors
60  const unsigned n_q_basis = this->nq_basis();
61 
62  // Loop over the basis vectors
63  for(unsigned l=0;l<n_q_basis;l++)
64  {
65  // Loop over the spatial components
66  for(unsigned i=0;i<DIM;i++)
67  {
68  // Zero the basis
69  q_basis(l,i) = 0.0;
70  }
71  }
72 
73  // Loop over the spatial components
74  for(unsigned i=0;i<DIM;i++)
75  {
76  // And again
77  for(unsigned j=0;j<DIM;j++)
78  {
79  // Get the element of the jacobian (must transpose it due to different
80  // conventions) and divide by the determinant
81  double jac_trans = jacobian(j,i)/det;
82 
83  // Loop over the computational basis vectors
84  for(unsigned l=0;l<n_q_basis;l++)
85  {
86  // Apply the appropriate div-conserving mapping
87  q_basis(l,i) += jac_trans*q_basis_local(l,j);
88  }
89  }
90  }
91 
92  // Scale the basis by the ratio of the length of the edge to the length of
93  // the corresponding edge on the reference element
94  scale_basis(q_basis);
95 
96  return det;
97  }
98 
99  //=====================================================================
100  /// \short Output FE representation of soln: x,y,q1,q2,div_q,p,q \cdot n
101  //=====================================================================
102  template<unsigned DIM>
104  std::ostream &outfile,
105  const unsigned &nplot,
106  const Vector<double> unit_normal)
107  {
108 
109  // Vector of local coordinates
110  Vector<double> s(DIM);
111 
112  // Tecplot header info
113  outfile << this->tecplot_zone_string(nplot);
114 
115  // Loop over plot points
116  unsigned num_plot_points=this->nplot_points(nplot);
117 
118  for(unsigned iplot=0;iplot<num_plot_points;iplot++)
119  {
120  // Get local coordinates of plot point
121  this->get_s_plot(iplot,nplot,s);
122 
123  // Output the components of the position
124  for(unsigned i=0;i<DIM;i++)
125  {
126  outfile << this->interpolated_x(s,i) << " ";
127  }
128 
129  // Output the components of the FE representation of q at s
130  for(unsigned i=0;i<DIM;i++)
131  {
132  outfile << this->interpolated_q(s,i) << " ";
133  }
134 
135  // Output FE representation of div q at s
136  outfile << this->interpolated_div_q(s) << " ";
137 
138  // Output FE representation of p at s
139  outfile << this->interpolated_p(s) << " ";
140 
141  // Fluxes projected into the direction of the face normal
142  double flux=0.0;
143  for (unsigned i=0;i<2;i++)
144  {
145  flux+=this->interpolated_q(s,i)*unit_normal[i];
146  }
147  outfile << flux << " ";
148 
149  outfile << std::endl;
150  }
151 
152  // Write tecplot footer (e.g. FE connectivity lists)
153  this->write_tecplot_zone_footer(outfile,nplot);
154 
155  }
156 
157 
158  //=====================================================================
159  /// \short Output FE representation of soln: x,y,q1,q2,div_q,p at
160  /// Nplot^DIM plot points
161  //=====================================================================
162  template<unsigned DIM>
163  void DarcyEquations<DIM>::output(std::ostream &outfile, const unsigned &nplot)
164  {
165  // Vector of local coordinates
166  Vector<double> s(DIM);
167 
168  // Tecplot header info
169  outfile << tecplot_zone_string(nplot);
170 
171  // Loop over plot points
172  unsigned num_plot_points=nplot_points(nplot);
173 
174  for(unsigned iplot=0;iplot<num_plot_points;iplot++)
175  {
176  // Get local coordinates of plot point
177  get_s_plot(iplot,nplot,s);
178 
179  // Output the components of the position
180  for(unsigned i=0;i<DIM;i++)
181  {
182  outfile << interpolated_x(s,i) << " ";
183  }
184 
185  // Output the components of the FE representation of q at s
186  for(unsigned i=0;i<DIM;i++)
187  {
188  outfile << interpolated_q(s,i) << " ";
189  }
190 
191  // Output FE representation of div q at s
192  outfile << interpolated_div_q(s) << " ";
193 
194  // Output FE representation of p at s
195  outfile << interpolated_p(s) << " ";
196 
197  outfile << std::endl;
198  }
199 
200  // Write tecplot footer (e.g. FE connectivity lists)
201  this->write_tecplot_zone_footer(outfile,nplot);
202  }
203 
204  //=====================================================================
205  /// \short Output FE representation of exact soln: x,y,q1,q2,div_q,p at
206  /// Nplot^DIM plot points
207  //=====================================================================
208  template<unsigned DIM>
210  std::ostream &outfile,
211  const unsigned &nplot,
213  {
214  // Vector of local coordinates
215  Vector<double> s(DIM);
216 
217  // Vector for coordintes
218  Vector<double> x(DIM);
219 
220  // Tecplot header info
221  outfile << this->tecplot_zone_string(nplot);
222 
223  // Exact solution Vector
224  Vector<double> exact_soln(DIM+2);
225 
226  // Loop over plot points
227  unsigned num_plot_points=this->nplot_points(nplot);
228 
229  for(unsigned iplot=0;iplot<num_plot_points;iplot++)
230  {
231  // Get local coordinates of plot point
232  this->get_s_plot(iplot,nplot,s);
233 
234  // Get x position as Vector
235  this->interpolated_x(s,x);
236 
237  // Get exact solution at this point
238  (*exact_soln_pt)(x,exact_soln);
239 
240  // Output x,y,q_exact,p_exact,div_q_exact
241  for(unsigned i=0;i<DIM;i++)
242  {
243  outfile << x[i] << " ";
244  }
245  for(unsigned i=0;i<DIM+2;i++)
246  {
247  outfile << exact_soln[i] << " ";
248  }
249  outfile << std::endl;
250  }
251 
252  // Write tecplot footer (e.g. FE connectivity lists)
253  this->write_tecplot_zone_footer(outfile,nplot);
254  }
255 
256  //=====================================================================
257  /// \short Compute the error between the FE solution and the exact solution
258  /// using the H(div) norm for q and L^2 norm for p
259  //=====================================================================
260  template<unsigned DIM>
262  std::ostream &outfile,
264  Vector<double>& error,
265  Vector<double>& norm)
266  {
267  for(unsigned i=0;i<2;i++)
268  {
269  error[i]=0.0;
270  norm[i]=0.0;
271  }
272 
273  // Vector of local coordinates
274  Vector<double> s(DIM);
275 
276  // Vector for coordinates
277  Vector<double> x(DIM);
278 
279  // Set the value of n_intpt
280  unsigned n_intpt = this->integral_pt()->nweight();
281 
282  outfile << "ZONE" << std::endl;
283 
284  // Exact solution Vector (u,v,[w])
285  Vector<double> exact_soln(DIM+2);
286 
287  // Loop over the integration points
288  for(unsigned ipt=0;ipt<n_intpt;ipt++)
289  {
290  // Assign values of s
291  for(unsigned i=0;i<DIM;i++)
292  {
293  s[i] = this->integral_pt()->knot(ipt,i);
294  }
295 
296  // Get the integral weight
297  double w = this->integral_pt()->weight(ipt);
298 
299  // Get jacobian of mapping
300  double J=this->J_eulerian(s);
301 
302  // Premultiply the weights and the Jacobian
303  double W = w*J;
304 
305  // Get x position as Vector
306  this->interpolated_x(s,x);
307 
308  // Get exact solution at this point
309  (*exact_soln_pt)(x,exact_soln);
310 
311  // Flux error
312  for(unsigned i=0;i<DIM;i++)
313  {
314  norm[0]+=exact_soln[i]*exact_soln[i]*W;
315  // Error due to q_i
316  error[0]+=(exact_soln[i]-this->interpolated_q(s,i))*
317  (exact_soln[i]-this->interpolated_q(s,i))*W;
318  }
319 
320  // // Flux divergence error
321  // norm[0]+=exact_soln[DIM]*exact_soln[DIM]*W;
322  // error[0]+=(exact_soln[DIM]-interpolated_div_q(s))*
323  // (exact_soln[DIM]-interpolated_div_q(s))*W;
324 
325  // Pressure error
326  norm[1]+=exact_soln[DIM+1]*exact_soln[DIM+1]*W;
327  error[1]+=(exact_soln[DIM+1]-this->interpolated_p(s))*
328  (exact_soln[DIM+1]-this->interpolated_p(s))*W;
329 
330  // Output x,y,[z]
331  for(unsigned i=0;i<DIM;i++)
332  {
333  outfile << x[i] << " ";
334  }
335 
336  // Output q_1_error,q_2_error,...
337  for(unsigned i=0;i<DIM;i++)
338  {
339  outfile << exact_soln[i]-this->interpolated_q(s,i) << " ";
340  }
341 
342  // Output p_error
343  outfile << exact_soln[DIM+1]-this->interpolated_p(s) << " ";
344 
345  outfile << std::endl;
346  }
347  }
348 
349  //=====================================================================
350  /// Fill in residuals and, if flag==true, jacobian
351  //=====================================================================
352  template<unsigned DIM>
354  Vector<double> &residuals,
355  DenseMatrix<double> &jacobian,
356  bool flag)
357  {
358  // Get the number of geometric nodes, total number of basis functions,
359  // and number of edges basis functions
360  const unsigned n_node = nnode();
361  const unsigned n_q_basis = nq_basis();
362  const unsigned n_q_basis_edge = nq_basis_edge();
363  const unsigned n_p_basis = np_basis();
364 
365  // Storage for the geometric and computational bases and the test functions
366  Shape psi(n_node), q_basis(n_q_basis,DIM), q_test(n_q_basis,DIM),
367  p_basis(n_p_basis), p_test(n_p_basis),
368  div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
369 
370  // Get the number of integration points
371  unsigned n_intpt = integral_pt()->nweight();
372 
373  // Storage for the local coordinates
374  Vector<double> s(DIM);
375 
376  // Storage for the source function
377  Vector<double> f(DIM);
378 
379  // Storage for the mass source function
380  double mass_source_local;
381 
382  // Local equation and unknown numbers
383  int local_eqn = 0;//, local_unknown = 0;
384 
385  // Loop over the integration points
386  for(unsigned ipt=0;ipt<n_intpt;ipt++)
387  {
388  // Find the local coordinates at the integration point
389  for(unsigned i=0;i<DIM;i++) { s[i] = integral_pt()->knot(ipt,i); }
390 
391  // Get the weight of the intetgration point
392  double w = integral_pt()->weight(ipt);
393 
394  // Call the basis functions and test functions and get the
395  // (geometric) jacobian of the current element
396  double J =
397  shape_basis_test_local_at_knot(
398  ipt,psi,q_basis,q_test,p_basis,p_test,div_q_basis_ds,div_q_test_ds);
399 
400  // Storage for interpolated values
401  Vector<double> interpolated_x(DIM,0.0);
402  Vector<double> interpolated_q(DIM,0.0);
403  double interpolated_div_q_ds=0.0;
404  double interpolated_p=0.0;
405 
406  // loop over the geometric basis functions to find interpolated x
407  for(unsigned l=0;l<n_node;l++)
408  {
409  for(unsigned i=0;i<DIM;i++)
410  {
411  interpolated_x[i] += nodal_position(l,i)*psi[l];
412  }
413  }
414 
415  // loop over the nodes and use the vector basis functions to find the
416  // interpolated flux
417  for(unsigned i=0;i<DIM;i++)
418  {
419  // Loop over the edge basis vectors
420  for(unsigned l=0;l<n_q_basis_edge;l++)
421  {
422  interpolated_q[i] += q_edge(l)*q_basis(l,i);
423  }
424  // Loop over the internal basis vectors
425  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
426  {
427  interpolated_q[i] += q_internal(l-n_q_basis_edge)*q_basis(l,i);
428  }
429  }
430 
431  // loop over the pressure basis and find the interpolated pressure
432  for(unsigned l=0;l<n_p_basis;l++)
433  {
434  interpolated_p += p_value(l)*p_basis(l);
435  }
436 
437  // loop over the q edge divergence basis and the q internal divergence
438  // basis to find interpolated div q
439  for(unsigned l=0;l<n_q_basis_edge;l++)
440  {
441  interpolated_div_q_ds += q_edge(l)*div_q_basis_ds(l);
442  }
443  for(unsigned l=n_q_basis_edge;l<n_q_basis;l++)
444  {
445  interpolated_div_q_ds += q_internal(l-n_q_basis_edge)*div_q_basis_ds(l);
446  }
447 
448  // Get the source function
449  this->source(interpolated_x,f);
450 
451  // Get the mass source function
452  this->mass_source(interpolated_x,mass_source_local);
453 
454  // Loop over the test functions
455  for(unsigned l=0;l<n_q_basis;l++)
456  {
457  if(l<n_q_basis_edge)
458  {
459  local_eqn = q_edge_local_eqn(l);
460  }
461  else // n_q_basis_edge <= l < n_basis
462  {
463  local_eqn = q_internal_local_eqn(l-n_q_basis_edge);
464  }
465 
466  // If it's not a boundary condition
467  if(local_eqn>=0)
468  {
469  for(unsigned i=0;i<DIM;i++)
470  {
471  residuals[local_eqn] +=
472  (interpolated_q[i]-f[i])*q_test(l,i)*w*J;
473  }
474 
475  // deliberately no jacobian factor in this integral
476  residuals[local_eqn] -= (interpolated_p*div_q_test_ds(l))*w;
477  }
478  } // End of loop over test functions
479 
480  // loop over pressure test functions
481  for(unsigned l=0;l<n_p_basis;l++)
482  {
483  // get the local equation number
484  local_eqn = p_local_eqn(l);
485 
486  // If it's not a boundary condition
487  if(local_eqn>=0)
488  {
489  // deliberately no jacobian factor in this integral
490  residuals[local_eqn] += interpolated_div_q_ds*p_test(l)*w;
491  residuals[local_eqn] -= mass_source_local*p_test(l)*w*J;
492  }
493  }
494  } // End of loop over integration points
495  }
496 
497 
498  //=====================================================================
499  // Force building of templates
500  //=====================================================================
501  template class DarcyEquations<2>;
502 
503 } // End of oomph namespace
504 
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 ...
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.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output FE representation of exact soln: x,y,q1,q2,div_q,p at Nplot^DIM plot points.
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile)
Output with default number of plot points.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot, const Vector< double > unit_normal)
Output incl. projection of fluxes into direction of the specified unit vector.
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...