Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONSISTENT_GRAVITY__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONSISTENT_GRAVITY__
35 :
36 : #include <vector>
37 :
38 : // other ug4 modules
39 : #include "common/common.h"
40 : #include "lib_disc/reference_element/reference_mapping.h"
41 :
42 : namespace ug{
43 :
44 : /// Class for the computation of the standard version ('Voss-Souza-type') of the consistent gravity
45 : /**
46 : * Density driven flow models involve the convection velocity depending on the
47 : * gradient of the pressure and the density-dependent gravity force, typically
48 : * of the form \f$\mathbf{v} = \mathbf{K} (- \nabla p + \rho \mathbf{g}) / \mu\f$,
49 : * where \f$\mathbf{K}\f$ is a matrix, \f$\mu\f$ a scalar, \f$\mathbf{g}\f$
50 : * gravity vector (all independent of \f$p\f$ or \f$\rho\f$), \f$p\f$ is
51 : * pressure, \f$rho\f$ is density (depending on the concentration, ...). (Both
52 : * the pressure and the concentration are unknown grid functions in the PDEs.)
53 : * Treating \f$p\f$ and \f$\rho\f$ as grid functions of the same class (like
54 : * the piecewise linear/bilinear interpolations of nodal values) leads to oscillatory
55 : * solutions because \f$\nabla p\f$ and \f$\rho \mathbf{g}\f$ belong to different
56 : * classes of grid functions. In particular, if \f$p\f$ is a piecewise linear
57 : * function, \f$\nabla p\f$ can never cancel the of the gravitaty force if
58 : * \f$\rho\f$ increases linearly from the top down, although this happens in the
59 : * analytical solution. The idea of the consistent gravity is to consider a
60 : * vector function \f$\mathbf{h} = (h_x, h_y, \dots)\f$ that is in some sence a
61 : * primitive function for \f$\rho \mathbf{g}\f$ so that the velocity can be
62 : * written as \f$\mathbf{v} = - \mathbf{K} (p_x - h_x, p_y - h_y, \dots) / \mu\f$.
63 : * In the discretization, the function \f$\mathbf{h}\f$ should belong to the
64 : * same class of the grid functions as \f$p\f$. In the present implementation,
65 : * \f$\mathbf{h}\f$ is a piecewise linear/bilinear grid function.
66 : *
67 : * Class method 'prepare' computes the values of \f$h\f$ at the corners
68 : * of an element. Using these values, method 'compute' computes the
69 : * consistent gravity force \f$\rho \mathbf{g}\f$ (not \f$\mathbf{h}\f$!) at any
70 : * given point.
71 : *
72 : * \remark (Consistent gravity in the computation of the Jacobian)
73 : * The consistent gravity computed by the implemented methods depends linearly
74 : * on the density. Thus, to compute the derivative of the consistent gravity
75 : * w.r.t. the density at one of the corners of the element, set the density to
76 : * 1 at that corner and to 0 at all other corners, then call the functions. Then
77 : * 'prepare' prepares the nodal values of the derivatives and 'compute' computes
78 : * the derivative itself.
79 : *
80 : * Alternatively, you can set the density at that corner to the derivative
81 : * of the density w.r.t. the concentration (instead of 1). Then you get
82 : * directly the derivative of the consistent gravity w.r.t. the concentration
83 : * at that node.
84 : *
85 : * \remark There is an enhanced version of the consistent gravity. \see StdLinConsistentGravityX
86 : *
87 : * References:
88 : * <ul>
89 : * <li> P. Frolkovic, P. Knabner, Consistent Velocity Approximations in Finite
90 : * Element or Volume Discretizations of Density Driven Flow, In: Computational
91 : * Methods in WaterResources XI, Vol. 1 (A.A. Aldama et al., eds.),
92 : * Computational Mechanics Publication, Southhampten, 1996, p. 93-100
93 : * </li>
94 : * </ul>
95 : *
96 : * \tparam refDim dimensionality of the reference element (e.g. 2 for triangles, 3 for tetrahedra)
97 : */
98 : template <int refDim>
99 : class StdLinConsistentGravity
100 : {
101 : private:
102 : // static constants
103 : static const size_t _X_ = 0;
104 : static const size_t _Y_ = 1;
105 : static const size_t _Z_ = 2;
106 :
107 : public:
108 :
109 : /// constructor (sets the 'not init.' flag)
110 0 : StdLinConsistentGravity () : m_nCo (0) {};
111 :
112 : /// computation of the primary function for the consistent gravity at corners, cf. the specializations
113 : template <int dim>
114 : inline void prepare
115 : (
116 : MathVector<refDim>* vConsGravity, ///< where to save the values (n_co vectors)
117 : const int n_co, ///< number of corners of the element
118 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
119 : const number* vDensity, ///< corner density (n_co scalar values)
120 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
121 : )
122 : {
123 : UG_THROW ("StdLinConsistentGravity: Combination of the world dim " << dim <<
124 : "and the reference element dim " << refDim << " is not implemented.");
125 : }
126 :
127 : /// computation of the consistent gravity at a given point
128 : template <int dim>
129 0 : inline void compute
130 : (
131 : MathVector<dim>& ConsistentGravity, ///< where to save the vector
132 : const MathVector<refDim>& LocalCoord, ///< local coordinates of the point
133 : const MathMatrix<dim, refDim>& JTInv, ///< inverse transposed Jacobian
134 : const MathVector<refDim>* vLocalGrad, ///< gradients of the shape functions at the given point
135 : const MathVector<refDim>* vConsGravity ///< primary function of the consistent gravity at corners
136 : )
137 : {
138 : UG_ASSERT (m_nCo > 0, "StdLinConsistentGravity: Object not initialized.");
139 :
140 : MathVector<refDim> LocalGravity;
141 : VecSet(LocalGravity, 0.0);
142 :
143 : // Loop shape functions
144 0 : for(size_t sh = 0; sh < (size_t) m_nCo; sh++)
145 0 : for(size_t d = 0; d < refDim; d++)
146 0 : LocalGravity[d] += vConsGravity[sh][d] * vLocalGrad[sh][d];
147 :
148 : // Multiply by JacobianTransposedInverse
149 : MatVecMult(ConsistentGravity, JTInv, LocalGravity);
150 0 : }
151 :
152 : protected:
153 :
154 : int m_nCo; ///< number of corners of the element for which the object is init. (0 if not init)
155 :
156 : /// computation of the primary function for the consistent gravity at corners of an edge
157 : template <int dim>
158 0 : inline void prepare_edge
159 : (
160 : MathVector<1>* vConsGravity, ///< where to save the values (2 vectors)
161 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (2 vectors)
162 : const number* vDensity, ///< corner density (2 scalar values)
163 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
164 : )
165 : {
166 : MathVector<1> LocalPoint;
167 : MathVector<1> LocalGravity;
168 : MathMatrix<1,dim> JT;
169 0 : static ReferenceMapping<ReferenceEdge, dim> EdgeMapping;
170 : EdgeMapping.update (vCorners);
171 :
172 : /* compute the local gravity */
173 : VecSet (LocalPoint, 0.0);
174 : EdgeMapping.jacobian_transposed (JT, LocalPoint);
175 : MatVecMult (LocalGravity, JT, PhysicalGravity);
176 :
177 0 : vConsGravity[0][_X_] = 0.0;
178 0 : vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
179 0 : }
180 :
181 : /// computation of the primary function for the consistent gravity at corners of a triangle
182 : template <int dim>
183 0 : inline void prepare_triangle
184 : (
185 : MathVector<2>* vConsGravity, ///< where to save the values (3 vectors)
186 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (3 vectors)
187 : const number* vDensity, ///< corner density (3 scalar values)
188 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
189 : )
190 : {
191 : MathVector<2> LocalPoint;
192 : MathVector<2> LocalGravity;
193 : MathMatrix<2,dim> JT;
194 0 : static ReferenceMapping<ReferenceTriangle, dim> TriangleMapping;
195 0 : TriangleMapping.update (vCorners);
196 :
197 : /* compute the local gravity */
198 : VecSet (LocalPoint, 0.0);
199 : TriangleMapping.jacobian_transposed (JT, LocalPoint);
200 : MatVecMult (LocalGravity, JT, PhysicalGravity);
201 :
202 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
203 0 : vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
204 :
205 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
206 0 : vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
207 0 : }
208 :
209 : /// computation of the primary function for the consistent gravity at corners of a quadrilateral
210 : template <int dim>
211 0 : inline void prepare_quadrilateral
212 : (
213 : MathVector<2>* vConsGravity, ///< where to save the values (4 vectors)
214 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (4 vectors)
215 : const number* vDensity, ///< corner density (4 scalar values)
216 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
217 : )
218 : {
219 : MathVector<2> LocalPoint;
220 : MathVector<2> LocalGravityAt000, LocalGravityAt110;
221 : MathMatrix<2,dim> JT;
222 0 : static ReferenceMapping<ReferenceQuadrilateral, dim> QuadMapping;
223 : QuadMapping.update (vCorners);
224 :
225 : /* compute the local gravity at local corner (0,0) */
226 : VecSet (LocalPoint, 0.0);
227 0 : QuadMapping.jacobian_transposed (JT, LocalPoint);
228 : MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
229 :
230 : /* compute the local gravity at local corner (1,1) */
231 : VecSet (LocalPoint, 1.0);
232 0 : QuadMapping.jacobian_transposed (JT, LocalPoint);
233 : MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
234 :
235 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
236 0 : vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
237 0 : vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
238 :
239 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
240 0 : vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
241 0 : vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
242 0 : }
243 :
244 : /// computation of the primary function for the consistent gravity at corners of a tetrahedron
245 : template <int dim>
246 0 : inline void prepare_tetrahedron
247 : (
248 : MathVector<3>* vConsGravity, ///< where to save the values (4 vectors)
249 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (4 vectors)
250 : const number* vDensity, ///< corner density (4 scalar values)
251 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
252 : )
253 : {
254 : MathVector<3> LocalPoint;
255 : MathVector<3> LocalGravity;
256 : MathMatrix<3,dim> JT;
257 0 : static ReferenceMapping<ReferenceTetrahedron, dim> TetMapping;
258 0 : TetMapping.update (vCorners);
259 :
260 : /* compute the local gravity */
261 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
262 0 : TetMapping.jacobian_transposed (JT, LocalPoint);
263 : MatVecMult (LocalGravity, JT, PhysicalGravity);
264 :
265 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
266 0 : vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
267 :
268 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[3][_Y_] = 0.0;
269 0 : vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
270 :
271 0 : vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
272 0 : vConsGravity[3][_Z_] = LocalGravity[_Z_]*(vDensity[0] + vDensity[3])*0.5;
273 0 : }
274 :
275 : /// computation of the primary function for the consistent gravity at corners of a pyramid
276 : /**
277 : * TODO: Verify this implementation! Cf. UG3.
278 : */
279 : template <int dim>
280 0 : inline void prepare_pyramid
281 : (
282 : MathVector<3>* vConsGravity, ///< where to save the values (5 vectors)
283 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (5 vectors)
284 : const number* vDensity, ///< corner density (5 scalar values)
285 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
286 : )
287 : {
288 : MathVector<3> LocalPoint;
289 : MathVector<3> LocalGravityAt000, LocalGravityAt110;
290 : MathMatrix<3,dim> JT;
291 0 : static ReferenceMapping<ReferencePyramid, dim> PyramidMapping;
292 : PyramidMapping.update (vCorners);
293 :
294 : /* compute the local gravity at local corner (0,0,0) */
295 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
296 0 : PyramidMapping.jacobian_transposed (JT, LocalPoint);
297 : MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
298 :
299 : /* compute the local gravity at local corner (1,1,0) */
300 0 : LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
301 0 : PyramidMapping.jacobian_transposed (JT, LocalPoint);
302 : MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
303 :
304 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0; vConsGravity[4][_X_] = 0.0;
305 0 : vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
306 0 : vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
307 :
308 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
309 0 : vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
310 0 : vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
311 :
312 0 : vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
313 0 : vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
314 0 : vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
315 0 : }
316 :
317 : /// computation of the primary function for the consistent gravity at corners of a prism
318 : template <int dim>
319 0 : inline void prepare_prism
320 : (
321 : MathVector<3>* vConsGravity, ///< where to save the values (6 vectors)
322 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (6 vectors)
323 : const number* vDensity, ///< corner density (6 scalar values)
324 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
325 : )
326 : {
327 : MathVector<3> LocalPoint;
328 : MathVector<3> LocalGravityAt000, LocalGravityAt101, LocalGravityAt011;
329 : MathMatrix<3,dim> JT;
330 0 : static ReferenceMapping<ReferencePrism, dim> PrismMapping;
331 : PrismMapping.update (vCorners);
332 :
333 : /* compute the local gravity at local corner (0,0,0) */
334 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
335 0 : PrismMapping.jacobian_transposed (JT, LocalPoint);
336 : MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
337 :
338 : /* compute the local gravity at local corner (1,0,1) */
339 0 : LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
340 0 : PrismMapping.jacobian_transposed (JT, LocalPoint);
341 : MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
342 :
343 : /* compute the local gravity at local corner (0,1,1) */
344 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
345 0 : PrismMapping.jacobian_transposed (JT, LocalPoint);
346 : MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
347 :
348 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
349 0 : vConsGravity[3][_X_] = 0.0; vConsGravity[5][_X_] = 0.0;
350 0 : vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
351 0 : vConsGravity[4][_X_] = LocalGravityAt011[_X_]*(vDensity[3] + vDensity[4])*0.5;
352 :
353 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
354 0 : vConsGravity[3][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
355 0 : vConsGravity[2][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[2])*0.5;
356 0 : vConsGravity[5][_Y_] = LocalGravityAt011[_Y_]*(vDensity[3] + vDensity[5])*0.5;
357 :
358 0 : vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
359 0 : vConsGravity[3][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[3])*0.5;
360 0 : vConsGravity[4][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[4])*0.5;
361 0 : vConsGravity[5][_Z_] = LocalGravityAt011[_Z_]*(vDensity[2] + vDensity[5])*0.5;
362 0 : }
363 :
364 : /// computation of the primary function for the consistent gravity at corners of a hexahedron
365 : template <int dim>
366 0 : inline void prepare_hexahedron
367 : (
368 : MathVector<3>* vConsGravity, ///< where to save the values (8 vectors)
369 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (8 vectors)
370 : const number* vDensity, ///< corner density (8 scalar values)
371 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
372 : )
373 : {
374 : MathVector<3> LocalPoint;
375 : MathVector<3> LocalGravityAt000, LocalGravityAt110, LocalGravityAt101, LocalGravityAt011;
376 : MathMatrix<3,dim> JT;
377 0 : static ReferenceMapping<ReferenceHexahedron, dim> HexMapping;
378 : HexMapping.update (vCorners);
379 :
380 : /* compute the local gravity at local corner (0,0,0) */
381 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
382 0 : HexMapping.jacobian_transposed (JT, LocalPoint);
383 : MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
384 :
385 : /* compute the local gravity at local corner (1,1,0) */
386 0 : LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
387 0 : HexMapping.jacobian_transposed (JT, LocalPoint);
388 : MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
389 :
390 : /* compute the local gravity at local corner (1,0,1) */
391 0 : LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
392 0 : HexMapping.jacobian_transposed (JT, LocalPoint);
393 : MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
394 :
395 : /* compute the local gravity at local corner (0,1,1) */
396 0 : LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
397 0 : HexMapping.jacobian_transposed (JT, LocalPoint);
398 : MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
399 :
400 0 : vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
401 0 : vConsGravity[4][_X_] = 0.0; vConsGravity[7][_X_] = 0.0;
402 0 : vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
403 0 : vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
404 0 : vConsGravity[5][_X_] = LocalGravityAt101[_X_]*(vDensity[4] + vDensity[5])*0.5;
405 0 : vConsGravity[6][_X_] = LocalGravityAt011[_X_]*(vDensity[6] + vDensity[7])*0.5;
406 :
407 0 : vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
408 0 : vConsGravity[4][_Y_] = 0.0; vConsGravity[5][_Y_] = 0.0;
409 0 : vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
410 0 : vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
411 0 : vConsGravity[6][_Y_] = LocalGravityAt101[_Y_]*(vDensity[5] + vDensity[6])*0.5;
412 0 : vConsGravity[7][_Y_] = LocalGravityAt011[_Y_]*(vDensity[4] + vDensity[7])*0.5;
413 :
414 0 : vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
415 0 : vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
416 0 : vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
417 0 : vConsGravity[5][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[5])*0.5;
418 0 : vConsGravity[6][_Z_] = LocalGravityAt110[_Z_]*(vDensity[2] + vDensity[6])*0.5;
419 0 : vConsGravity[7][_Z_] = LocalGravityAt011[_Z_]*(vDensity[3] + vDensity[7])*0.5;
420 0 : }
421 : };
422 :
423 : /// spacialization of the method for edges (reference dimension 1)
424 : template <>
425 : template <int dim>
426 : void StdLinConsistentGravity<1>::prepare
427 : (
428 : MathVector<1>* vConsGravity, ///< where to save the values (n_co vectors)
429 : const int n_co, ///< number of corners of the element (should be 2)
430 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
431 : const number* vDensity, ///< corner density (n_co scalar values)
432 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
433 : )
434 : {
435 : UG_ASSERT (n_co == 2, "StdLinConsistentGravity: Illegal number of corners of an edge.");
436 0 : m_nCo = n_co;
437 0 : this->template prepare_edge<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
438 : }
439 :
440 : /// spacialization of the method for faces (reference dimension 2)
441 : template <>
442 : template <int dim>
443 0 : void StdLinConsistentGravity<2>::prepare
444 : (
445 : MathVector<2>* vConsGravity, ///< where to save the values (n_co vectors)
446 : const int n_co, ///< number of corners of the element
447 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
448 : const number* vDensity, ///< corner density (n_co scalar values)
449 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
450 : )
451 : {
452 0 : switch (m_nCo = n_co)
453 : {
454 0 : case 3:
455 0 : this->template prepare_triangle<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
456 0 : break;
457 0 : case 4:
458 0 : this->template prepare_quadrilateral<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
459 0 : break;
460 0 : default:
461 0 : UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
462 : << n_co << ") of an element with reference dimension 2.");
463 : }
464 0 : }
465 :
466 : /// spacialization of the method for volumes (reference dimension 3)
467 : template <>
468 : template <int dim>
469 0 : void StdLinConsistentGravity<3>::prepare
470 : (
471 : MathVector<3>* vConsGravity, ///< where to save the values (n_co vectors)
472 : const int n_co, ///< number of corners of the element
473 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
474 : const number* vDensity, ///< corner density (n_co scalar values)
475 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
476 : )
477 : {
478 0 : switch (m_nCo = n_co)
479 : {
480 0 : case 4:
481 0 : this->template prepare_tetrahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
482 0 : break;
483 0 : case 5:
484 0 : this->template prepare_pyramid<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
485 0 : break;
486 0 : case 6:
487 0 : this->template prepare_prism<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
488 0 : break;
489 0 : case 8:
490 0 : this->template prepare_hexahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
491 0 : break;
492 0 : default:
493 0 : UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
494 : << n_co << ") of an element with reference dimension 3.");
495 : }
496 0 : }
497 :
498 : /// Class for the computation of the enhanced version ('Frolkovic-type') of the consistent gravity
499 : /**
500 : * Implementation of the enhanced ('Frolkovic-type') version of the consistent
501 : * gravity for simplices (triangles and tetrahedra) in the case of the gravity
502 : * force parallel to the z-axis. For all other elements and other gravities,
503 : * the same method as in StdLinConsistentGravity is used.
504 : * \see StdLinConsistentGravity
505 : *
506 : * References:
507 : * <ul>
508 : * <li> P. Frolkovic, Consistent velocity approximation for density driven
509 : * flow and transport, In: Advanced Computational Methods in Engineering,
510 : * Part 2: Contributed papers; (R. Van Keer at al., eds.), Shaker Publishing,
511 : * Maastricht, 1998, p. 603-611
512 : * </li>
513 : * </ul>
514 : *
515 : * \tparam refDim dimensionality of the reference element (e.g. 2 for triangles, 3 for tetrahedra)
516 : */
517 : template <int refDim>
518 : class StdLinConsistentGravityX : public StdLinConsistentGravity<refDim>
519 : {
520 : typedef StdLinConsistentGravity<refDim> base_type;
521 :
522 : public:
523 :
524 : /// constructor
525 : StdLinConsistentGravityX () {}
526 :
527 : /// computation of the primary function for the consistent gravity at corners, cf. the specializations
528 : template <int dim>
529 : inline void prepare
530 : (
531 : MathVector<refDim>* vConsGravity, ///< where to save the values (n_co vectors)
532 : const int n_co, ///< number of corners of the element
533 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
534 : const number* vDensity, ///< corner density (n_co scalar values)
535 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
536 : )
537 : {
538 : base_type::template prepare<dim>
539 : (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
540 : // C.f. the specializations below for the enhanced version
541 : }
542 :
543 : /// computation of the consistent gravity at a given point
544 : template <int dim>
545 : inline void compute
546 : (
547 : MathVector<dim>& ConsistentGravity, ///< where to save the vector
548 : const MathVector<refDim>& LocalCoord, ///< local coordinates of the point
549 : const MathMatrix<dim, refDim>& JTInv, ///< inverse transposed Jacobian
550 : const MathVector<refDim>* vLocalGrad, ///< gradients of the shape functions at the given point
551 : const MathVector<refDim>* vConsGravity ///< primary function of the consistent gravity at corners
552 : )
553 : {
554 : if (base_type::m_nCo > 0) // use the standard version
555 : base_type::template compute<dim> (ConsistentGravity, LocalCoord, JTInv, vLocalGrad, vConsGravity);
556 : else
557 : {
558 : // special method for triangles and tetrahedra, currently only in the full dimension
559 : UG_ASSERT (dim == refDim && base_type::m_nCo == -(refDim+1), "StdLinConsistentGravityX: Illegal initialization of the object.");
560 :
561 : MathVector<refDim> LocalGravity;
562 : VecSet (LocalGravity, 0.0);
563 :
564 : for (size_t d = 0; d < refDim; d++)
565 : LocalGravity[d] = vConsGravity[d+1][dim-1];
566 :
567 : // multiply by JacobianTransposedInverse
568 : MatVecMult(ConsistentGravity, JTInv, LocalGravity);
569 :
570 : // correct the first coordinates
571 : for (size_t d = 0; d < dim-1; d++)
572 : for (size_t sh = 1; sh <= refDim; sh++)
573 : ConsistentGravity[d] -= vConsGravity[sh][d] * LocalCoord[sh-1];
574 : }
575 : }
576 :
577 : protected:
578 :
579 : /// computation of the extended version of the corner consistent gravity for simplices (only in full dimension!)
580 : template <typename refElem, int dim>
581 : inline void prepare_simplex
582 : (
583 : MathVector<refDim>* vConsGravity, ///< where to save the values (dim+1 vectors)
584 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (dim+1 vectors)
585 : const number* vDensity, ///< corner density (dim+1 scalar values)
586 : const MathVector<dim>& PhysicalGravity ///< the gravity vector (MUST BE (0, ..., 0, g))
587 : )
588 : {
589 : number DensityIP, Diff, gradient;
590 : MathVector<refDim> LocalPoint;
591 : MathVector<dim> ShiftedGlobalPoint;
592 : MathMatrix<refDim,dim> JT;
593 : static ReferenceMapping<refElem, dim> Mapping;
594 : Mapping.update (vCorners);
595 : Mapping.jacobian_transposed_inverse (JT, vCorners[0]); // the 2. argument is dummy
596 :
597 : for (size_t i = 1; i < dim+1; i++)
598 : {
599 : VecSubtract (ShiftedGlobalPoint, vCorners[i], vCorners[0]);
600 : Diff = ShiftedGlobalPoint [dim-1];
601 :
602 : ShiftedGlobalPoint [dim-1] = 0.0;
603 : TransposedMatVecMult (LocalPoint, JT, ShiftedGlobalPoint);
604 :
605 : DensityIP = 1;
606 : for (size_t j = 0; j < dim; j++) DensityIP -= LocalPoint[j]; //GN(0)=1-local(0)-local(1)-...
607 : DensityIP *= vDensity[0];
608 : for (size_t j = 1; j < dim+1; j++)
609 : DensityIP += vDensity[j] * LocalPoint[j-1]; //GN(j)=local(j-1) for j=1,2,...
610 :
611 : for(size_t j = 0; j < dim-1; ++j)
612 : {
613 : gradient = 0.0;
614 : for (size_t k = 0; k < dim; k++)
615 : gradient += JT[j][k] * (vDensity[k+1] - vDensity[0]);
616 : vConsGravity[i][j] = Diff * PhysicalGravity[dim-1] * gradient;
617 : }
618 : vConsGravity[i][dim-1] = Diff * PhysicalGravity[dim-1] * (vDensity[i] + DensityIP)*0.5;
619 : }
620 : }
621 : };
622 :
623 : /// spacialization of the method for faces (reference dimension 2)
624 : template <>
625 : template <int dim>
626 : void StdLinConsistentGravityX<2>::prepare
627 : (
628 : MathVector<2>* vConsGravity, ///< where to save the values (n_co vectors)
629 : const int n_co, ///< number of corners of the element
630 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
631 : const number* vDensity, ///< corner density (n_co scalar values)
632 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
633 : )
634 : {
635 : if (dim == 2 && n_co == 3 && PhysicalGravity[0] == 0.0) // use the enhanced version
636 : {
637 : m_nCo = -3;
638 : this->template prepare_simplex<ReferenceTriangle, dim>
639 : (vConsGravity, vCorners, vDensity, PhysicalGravity);
640 : }
641 : else // use the standard version
642 : base_type::template prepare<dim>
643 : (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
644 : }
645 :
646 : /// spacialization of the method for volumes (reference dimension 3)
647 : template <>
648 : template <int dim>
649 : void StdLinConsistentGravityX<3>::prepare
650 : (
651 : MathVector<3>* vConsGravity, ///< where to save the values (n_co vectors)
652 : const int n_co, ///< number of corners of the element
653 : const MathVector<dim>* vCorners, ///< (global) coordinates of the corners (n_co vectors)
654 : const number* vDensity, ///< corner density (n_co scalar values)
655 : const MathVector<dim>& PhysicalGravity ///< the gravity vector
656 : )
657 : {
658 : if (dim == 3 && n_co == 4 && PhysicalGravity[0] == 0.0 && PhysicalGravity[1] == 0.0) // use the enhanced version
659 : {
660 : m_nCo = -4;
661 : this->template prepare_simplex<ReferenceTetrahedron, dim>
662 : (vConsGravity, vCorners, vDensity, PhysicalGravity);
663 : }
664 : else // use the standard version
665 : base_type::template prepare<dim>
666 : (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
667 : }
668 :
669 : } // end namespace ug
670 :
671 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONSISTENT_GRAVITY__ */
|
