Tdarcy_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-2023 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//====================================================================
26 #include "Tdarcy_elements.h"
27 
28 namespace oomph
29 {
30  /// ///////////////////////////////////////////////////////////
31  /// ///////////////////////////////////////////////////////////
32  /// Lowest-order Raviart-Thomas based Darcy equation element
33  /// ///////////////////////////////////////////////////////////
34  /// ///////////////////////////////////////////////////////////
35 
36 
37  //===========================================================================
38  /// Constructor. Order 0 elements have 1 pressure dof and no internal
39  /// velocity dofs
40  //===========================================================================
41  template<>
43  : TElement<2, 3>(), DarcyEquations<2>(), Sign_edge(3, 1)
44  {
46  }
47 
48  //===========================================================================
49  /// Destructor
50  //===========================================================================
51  template<>
53  {
54  }
55 
56 
57  //===========================================================================
58  /// Return the number of edge basis functions for q
59  //===========================================================================
60  template<>
62  {
63  return 3;
64  }
65 
66 
67  //===========================================================================
68  /// Return the number of internal basis functions for q
69  //===========================================================================
70  template<>
72  {
73  return 0;
74  }
75 
76  //===========================================================================
77  /// Compute the local form of the q basis at local coordinate s
78  //===========================================================================
79  template<>
81  Shape& q_basis) const
82  {
83  // RT_0 basis functions
84  q_basis(0, 0) = Sign_edge[0] * std::sqrt(2) * s[0];
85  q_basis(0, 1) = Sign_edge[0] * std::sqrt(2) * s[1];
86 
87  q_basis(1, 0) = Sign_edge[1] * (s[0] - 1);
88  q_basis(1, 1) = Sign_edge[1] * s[1];
89 
90  q_basis(2, 0) = Sign_edge[2] * s[0];
91  q_basis(2, 1) = Sign_edge[2] * (s[1] - 1);
92  }
93 
94  //===========================================================================
95  /// Compute the local form of the q basis and dbasis/ds at local coordinate s
96  //===========================================================================
97  template<>
99  const Vector<double>& s, Shape& div_q_basis_ds) const
100  {
101  div_q_basis_ds(0) = Sign_edge[0] * 2 * std::sqrt(2);
102  div_q_basis_ds(1) = Sign_edge[1] * 2;
103  div_q_basis_ds(2) = Sign_edge[2] * 2;
104 
105  // Scale the basis by the ratio of the length of the edge to the length of
106  // the corresponding edge on the reference element
107  // hierher explain please
108  scale_basis(div_q_basis_ds);
109  }
110 
111  //===========================================================================
112  /// Return the number of flux interpolation points along each edge
113  /// of the element
114  //===========================================================================
115  template<>
117  {
118  return 1;
119  }
120 
121  //===========================================================================
122  /// Return the equation number of the n-th pressure degree of freedom
123  //===========================================================================
124  template<>
125  int TRaviartThomasDarcyElement<0>::p_local_eqn(const unsigned& n) const
126  {
127  return this->internal_local_eqn(P_internal_data_index, n);
128  }
129 
130  //===========================================================================
131  /// Return the nth pressure value
132  //===========================================================================
133  template<>
134  double TRaviartThomasDarcyElement<0>::p_value(const unsigned& n) const
135  {
136  return this->internal_data_pt(P_internal_data_index)->value(n);
137  }
138 
139  //===========================================================================
140  /// Return the total number of pressure basis functions
141  //===========================================================================
142  template<>
144  {
145  return 1;
146  }
147 
148  //===========================================================================
149  /// Compute the pressure basis
150  //===========================================================================
151  template<>
153  Shape& p_basis) const
154  {
155  p_basis(0) = 1.0;
156  }
157 
158 
159  //===========================================================================
160  /// Recovery order for Z2 error estimator
161  //===========================================================================
162  template<>
164  {
165  return 2; // Need to experiment with this...
166  }
167 
168  //===========================================================================
169  /// The number of values stored at each node: A single flux dof is
170  /// stored at each midside node
171  //===========================================================================
172  template<>
174  0, 0, 0, 1, 1, 1};
175 
176 
177  //===========================================================================
178  /// Face index associated with edge flux degree of freedom
179  //===========================================================================
180  template<>
182  2, 0, 1};
183 
184 
185  //===========================================================================
186  /// Conversion scheme from an edge degree of freedom to the node it's stored
187  /// at
188  /// Fluxes are stored at mid-side nodes
189  //===========================================================================
190  template<>
191  const unsigned TRaviartThomasDarcyElement<0>::Q_edge_conv[3] = {3, 4, 5};
192 
193 
194  //===========================================================================
195  /// The points along each edge where the fluxes are taken to be
196  // hierher explain please
197  //===========================================================================
198  template<>
200  0.5};
201 
202 
203  /// ///////////////////////////////////////////////////////////
204  /// ///////////////////////////////////////////////////////////
205  /// Second-order Raviart-Thomas based Darcy equation element
206  /// ///////////////////////////////////////////////////////////
207  /// ///////////////////////////////////////////////////////////
208 
209 
210  //===========================================================================
211  /// Constructor. Order 1 elements have 3 pressure dofs and 2 internal
212  /// velocity dofs
213  //===========================================================================
214  template<>
216  : TElement<2, 3>(), DarcyEquations<2>(), Sign_edge(3, 1)
217  {
218  // RT_1 elements have 2 internal degrees of freedom for q, and 3 for p
221  }
222 
223  //===========================================================================
224  /// Destructor
225  //===========================================================================
226  template<>
228  {
229  }
230 
231 
232  //===========================================================================
233  /// Return the number of edge basis functions for q
234  //===========================================================================
235  template<>
237  {
238  return 6;
239  }
240 
241 
242  //===========================================================================
243  /// Return the number of internal basis functions for q
244  //===========================================================================
245  template<>
247  {
248  return 2;
249  }
250 
251  //===========================================================================
252  /// Compute the local form of the q basis at local coordinate s
253  //===========================================================================
254  template<>
256  Shape& q_basis) const
257  {
258  // RT_1 basis functions
259 
260  Vector<double> g1_vect = edge_flux_interpolation_point(0, 0);
261  Vector<double> g2_vect = edge_flux_interpolation_point(0, 1);
262  double g1 = g1_vect[0];
263  double g2 = g2_vect[0];
264  q_basis(0, 0) =
265  Sign_edge[0] * std::sqrt(2.0) * s[0] * (s[1] - g2) / (g1 - g2);
266  q_basis(0, 1) =
267  Sign_edge[0] * std::sqrt(2.0) * s[1] * (s[1] - g2) / (g1 - g2);
268 
269  q_basis(1, 0) =
270  Sign_edge[0] * std::sqrt(2.0) * s[0] * (s[1] - g1) / (g2 - g1);
271  q_basis(1, 1) =
272  Sign_edge[0] * std::sqrt(2.0) * s[1] * (s[1] - g1) / (g2 - g1);
273 
274 
275  g1_vect = edge_flux_interpolation_point(1, 0);
276  g2_vect = edge_flux_interpolation_point(1, 1);
277  g1 = g1_vect[0];
278  g2 = g2_vect[0];
279  q_basis(2, 0) = Sign_edge[1] * (s[0] - 1.0) * (s[1] - g1) / (g2 - g1);
280  q_basis(2, 1) = Sign_edge[1] * s[1] * (s[1] - g1) / (g2 - g1);
281 
282  q_basis(3, 0) = Sign_edge[1] * (s[0] - 1.0) * (s[1] - g2) / (g1 - g2);
283  q_basis(3, 1) = Sign_edge[1] * s[1] * (s[1] - g2) / (g1 - g2);
284 
285 
286  g1_vect = edge_flux_interpolation_point(2, 0);
287  g2_vect = edge_flux_interpolation_point(2, 1);
288  g1 = g1_vect[0];
289  g2 = g2_vect[0];
290  q_basis(4, 0) = Sign_edge[2] * s[0] * (s[0] - g2) / (g1 - g2);
291  q_basis(4, 1) = Sign_edge[2] * (s[1] - 1.0) * (s[0] - g2) / (g1 - g2);
292 
293  q_basis(5, 0) = Sign_edge[2] * s[0] * (s[0] - g1) / (g2 - g1);
294  q_basis(5, 1) = Sign_edge[2] * (s[1] - 1.0) * (s[0] - g1) / (g2 - g1);
295 
296  q_basis(6, 0) = s[1] * s[0];
297  q_basis(6, 1) = s[1] * (s[1] - 1.0);
298 
299  q_basis(7, 0) = s[0] * (s[0] - 1.0);
300  q_basis(7, 1) = s[0] * s[1];
301  }
302 
303 
304  //===========================================================================
305  /// Compute the local form of the q basis and dbasis/ds at local coordinate s
306  //===========================================================================
307  template<>
309  const Vector<double>& s, Shape& div_q_basis_ds) const
310  {
311  Vector<double> g1_vect = edge_flux_interpolation_point(0, 0);
312  Vector<double> g2_vect = edge_flux_interpolation_point(0, 1);
313  double g1 = g1_vect[0];
314  double g2 = g2_vect[0];
315  div_q_basis_ds(0) =
316  Sign_edge[0] * std::sqrt(2.0) * (3.0 * s[1] - 2.0 * g2) / (g1 - g2);
317  div_q_basis_ds(1) =
318  Sign_edge[0] * std::sqrt(2.0) * (2.0 * g1 - 3.0 * s[1]) / (g1 - g2);
319 
320 
321  g1_vect = edge_flux_interpolation_point(1, 0);
322  g2_vect = edge_flux_interpolation_point(1, 1);
323  g1 = g1_vect[0];
324  g2 = g2_vect[0];
325  div_q_basis_ds(2) = Sign_edge[1] * (2.0 * g1 - 3.0 * s[1]) / (g1 - g2);
326  div_q_basis_ds(3) = Sign_edge[1] * (3.0 * s[1] - 2.0 * g2) / (g1 - g2);
327 
328 
329  g1_vect = edge_flux_interpolation_point(2, 0);
330  g2_vect = edge_flux_interpolation_point(2, 1);
331  g1 = g1_vect[0];
332  g2 = g2_vect[0];
333  div_q_basis_ds(4) = Sign_edge[2] * (3.0 * s[0] - 2.0 * g2) / (g1 - g2);
334  div_q_basis_ds(5) = Sign_edge[2] * (2.0 * g1 - 3.0 * s[0]) / (g1 - g2);
335 
336  div_q_basis_ds(6) = 3.0 * s[1] - 1.0;
337  div_q_basis_ds(7) = 3.0 * s[0] - 1.0;
338 
339  // Scale the basis by the ratio of the length of the edge to the length of
340  // the corresponding edge on the reference element
341  scale_basis(div_q_basis_ds);
342  }
343 
344  //===========================================================================
345  /// Return the number of flux_interpolation points along each edge of
346  /// the element
347  //===========================================================================
348  template<>
350  {
351  return 2;
352  }
353 
354  //===========================================================================
355  /// Return the equation number of the n-th pressure degree of freedom
356  //===========================================================================
357  template<>
358  int TRaviartThomasDarcyElement<1>::p_local_eqn(const unsigned& n) const
359  {
360  return this->internal_local_eqn(P_internal_data_index, n);
361  }
362 
363  //===========================================================================
364  /// Return the nth pressure value
365  //===========================================================================
366  template<>
367  double TRaviartThomasDarcyElement<1>::p_value(const unsigned& n) const
368  {
369  return this->internal_data_pt(P_internal_data_index)->value(n);
370  }
371 
372  //===========================================================================
373  /// Return the total number of pressure basis functions
374  //===========================================================================
375  template<>
377  {
378  return 3;
379  }
380 
381  //===========================================================================
382  /// Compute the pressure basis
383  //===========================================================================
384  template<>
386  Shape& p_basis) const
387  {
388  p_basis(0) = 1.0;
389  p_basis(1) = s[0];
390  p_basis(2) = s[1];
391  }
392 
393 
394  //===========================================================================
395  /// Recovery order for Z2 error estimator
396  //===========================================================================
397  template<>
399  {
400  return 2; // Need to experiment with this...
401  }
402 
403  //===========================================================================
404  /// The number of values stored at each node. Two flux values are stored
405  /// at each midside node.
406  //===========================================================================
407  template<>
409  0, 0, 0, 2, 2, 2};
410 
411 
412  //===========================================================================
413  /// Face index associated with edge flux degree of freedom
414  //===========================================================================
415  template<>
417  2, 2, 0, 0, 1, 1};
418 
419  //===========================================================================
420  /// Conversion scheme from an edge degree of freedom to the node it's stored
421  /// at
422  /// Fluxes are stored at midside nodes
423  //===========================================================================
424  template<>
425  const unsigned TRaviartThomasDarcyElement<1>::Q_edge_conv[3] = {3, 4, 5};
426 
427 
428  //===========================================================================
429  /// The points along each edge where the fluxes are taken to be
430  //===========================================================================
431  template<>
433  0.5 - std::sqrt(3.0) / 6.0, 0.5 + std::sqrt(3.0) / 6.0};
434 
435 
436  //===========================================================================
437  // Force building of templates
438  //===========================================================================
439  template class TRaviartThomasDarcyElement<0>;
440  template class TRaviartThomasDarcyElement<1>;
441 
442 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
Class implementing the generic maths of the Darcy equations using Raviart-Thomas elements with both e...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
General TElement class.
Definition: Telements.h:1208
Element which solves the Darcy equations using TElements. Geometrically the element is always a six n...
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.
int p_local_eqn(const unsigned &n) const
Return the equation number of the n-th pressure degree of freedom.
double p_value(const unsigned &n) const
Return the nth pressure value.
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 np_basis() const
Return the total number of pressure basis functions.
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 nrecovery_order()
Recovery order for Z2 error estimator.
unsigned Q_internal_data_index
The internal data index where the internal q degrees of freedom are stored.
unsigned nedge_flux_interpolation_point() const
Return the number of flux interpolation points along each edge of the element.
unsigned nq_basis_internal() const
Return the number of internal basis functions for q.
unsigned nq_basis_edge() const
Return the number of edge basis functions for q.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...