Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
conv_shape.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2011-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__CONV_SHAPE__
34#define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE__
35
36#include <boost/mpl/range_c.hpp>
37#include <boost/mpl/for_each.hpp>
38
39// ug4 headers
43
44namespace ug{
45
47// No Upwind
49
50template <int TDim>
52 : public IConvectionShapes<TDim>
53{
54 public:
57
60
62 static const int dim = TDim;
63
64 protected:
65 // explicitly forward some function
68 using base_type::D_vel;
72
73 public:
76 {
77 // the shapes do not depend on the DiffDisp. Thus, we can set the
78 // derivative to be always zero w.r.t. the DiffDisp for all shapes
80
81 // register evaluation function
82 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
83 boost::mpl::for_each< boost::mpl::range_c<int, 1, dim+1> >(RegisterRefDimFunc(this));
84 }
85
87 template <typename TFVGeom>
88 bool update(const TFVGeom* geo,
89 const MathVector<dim>* Velocity,
90 const MathMatrix<dim, dim>* DiffDisp,
91 bool computeDeriv);
92
93 private:
94
97 {
99 RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
100 template<typename TElem> void operator() (TElem &)
101 {
102 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
103 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
104 m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
105 }
106 };
107
109 template <typename TElem, typename TFVGeom>
111 {
112 typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
113 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
114 }
115
118 {
121 template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
122 };
123
125 template <int refDim>
127 {
128 typedef DimFV1Geometry<refDim, dim> TGeom;
129 typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
130 base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
131 }
132};
133
134template <int TDim>
135template <typename TFVGeom>
136bool
138update(const TFVGeom* geo,
139 const MathVector<dim>* Velocity,
140 const MathMatrix<dim, dim>* DiffDisp,
141 bool computeDeriv)
142{
143 UG_ASSERT(geo != nullptr, "Null pointer");
144 UG_ASSERT(Velocity != nullptr, "Null pointer");
145
146// \todo: think about: this should be something like scvf.num_sh()
147 const size_t numSH = geo->num_sh();
148
149// loop subcontrol volume faces
150 for(size_t ip = 0; ip < geo->num_scvf(); ++ip)
151 {
152 // get subcontrol volume face
153 const typename TFVGeom::SCVF& scvf = geo->scvf(ip);
154
155 // Compute flux
156 const number flux = VecDot(scvf.normal(), Velocity[ip]);
157
158 // Write Shapes
159 for(size_t sh = 0; sh < scvf.num_sh(); sh++)
160 conv_shape(ip, sh) = flux * scvf.shape(sh);
161
162 // this is introduced here, hopefully temporarily: The problem is, that
163 // for hanging nodes the number of shape function is not the number of
164 // corners, but scvf.num_sh() currently returns the number of corners.
165 // this is actually enough to interpolate the function, but still we
166 // should reset the interpolation adding for hanging dofs to zero
167 for(size_t sh = scvf.num_sh(); sh < numSH; sh++)
168 conv_shape(ip, sh) = 0.0;
169
170 // Write Derivatives if wanted
171 if(computeDeriv){
172 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
173 VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
174
175 // temporary, see above
176 for(size_t sh = scvf.num_sh(); sh < numSH; sh++)
177 VecSet(D_vel(ip, sh), 0.0);
178 }
179
180 // The shapes do not depend of the diffusion tensor
181 }
182
183// we're done
184 return true;
185}
186
188// Full Upwind
190
191template <int TDim>
193 : public IConvectionShapes<TDim>
194{
195 public:
198
201
203 static const int dim = TDim;
204
205 protected:
206 // explicitly forward some function
209 using base_type::D_vel;
213
214 public:
217 {
218 // the shapes do not depend on the DiffDisp. Thus, we can set the
219 // derivative to be always zero w.r.t. the DiffDisp for all shapes
221
222 // register evaluation function
223 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
224 boost::mpl::for_each< boost::mpl::range_c<int, 1, dim+1> >(RegisterRefDimFunc(this));
225 }
226
228 template <typename TFVGeom>
229 bool update(const TFVGeom* geo,
230 const MathVector<dim>* Velocity,
231 const MathMatrix<dim, dim>* DiffDisp,
232 bool computeDeriv);
233
234 private:
235
238 {
240 RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
241 template<typename TElem> void operator() (TElem &)
242 {
243 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
244 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
245 m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
246 }
247 };
248
250 template <typename TElem, typename TFVGeom>
252 {
253 typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
254 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
255 }
256
259 {
262 template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
263 };
264
266 template <int refDim>
268 {
269 typedef DimFV1Geometry<refDim, dim> TGeom;
270 typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
271 base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
272 }
273};
274
275template <int TDim>
276template <typename TFVGeom>
278update(const TFVGeom* geo,
279 const MathVector<dim>* Velocity,
280 const MathMatrix<dim, dim>* DiffDisp,
281 bool computeDeriv)
282{
283 UG_ASSERT(geo != nullptr, "Null pointer");
284 UG_ASSERT(Velocity != nullptr, "Null pointer");
285
286// \todo: think about: this should be something like scvf.num_sh()
287 const size_t numSH = geo->num_sh();
288
289// loop subcontrol volume faces
290 for(size_t ip = 0; ip < geo->num_scvf(); ++ip)
291 {
292 // get subcontrol volume face
293 const typename TFVGeom::SCVF& scvf = geo->scvf(ip);
294
295 // Compute flux
296 const number flux = VecDot(scvf.normal(), Velocity[ip]);
297
298 size_t from = scvf.from();
299 size_t to = scvf.to();
300
301 // if e.g. hanging nodes are involved, no upwind can be performed...
302 if((from >= scvf.num_sh()) || (to >= scvf.num_sh())){
303 // No upwind...
304 // Write Shapes
305 for(size_t sh = 0; sh < scvf.num_sh(); sh++)
306 conv_shape(ip, sh) = flux * scvf.shape(sh);
307
308 UG_ASSERT(scvf.num_sh() == numSH, "sh's have to match!");
309
310 // Write Derivatives if wanted
311 if(computeDeriv){
312 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
313 VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
314 }
315
316 continue;
317 }
318
319 // Full upwind below...
320 // Choose Upwind corner
321 size_t up = (flux >= 0) ? from : to;
322
323 // Write Shapes
324 for(size_t sh = 0; sh < scvf.num_sh(); ++sh) conv_shape(ip, sh) = 0.0;
325
326 conv_shape(ip, up) = flux;
327
328 // Write Derivatives if wanted
329 if(computeDeriv)
330 {
331 for(size_t sh = 0; sh < numSH; ++sh) VecSet(D_vel(ip, sh), 0.0);
332 D_vel(ip, up) = scvf.normal();
333 }
334
335 // The shapes do not depend of the diffusion tensor
336 }
337
338// we're done
339 return true;
340}
341
343// Weighted Upwind
345
346template <int TDim>
348 : public IConvectionShapes<TDim>
349{
350 public:
353
356
358 static const int dim = TDim;
359
360 protected:
361 // explicitly forward some function
364 using base_type::D_vel;
368
369 public:
372 {
373 // the shapes do not depend on the DiffDisp. Thus, we can set the
374 // derivative to be always zero w.r.t. the DiffDisp for all shapes
376
377 // register evaluation function
378 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
379 boost::mpl::for_each< boost::mpl::range_c<int, 1, dim+1> >(RegisterRefDimFunc(this));
380 }
381
384 {
385 set_weight(weight);
386
387 // the shapes do not depend on the DiffDisp. Thus, we can set the
388 // derivative to be always zero w.r.t. the DiffDisp for all shapes
390
391 // register evaluation function
392 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
393 boost::mpl::for_each< boost::mpl::range_c<int, 1, dim+1> >(RegisterRefDimFunc(this));
394 }
395
397 void set_weight(number weight) {m_weight = weight;}
398
400 template <typename TFVGeom>
401 bool update(const TFVGeom* geo,
402 const MathVector<dim>* Velocity,
403 const MathMatrix<dim, dim>* DiffDisp,
404 bool computeDeriv);
405
406 private:
407 // weight between no and full upwind (1.0 -> full upwind, 0.0 -> no upwind)
409
410 private:
411
414 {
416 RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
417 template<typename TElem> void operator() (TElem &)
418 {
419 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
420 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
421 m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
422 }
423 };
424
426 template <typename TElem, typename TFVGeom>
428 {
429 typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
430 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
431 }
432
435 {
438 template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
439 };
440
442 template <int refDim>
444 {
445 typedef DimFV1Geometry<refDim, dim> TGeom;
446 typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
447 base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
448 }
449};
450
451template <int TDim>
452template <typename TFVGeom>
454update(const TFVGeom* geo,
455 const MathVector<dim>* Velocity,
456 const MathMatrix<dim, dim>* DiffDisp,
457 bool computeDeriv)
458{
459 UG_ASSERT(geo != nullptr, "Null pointer");
460 UG_ASSERT(Velocity != nullptr, "Null pointer");
461
462// \todo: think about: this should be something like scvf.num_sh()
463 const size_t numSH = geo->num_sh();
464
465// loop subcontrol volume faces
466 for(size_t ip = 0; ip < geo->num_scvf(); ++ip)
467 {
468 // get subcontrol volume face
469 const typename TFVGeom::SCVF& scvf = geo->scvf(ip);
470
471 // Compute flux
472 const number flux = VecDot(scvf.normal(), Velocity[ip]);
473
474 size_t from = scvf.from();
475 size_t to = scvf.to();
476
477 // if e.g. hanging nodes are involved, no upwind can be performed...
478 if((from >= scvf.num_sh()) || (to >= scvf.num_sh())){
479 // No upwind...
480 // Write Shapes
481 for(size_t sh = 0; sh < scvf.num_sh(); sh++)
482 conv_shape(ip, sh) = flux * scvf.shape(sh);
483
484 UG_ASSERT(scvf.num_sh() == numSH, "sh's have to match!");
485
486 // Write Derivatives if wanted
487 if(computeDeriv){
488 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
489 VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
490 }
491
492 continue;
493 }
494
495 // Choose Upwind corner
496 size_t up = (flux >= 0) ? from : to;
497
498 // write no upwind part of shapes
499 const number noUpFlux = (1.-m_weight)*flux;
500 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
501 conv_shape(ip, sh) = noUpFlux * scvf.shape(sh);
502
503 // add full upwind part of shapes
504 conv_shape(ip, up) += m_weight * flux;
505
506 // Write Derivatives if wanted
507 if(computeDeriv)
508 {
509 // write no upwind part of derivatives
510 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
511 VecScale(D_vel(ip, sh), scvf.normal(),
512 (1.-m_weight)*scvf.shape(sh));
513 // see comment above
514 for(size_t sh = scvf.num_sh(); sh < numSH; sh++)
515 VecSet(D_vel(ip, sh), 0.0);
516
517 // add full upwind part of derivatives
518 VecScaleAppend(D_vel(ip, up), m_weight, scvf.normal());
519 }
520
521 // The shapes do not depend of the diffusion tensor
522 }
523
524// we're done
525 return true;
526}
527
528
530// Partial Upwind
532
533template <int TDim>
535 : public IConvectionShapes<TDim>
536{
537 public:
540
543
545 static const int dim = TDim;
546
547 protected:
548 // explicitly forward some function
551 using base_type::D_vel;
555
556 public:
559 {
560 // register evaluation function
561 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
562 boost::mpl::for_each< boost::mpl::range_c<int, 1, dim+1> >(RegisterRefDimFunc(this));
563 }
564
566 template <typename TFVGeom>
567 bool update(const TFVGeom* geo,
568 const MathVector<dim>* Velocity,
569 const MathMatrix<dim, dim>* DiffDisp,
570 bool computeDeriv);
571
572 private:
573
576 {
578 RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
579 template<typename TElem> void operator() (TElem &)
580 {
581 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
582 m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
583 }
584 };
585
587 template <typename TElem, typename TFVGeom>
589 {
590 typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
591 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
592 }
593
596 {
599 template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
600 };
601
603 template <int refDim>
605 {
606 typedef DimFV1Geometry<refDim, dim> TGeom;
607 typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
608 base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
609 }
610};
611
612template <int TDim>
613template <typename TFVGeom>
614bool
616update(const TFVGeom* geo,
617 const MathVector<dim>* Velocity,
618 const MathMatrix<dim, dim>* DiffDisp,
619 bool computeDeriv)
620{
621 UG_ASSERT(geo != nullptr, "Null pointer");
622 UG_ASSERT(Velocity != nullptr, "Null pointer");
623// UG_ASSERT(DiffDisp != nullptr, "Null pointer");
624
625// Compute Volume of Element
626// typedef typename TFVGeom::ref_elem_type ref_elem_type;
627 const number vol = ElementSize<dim>(geo->roid(), geo->corners());
628
629// loop subcontrol volume faces
630 for(size_t i = 0; i < geo->num_scvf(); ++i)
631 {
632 // get subcontrol volume face
633 const typename TFVGeom::SCVF& scvf = geo->scvf(i);
634
635 // get corners
636 const size_t from = scvf.from();
637 const size_t to = scvf.to();
638
639 // get gradients
640 const MathVector<dim>& gradTo = scvf.global_grad(to);
641 const MathVector<dim>& gradFrom = scvf.global_grad(from);
642
643 // set lambda negative as default
644 number lambda = -1;
645
646 // if DiffDisp-Tensor passed, compute lambda
647 if(DiffDisp != nullptr)
648 {
649 // Get Gradients
650 MathVector<dim> DiffGrad;
651
652 // Compute DiffGrad = D * Grad Phi_to
653 MatVecMult(DiffGrad, DiffDisp[i], gradTo);
654
655 // Compute GradDiffGrad = < Grad Phi_from, DiffGrad >
656 const number GradDiffGrad = VecDot(DiffGrad, gradFrom);
657
658 // Set lambda
659 lambda = - GradDiffGrad * vol;
660 }
661
662 // Compute flux
663 const number flux = VecDot(scvf.normal(), Velocity[i]);
664
665 // reset values
666 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
667 conv_shape(i, sh) = 0.0;
668 if(computeDeriv)
669 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
670 VecSet(D_vel(i, sh), 0.0);
671
672
673 // check: Currently hanging nodes not supported.
674 // /todo: add support for hanging nodes
675 if(from >= scvf.num_sh() || to >= scvf.num_sh())
676 UG_THROW("PartialUpwind: Currently not implemented for hanging nodes.")
677
678
679 // Case 1:
680 // full upwind is used
682 if(lambda <= 0 || DiffDisp == nullptr)
683 {
684 // Choose Upwind corner
685 const size_t up = (flux >= 0) ? scvf.from() : scvf.to();
686
687 // Write Shapes
688 conv_shape(i, up) = flux;
689
690 // Write Derivatives if wanted
691 if(computeDeriv)
692 {
693 // set derivative
694 D_vel(i, up) = scvf.normal();
695
696 // does not depend on diffusion
697 set_non_zero_deriv_diffusion_flag(false);
698 }
699
700 // everything done
701 continue;
702 }
703
705 // Case 2:
706 // The case of the diffusion dominance (central differences)
708 if (2 * lambda > fabs(flux))
709 {
710 conv_shape(i, from) = flux / 2.0;
711 conv_shape(i, to) = flux / 2.0;
712
713 if(computeDeriv)
714 {
715 set_non_zero_deriv_diffusion_flag(false);
716
717 VecScale(D_vel(i,from), scvf.normal(), 1.0/2.0);
718 VecScale(D_vel(i, to), scvf.normal(), 1.0/2.0);
719 }
720
721 // everything done
722 continue;
723 }
724
726 // Case 3:
727 // The cases of the convection dominance
729 set_non_zero_deriv_diffusion_flag(true);
730 if (flux >= 0)
731 {
732 conv_shape(i, from) = flux - lambda;
733 conv_shape(i, to) = lambda;
734
735 if(computeDeriv)
736 D_vel(i,from) = scvf.normal();
737 }
738 else
739 {
740 conv_shape(i, from) = - lambda;
741 conv_shape(i, to) = flux + lambda;
742
743 if(computeDeriv)
744 D_vel(i,to) = scvf.normal();
745 }
746
747 if (computeDeriv)
748 {
749 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
750 MatSet(conv_shape_diffusion(i, sh), 0.0);
751
752 for (size_t k = 0; k < (size_t)dim; k++)
753 for (size_t l = 0; l < (size_t)dim; l++)
754 {
755 conv_shape_diffusion(i, from)(k,l) = gradFrom[k]*gradTo[l]*vol;
756 conv_shape_diffusion(i, to)(k,l) = - gradFrom[k]*gradTo[l]*vol;
757 }
758 }
759 }
760
761// we're done
762 return true;
763}
764
766// Skewed Upwind
768
769template <int TDim>
771 : public IConvectionShapes<TDim>
772{
773public:
776
779
781 static const int dim = TDim;
782
783protected:
784 // explicitly forward some function
787 using base_type::D_vel;
791
792public:
795 {
796 // the shapes do not depend on the DiffDisp. Thus, we can set the
797 // derivative to be always zero w.r.t. the DiffDisp for all shapes
799
800 // register evaluation function
801 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
802 boost::mpl::for_each<boost::mpl::range_c<int, 1, dim + 1>>(RegisterRefDimFunc(this));
803 }
804
806 template <typename TFVGeom>
807 bool update(const TFVGeom *geo,
808 const MathVector<dim> *Velocity,
809 const MathMatrix<dim, dim> *DiffDisp,
810 bool computeDeriv);
811
812private:
815 {
818 template <typename TElem>
819 void operator()(TElem &)
820 {
821 m_pThis->template register_func_for_elem_fvgeom<TElem, FV1Geometry<TElem, dim>>();
822 //m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
823 //m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
824 }
825 };
826
828 template <typename TElem, typename TFVGeom>
830 {
831 typedef bool (this_type::*TFunc)(const TFVGeom *, const MathVector<dim> *, const MathMatrix<dim, dim> *, bool);
832 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
833 }
834
837 {
840 template <typename TRefDim>
841 void operator()(TRefDim &) { m_pThis->register_func_for_refDim<TRefDim::value>(); }
842 };
843
845 template <int refDim>
847 {
848 //typedef DimFV1Geometry<refDim, dim> TGeom;
849 //typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
850 //base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
851 }
852};
853
855template <typename TRefElem, int TWorldDim>
856void GetNodeNextToCut(size_t &coOut,
857 const MathVector<TWorldDim> &IP,
858 const MathVector<TWorldDim> &IPVel,
859 const MathVector<TWorldDim> *vCornerCoords)
860{
861 // help variables
862 size_t side = 0;
863 MathVector<TWorldDim> globalIntersection;
864 MathVector<TRefElem::dim> localIntersection;
865
866 // compute intersection of ray in direction of ip velocity with elem side
867 // we search the ray only in upwind direction
868 if (!ElementSideRayIntersection<TRefElem, TWorldDim>(side, globalIntersection, localIntersection,
869 IP, IPVel, false /* i.e. search upwind */, vCornerCoords))
870 UG_THROW("GetNodeNextToCut: ElementSideRayIntersection failed");
871
872 // get reference element
873 static const TRefElem &rRefElem = Provider<TRefElem>::get();
874 const int dim = TRefElem::dim;
875
876 // reset minimum
877 number min = std::numeric_limits<number>::max();
878
879 // loop corners of side
880 for (size_t i = 0; i < rRefElem.num(dim - 1, side, 0); ++i)
881 {
882 // get corner
883 const size_t co = rRefElem.id(dim - 1, side, 0, i);
884
885 // Compute Distance to intersection
886 number dist = VecDistanceSq(globalIntersection, vCornerCoords[co]);
887
888 // if distance is smaller, choose this node
889 if (dist < min)
890 {
891 min = dist;
892 coOut = co;
893 }
894 }
895}
896
897template <int TDim>
898template <typename TFVGeom>
900 update(const TFVGeom *geo,
901 const MathVector<dim> *Velocity,
902 const MathMatrix<dim, dim> *DiffDisp,
903 bool computeDeriv)
904{
905 UG_ASSERT(geo != nullptr, "Null pointer");
906 UG_ASSERT(Velocity != nullptr, "Null pointer");
907
908 // \todo: think about: this should be something like scvf.num_sh()
909 const size_t numSH = geo->num_sh();
910
911 // loop subcontrol volume faces
912 for (size_t ip = 0; ip < geo->num_scvf(); ++ip)
913 {
914 // get subcontrol volume face
915 const typename TFVGeom::SCVF &scvf = geo->scvf(ip);
916
917 // Compute flux
918 const number flux = VecDot(scvf.normal(), Velocity[ip]);
919
920 // Switch to no upwind in case of small velocity
921 if (VecTwoNorm(Velocity[ip]) < 1e-14)
922 {
923 // No upwind!
924 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
925 conv_shape(ip, sh) = flux * scvf.shape(sh);
926 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
927 conv_shape(ip, sh) = 0.0;
928 if (computeDeriv)
929 {
930 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
931 VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
932 // temporary: case of hanging nodes
933 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
934 VecSet(D_vel(ip, sh), 0.0);
935 }
936 continue;
937 }
938
939 // initialize shapes with zeroes
940 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
941 conv_shape(ip, sh) = 0.0;
942
943 // Hanging nodes
944 // this is introduced here, hopefully temporarily: The problem is, that
945 // for hanging nodes the number of shape function is not the number of
946 // corners, but scvf.num_sh() currently returns the number of corners.
947 // this is actually enough to interpolate the function, but still we
948 // should reset the interpolation adding for hanging dofs to zero
949 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
950 conv_shape(ip, sh) = 0.0;
951
952 const MathVector<dim> *vCornerCoords = geo->corners();
953 size_t up = 0;
954 try
955 {
956 GetNodeNextToCut<typename TFVGeom::ref_elem_type, dim>(up, scvf.global_ip(), Velocity[ip], vCornerCoords);
957 }
958 UG_CATCH_THROW("Skewed upwind: Cannot find node next to cut");
959
960 conv_shape(ip, up) = flux;
961
962 // Write Derivatives if wanted
963 if (computeDeriv)
964 {
965 for (size_t sh = 0; sh < numSH; ++sh)
966 VecSet(D_vel(ip, sh), 0.0);
967 D_vel(ip, up) = scvf.normal();
968 }
969
970 // The shapes do not depend of the diffusion tensor
971 }
972
973 // we're done
974 return true;
975}
976
978// Linear Profile Skewed Upwind
980
981template <int TDim>
983 : public IConvectionShapes<TDim>
984{
985public:
988
991
993 static const int dim = TDim;
994
995protected:
996 // explicitly forward some function
999 using base_type::D_vel;
1003
1004public:
1007 {
1008 // the shapes do not depend on the DiffDisp. Thus, we can set the
1009 // derivative to be always zero w.r.t. the DiffDisp for all shapes
1011
1012 // register evaluation function
1013 boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
1014 boost::mpl::for_each<boost::mpl::range_c<int, 1, dim + 1>>(RegisterRefDimFunc(this));
1015 }
1016
1018 template <typename TFVGeom>
1019 bool update(const TFVGeom *geo,
1020 const MathVector<dim> *Velocity,
1021 const MathMatrix<dim, dim> *DiffDisp,
1022 bool computeDeriv);
1023
1024private:
1027 {
1030 template <typename TElem>
1031 void operator()(TElem &)
1032 {
1033 m_pThis->template register_func_for_elem_fvgeom<TElem, FV1Geometry<TElem, dim>>();
1034 //m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
1035 //m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
1036 }
1037 };
1038
1040 template <typename TElem, typename TFVGeom>
1042 {
1043 typedef bool (this_type::*TFunc)(const TFVGeom *, const MathVector<dim> *, const MathMatrix<dim, dim> *, bool);
1044 base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
1045 }
1046
1049 {
1052 template <typename TRefDim>
1053 void operator()(TRefDim &) { m_pThis->register_func_for_refDim<TRefDim::value>(); }
1054 };
1055
1057 template <int refDim>
1059 {
1060 //typedef DimFV1Geometry<refDim, dim> TGeom;
1061 //typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
1062 //base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
1063 }
1064};
1065
1066template <int TDim>
1067template <typename TFVGeom>
1069 update(const TFVGeom *geo,
1070 const MathVector<dim> *Velocity,
1071 const MathMatrix<dim, dim> *DiffDisp,
1072 bool computeDeriv)
1073{
1074 UG_ASSERT(geo != nullptr, "Null pointer");
1075 UG_ASSERT(Velocity != nullptr, "Null pointer");
1076
1077 // \todo: think about: this should be something like scvf.num_sh()
1078 const size_t numSH = geo->num_sh();
1079
1080 // loop subcontrol volume faces
1081 for (size_t ip = 0; ip < geo->num_scvf(); ++ip)
1082 {
1083 // get subcontrol volume face
1084 const typename TFVGeom::SCVF &scvf = geo->scvf(ip);
1085
1086 // Compute flux
1087 const number flux = VecDot(scvf.normal(), Velocity[ip]);
1088
1089 // Switch to no upwind in case of small velocity
1090 if (VecTwoNorm(Velocity[ip]) < 1e-14)
1091 {
1092 // No upwind!
1093 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1094 conv_shape(ip, sh) = flux * scvf.shape(sh);
1095 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1096 conv_shape(ip, sh) = 0.0;
1097 if (computeDeriv)
1098 {
1099 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1100 VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
1101 // temporary, see above
1102 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1103 VecSet(D_vel(ip, sh), 0.0);
1104 }
1105 continue;
1106 }
1107
1108 // initialize shapes with zeroes
1109 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1110 conv_shape(ip, sh) = 0.0;
1111
1112 // Hanging nodes
1113 // this is introduced here, hopefully temporarily: The problem is, that
1114 // for hanging nodes the number of shape function is not the number of
1115 // corners, but scvf.num_sh() currently returns the number of corners.
1116 // this is actually enough to interpolate the function, but still we
1117 // should reset the interpolation adding for hanging dofs to zero
1118 for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1119 conv_shape(ip, sh) = 0.0;
1120
1121 const MathVector<dim> *vCornerCoords = geo->corners();
1122 // Reference element type
1123 typedef typename TFVGeom::ref_elem_type TRefElem;
1124 size_t side = 0;
1125 MathVector<dim> globalIntersection;
1126 MathVector<TRefElem::dim> localIntersection;
1127
1128 try
1129 {
1130 ElementSideRayIntersection<TRefElem, dim>(side, globalIntersection, localIntersection,
1131 scvf.global_ip(), Velocity[ip], false /* search upwind */, vCornerCoords);
1132 }
1133 UG_CATCH_THROW("GetLinearProfileSkewedUpwindShapes: Cannot find cut side.");
1134
1135 // get linear trial space
1136 static const ReferenceObjectID roid = TRefElem::REFERENCE_OBJECT_ID;
1137 const LocalShapeFunctionSet<TRefElem::dim> &TrialSpace =
1138 LocalFiniteElementProvider::get<TRefElem::dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
1139
1140 // get Reference Element
1141 static const TRefElem &rRefElem = Provider<TRefElem>::get();
1142
1143 // Shape function values
1144 std::vector<number> vShape;
1145 // Values of shape function gradients
1146 std::vector<MathVector<TRefElem::dim>> vShapeGrad;
1147 // Get values at the intersection point
1148 TrialSpace.shapes(vShape, localIntersection);
1149 TrialSpace.grads(vShapeGrad, localIntersection);
1150
1151 size_t num_corners_inters_side = rRefElem.num(dim - 1, side, 0);
1152
1153 // loop corners of side
1154 for (size_t j = 0; j < num_corners_inters_side; ++j)
1155 {
1156 // get corner
1157 const size_t co = rRefElem.id(dim-1, side, 0, j);
1158
1159 // write shape
1160 conv_shape(ip, co) = flux * vShape[co];
1161 }
1162
1163 // Write Derivatives if wanted
1164 if (computeDeriv)
1165 {
1166 // No derivatives yet!
1167 }
1168
1169 // The shapes do not depend of the diffusion tensor
1170 }
1171
1172 // we're done
1173 return true;
1174}
1175
1176
1177} // end namespace ug
1178
1179#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE__ */
Definition conv_shape.h:194
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:197
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:251
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:267
static const int dim
Dimension.
Definition conv_shape.h:203
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV1Geometry
Definition conv_shape.h:278
ConvectionShapesFullUpwind()
constructor
Definition conv_shape.h:216
ConvectionShapesFullUpwind< TDim > this_type
This class.
Definition conv_shape.h:200
ConvectionShapesLinearProfileSkewedUpwind< TDim > this_type
This class.
Definition conv_shape.h:990
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:1058
ConvectionShapesLinearProfileSkewedUpwind()
constructor
Definition conv_shape.h:1006
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:1041
static const int dim
Dimension.
Definition conv_shape.h:993
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:987
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV geometry
Definition conv_shape.h:1069
Definition conv_shape.h:53
static const int dim
Dimension.
Definition conv_shape.h:62
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:110
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV geometry
Definition conv_shape.h:138
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:56
ConvectionShapesNoUpwind()
constructor
Definition conv_shape.h:75
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:126
ConvectionShapesNoUpwind< TDim > this_type
This class.
Definition conv_shape.h:59
Definition conv_shape.h:536
ConvectionShapesPartialUpwind< TDim > this_type
This class.
Definition conv_shape.h:542
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:539
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:588
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:604
ConvectionShapesPartialUpwind()
constructor
Definition conv_shape.h:558
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV1Geometry
Definition conv_shape.h:616
static const int dim
Dimension.
Definition conv_shape.h:545
Definition conv_shape.h:772
static const int dim
Dimension.
Definition conv_shape.h:781
ConvectionShapesSkewedUpwind< TDim > this_type
This class.
Definition conv_shape.h:778
ConvectionShapesSkewedUpwind()
constructor
Definition conv_shape.h:794
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:846
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV geometry
Definition conv_shape.h:900
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:829
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:775
Definition conv_shape.h:349
void register_func_for_elem_fvgeom()
registers the update function for an element type and a FV geometry
Definition conv_shape.h:427
void register_func_for_refDim()
registers the update function for a reference dimension
Definition conv_shape.h:443
number m_weight
Definition conv_shape.h:408
ConvectionShapesWeightedUpwind< TDim > this_type
This class.
Definition conv_shape.h:355
ConvectionShapesWeightedUpwind(number weight)
constructor
Definition conv_shape.h:383
IConvectionShapes< TDim > base_type
Base class.
Definition conv_shape.h:352
bool update(const TFVGeom *geo, const MathVector< dim > *Velocity, const MathMatrix< dim, dim > *DiffDisp, bool computeDeriv)
update of values for FV1Geometry
Definition conv_shape.h:454
void set_weight(number weight)
set weighting between full upwind (1.0) and no upwind (0.0)
Definition conv_shape.h:397
ConvectionShapesWeightedUpwind()
constructor
Definition conv_shape.h:371
static const int dim
Dimension.
Definition conv_shape.h:358
Geometry and shape functions for 1st order Vertex-Centered Finite Volume.
Definition fv1_geom.h:665
Interface class for upwind methods.
Definition conv_shape_interface.h:100
void register_update_func(TAssFunc func)
register a update function for a Geometry
Definition conv_shape_interface.h:249
void set_non_zero_deriv_diffusion_flag(bool flag)
sets the shape ip flag
Definition conv_shape_interface.h:174
MathMatrix< dim, dim > & conv_shape_diffusion(size_t scvf, size_t sh)
non-const access to upwind shapes for ip vel
Definition conv_shape_interface.h:193
const MathVector< dim > & D_vel(size_t scvf, size_t sh) const
upwind shape for corner vel
Definition conv_shape_interface.h:141
number & conv_shape(size_t scvf, size_t sh)
non-const access to ip velocity (i.e. interpolated velocity at ip)
Definition conv_shape_interface.h:177
bool non_zero_deriv_diffusion() const
returns if upwind shape w.r.t. ip vel is non-zero
Definition conv_shape_interface.h:149
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
@ LAGRANGE
Definition local_finite_element_id.h:104
virtual base class for local shape function sets
Definition local_shape_function_set.h:70
virtual void grads(grad_type *vGrad, const MathVector< dim > &x) const =0
returns all gradients evaluated at a point
virtual void shapes(shape_type *vShape, const MathVector< dim > &x) const =0
returns all shape functions evaluated at a point
A class for fixed size, dense matrices.
Definition math_matrix.h:63
a mathematical Vector with N entries.
Definition math_vector.h:97
static TClass & get()
returns a singleton based on the identifier
Definition provider.h:54
void MatSet(matrix_t &mInOut, typename matrix_t::value_type s)
Set each matrix entry to a scalar (componentwise)
Definition math_matrix_functions_common_impl.hpp:794
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
double number
Definition types.h:124
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
void VecScaleAppend(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1)
Scales a Vector and adds it to a second vector.
Definition math_vector_functions_common_impl.hpp:126
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
vector_t::value_type VecDistanceSq(const vector_t &v1, const vector_t &v2)
returns the squared distance of two vector_ts.
Definition math_vector_functions_common_impl.hpp:351
void VecScale(vector_t &vOut, const vector_t &v, typename vector_t::value_type s)
scales a MathVector<N>
Definition math_vector_functions_common_impl.hpp:252
vector_t::value_type VecTwoNorm(const vector_t &v)
Definition math_vector_functions_common_impl.hpp:577
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
returns the dot-product of two vector_ts
Definition math_vector_functions_common_impl.hpp:385
the ug namespace
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition grid_base_objects.h:74
void GetNodeNextToCut(size_t &coOut, const MathVector< TWorldDim > &IP, const MathVector< TWorldDim > &IPVel, const MathVector< TWorldDim > *vCornerCoords)
computes the closest node to a elem side ray intersection
Definition conv_shape.h:856
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:238
void operator()(TElem &)
Definition conv_shape.h:241
this_type * m_pThis
Definition conv_shape.h:239
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:240
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:259
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:261
this_type * m_pThis
Definition conv_shape.h:260
void operator()(TRefDim &)
Definition conv_shape.h:262
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:1027
void operator()(TElem &)
Definition conv_shape.h:1031
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:1029
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:1049
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:1051
void operator()(TRefDim &)
Definition conv_shape.h:1053
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:97
this_type * m_pThis
Definition conv_shape.h:98
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:99
void operator()(TElem &)
Definition conv_shape.h:100
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:118
void operator()(TRefDim &)
Definition conv_shape.h:121
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:120
this_type * m_pThis
Definition conv_shape.h:119
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:576
void operator()(TElem &)
Definition conv_shape.h:579
this_type * m_pThis
Definition conv_shape.h:577
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:578
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:596
this_type * m_pThis
Definition conv_shape.h:597
void operator()(TRefDim &)
Definition conv_shape.h:599
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:598
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:815
this_type * m_pThis
Definition conv_shape.h:816
void operator()(TElem &)
Definition conv_shape.h:819
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:817
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:837
this_type * m_pThis
Definition conv_shape.h:838
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:839
void operator()(TRefDim &)
Definition conv_shape.h:841
functor for registering the shapes for the element-templated FV geometries
Definition conv_shape.h:414
RegisterElemFunc(this_type *pThis)
Definition conv_shape.h:416
void operator()(TElem &)
Definition conv_shape.h:417
this_type * m_pThis
Definition conv_shape.h:415
functor for registering the shapes for the reference-dimension-templated FV geometries
Definition conv_shape.h:435
void operator()(TRefDim &)
Definition conv_shape.h:438
RegisterRefDimFunc(this_type *pThis)
Definition conv_shape.h:437
this_type * m_pThis
Definition conv_shape.h:436