Tdarcy_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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//====================================================================
30 #ifndef OOMPH_TRAVIART_THOMAS_DARCY_HEADER
31 #define OOMPH_TRAVIART_THOMAS_DARCY_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 #include "darcy_elements.h"
39 #include "../generic/Telements.h"
40 
41 namespace oomph
42 {
43 
44 
45  //================================================================
46  /// Element which solves the Darcy equations using TElements.
47  /// Geometrically the element is always a six noded triangle.
48  /// We use the mid-side nodes to store edge-based flux degrees of
49  /// freedom and internal data for the discontinuous pressure
50  /// and internal flux dofs.
51  //================================================================
52  template<unsigned ORDER>
53  class TRaviartThomasDarcyElement : public TElement<2,3>,
54  public DarcyEquations<2>
55  {
56  private:
57 
58  /// Face index associated with edge flux degree of freedom
59  static const unsigned Face_index_of_edge_flux[];
60 
61  /// \short Conversion scheme from an edge degree of freedom to the node it's
62  /// stored at
63  static const unsigned Q_edge_conv[];
64 
65  /// The points along each edge where the fluxes are interpolated
66  static const double Flux_interpolation_point[];
67 
68  /// The internal data index where the internal q degrees of freedom are stored
70 
71  /// The internal data index where the p degrees of freedom are stored
73 
74  /// \short Unit normal signs associated with each edge to ensure inter-element
75  /// continuity of the flux
76  std::vector<short> Sign_edge;
77 
78  public:
79 
80  /// Constructor
82 
83  /// Destructor
85 
86  /// Number of values required at node n
87  unsigned required_nvalue(const unsigned &n) const
88  {
89  return Initial_Nvalue[n];
90  }
91 
92  /// Return the face index associated with specified edge
93  unsigned face_index_of_edge(const unsigned& j) const
94  {
95  return (j+2)%3;
96  }
97 
98  /// \short Compute the face element coordinates of the nth flux interpolation
99  /// point along specified edge
101  const unsigned &n,
102  Vector<double> &s)
103  const
104  {
105  // Get the location of the n-th flux interpolation point along
106  // the edge in terms of the distance along the edge itself
107  Vector<double> flux_interpolation_point=
109 
110  // Convert the edge number to the number of the mid-edge node along that
111  // edge
112  unsigned node_number=Q_edge_conv[edge];
113 
114  // The edge basis functions are defined in a clockwise manner, so we have
115  // to effectively "flip" some coordinates
116  switch(node_number)
117  {
118  case 3:
119  s[0]=flux_interpolation_point[0];
120  break;
121  case 4:
122  s[0]=1.0-flux_interpolation_point[0];
123  break;
124  case 5:
125  s[0]=flux_interpolation_point[0];
126  break;
127  }
128 
129  }
130 
131  /// Return the face index associated with edge flux degree of freedom
132  unsigned face_index_of_q_edge_basis_fct(const unsigned& j) const
133  {
134  return Face_index_of_edge_flux[j];
135  }
136 
137  /// Return the equation nunmber of the n-th edge (flux) degree of freedom
138  int q_edge_local_eqn(const unsigned &n) const
139  {
140  return this->nodal_local_eqn(q_edge_node_number(n),q_edge_index(n));
141  }
142 
143  /// Return the equation number of the n-th internal degree of freedom
144  int q_internal_local_eqn(const unsigned &n) const
145  {
147  }
148 
149  /// \short Return vector of pointers to the Data objects that store the
150  /// edge flux values
152  {
153  // It's the mid-side nodes:
154  Vector<Data*> data_pt(3);
155  data_pt[0]=node_pt(3);
156  data_pt[1]=node_pt(4);
157  data_pt[2]=node_pt(5);
158  return data_pt;
159  }
160 
161  /// Return pointer to the Data object that stores the internal flux values
163  {
165  }
166 
167  /// \short Return the index of the internal data where the q_internal degrees
168  /// of freedom are stored
169  unsigned q_internal_index() const
170  {
171  return Q_internal_data_index;
172  }
173 
174  /// Return the nodal index at which the nth edge unknown is stored
175  unsigned q_edge_index(const unsigned &n) const
176  {
177  return n%(ORDER+1);
178  }
179 
180  /// \short Return the local node number of the node where the nth edge
181  /// unknown is stored
182  unsigned q_edge_node_number(const unsigned &n) const
183  {
184  return Q_edge_conv[n/(ORDER+1)];
185  }
186 
187  /// Get pointer to node that stores the edge flux dofs for specified edge
188  Node* edge_flux_node_pt(const unsigned& edge)
189  {
190  return node_pt(Q_edge_conv[edge]);
191  }
192 
193  /// Return the values of the edge (flux) degree of freedom
194  double q_edge(const unsigned &n) const
195  {
197  }
198 
199  /// Return the values of the internal degree of freedom
200  double q_internal(const unsigned &n) const
201  {
202  return this->internal_data_pt(q_internal_index())->value(n);
203  }
204 
205  /// Set the values of the edge (flux) degree of freedom
206  void set_q_edge(const unsigned &n, const double& value)
207  {
209  }
210 
211  /// Set the values of the internal degree of freedom
212  void set_q_internal(const unsigned &n, const double& value)
213  {
214  this->internal_data_pt(q_internal_index())->set_value(n,value);
215  }
216 
217  /// Return the number of edge basis functions for q
218  unsigned nq_basis_edge() const;
219 
220  /// Return the number of internal basis functions for q
221  unsigned nq_basis_internal() const;
222 
223  /// Return the local form of the q basis at local coordinate s
224  void get_q_basis_local(const Vector<double> &s,
225  Shape &q_basis) const;
226 
227  /// Return the local form of the q basis and dbasis/ds at local coordinate s
229  Shape &div_q_basis_ds) const;
230 
231  /// \short Return the number of flux interpolation points along each
232  /// edge of the element
233  unsigned nedge_flux_interpolation_point() const;
234 
235  /// Return the local coordinate of the nth flux interpolation
236  /// point along an edge
238  const unsigned &n) const
239  {
240  Vector<double> coord(1);
241  coord[0]=(1.0-sign_edge(edge))/2.0+sign_edge(edge)*
243  return coord;
244  }
245 
246 
247  /// \short Compute the global coordinates of the flux interpolation
248  /// point associated with the j-th edge-based q basis fct
249  void edge_flux_interpolation_point_global(const unsigned &j,
250  Vector<double> &x) const
251  {
252  unsigned n=j%nedge_flux_interpolation_point();
253  unsigned edge=(j-n)/nedge_flux_interpolation_point();
255  }
256 
257 
258  /// \short Compute the global coordinates of the nth flux interpolation
259  /// point along an edge
260  void edge_flux_interpolation_point_global(const unsigned &edge,
261  const unsigned &n,
262  Vector<double> &x) const
263  {
264  // Get the location of the n-th flux interpolation point along
265  // the edge in terms of the distance along the edge itself
266  Vector<double> flux_interpolation_point=
268 
269  // Convert the edge number to the number of the mid-edge node along that
270  // edge
271  unsigned node_number=Q_edge_conv[edge];
272 
273  // Storage for the local coords of the flux_interpolation point
274  Vector<double> s_flux_interpolation(2,0.0);
275 
276  // The edge basis functions are defined in a clockwise manner, so we have
277  // to effectively "flip" the coordinates along edges 0 and 1 to match this
278  switch(node_number)
279  {
280  case 3:
281  s_flux_interpolation[0]=1.0-flux_interpolation_point[0];
282  s_flux_interpolation[1]=flux_interpolation_point[0];
283  break;
284  case 4:
285  s_flux_interpolation[0]=0.0;
286  s_flux_interpolation[1]=1.0-flux_interpolation_point[0];
287  break;
288  case 5:
289  s_flux_interpolation[0]=flux_interpolation_point[0];
290  s_flux_interpolation[1]=0.0;
291  break;
292  }
293 
294  // Calculate the global coordinates from the local ones
295  interpolated_x(s_flux_interpolation,x);
296  }
297 
298  /// Pin the nth internal q value
299  void pin_q_internal_value(const unsigned &n)
300  {
302  }
303 
304  /// Return the equation number of the n-th pressure degree of freedom
305  int p_local_eqn(const unsigned &n) const;
306 
307  /// Return the nth pressure value
308  double p_value(const unsigned &n) const;
309 
310  /// Return the total number of pressure basis functions
311  unsigned np_basis() const;
312 
313  /// Compute the pressure basis
314  void get_p_basis(const Vector<double> &s,
315  Shape &p_basis) const;
316 
317  /// Pin the nth pressure value
318  void pin_p_value(const unsigned &n)
319  {
321  }
322 
323  /// Return pointer to the Data object that stores the pressure values
324  Data* p_data_pt() const
325  {
327  }
328 
329  /// Set the nth pressure value
330  void set_p_value(const unsigned &n, const double& value)
331  {
333  }
334 
335  /// Scale the edge basis to allow arbitrary edge mappings
336  // hierher explain please
337  void scale_basis(Shape& basis) const
338  {
339 
340  // Storage for the lengths of the edges of the element
341  Vector<double> length(3,0.0);
342 
343  // Temporary storage for the vertex positions
344  double x0,y0,x1,y1;
345 
346  // loop over the edges of the element and calculate their lengths (in x-y
347  // space)
348  for(unsigned i=0;i<3;i++)
349  {
350  x0=this->node_pt(i)->x(0);
351  y0=this->node_pt(i)->x(1);
352  x1=this->node_pt((i+1)%3)->x(0);
353  y1=this->node_pt((i+1)%3)->x(1);
354 
355  length[i]=std::sqrt(std::pow(y1-y0,2)+std::pow(x1-x0,2));
356  }
357 
358  // lengths of the sides of the reference element (in the same order as the
359  // basis functions)
360  const double ref_length[3]={std::sqrt(2.0),1,1};
361 
362  // get the number of basis functions associated with the edges
363  unsigned n_q_basis_edge=nq_basis_edge();
364 
365  // rescale the edge basis functions to allow arbitrary edge mappings from
366  // element to ref. element
367  const unsigned n_index2 = basis.nindex2();
368  for(unsigned i=0;i<n_index2;i++)
369  {
370  for(unsigned l=0;l<n_q_basis_edge;l++)
371  {
372  basis(l,i)*=(length[l/(ORDER+1)]/ref_length[l/(ORDER+1)]);
373  }
374  }
375  }
376 
377  /// Accessor for the unit normal sign of edge n (const version)
378  const short &sign_edge(const unsigned &n) const
379  {
380  return Sign_edge[n];
381  }
382 
383  /// Accessor for the unit normal sign of edge n
384  short &sign_edge(const unsigned &n)
385  {
386  return Sign_edge[n];
387  }
388 
389  /// Output with default number of plot points
390  void output(std::ostream &outfile)
391  {
392  DarcyEquations<2>::output(outfile);
393  }
394 
395  /// \short Output FE representation of soln: x,y,u1,u2,div_q,p at
396  /// Nplot^DIM plot points
397  void output(std::ostream &outfile, const unsigned &Nplot)
398  {
399  DarcyEquations<2>::output(outfile,Nplot);
400  }
401 
402  /// \short Number of vertex nodes in the element
403  unsigned nvertex_node() const
404  {return TElement<2,3>::nvertex_node();}
405 
406  /// \short Pointer to the j-th vertex node in the element
407  Node* vertex_node_pt(const unsigned& j) const
408  {return TElement<2,3>::vertex_node_pt(j);}
409 
410 
411 /// Recovery order for Z2 error estimator
412  unsigned nrecovery_order();
413 
414  protected:
415 
416  /// \short Return the geometric basis, and the q, p and divergence basis
417  /// functions and test functions at local coordinate s
419  Shape &psi,
420  Shape &q_basis,
421  Shape &q_test,
422  Shape &p_basis,
423  Shape &p_test,
424  Shape &div_q_basis_ds,
425  Shape &div_q_test_ds) const
426  {
427  const unsigned n_q_basis = this->nq_basis();
428 
429  Shape q_basis_local(n_q_basis,2);
430  this->get_q_basis_local(s,q_basis_local);
431  this->get_p_basis(s,p_basis);
432  this->get_div_q_basis_local(s,div_q_basis_ds);
433 
434  double J = this->transform_basis(s,q_basis_local,psi,q_basis);
435 
436  q_test = q_basis;
437  p_test = p_basis;
438  div_q_test_ds = div_q_basis_ds;
439 
440  return J;
441  }
442 
443  /// \short Return the geometric basis, and the q, p and divergence basis
444  /// functions and test functions at integration point ipt
445  double shape_basis_test_local_at_knot(const unsigned &ipt,
446  Shape &psi,
447  Shape &q_basis,
448  Shape &q_test,
449  Shape &p_basis,
450  Shape &p_test,
451  Shape &div_q_basis_ds,
452  Shape &div_q_test_ds) const
453  {
454  Vector<double> s(2);
455  for(unsigned i=0;i<2;i++) { s[i] = this->integral_pt()->knot(ipt,i); }
456 
457  return shape_basis_test_local(
458  s,psi,q_basis,q_test,p_basis,p_test,div_q_basis_ds,div_q_test_ds);
459  }
460 
461  /// The number of values stored at each node
462  static const unsigned Initial_Nvalue[];
463 
464  };
465 
466 
467 
468 
469  //============================================================
470  /// Face geometry for TRaviartThomasDarcyElement<0>
471  //============================================================
472  template<>
474  public virtual TElement<1,3>
475  {
476  public:
477 
478  /// Constructor: Call constructor of base
479  FaceGeometry() : TElement<1,3>() {}
480  };
481 
482 
483  //============================================================
484  /// Face geometry for TRaviartThomasDarcyElement<1>
485  //============================================================
486  template<>
488  public virtual TElement<1,3>
489  {
490  public:
491 
492  /// Constructor: Call constructor of base class
493  FaceGeometry() : TElement<1,3>() {}
494  };
495 
496 
497 } // End of oomph namespace
498 
499 #endif
500 
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
Class implementing the generic maths of the Darcy equations using Raviart-Thomas elements with both e...
unsigned nq_basis_edge() const
Return the number of edge basis functions for q.
short & sign_edge(const unsigned &n)
Accessor for the unit normal sign of edge n.
unsigned nrecovery_order()
Recovery order for Z2 error estimator.
int p_local_eqn(const unsigned &n) const
Return the equation number of the n-th pressure degree of freedom.
double q_internal(const unsigned &n) const
Return the values of the internal degree of freedom.
TRaviartThomasDarcyElement()
Constructor.
cstr elem_len * i
Definition: cfortran.h:607
unsigned Q_internal_data_index
The internal data index where the internal q degrees of freedom are stored.
double shape_basis_test_local(const Vector< double > &s, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const
Return the geometric basis, and the q, p and divergence basis functions and test functions at local c...
double p_value(const unsigned &n) const
Return the nth pressure value.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
Vector< Data * > q_edge_data_pt() const
Return vector of pointers to the Data objects that store the edge flux values.
unsigned q_edge_node_number(const unsigned &n) const
Return the local node number of the node where the nth edge unknown is stored.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void set_q_internal(const unsigned &n, const double &value)
Set the values of the internal degree of freedom.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
double shape_basis_test_local_at_knot(const unsigned &ipt, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const
Return the geometric basis, and the q, p and divergence basis functions and test functions at integra...
Vector< double > edge_flux_interpolation_point(const unsigned &edge, const unsigned &n) const
void set_p_value(const unsigned &n, const double &value)
Set the nth pressure value.
int q_edge_local_eqn(const unsigned &n) const
Return the equation nunmber of the n-th edge (flux) degree of freedom.
FaceGeometry()
Constructor: Call constructor of base.
void output(std::ostream &outfile, const unsigned &Nplot)
Output FE representation of soln: x,y,u1,u2,div_q,p at Nplot^DIM plot points.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
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:3841
double q_edge(const unsigned &n) const
Return the values of the edge (flux) degree of freedom.
void scale_basis(Shape &basis) const
Scale the edge basis to allow arbitrary edge mappings.
void face_local_coordinate_of_flux_interpolation_point(const unsigned &edge, const unsigned &n, Vector< double > &s) const
Compute the face element coordinates of the nth flux interpolation point along specified edge...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned q_edge_index(const unsigned &n) const
Return the nodal index at which the nth edge unknown is stored.
FaceGeometry()
Constructor: Call constructor of base class.
unsigned nq_basis() const
Return the total number of computational basis functions for q.
static char t char * s
Definition: cfortran.h:572
Data * p_data_pt() const
Return pointer to the Data object that stores the pressure values.
void set_q_edge(const unsigned &n, const double &value)
Set the values of the edge (flux) degree of freedom.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void pin_p_value(const unsigned &n)
Pin the nth pressure value.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void output(std::ostream &outfile)
Output with default number of plot points.
void pin_q_internal_value(const unsigned &n)
Pin the nth internal q value.
unsigned nedge_flux_interpolation_point() const
Return the number of flux interpolation points along each edge of the element.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void edge_flux_interpolation_point_global(const unsigned &j, Vector< double > &x) const
Compute the global coordinates of the flux interpolation point associated with the j-th edge-based q ...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void get_p_basis(const Vector< double > &s, Shape &p_basis) const
Compute the pressure basis.
unsigned P_internal_data_index
The internal data index where the p degrees of freedom are stored.
unsigned face_index_of_edge(const unsigned &j) const
Return the face index associated with specified edge.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
unsigned nvertex_node() const
Number of vertex nodes in the element.
std::vector< short > Sign_edge
Unit normal signs associated with each edge to ensure inter-element continuity of the flux...
unsigned required_nvalue(const unsigned &n) const
Number of values required at node n.
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:1383
unsigned q_internal_index() const
Return the index of the internal data where the q_internal degrees of freedom are stored...
int q_internal_local_eqn(const unsigned &n) const
Return the equation number of the n-th internal degree of freedom.
unsigned np_basis() const
Return the total number of pressure basis functions.
const short & sign_edge(const unsigned &n) const
Accessor for the unit normal sign of edge n (const version)
unsigned nq_basis_internal() const
Return the number of internal basis functions for q.
Node * edge_flux_node_pt(const unsigned &edge)
Get pointer to node that stores the edge flux dofs for specified edge.
void edge_flux_interpolation_point_global(const unsigned &edge, const unsigned &n, Vector< double > &x) const
Compute the global coordinates of the nth flux interpolation point along an edge. ...
Data * q_internal_data_pt() const
Return pointer to the Data object that stores the internal flux values.
void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const
Return the local form of the q basis and dbasis/ds at local coordinate s.
unsigned face_index_of_q_edge_basis_fct(const unsigned &j) const
Return the face index associated with edge flux degree of freedom.
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
void output(std::ostream &outfile)
Output with default number of plot points.
void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const
Return the local form of the q basis at local coordinate s.
unsigned nindex2() const
Return the range of index 2 of the shape function object.
Definition: shape.h:265
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 ...