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#include <vector>
37
38// other ug4 modules
39#include "common/common.h"
40#include "lib_disc/reference_element/reference_mapping.h"
41
42namespace ug{
43
45
98template <int refDim>
99class StdLinConsistentGravity
100{
101private:
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
107public:
108
110 StdLinConsistentGravity () : m_nCo (0) {};
111
113 template <int dim>
114 inline void prepare
115 (
116 MathVector<refDim>* vConsGravity,
117 const int n_co,
118 const MathVector<dim>* vCorners,
119 const number* vDensity,
120 const MathVector<dim>& PhysicalGravity
121 )
122 {
123 UG_THROW ("StdLinConsistentGravity: Combination of the world dim " << dim <<
124 "and the reference element dim " << refDim << " is not implemented.");
125 }
126
128 template <int dim>
129 inline void compute
130 (
131 MathVector<dim>& ConsistentGravity,
132 const MathVector<refDim>& LocalCoord,
133 const MathMatrix<dim, refDim>& JTInv,
134 const MathVector<refDim>* vLocalGrad,
135 const MathVector<refDim>* vConsGravity
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 for(size_t sh = 0; sh < (size_t) m_nCo; sh++)
145 for(size_t d = 0; d < refDim; d++)
146 LocalGravity[d] += vConsGravity[sh][d] * vLocalGrad[sh][d];
147
148 // Multiply by JacobianTransposedInverse
149 MatVecMult(ConsistentGravity, JTInv, LocalGravity);
150 }
151
152protected:
153
154 int m_nCo;
155
157 template <int dim>
158 inline void prepare_edge
159 (
160 MathVector<1>* vConsGravity,
161 const MathVector<dim>* vCorners,
162 const number* vDensity,
163 const MathVector<dim>& PhysicalGravity
164 )
165 {
166 MathVector<1> LocalPoint;
167 MathVector<1> LocalGravity;
168 MathMatrix<1,dim> JT;
169 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 vConsGravity[0][_X_] = 0.0;
178 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
179 }
180
182 template <int dim>
183 inline void prepare_triangle
184 (
185 MathVector<2>* vConsGravity,
186 const MathVector<dim>* vCorners,
187 const number* vDensity,
188 const MathVector<dim>& PhysicalGravity
189 )
190 {
191 MathVector<2> LocalPoint;
192 MathVector<2> LocalGravity;
193 MathMatrix<2,dim> JT;
194 static ReferenceMapping<ReferenceTriangle, dim> TriangleMapping;
195 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 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
203 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
204
205 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
206 vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
207 }
208
210 template <int dim>
211 inline void prepare_quadrilateral
212 (
213 MathVector<2>* vConsGravity,
214 const MathVector<dim>* vCorners,
215 const number* vDensity,
216 const MathVector<dim>& PhysicalGravity
217 )
218 {
219 MathVector<2> LocalPoint;
220 MathVector<2> LocalGravityAt000, LocalGravityAt110;
221 MathMatrix<2,dim> JT;
222 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 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 QuadMapping.jacobian_transposed (JT, LocalPoint);
233 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
234
235 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
236 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
237 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
238
239 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
240 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
241 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
242 }
243
245 template <int dim>
246 inline void prepare_tetrahedron
247 (
248 MathVector<3>* vConsGravity,
249 const MathVector<dim>* vCorners,
250 const number* vDensity,
251 const MathVector<dim>& PhysicalGravity
252 )
253 {
254 MathVector<3> LocalPoint;
255 MathVector<3> LocalGravity;
256 MathMatrix<3,dim> JT;
257 static ReferenceMapping<ReferenceTetrahedron, dim> TetMapping;
258 TetMapping.update (vCorners);
259
260 /* compute the local gravity */
261 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
262 TetMapping.jacobian_transposed (JT, LocalPoint);
263 MatVecMult (LocalGravity, JT, PhysicalGravity);
264
265 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
266 vConsGravity[1][_X_] = LocalGravity[_X_]*(vDensity[0] + vDensity[1])*0.5;
267
268 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[3][_Y_] = 0.0;
269 vConsGravity[2][_Y_] = LocalGravity[_Y_]*(vDensity[0] + vDensity[2])*0.5;
270
271 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
272 vConsGravity[3][_Z_] = LocalGravity[_Z_]*(vDensity[0] + vDensity[3])*0.5;
273 }
274
276
279 template <int dim>
280 inline void prepare_pyramid
281 (
282 MathVector<3>* vConsGravity,
283 const MathVector<dim>* vCorners,
284 const number* vDensity,
285 const MathVector<dim>& PhysicalGravity
286 )
287 {
288 MathVector<3> LocalPoint;
289 MathVector<3> LocalGravityAt000, LocalGravityAt110;
290 MathMatrix<3,dim> JT;
291 static ReferenceMapping<ReferencePyramid, dim> PyramidMapping;
292 PyramidMapping.update (vCorners);
293
294 /* compute the local gravity at local corner (0,0,0) */
295 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
296 PyramidMapping.jacobian_transposed (JT, LocalPoint);
297 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
298
299 /* compute the local gravity at local corner (1,1,0) */
300 LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
301 PyramidMapping.jacobian_transposed (JT, LocalPoint);
302 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
303
304 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0; vConsGravity[4][_X_] = 0.0;
305 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
306 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
307
308 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
309 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
310 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
311
312 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
313 vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
314 vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
315 }
316
318 template <int dim>
319 inline void prepare_prism
320 (
321 MathVector<3>* vConsGravity,
322 const MathVector<dim>* vCorners,
323 const number* vDensity,
324 const MathVector<dim>& PhysicalGravity
325 )
326 {
327 MathVector<3> LocalPoint;
328 MathVector<3> LocalGravityAt000, LocalGravityAt101, LocalGravityAt011;
329 MathMatrix<3,dim> JT;
330 static ReferenceMapping<ReferencePrism, dim> PrismMapping;
331 PrismMapping.update (vCorners);
332
333 /* compute the local gravity at local corner (0,0,0) */
334 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
335 PrismMapping.jacobian_transposed (JT, LocalPoint);
336 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
337
338 /* compute the local gravity at local corner (1,0,1) */
339 LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
340 PrismMapping.jacobian_transposed (JT, LocalPoint);
341 MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
342
343 /* compute the local gravity at local corner (0,1,1) */
344 LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
345 PrismMapping.jacobian_transposed (JT, LocalPoint);
346 MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
347
348 vConsGravity[0][_X_] = 0.0; vConsGravity[2][_X_] = 0.0;
349 vConsGravity[3][_X_] = 0.0; vConsGravity[5][_X_] = 0.0;
350 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
351 vConsGravity[4][_X_] = LocalGravityAt011[_X_]*(vDensity[3] + vDensity[4])*0.5;
352
353 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
354 vConsGravity[3][_Y_] = 0.0; vConsGravity[4][_Y_] = 0.0;
355 vConsGravity[2][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[2])*0.5;
356 vConsGravity[5][_Y_] = LocalGravityAt011[_Y_]*(vDensity[3] + vDensity[5])*0.5;
357
358 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0; vConsGravity[2][_Z_] = 0.0;
359 vConsGravity[3][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[3])*0.5;
360 vConsGravity[4][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[4])*0.5;
361 vConsGravity[5][_Z_] = LocalGravityAt011[_Z_]*(vDensity[2] + vDensity[5])*0.5;
362 }
363
365 template <int dim>
366 inline void prepare_hexahedron
367 (
368 MathVector<3>* vConsGravity,
369 const MathVector<dim>* vCorners,
370 const number* vDensity,
371 const MathVector<dim>& PhysicalGravity
372 )
373 {
374 MathVector<3> LocalPoint;
375 MathVector<3> LocalGravityAt000, LocalGravityAt110, LocalGravityAt101, LocalGravityAt011;
376 MathMatrix<3,dim> JT;
377 static ReferenceMapping<ReferenceHexahedron, dim> HexMapping;
378 HexMapping.update (vCorners);
379
380 /* compute the local gravity at local corner (0,0,0) */
381 LocalPoint.x() = 0.0; LocalPoint.y() = 0.0; LocalPoint.z() = 0.0;
382 HexMapping.jacobian_transposed (JT, LocalPoint);
383 MatVecMult (LocalGravityAt000, JT, PhysicalGravity);
384
385 /* compute the local gravity at local corner (1,1,0) */
386 LocalPoint.x() = 1.0; LocalPoint.y() = 1.0; LocalPoint.z() = 0.0;
387 HexMapping.jacobian_transposed (JT, LocalPoint);
388 MatVecMult (LocalGravityAt110, JT, PhysicalGravity);
389
390 /* compute the local gravity at local corner (1,0,1) */
391 LocalPoint.x() = 1.0; LocalPoint.y() = 0.0; LocalPoint.z() = 1.0;
392 HexMapping.jacobian_transposed (JT, LocalPoint);
393 MatVecMult (LocalGravityAt101, JT, PhysicalGravity);
394
395 /* compute the local gravity at local corner (0,1,1) */
396 LocalPoint.x() = 0.0; LocalPoint.y() = 1.0; LocalPoint.z() = 1.0;
397 HexMapping.jacobian_transposed (JT, LocalPoint);
398 MatVecMult (LocalGravityAt011, JT, PhysicalGravity);
399
400 vConsGravity[0][_X_] = 0.0; vConsGravity[3][_X_] = 0.0;
401 vConsGravity[4][_X_] = 0.0; vConsGravity[7][_X_] = 0.0;
402 vConsGravity[1][_X_] = LocalGravityAt000[_X_]*(vDensity[0] + vDensity[1])*0.5;
403 vConsGravity[2][_X_] = LocalGravityAt110[_X_]*(vDensity[2] + vDensity[3])*0.5;
404 vConsGravity[5][_X_] = LocalGravityAt101[_X_]*(vDensity[4] + vDensity[5])*0.5;
405 vConsGravity[6][_X_] = LocalGravityAt011[_X_]*(vDensity[6] + vDensity[7])*0.5;
406
407 vConsGravity[0][_Y_] = 0.0; vConsGravity[1][_Y_] = 0.0;
408 vConsGravity[4][_Y_] = 0.0; vConsGravity[5][_Y_] = 0.0;
409 vConsGravity[2][_Y_] = LocalGravityAt110[_Y_]*(vDensity[1] + vDensity[2])*0.5;
410 vConsGravity[3][_Y_] = LocalGravityAt000[_Y_]*(vDensity[0] + vDensity[3])*0.5;
411 vConsGravity[6][_Y_] = LocalGravityAt101[_Y_]*(vDensity[5] + vDensity[6])*0.5;
412 vConsGravity[7][_Y_] = LocalGravityAt011[_Y_]*(vDensity[4] + vDensity[7])*0.5;
413
414 vConsGravity[0][_Z_] = 0.0; vConsGravity[1][_Z_] = 0.0;
415 vConsGravity[2][_Z_] = 0.0; vConsGravity[3][_Z_] = 0.0;
416 vConsGravity[4][_Z_] = LocalGravityAt000[_Z_]*(vDensity[0] + vDensity[4])*0.5;
417 vConsGravity[5][_Z_] = LocalGravityAt101[_Z_]*(vDensity[1] + vDensity[5])*0.5;
418 vConsGravity[6][_Z_] = LocalGravityAt110[_Z_]*(vDensity[2] + vDensity[6])*0.5;
419 vConsGravity[7][_Z_] = LocalGravityAt011[_Z_]*(vDensity[3] + vDensity[7])*0.5;
420 }
421};
422
424template <>
425template <int dim>
426void StdLinConsistentGravity<1>::prepare
427(
428 MathVector<1>* vConsGravity,
429 const int n_co,
430 const MathVector<dim>* vCorners,
431 const number* vDensity,
432 const MathVector<dim>& PhysicalGravity
433)
434{
435 UG_ASSERT (n_co == 2, "StdLinConsistentGravity: Illegal number of corners of an edge.");
436 m_nCo = n_co;
437 this->template prepare_edge<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
438}
439
441template <>
442template <int dim>
443void StdLinConsistentGravity<2>::prepare
444(
445 MathVector<2>* vConsGravity,
446 const int n_co,
447 const MathVector<dim>* vCorners,
448 const number* vDensity,
449 const MathVector<dim>& PhysicalGravity
450)
451{
452 switch (m_nCo = n_co)
453 {
454 case 3:
455 this->template prepare_triangle<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
456 break;
457 case 4:
458 this->template prepare_quadrilateral<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
459 break;
460 default:
461 UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
462 << n_co << ") of an element with reference dimension 2.");
463 }
464}
465
467template <>
468template <int dim>
469void StdLinConsistentGravity<3>::prepare
470(
471 MathVector<3>* vConsGravity,
472 const int n_co,
473 const MathVector<dim>* vCorners,
474 const number* vDensity,
475 const MathVector<dim>& PhysicalGravity
476)
477{
478 switch (m_nCo = n_co)
479 {
480 case 4:
481 this->template prepare_tetrahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
482 break;
483 case 5:
484 this->template prepare_pyramid<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
485 break;
486 case 6:
487 this->template prepare_prism<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
488 break;
489 case 8:
490 this->template prepare_hexahedron<dim> (vConsGravity, vCorners, vDensity, PhysicalGravity);
491 break;
492 default:
493 UG_THROW ("StdLinConsistentGravity: Illegal number of corners ("
494 << n_co << ") of an element with reference dimension 3.");
495 }
496}
497
499
517template <int refDim>
518class StdLinConsistentGravityX : public StdLinConsistentGravity<refDim>
519{
520 typedef StdLinConsistentGravity<refDim> base_type;
521
522public:
523
525 StdLinConsistentGravityX () {}
526
528 template <int dim>
529 inline void prepare
530 (
531 MathVector<refDim>* vConsGravity,
532 const int n_co,
533 const MathVector<dim>* vCorners,
534 const number* vDensity,
535 const MathVector<dim>& PhysicalGravity
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
544 template <int dim>
545 inline void compute
546 (
547 MathVector<dim>& ConsistentGravity,
548 const MathVector<refDim>& LocalCoord,
549 const MathMatrix<dim, refDim>& JTInv,
550 const MathVector<refDim>* vLocalGrad,
551 const MathVector<refDim>* vConsGravity
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
577protected:
578
580 template <typename refElem, int dim>
581 inline void prepare_simplex
582 (
583 MathVector<refDim>* vConsGravity,
584 const MathVector<dim>* vCorners,
585 const number* vDensity,
586 const MathVector<dim>& PhysicalGravity
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
624template <>
625template <int dim>
626void StdLinConsistentGravityX<2>::prepare
627(
628 MathVector<2>* vConsGravity,
629 const int n_co,
630 const MathVector<dim>* vCorners,
631 const number* vDensity,
632 const MathVector<dim>& PhysicalGravity
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
647template <>
648template <int dim>
649void StdLinConsistentGravityX<3>::prepare
650(
651 MathVector<3>* vConsGravity,
652 const int n_co,
653 const MathVector<dim>* vCorners,
654 const number* vDensity,
655 const MathVector<dim>& PhysicalGravity
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__ */
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:100
ug::StdLinConsistentGravity::_Y_
static const size_t _Y_
Definition consistent_gravity.h:104
ug::StdLinConsistentGravity::_X_
static const size_t _X_
Definition consistent_gravity.h:103
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:184
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:115
ug::StdLinConsistentGravity::StdLinConsistentGravity
StdLinConsistentGravity()
constructor (sets the 'not init.' flag)
Definition consistent_gravity.h:110
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:212
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:159
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:281
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:130
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:154
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:247
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:320
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:367
ug::StdLinConsistentGravity::_Z_
static const size_t _Z_
Definition consistent_gravity.h:105
ug::StdLinConsistentGravityX
Class for the computation of the enhanced version ('Frolkovic-type') of the consistent gravity.
Definition consistent_gravity.h:519
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:530
ug::StdLinConsistentGravityX::StdLinConsistentGravityX
StdLinConsistentGravityX()
constructor
Definition consistent_gravity.h:525
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:546
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:582
ug::StdLinConsistentGravityX::base_type
StdLinConsistentGravity< refDim > base_type
Definition consistent_gravity.h:520
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
reference_mapping.h