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-2025 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
28namespace 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<>
47
48 //===========================================================================
49 /// Destructor
50 //===========================================================================
51 template<>
55
56
57 //===========================================================================
58 /// Return the number of edge basis functions for q
59 //===========================================================================
60 template<>
61 unsigned TRaviartThomasDarcyElement<0>::nq_basis_edge() const
62 {
63 return 3;
64 }
65
66
67 //===========================================================================
68 /// Return the number of internal basis functions for q
69 //===========================================================================
70 template<>
71 unsigned TRaviartThomasDarcyElement<0>::nq_basis_internal() const
72 {
73 return 0;
74 }
75
76 //===========================================================================
77 /// Compute the local form of the q basis at local coordinate s
78 //===========================================================================
79 template<>
80 void TRaviartThomasDarcyElement<0>::get_q_basis_local(const Vector<double>& s,
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<>
98 void TRaviartThomasDarcyElement<0>::get_div_q_basis_local(
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<>
116 unsigned TRaviartThomasDarcyElement<0>::nedge_flux_interpolation_point() const
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<>
143 unsigned TRaviartThomasDarcyElement<0>::np_basis() const
144 {
145 return 1;
146 }
147
148 //===========================================================================
149 /// Compute the pressure basis
150 //===========================================================================
151 template<>
152 void TRaviartThomasDarcyElement<0>::get_p_basis(const Vector<double>& s,
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<>
163 unsigned TRaviartThomasDarcyElement<0>::nrecovery_order()
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<>
173 const unsigned TRaviartThomasDarcyElement<0>::Initial_Nvalue[6] = {
174 0, 0, 0, 1, 1, 1};
175
176
177 //===========================================================================
178 /// Face index associated with edge flux degree of freedom
179 //===========================================================================
180 template<>
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<>
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<>
215 TRaviartThomasDarcyElement<1>::TRaviartThomasDarcyElement()
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
219 Q_internal_data_index = this->add_internal_data(new Data(2));
220 P_internal_data_index = this->add_internal_data(new Data(3));
221 }
222
223 //===========================================================================
224 /// Destructor
225 //===========================================================================
226 template<>
230
231
232 //===========================================================================
233 /// Return the number of edge basis functions for q
234 //===========================================================================
235 template<>
236 unsigned TRaviartThomasDarcyElement<1>::nq_basis_edge() const
237 {
238 return 6;
239 }
240
241
242 //===========================================================================
243 /// Return the number of internal basis functions for q
244 //===========================================================================
245 template<>
246 unsigned TRaviartThomasDarcyElement<1>::nq_basis_internal() const
247 {
248 return 2;
249 }
250
251 //===========================================================================
252 /// Compute the local form of the q basis at local coordinate s
253 //===========================================================================
254 template<>
255 void TRaviartThomasDarcyElement<1>::get_q_basis_local(const Vector<double>& s,
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<>
308 void TRaviartThomasDarcyElement<1>::get_div_q_basis_local(
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<>
349 unsigned TRaviartThomasDarcyElement<1>::nedge_flux_interpolation_point() const
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<>
376 unsigned TRaviartThomasDarcyElement<1>::np_basis() const
377 {
378 return 3;
379 }
380
381 //===========================================================================
382 /// Compute the pressure basis
383 //===========================================================================
384 template<>
385 void TRaviartThomasDarcyElement<1>::get_p_basis(const Vector<double>& s,
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<>
398 unsigned TRaviartThomasDarcyElement<1>::nrecovery_order()
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<>
408 const unsigned TRaviartThomasDarcyElement<1>::Initial_Nvalue[6] = {
409 0, 0, 0, 2, 2, 2};
410
411
412 //===========================================================================
413 /// Face index associated with edge flux degree of freedom
414 //===========================================================================
415 template<>
416 const unsigned TRaviartThomasDarcyElement<1>::Face_index_of_edge_flux[6] = {
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<>
432 const double TRaviartThomasDarcyElement<1>::Flux_interpolation_point[2] = {
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
Tdarcy_elements.h
s
static char t char * s
Definition cfortran.h:568
oomph::DarcyEquations
Class implementing the generic maths of the Darcy equations using Raviart-Thomas elements with both e...
Definition darcy_elements.h:49
oomph::Data
A class that represents a collection of data; each Data object may contain many different individual ...
Definition nodes.h:86
oomph::Data::value
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:293
oomph::GeneralisedElement::add_internal_data
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:67
oomph::GeneralisedElement::internal_data_pt
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition elements.h:605
oomph::GeneralisedElement::internal_local_eqn
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:267
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition shape.h:76
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TAdvectionDiffusionReactionElement::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Definition Tadvection_diffusion_reaction_elements.h:70
oomph::TElement
General TElement class.
Definition Telements.h:1208
oomph::TRaviartThomasDarcyElement
Element which solves the Darcy equations using TElements. Geometrically the element is always a six n...
Definition Tdarcy_elements.h:49
oomph::TRaviartThomasDarcyElement::get_q_basis_local
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.
oomph::TRaviartThomasDarcyElement::p_local_eqn
int p_local_eqn(const unsigned &n) const
Return the equation number of the n-th pressure degree of freedom.
oomph::TRaviartThomasDarcyElement::p_value
double p_value(const unsigned &n) const
Return the nth pressure value.
oomph::TRaviartThomasDarcyElement::get_div_q_basis_local
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.
oomph::TRaviartThomasDarcyElement::np_basis
unsigned np_basis() const
Return the total number of pressure basis functions.
oomph::TRaviartThomasDarcyElement::get_p_basis
void get_p_basis(const Vector< double > &s, Shape &p_basis) const
Compute the pressure basis.
oomph::TRaviartThomasDarcyElement::P_internal_data_index
unsigned P_internal_data_index
The internal data index where the p degrees of freedom are stored.
Definition Tdarcy_elements.h:66
oomph::TRaviartThomasDarcyElement::nrecovery_order
unsigned nrecovery_order()
Recovery order for Z2 error estimator.
oomph::TRaviartThomasDarcyElement::Q_internal_data_index
unsigned Q_internal_data_index
The internal data index where the internal q degrees of freedom are stored.
Definition Tdarcy_elements.h:63
oomph::TRaviartThomasDarcyElement::nedge_flux_interpolation_point
unsigned nedge_flux_interpolation_point() const
Return the number of flux interpolation points along each edge of the element.
oomph::TRaviartThomasDarcyElement::nq_basis_internal
unsigned nq_basis_internal() const
Return the number of internal basis functions for q.
oomph::TRaviartThomasDarcyElement::nq_basis_edge
unsigned nq_basis_edge() const
Return the number of edge basis functions for q.
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30