ug4
Loading...
Searching...
No Matches
consistent_gravity.h
Go to the documentation of this file.
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// other ug4 modules
37#include "common/common.h"
38#include <vector>
39
40namespace ug{
41
43
96template <int refDim>
97class StdLinConsistentGravity
98{
99private:
100// static constants
101 static const size_t _X_ = 0;
102 static const size_t _Y_ = 1;
103 static const size_t _Z_ = 2;
104
105public:
106
108 StdLinConsistentGravity () : m_nCo (0) {};
109
111 template <int dim>
112 inline void prepare
113 (
114 MathVector<refDim>* vConsGravity,
115 const int n_co,
116 const MathVector<dim>* vCorners,
117 const number* vDensity,
118 const MathVector<dim>& PhysicalGravity
119 )
120 {
121 UG_THROW ("StdLinConsistentGravity: Combination of the world dim " << dim <<
122 "and the reference element dim " << refDim << " is not implemented.");
123 }
124
126 template <int dim>
127 inline void compute
128 (
129 MathVector<dim>& ConsistentGravity,
130 const MathVector<refDim>& LocalCoord,
131 const MathMatrix<dim, refDim>& JTInv,
132 const MathVector<refDim>* vLocalGrad,
133 const MathVector<refDim>* vConsGravity
134 )
135 {
136 UG_ASSERT (m_nCo > 0, "StdLinConsistentGravity: Object not initialized.");
137
138 MathVector<refDim> LocalGravity;
139 VecSet(LocalGravity, 0.0);
140
141 // Loop shape functions
142 for(size_t sh = 0; sh < (size_t) m_nCo; sh++)
143 for(size_t d = 0; d < refDim; d++)
144 LocalGravity[d] += vConsGravity[sh][d] * vLocalGrad[sh][d];
145
146 // Multiply by JacobianTransposedInverse
147 MatVecMult(ConsistentGravity, JTInv, LocalGravity);
148 }
149
150protected:
151
152 int m_nCo;
153
155 template <int dim>
156 inline void prepare_edge
157 (
158 MathVector<1>* vConsGravity,
159 const MathVector<dim>* vCorners,
160 const number* vDensity,
161 const MathVector<dim>& PhysicalGravity
162 )
163 {
164 MathVector<1> LocalPoint;
165 MathVector<1> LocalGravity;
166 MathMatrix<1,dim> JT;
167 static ReferenceMapping<ReferenceEdge, dim> EdgeMapping;
168 EdgeMapping.update (vCorners);
169
170 /* compute the local gravity */
171 VecSet (LocalPoint, 0.0);
172 EdgeMapping.jacobian_transposed (JT, LocalPoint);
173 MatVecMult (LocalGravity, JT, PhysicalGravity);
174
175 vConsGravity[0][_X_] = 0.0;
176 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
177 }
178
180 template <int dim>
181 inline void prepare_triangle
182 (
183 MathVector<2>* vConsGravity,
184 const MathVector<dim>* vCorners,
185 const number* vDensity,
186 const MathVector<dim>& PhysicalGravity
187 )
188 {
189 MathVector<2> LocalPoint;
190 MathVector<2> LocalGravity;
191 MathMatrix<2,dim> JT;
192 static ReferenceMapping<ReferenceTriangle, dim> TriangleMapping;
193 TriangleMapping.update (vCorners);
194
195 /* compute the local gravity */
196 VecSet (LocalPoint, 0.0);
197 TriangleMapping.jacobian_transposed (JT, LocalPoint);
198 MatVecMult (LocalGravity, JT, PhysicalGravity);
199
200 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
201 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
202
203 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
204 vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
205 }
206
208 template <int dim>
209 inline void prepare_quadrilateral
210 (
211 MathVector<2>* vConsGravity,
212 const MathVector<dim>* vCorners,
213 const number* vDensity,
214 const MathVector<dim>& PhysicalGravity
215 )
216 {
217 MathVector<2> LocalPoint;
218 MathVector<2> LocalGravityAt000, LocalGravityAt110;
219 MathMatrix<2,dim> JT;
220 static ReferenceMapping<ReferenceQuadrilateral, dim> QuadMapping;
221 QuadMapping.update (vCorners);
222
223 /* compute the local gravity at local corner (0,0) */
224 VecSet (LocalPoint, 0.0);
225 QuadMapping.jacobian_transposed (JT, LocalPoint);
226 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
227
228 /* compute the local gravity at local corner (1,1) */
229 VecSet (LocalPoint, 1.0);
230 QuadMapping.jacobian_transposed (JT, LocalPoint);
231 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
232
233 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
234 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
235 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
236
237 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
238 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
239 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
240 }
241
243 template <int dim>
244 inline void prepare_tetrahedron
245 (
246 MathVector<3>* vConsGravity,
247 const MathVector<dim>* vCorners,
248 const number* vDensity,
249 const MathVector<dim>& PhysicalGravity
250 )
251 {
252 MathVector<3> LocalPoint;
253 MathVector<3> LocalGravity;
254 MathMatrix<3,dim> JT;
255 static ReferenceMapping<ReferenceTetrahedron, dim> TetMapping;
256 TetMapping.update (vCorners);
257
258 /* compute the local gravity */
259 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
260 TetMapping.jacobian_transposed (JT, LocalPoint);
261 MatVecMult (LocalGravity, JT, PhysicalGravity);
262
263 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
264 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
265
266 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[3][_Y_] = 0.0;
267 vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
268
269 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
270 vConsGravity[3][_Z_] = LocalGravity[_Z_]*(vDensity[0] + vDensity[3])*0.5;
271 }
272
274
277 template <int dim>
278 inline void prepare_pyramid
279 (
280 MathVector<3>* vConsGravity,
281 const MathVector<dim>* vCorners,
282 const number* vDensity,
283 const MathVector<dim>& PhysicalGravity
284 )
285 {
286 MathVector<3> LocalPoint;
287 MathVector<3> LocalGravityAt000, LocalGravityAt110;
288 MathMatrix<3,dim> JT;
289 static ReferenceMapping<ReferencePyramid, dim> PyramidMapping;
290 PyramidMapping.update (vCorners);
291
292 /* compute the local gravity at local corner (0,0,0) */
293 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
294 PyramidMapping.jacobian_transposed (JT, LocalPoint);
295 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
296
297 /* compute the local gravity at local corner (1,1,0) */
298 LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
299 PyramidMapping.jacobian_transposed (JT, LocalPoint);
300 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
301
302 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0; vConsGravity[4][_X_] = 0.0;
303 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
304 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
305
306 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
307 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
308 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
309
310 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
311 vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
312 vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
313 }
314
316 template <int dim>
317 inline void prepare_prism
318 (
319 MathVector<3>* vConsGravity,
320 const MathVector<dim>* vCorners,
321 const number* vDensity,
322 const MathVector<dim>& PhysicalGravity
323 )
324 {
325 MathVector<3> LocalPoint;
326 MathVector<3> LocalGravityAt000, LocalGravityAt101, LocalGravityAt011;
327 MathMatrix<3,dim> JT;
328 static ReferenceMapping<ReferencePrism, dim> PrismMapping;
329 PrismMapping.update (vCorners);
330
331 /* compute the local gravity at local corner (0,0,0) */
332 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
333 PrismMapping.jacobian_transposed (JT, LocalPoint);
334 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
335
336 /* compute the local gravity at local corner (1,0,1) */
337 LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
338 PrismMapping.jacobian_transposed (JT, LocalPoint);
339 MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
340
341 /* compute the local gravity at local corner (0,1,1) */
342 LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
343 PrismMapping.jacobian_transposed (JT, LocalPoint);
344 MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
345
346 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
347 vConsGravity[3][_X_] = 0.0; vConsGravity[5][_X_] = 0.0;
348 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
349 vConsGravity[4][_X_] = LocalGravityAt011[_X_]*(vDensity[3] + vDensity[4])*0.5;
350
351 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
352 vConsGravity[3][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
353 vConsGravity[2][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[2])*0.5;
354 vConsGravity[5][_Y_] = LocalGravityAt011[_Y_]*(vDensity[3] + vDensity[5])*0.5;
355
356 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
357 vConsGravity[3][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[3])*0.5;
358 vConsGravity[4][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[4])*0.5;
359 vConsGravity[5][_Z_] = LocalGravityAt011[_Z_]*(vDensity[2] + vDensity[5])*0.5;
360 }
361
363 template <int dim>
364 inline void prepare_hexahedron
365 (
366 MathVector<3>* vConsGravity,
367 const MathVector<dim>* vCorners,
368 const number* vDensity,
369 const MathVector<dim>& PhysicalGravity
370 )
371 {
372 MathVector<3> LocalPoint;
373 MathVector<3> LocalGravityAt000, LocalGravityAt110, LocalGravityAt101, LocalGravityAt011;
374 MathMatrix<3,dim> JT;
375 static ReferenceMapping<ReferenceHexahedron, dim> HexMapping;
376 HexMapping.update (vCorners);
377
378 /* compute the local gravity at local corner (0,0,0) */
379 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
380 HexMapping.jacobian_transposed (JT, LocalPoint);
381 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
382
383 /* compute the local gravity at local corner (1,1,0) */
384 LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
385 HexMapping.jacobian_transposed (JT, LocalPoint);
386 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
387
388 /* compute the local gravity at local corner (1,0,1) */
389 LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
390 HexMapping.jacobian_transposed (JT, LocalPoint);
391 MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
392
393 /* compute the local gravity at local corner (0,1,1) */
394 LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
395 HexMapping.jacobian_transposed (JT, LocalPoint);
396 MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
397
398 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
399 vConsGravity[4][_X_] = 0.0; vConsGravity[7][_X_] = 0.0;
400 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
401 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
402 vConsGravity[5][_X_] = LocalGravityAt101[_X_]*(vDensity[4] + vDensity[5])*0.5;
403 vConsGravity[6][_X_] = LocalGravityAt011[_X_]*(vDensity[6] + vDensity[7])*0.5;
404
405 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
406 vConsGravity[4][_Y_] = 0.0; vConsGravity[5][_Y_] = 0.0;
407 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
408 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
409 vConsGravity[6][_Y_] = LocalGravityAt101[_Y_]*(vDensity[5] + vDensity[6])*0.5;
410 vConsGravity[7][_Y_] = LocalGravityAt011[_Y_]*(vDensity[4] + vDensity[7])*0.5;
411
412 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
413 vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
414 vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
415 vConsGravity[5][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[5])*0.5;
416 vConsGravity[6][_Z_] = LocalGravityAt110[_Z_]*(vDensity[2] + vDensity[6])*0.5;
417 vConsGravity[7][_Z_] = LocalGravityAt011[_Z_]*(vDensity[3] + vDensity[7])*0.5;
418 }
419};
420
422template <>
423template <int dim>
424void StdLinConsistentGravity<1>::prepare
425(
426 MathVector<1>* vConsGravity,
427 const int n_co,
428 const MathVector<dim>* vCorners,
429 const number* vDensity,
430 const MathVector<dim>& PhysicalGravity
431)
432{
433 UG_ASSERT (n_co == 2, "StdLinConsistentGravity: Illegal number of corners of an edge.");
434 m_nCo = n_co;
435 this->template prepare_edge<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
436}
437
439template <>
440template <int dim>
441void StdLinConsistentGravity<2>::prepare
442(
443 MathVector<2>* vConsGravity,
444 const int n_co,
445 const MathVector<dim>* vCorners,
446 const number* vDensity,
447 const MathVector<dim>& PhysicalGravity
448)
449{
450 switch (m_nCo = n_co)
451 {
452 case 3:
453 this->template prepare_triangle<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
454 break;
455 case 4:
456 this->template prepare_quadrilateral<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
457 break;
458 default:
459 UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
460 << n_co << ") of an element with reference dimension 2.");
461 }
462}
463
465template <>
466template <int dim>
467void StdLinConsistentGravity<3>::prepare
468(
469 MathVector<3>* vConsGravity,
470 const int n_co,
471 const MathVector<dim>* vCorners,
472 const number* vDensity,
473 const MathVector<dim>& PhysicalGravity
474)
475{
476 switch (m_nCo = n_co)
477 {
478 case 4:
479 this->template prepare_tetrahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
480 break;
481 case 5:
482 this->template prepare_pyramid<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
483 break;
484 case 6:
485 this->template prepare_prism<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
486 break;
487 case 8:
488 this->template prepare_hexahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
489 break;
490 default:
491 UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
492 << n_co << ") of an element with reference dimension 3.");
493 }
494}
495
497
515template <int refDim>
516class StdLinConsistentGravityX : public StdLinConsistentGravity<refDim>
517{
518 typedef StdLinConsistentGravity<refDim> base_type;
519
520public:
521
523 StdLinConsistentGravityX () {}
524
526 template <int dim>
527 inline void prepare
528 (
529 MathVector<refDim>* vConsGravity,
530 const int n_co,
531 const MathVector<dim>* vCorners,
532 const number* vDensity,
533 const MathVector<dim>& PhysicalGravity
534 )
535 {
536 base_type::template prepare<dim>
537 (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
538 // C.f. the specializations below for the enhanced version
539 }
540
542 template <int dim>
543 inline void compute
544 (
545 MathVector<dim>& ConsistentGravity,
546 const MathVector<refDim>& LocalCoord,
547 const MathMatrix<dim, refDim>& JTInv,
548 const MathVector<refDim>* vLocalGrad,
549 const MathVector<refDim>* vConsGravity
550 )
551 {
552 if (base_type::m_nCo > 0) // use the standard version
553 base_type::template compute<dim> (ConsistentGravity, LocalCoord, JTInv, vLocalGrad, vConsGravity);
554 else
555 {
556 // special method for triangles and tetrahedra, currently only in the full dimension
557 UG_ASSERT (dim == refDim && base_type::m_nCo == -(refDim+1), "StdLinConsistentGravityX: Illegal initialization of the object.");
558
559 MathVector<refDim> LocalGravity;
560 VecSet (LocalGravity, 0.0);
561
562 for (size_t d = 0; d < refDim; d++)
563 LocalGravity[d] = vConsGravity[d+1][dim-1];
564
565 // multiply by JacobianTransposedInverse
566 MatVecMult(ConsistentGravity, JTInv, LocalGravity);
567
568 // correct the first coordinates
569 for (size_t d = 0; d < dim-1; d++)
570 for (size_t sh = 1; sh <= refDim; sh++)
571 ConsistentGravity[d] -= vConsGravity[sh][d] * LocalCoord[sh-1];
572 }
573 }
574
575protected:
576
578 template <typename refElem, int dim>
579 inline void prepare_simplex
580 (
581 MathVector<refDim>* vConsGravity,
582 const MathVector<dim>* vCorners,
583 const number* vDensity,
584 const MathVector<dim>& PhysicalGravity
585 )
586 {
587 number DensityIP, Diff, gradient;
588 MathVector<refDim> LocalPoint;
589 MathVector<dim> ShiftedGlobalPoint;
590 MathMatrix<refDim,dim> JT;
591 static ReferenceMapping<refElem, dim> Mapping;
592 Mapping.update (vCorners);
593 Mapping.jacobian_transposed_inverse (JT, vCorners[0]); // the 2. argument is dummy
594
595 for (size_t i = 1; i < dim+1; i++)
596 {
597 VecSubtract (ShiftedGlobalPoint, vCorners[i], vCorners[0]);
598 Diff = ShiftedGlobalPoint [dim-1];
599
600 ShiftedGlobalPoint [dim-1] = 0.0;
601 TransposedMatVecMult (LocalPoint, JT, ShiftedGlobalPoint);
602
603 DensityIP = 1;
604 for (size_t j = 0; j < dim; j++) DensityIP -= LocalPoint[j]; //GN(0)=1-local(0)-local(1)-...
605 DensityIP *= vDensity[0];
606 for (size_t j = 1; j < dim+1; j++)
607 DensityIP += vDensity[j] * LocalPoint[j-1]; //GN(j)=local(j-1) for j=1,2,...
608
609 for(size_t j = 0; j < dim-1; ++j)
610 {
611 gradient = 0.0;
612 for (size_t k = 0; k < dim; k++)
613 gradient += JT[j][k] * (vDensity[k+1] - vDensity[0]);
614 vConsGravity[i][j] = Diff * PhysicalGravity[dim-1] * gradient;
615 }
616 vConsGravity[i][dim-1] = Diff * PhysicalGravity[dim-1] * (vDensity[i] + DensityIP)*0.5;
617 }
618 }
619};
620
622template <>
623template <int dim>
624void StdLinConsistentGravityX<2>::prepare
625(
626 MathVector<2>* vConsGravity,
627 const int n_co,
628 const MathVector<dim>* vCorners,
629 const number* vDensity,
630 const MathVector<dim>& PhysicalGravity
631)
632{
633 if (dim == 2 && n_co == 3 && PhysicalGravity[0] == 0.0) // use the enhanced version
634 {
635 m_nCo = -3;
636 this->template prepare_simplex<ReferenceTriangle, dim>
637 (vConsGravity, vCorners, vDensity, PhysicalGravity);
638 }
639 else // use the standard version
640 base_type::template prepare<dim>
641 (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
642}
643
645template <>
646template <int dim>
647void StdLinConsistentGravityX<3>::prepare
648(
649 MathVector<3>* vConsGravity,
650 const int n_co,
651 const MathVector<dim>* vCorners,
652 const number* vDensity,
653 const MathVector<dim>& PhysicalGravity
654)
655{
656 if (dim == 3 && n_co == 4 && PhysicalGravity[0] == 0.0 && PhysicalGravity[1] == 0.0) // use the enhanced version
657 {
658 m_nCo = -4;
659 this->template prepare_simplex<ReferenceTetrahedron, dim>
660 (vConsGravity, vCorners, vDensity, PhysicalGravity);
661 }
662 else // use the standard version
663 base_type::template prepare<dim>
664 (vConsGravity, n_co, vCorners, vDensity, PhysicalGravity);
665}
666
667} // end namespace ug
668
669#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONSISTENT_GRAVITY__ */
ug::MathMatrix
A class for fixed size, dense matrices.
Definition math_matrix.h:63
ug::MathVector
a mathematical Vector with N entries.
Definition math_vector.h:97
ug::ReferenceMapping
Definition reference_mapping.h:65
ug::ReferenceMapping::jacobian_transposed_inverse
number jacobian_transposed_inverse(MathMatrix< worldDim, dim > &JTInv, const MathVector< dim > &locPos) const
returns transposed of the inverse of the jacobian and sqrt of gram determinante
ug::ReferenceMapping::update
void update(const MathVector< worldDim > *vCornerCoord)
refresh mapping for new set of corners
ug::ReferenceMapping::jacobian_transposed
void jacobian_transposed(MathMatrix< dim, worldDim > &JT, const MathVector< dim > &locPos) const
returns transposed of jacobian
ug::StdLinConsistentGravity
Class for the computation of the standard version ('Voss-Souza-type') of the consistent gravity.
Definition consistent_gravity.h:98
ug::StdLinConsistentGravity::_Y_
static const size_t _Y_
Definition consistent_gravity.h:102
ug::StdLinConsistentGravity::_X_
static const size_t _X_
Definition consistent_gravity.h:101
ug::StdLinConsistentGravity::prepare_triangle
void prepare_triangle(MathVector< 2 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a triangle
Definition consistent_gravity.h:182
ug::StdLinConsistentGravity::prepare
void prepare(MathVector< refDim > *vConsGravity, const int n_co, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners, cf. the specializations
Definition consistent_gravity.h:113
ug::StdLinConsistentGravity::StdLinConsistentGravity
StdLinConsistentGravity()
constructor (sets the 'not init.' flag)
Definition consistent_gravity.h:108
ug::StdLinConsistentGravity::prepare_quadrilateral
void prepare_quadrilateral(MathVector< 2 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a quadrilateral
Definition consistent_gravity.h:210
ug::StdLinConsistentGravity::prepare_edge
void prepare_edge(MathVector< 1 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of an edge
Definition consistent_gravity.h:157
ug::StdLinConsistentGravity::prepare_pyramid
void prepare_pyramid(MathVector< 3 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a pyramid
Definition consistent_gravity.h:279
ug::StdLinConsistentGravity::compute
void compute(MathVector< dim > &ConsistentGravity, const MathVector< refDim > &LocalCoord, const MathMatrix< dim, refDim > &JTInv, const MathVector< refDim > *vLocalGrad, const MathVector< refDim > *vConsGravity)
computation of the consistent gravity at a given point
Definition consistent_gravity.h:128
ug::StdLinConsistentGravity::m_nCo
int m_nCo
number of corners of the element for which the object is init. (0 if not init)
Definition consistent_gravity.h:152
ug::StdLinConsistentGravity::prepare_tetrahedron
void prepare_tetrahedron(MathVector< 3 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a tetrahedron
Definition consistent_gravity.h:245
ug::StdLinConsistentGravity::prepare_prism
void prepare_prism(MathVector< 3 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a prism
Definition consistent_gravity.h:318
ug::StdLinConsistentGravity::prepare_hexahedron
void prepare_hexahedron(MathVector< 3 > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners of a hexahedron
Definition consistent_gravity.h:365
ug::StdLinConsistentGravity::_Z_
static const size_t _Z_
Definition consistent_gravity.h:103
ug::StdLinConsistentGravityX
Class for the computation of the enhanced version ('Frolkovic-type') of the consistent gravity.
Definition consistent_gravity.h:517
ug::StdLinConsistentGravityX::prepare
void prepare(MathVector< refDim > *vConsGravity, const int n_co, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the primary function for the consistent gravity at corners, cf. the specializations
Definition consistent_gravity.h:528
ug::StdLinConsistentGravityX::StdLinConsistentGravityX
StdLinConsistentGravityX()
constructor
Definition consistent_gravity.h:523
ug::StdLinConsistentGravityX::compute
void compute(MathVector< dim > &ConsistentGravity, const MathVector< refDim > &LocalCoord, const MathMatrix< dim, refDim > &JTInv, const MathVector< refDim > *vLocalGrad, const MathVector< refDim > *vConsGravity)
computation of the consistent gravity at a given point
Definition consistent_gravity.h:544
ug::StdLinConsistentGravityX::prepare_simplex
void prepare_simplex(MathVector< refDim > *vConsGravity, const MathVector< dim > *vCorners, const number *vDensity, const MathVector< dim > &PhysicalGravity)
computation of the extended version of the corner consistent gravity for simplices (only in full dime...
Definition consistent_gravity.h:580
ug::StdLinConsistentGravityX::base_type
StdLinConsistentGravity< refDim > base_type
Definition consistent_gravity.h:518
UG_ASSERT
#define UG_ASSERT(expr, msg)
Definition assert.h:70
UG_THROW
#define UG_THROW(msg)
Definition error.h:57
number
double number
Definition types.h:124
ug::TransposedMatVecMult
void TransposedMatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Transposed Matrix - Vector Muliplication.
Definition math_matrix_vector_functions_common_impl.hpp:111
ug::MatVecMult
void MatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication.
Definition math_matrix_vector_functions_common_impl.hpp:49
ug::VecSet
void VecSet(vector_t &vInOut, typename vector_t::value_type s)
Set each vector component to scalar (componentwise)
Definition math_vector_functions_common_impl.hpp:539
ug::VecSubtract
void VecSubtract(vector_t &vOut, const vector_t &v1, const vector_t &v2)
subtracts v2 from v1 and stores the result in a vOut
Definition math_vector_functions_common_impl.hpp:226
ug
the ug namespace