ug4
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 
40 namespace ug{
41 
43 
96 template <int refDim>
97 class StdLinConsistentGravity
98 {
99 private:
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 
105 public:
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 
150 protected:
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 
422 template <>
423 template <int dim>
424 void 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 
439 template <>
440 template <int dim>
441 void 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 
465 template <>
466 template <int dim>
467 void 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 
515 template <int refDim>
516 class StdLinConsistentGravityX : public StdLinConsistentGravity<refDim>
517 {
518  typedef StdLinConsistentGravity<refDim> base_type;
519 
520 public:
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 
575 protected:
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 
622 template <>
623 template <int dim>
624 void 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 
645 template <>
646 template <int dim>
647 void 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:52
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
dim
static const int dim
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