Line data Source code
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
40 : #include "conv_shape_interface.h"
41 : #include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
42 : #include "lib_disc/spatial_disc/disc_util/hfv1_geom.h"
43 :
44 : namespace ug{
45 :
46 : /////////////////////////////////////////////////////////////////////////////
47 : // No Upwind
48 : /////////////////////////////////////////////////////////////////////////////
49 :
50 : template <int TDim>
51 : class ConvectionShapesNoUpwind
52 : : public IConvectionShapes<TDim>
53 : {
54 : public:
55 : /// Base class
56 : typedef IConvectionShapes<TDim> base_type;
57 :
58 : /// This class
59 : typedef ConvectionShapesNoUpwind<TDim> this_type;
60 :
61 : /// Dimension
62 : static const int dim = TDim;
63 :
64 : protected:
65 : // explicitly forward some function
66 : using base_type::set_non_zero_deriv_diffusion_flag;
67 : using base_type::conv_shape;
68 : using base_type::D_vel;
69 : using base_type::conv_shape_diffusion;
70 : using base_type::non_zero_deriv_diffusion;
71 : using base_type::register_update_func;
72 :
73 : public:
74 : /// constructor
75 0 : ConvectionShapesNoUpwind()
76 0 : {
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
79 : set_non_zero_deriv_diffusion_flag(false);
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 0 : }
85 :
86 : /// update of values for FV geometry
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 :
95 : /// functor for registering the shapes for the element-templated FV geometries
96 : struct RegisterElemFunc
97 : {
98 : this_type* m_pThis;
99 : RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
100 0 : template<typename TElem> void operator() (TElem &)
101 : {
102 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
103 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
104 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
105 0 : }
106 : };
107 :
108 : /// registers the update function for an element type and a FV geometry
109 : template <typename TElem, typename TFVGeom>
110 : void register_func_for_elem_fvgeom()
111 : {
112 : typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
113 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
114 : }
115 :
116 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
117 : struct RegisterRefDimFunc
118 : {
119 : this_type* m_pThis;
120 : RegisterRefDimFunc(this_type * pThis) : m_pThis(pThis) {}
121 0 : template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
122 : };
123 :
124 : /// registers the update function for a reference dimension
125 : template <int refDim>
126 : void register_func_for_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 0 : base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
131 : }
132 : };
133 :
134 : template <int TDim>
135 : template <typename TFVGeom>
136 : bool
137 0 : ConvectionShapesNoUpwind<TDim>::
138 : update(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 0 : 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 0 : const number flux = VecDot(scvf.normal(), Velocity[ip]);
157 :
158 : // Write Shapes
159 0 : for(size_t sh = 0; sh < scvf.num_sh(); sh++)
160 0 : 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 0 : for(size_t sh = scvf.num_sh(); sh < numSH; sh++)
168 0 : conv_shape(ip, sh) = 0.0;
169 :
170 : // Write Derivatives if wanted
171 0 : if(computeDeriv){
172 0 : 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 0 : 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 0 : return true;
185 : }
186 :
187 : /////////////////////////////////////////////////////////////////////////////
188 : // Full Upwind
189 : /////////////////////////////////////////////////////////////////////////////
190 :
191 : template <int TDim>
192 : class ConvectionShapesFullUpwind
193 : : public IConvectionShapes<TDim>
194 : {
195 : public:
196 : /// Base class
197 : typedef IConvectionShapes<TDim> base_type;
198 :
199 : /// This class
200 : typedef ConvectionShapesFullUpwind<TDim> this_type;
201 :
202 : /// Dimension
203 : static const int dim = TDim;
204 :
205 : protected:
206 : // explicitly forward some function
207 : using base_type::set_non_zero_deriv_diffusion_flag;
208 : using base_type::conv_shape;
209 : using base_type::D_vel;
210 : using base_type::conv_shape_diffusion;
211 : using base_type::non_zero_deriv_diffusion;
212 : using base_type::register_update_func;
213 :
214 : public:
215 : /// constructor
216 0 : ConvectionShapesFullUpwind()
217 0 : {
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
220 : set_non_zero_deriv_diffusion_flag(false);
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 0 : }
226 :
227 : /// update of values for FV1Geometry
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 :
236 : /// functor for registering the shapes for the element-templated FV geometries
237 : struct RegisterElemFunc
238 : {
239 : this_type* m_pThis;
240 : RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
241 0 : template<typename TElem> void operator() (TElem &)
242 : {
243 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
244 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
245 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
246 0 : }
247 : };
248 :
249 : /// registers the update function for an element type and a FV geometry
250 : template <typename TElem, typename TFVGeom>
251 : void register_func_for_elem_fvgeom()
252 : {
253 : typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
254 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
255 : }
256 :
257 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
258 : struct RegisterRefDimFunc
259 : {
260 : this_type* m_pThis;
261 : RegisterRefDimFunc(this_type * pThis) : m_pThis(pThis) {}
262 0 : template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
263 : };
264 :
265 : /// registers the update function for a reference dimension
266 : template <int refDim>
267 : void register_func_for_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 0 : base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
272 : }
273 : };
274 :
275 : template <int TDim>
276 : template <typename TFVGeom>
277 0 : bool ConvectionShapesFullUpwind<TDim>::
278 : update(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 0 : 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 0 : 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 0 : if((from >= scvf.num_sh()) || (to >= scvf.num_sh())){
303 : // No upwind...
304 : // Write Shapes
305 0 : for(size_t sh = 0; sh < scvf.num_sh(); sh++)
306 0 : 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 0 : if(computeDeriv){
312 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
313 : VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
314 : }
315 :
316 0 : continue;
317 0 : }
318 :
319 : // Full upwind below...
320 : // Choose Upwind corner
321 0 : size_t up = (flux >= 0) ? from : to;
322 :
323 : // Write Shapes
324 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh) conv_shape(ip, sh) = 0.0;
325 :
326 0 : conv_shape(ip, up) = flux;
327 :
328 : // Write Derivatives if wanted
329 0 : if(computeDeriv)
330 : {
331 0 : 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 0 : return true;
340 : }
341 :
342 : /////////////////////////////////////////////////////////////////////////////
343 : // Weighted Upwind
344 : /////////////////////////////////////////////////////////////////////////////
345 :
346 : template <int TDim>
347 : class ConvectionShapesWeightedUpwind
348 : : public IConvectionShapes<TDim>
349 : {
350 : public:
351 : /// Base class
352 : typedef IConvectionShapes<TDim> base_type;
353 :
354 : /// This class
355 : typedef ConvectionShapesWeightedUpwind<TDim> this_type;
356 :
357 : /// Dimension
358 : static const int dim = TDim;
359 :
360 : protected:
361 : // explicitly forward some function
362 : using base_type::set_non_zero_deriv_diffusion_flag;
363 : using base_type::conv_shape;
364 : using base_type::D_vel;
365 : using base_type::conv_shape_diffusion;
366 : using base_type::non_zero_deriv_diffusion;
367 : using base_type::register_update_func;
368 :
369 : public:
370 : /// constructor
371 0 : ConvectionShapesWeightedUpwind() : m_weight(0.5)
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
375 : set_non_zero_deriv_diffusion_flag(false);
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 0 : }
381 :
382 : /// constructor
383 0 : ConvectionShapesWeightedUpwind(number weight)
384 0 : {
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
389 : set_non_zero_deriv_diffusion_flag(false);
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 0 : }
395 :
396 : /// set weighting between full upwind (1.0) and no upwind (0.0)
397 0 : void set_weight(number weight) {m_weight = weight;}
398 :
399 : /// update of values for FV1Geometry
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)
408 : number m_weight;
409 :
410 : private:
411 :
412 : /// functor for registering the shapes for the element-templated FV geometries
413 : struct RegisterElemFunc
414 : {
415 : this_type* m_pThis;
416 : RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
417 0 : template<typename TElem> void operator() (TElem &)
418 : {
419 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
420 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
421 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
422 0 : }
423 : };
424 :
425 : /// registers the update function for an element type and a FV geometry
426 : template <typename TElem, typename TFVGeom>
427 : void register_func_for_elem_fvgeom()
428 : {
429 : typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
430 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
431 : }
432 :
433 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
434 : struct RegisterRefDimFunc
435 : {
436 : this_type* m_pThis;
437 : RegisterRefDimFunc(this_type * pThis) : m_pThis(pThis) {}
438 0 : template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
439 : };
440 :
441 : /// registers the update function for a reference dimension
442 : template <int refDim>
443 : void register_func_for_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 0 : base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
448 : }
449 : };
450 :
451 : template <int TDim>
452 : template <typename TFVGeom>
453 0 : bool ConvectionShapesWeightedUpwind<TDim>::
454 : update(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 0 : 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 0 : 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 0 : if((from >= scvf.num_sh()) || (to >= scvf.num_sh())){
479 : // No upwind...
480 : // Write Shapes
481 0 : for(size_t sh = 0; sh < scvf.num_sh(); sh++)
482 0 : 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 0 : if(computeDeriv){
488 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
489 : VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
490 : }
491 :
492 0 : continue;
493 0 : }
494 :
495 : // Choose Upwind corner
496 0 : size_t up = (flux >= 0) ? from : to;
497 :
498 : // write no upwind part of shapes
499 0 : const number noUpFlux = (1.-m_weight)*flux;
500 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
501 0 : conv_shape(ip, sh) = noUpFlux * scvf.shape(sh);
502 :
503 : // add full upwind part of shapes
504 0 : conv_shape(ip, up) += m_weight * flux;
505 :
506 : // Write Derivatives if wanted
507 0 : if(computeDeriv)
508 : {
509 : // write no upwind part of derivatives
510 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
511 0 : VecScale(D_vel(ip, sh), scvf.normal(),
512 0 : (1.-m_weight)*scvf.shape(sh));
513 : // see comment above
514 0 : 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 0 : 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 0 : return true;
526 : }
527 :
528 :
529 : /////////////////////////////////////////////////////////////////////////////
530 : // Partial Upwind
531 : /////////////////////////////////////////////////////////////////////////////
532 :
533 : template <int TDim>
534 : class ConvectionShapesPartialUpwind
535 : : public IConvectionShapes<TDim>
536 : {
537 : public:
538 : /// Base class
539 : typedef IConvectionShapes<TDim> base_type;
540 :
541 : /// This class
542 : typedef ConvectionShapesPartialUpwind<TDim> this_type;
543 :
544 : /// Dimension
545 : static const int dim = TDim;
546 :
547 : protected:
548 : // explicitly forward some function
549 : using base_type::set_non_zero_deriv_diffusion_flag;
550 : using base_type::conv_shape;
551 : using base_type::D_vel;
552 : using base_type::conv_shape_diffusion;
553 : using base_type::non_zero_deriv_diffusion;
554 : using base_type::register_update_func;
555 :
556 : public:
557 : /// constructor
558 0 : ConvectionShapesPartialUpwind()
559 0 : {
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 0 : }
564 :
565 : /// update of values for FV1Geometry
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 :
574 : /// functor for registering the shapes for the element-templated FV geometries
575 : struct RegisterElemFunc
576 : {
577 : this_type* m_pThis;
578 : RegisterElemFunc(this_type * pThis) : m_pThis(pThis) {}
579 0 : template<typename TElem> void operator() (TElem &)
580 : {
581 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1Geometry<TElem, dim> >();
582 0 : m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
583 0 : }
584 : };
585 :
586 : /// registers the update function for an element type and a FV geometry
587 : template <typename TElem, typename TFVGeom>
588 : void register_func_for_elem_fvgeom()
589 : {
590 : typedef bool (this_type::*TFunc) (const TFVGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
591 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
592 : }
593 :
594 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
595 : struct RegisterRefDimFunc
596 : {
597 : this_type* m_pThis;
598 : RegisterRefDimFunc(this_type * pThis) : m_pThis(pThis) {}
599 0 : template<typename TRefDim> void operator() (TRefDim &) {m_pThis->register_func_for_refDim<TRefDim::value> ();}
600 : };
601 :
602 : /// registers the update function for a reference dimension
603 : template <int refDim>
604 : void register_func_for_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 0 : base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
609 : }
610 : };
611 :
612 : template <int TDim>
613 : template <typename TFVGeom>
614 : bool
615 0 : ConvectionShapesPartialUpwind<TDim>::
616 : update(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 0 : const number vol = ElementSize<dim>(geo->roid(), geo->corners());
628 :
629 : // loop subcontrol volume faces
630 0 : 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 0 : if(DiffDisp != nullptr)
648 : {
649 : // Get Gradients
650 : MathVector<dim> DiffGrad;
651 :
652 : // Compute DiffGrad = D * Grad Phi_to
653 0 : 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 0 : lambda = - GradDiffGrad * vol;
660 : }
661 :
662 : // Compute flux
663 0 : const number flux = VecDot(scvf.normal(), Velocity[i]);
664 :
665 : // reset values
666 0 : for (size_t sh = 0; sh < scvf.num_sh(); ++sh)
667 0 : conv_shape(i, sh) = 0.0;
668 0 : if (computeDeriv)
669 0 : for (size_t sh = 0; sh < scvf.num_sh(); ++sh) {
670 : VecSet(D_vel(i, sh), 0.0);
671 : MatSet(conv_shape_diffusion(i, sh), 0.0);
672 : }
673 :
674 :
675 : // check: Currently hanging nodes not supported.
676 : // /todo: add support for hanging nodes
677 0 : if (from >= scvf.num_sh() || to >= scvf.num_sh())
678 0 : UG_THROW("PartialUpwind: Currently not implemented for hanging nodes.")
679 :
680 : ///////////////////////////////////////////////////////////////////
681 : // Case 1:
682 : // full upwind is used
683 : ///////////////////////////////////////////////////////////////////
684 0 : if(lambda <= 0 || DiffDisp == nullptr)
685 : {
686 : // Choose Upwind corner
687 0 : const size_t up = (flux >= 0) ? scvf.from() : scvf.to();
688 :
689 : // Write Shapes
690 0 : conv_shape(i, up) = flux;
691 :
692 : // Write Derivatives if wanted
693 0 : if(computeDeriv)
694 : {
695 : // set derivative
696 : D_vel(i, up) = scvf.normal();
697 :
698 : // does not depend on diffusion
699 : set_non_zero_deriv_diffusion_flag(false);
700 : }
701 :
702 : // everything done
703 0 : continue;
704 0 : }
705 :
706 : ///////////////////////////////////////////////////////////////////
707 : // Case 2:
708 : // The case of the diffusion dominance (central differences)
709 : ///////////////////////////////////////////////////////////////////
710 0 : if (2 * lambda > fabs(flux))
711 : {
712 0 : conv_shape(i, from) = flux / 2.0;
713 0 : conv_shape(i, to) = flux / 2.0;
714 :
715 0 : if(computeDeriv)
716 : {
717 : set_non_zero_deriv_diffusion_flag(false);
718 :
719 : VecScale(D_vel(i,from), scvf.normal(), 1.0/2.0);
720 : VecScale(D_vel(i, to), scvf.normal(), 1.0/2.0);
721 : }
722 :
723 : // everything done
724 0 : continue;
725 : }
726 :
727 : ///////////////////////////////////////////////////////////////////
728 : // Case 3:
729 : // The cases of the convection dominance
730 : ///////////////////////////////////////////////////////////////////
731 : set_non_zero_deriv_diffusion_flag(true);
732 0 : if (flux >= 0)
733 : {
734 0 : conv_shape(i, from) = flux - lambda;
735 0 : conv_shape(i, to) = lambda;
736 :
737 0 : if(computeDeriv)
738 : D_vel(i,from) = scvf.normal();
739 : }
740 : else
741 : {
742 0 : conv_shape(i, from) = - lambda;
743 0 : conv_shape(i, to) = flux + lambda;
744 :
745 0 : if(computeDeriv)
746 : D_vel(i,to) = scvf.normal();
747 : }
748 :
749 0 : if (computeDeriv)
750 : {
751 0 : for (size_t k = 0; k < (size_t)dim; k++)
752 0 : for (size_t l = 0; l < (size_t)dim; l++)
753 : {
754 0 : conv_shape_diffusion(i, from)(k,l) = gradFrom[k]*gradTo[l]*vol;
755 0 : conv_shape_diffusion(i, to)(k,l) = - gradFrom[k]*gradTo[l]*vol;
756 : }
757 : }
758 : }
759 :
760 : // we're done
761 0 : return true;
762 : }
763 :
764 : /////////////////////////////////////////////////////////////////////////////
765 : // Skewed Upwind
766 : /////////////////////////////////////////////////////////////////////////////
767 :
768 : template <int TDim>
769 : class ConvectionShapesSkewedUpwind
770 : : public IConvectionShapes<TDim>
771 : {
772 : public:
773 : /// Base class
774 : typedef IConvectionShapes<TDim> base_type;
775 :
776 : /// This class
777 : typedef ConvectionShapesSkewedUpwind<TDim> this_type;
778 :
779 : /// Dimension
780 : static const int dim = TDim;
781 :
782 : protected:
783 : // explicitly forward some function
784 : using base_type::conv_shape;
785 : using base_type::conv_shape_diffusion;
786 : using base_type::D_vel;
787 : using base_type::non_zero_deriv_diffusion;
788 : using base_type::register_update_func;
789 : using base_type::set_non_zero_deriv_diffusion_flag;
790 :
791 : public:
792 : /// constructor
793 0 : ConvectionShapesSkewedUpwind()
794 0 : {
795 : // the shapes do not depend on the DiffDisp. Thus, we can set the
796 : // derivative to be always zero w.r.t. the DiffDisp for all shapes
797 : set_non_zero_deriv_diffusion_flag(false);
798 :
799 : // register evaluation function
800 : boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
801 : boost::mpl::for_each<boost::mpl::range_c<int, 1, dim + 1>>(RegisterRefDimFunc(this));
802 0 : }
803 :
804 : /// update of values for FV geometry
805 : template <typename TFVGeom>
806 : bool update(const TFVGeom *geo,
807 : const MathVector<dim> *Velocity,
808 : const MathMatrix<dim, dim> *DiffDisp,
809 : bool computeDeriv);
810 :
811 : private:
812 : /// functor for registering the shapes for the element-templated FV geometries
813 : struct RegisterElemFunc
814 : {
815 : this_type *m_pThis;
816 : RegisterElemFunc(this_type *pThis) : m_pThis(pThis) {}
817 : template <typename TElem>
818 : void operator()(TElem &)
819 : {
820 0 : m_pThis->template register_func_for_elem_fvgeom<TElem, FV1Geometry<TElem, dim>>();
821 : //m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
822 : //m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
823 0 : }
824 : };
825 :
826 : /// registers the update function for an element type and a FV geometry
827 : template <typename TElem, typename TFVGeom>
828 : void register_func_for_elem_fvgeom()
829 : {
830 : typedef bool (this_type::*TFunc)(const TFVGeom *, const MathVector<dim> *, const MathMatrix<dim, dim> *, bool);
831 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
832 : }
833 :
834 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
835 : struct RegisterRefDimFunc
836 : {
837 : this_type *m_pThis;
838 : RegisterRefDimFunc(this_type *pThis) : m_pThis(pThis) {}
839 : template <typename TRefDim>
840 : void operator()(TRefDim &) { m_pThis->register_func_for_refDim<TRefDim::value>(); }
841 : };
842 :
843 : /// registers the update function for a reference dimension
844 : template <int refDim>
845 : void register_func_for_refDim()
846 : {
847 : //typedef DimFV1Geometry<refDim, dim> TGeom;
848 : //typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
849 : //base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
850 : }
851 : };
852 :
853 : /// computes the closest node to a elem side ray intersection
854 : template <typename TRefElem, int TWorldDim>
855 0 : void GetNodeNextToCut(size_t &coOut,
856 : const MathVector<TWorldDim> &IP,
857 : const MathVector<TWorldDim> &IPVel,
858 : const MathVector<TWorldDim> *vCornerCoords)
859 : {
860 : // help variables
861 0 : size_t side = 0;
862 : MathVector<TWorldDim> globalIntersection;
863 : MathVector<TRefElem::dim> localIntersection;
864 :
865 : // compute intersection of ray in direction of ip velocity with elem side
866 : // we search the ray only in upwind direction
867 0 : if (!ElementSideRayIntersection<TRefElem, TWorldDim>(side, globalIntersection, localIntersection,
868 : IP, IPVel, false /* i.e. search upwind */, vCornerCoords))
869 0 : UG_THROW("GetNodeNextToCut: ElementSideRayIntersection failed");
870 :
871 : // get reference element
872 0 : static const TRefElem &rRefElem = Provider<TRefElem>::get();
873 : const int dim = TRefElem::dim;
874 :
875 : // reset minimum
876 : number min = std::numeric_limits<number>::max();
877 :
878 : // loop corners of side
879 0 : for (size_t i = 0; i < rRefElem.num(dim - 1, side, 0); ++i)
880 : {
881 : // get corner
882 0 : const size_t co = rRefElem.id(dim - 1, side, 0, i);
883 :
884 : // Compute Distance to intersection
885 0 : number dist = VecDistanceSq(globalIntersection, vCornerCoords[co]);
886 :
887 : // if distance is smaller, choose this node
888 0 : if (dist < min)
889 : {
890 : min = dist;
891 0 : coOut = co;
892 : }
893 : }
894 0 : }
895 :
896 : template <int TDim>
897 : template <typename TFVGeom>
898 0 : bool ConvectionShapesSkewedUpwind<TDim>::
899 : update(const TFVGeom *geo,
900 : const MathVector<dim> *Velocity,
901 : const MathMatrix<dim, dim> *DiffDisp,
902 : bool computeDeriv)
903 : {
904 : UG_ASSERT(geo != nullptr, "Null pointer");
905 : UG_ASSERT(Velocity != nullptr, "Null pointer");
906 :
907 : // \todo: think about: this should be something like scvf.num_sh()
908 : const size_t numSH = geo->num_sh();
909 :
910 : // loop subcontrol volume faces
911 0 : for (size_t ip = 0; ip < geo->num_scvf(); ++ip)
912 : {
913 : // get subcontrol volume face
914 : const typename TFVGeom::SCVF &scvf = geo->scvf(ip);
915 :
916 : // Compute flux
917 0 : const number flux = VecDot(scvf.normal(), Velocity[ip]);
918 :
919 : // Switch to no upwind in case of small velocity
920 0 : if (VecTwoNorm(Velocity[ip]) < 1e-14)
921 : {
922 : // No upwind!
923 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
924 0 : conv_shape(ip, sh) = flux * scvf.shape(sh);
925 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
926 : conv_shape(ip, sh) = 0.0;
927 0 : if (computeDeriv)
928 : {
929 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
930 : VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
931 : // temporary: case of hanging nodes
932 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
933 : VecSet(D_vel(ip, sh), 0.0);
934 : }
935 0 : continue;
936 0 : }
937 :
938 : // initialize shapes with zeroes
939 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
940 0 : conv_shape(ip, sh) = 0.0;
941 :
942 : // Hanging nodes
943 : // this is introduced here, hopefully temporarily: The problem is, that
944 : // for hanging nodes the number of shape function is not the number of
945 : // corners, but scvf.num_sh() currently returns the number of corners.
946 : // this is actually enough to interpolate the function, but still we
947 : // should reset the interpolation adding for hanging dofs to zero
948 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
949 : conv_shape(ip, sh) = 0.0;
950 :
951 : const MathVector<dim> *vCornerCoords = geo->corners();
952 0 : size_t up = 0;
953 : try
954 : {
955 0 : GetNodeNextToCut<typename TFVGeom::ref_elem_type, dim>(up, scvf.global_ip(), Velocity[ip], vCornerCoords);
956 : }
957 0 : UG_CATCH_THROW("Skewed upwind: Cannot find node next to cut");
958 :
959 0 : conv_shape(ip, up) = flux;
960 :
961 : // Write Derivatives if wanted
962 0 : if (computeDeriv)
963 : {
964 0 : for (size_t sh = 0; sh < numSH; ++sh)
965 : VecSet(D_vel(ip, sh), 0.0);
966 : D_vel(ip, up) = scvf.normal();
967 : }
968 :
969 : // The shapes do not depend of the diffusion tensor
970 : }
971 :
972 : // we're done
973 0 : return true;
974 : }
975 :
976 : /////////////////////////////////////////////////////////////////////////////
977 : // Linear Profile Skewed Upwind
978 : /////////////////////////////////////////////////////////////////////////////
979 :
980 : template <int TDim>
981 : class ConvectionShapesLinearProfileSkewedUpwind
982 : : public IConvectionShapes<TDim>
983 : {
984 : public:
985 : /// Base class
986 : typedef IConvectionShapes<TDim> base_type;
987 :
988 : /// This class
989 : typedef ConvectionShapesLinearProfileSkewedUpwind<TDim> this_type;
990 :
991 : /// Dimension
992 : static const int dim = TDim;
993 :
994 : protected:
995 : // explicitly forward some function
996 : using base_type::conv_shape;
997 : using base_type::conv_shape_diffusion;
998 : using base_type::D_vel;
999 : using base_type::non_zero_deriv_diffusion;
1000 : using base_type::register_update_func;
1001 : using base_type::set_non_zero_deriv_diffusion_flag;
1002 :
1003 : public:
1004 : /// constructor
1005 0 : ConvectionShapesLinearProfileSkewedUpwind()
1006 0 : {
1007 : // the shapes do not depend on the DiffDisp. Thus, we can set the
1008 : // derivative to be always zero w.r.t. the DiffDisp for all shapes
1009 : set_non_zero_deriv_diffusion_flag(false);
1010 :
1011 : // register evaluation function
1012 : boost::mpl::for_each<typename domain_traits<dim>::AllElemList>(RegisterElemFunc(this));
1013 : boost::mpl::for_each<boost::mpl::range_c<int, 1, dim + 1>>(RegisterRefDimFunc(this));
1014 0 : }
1015 :
1016 : /// update of values for FV geometry
1017 : template <typename TFVGeom>
1018 : bool update(const TFVGeom *geo,
1019 : const MathVector<dim> *Velocity,
1020 : const MathMatrix<dim, dim> *DiffDisp,
1021 : bool computeDeriv);
1022 :
1023 : private:
1024 : /// functor for registering the shapes for the element-templated FV geometries
1025 : struct RegisterElemFunc
1026 : {
1027 : this_type *m_pThis;
1028 : RegisterElemFunc(this_type *pThis) : m_pThis(pThis) {}
1029 : template <typename TElem>
1030 : void operator()(TElem &)
1031 : {
1032 0 : m_pThis->template register_func_for_elem_fvgeom<TElem, FV1Geometry<TElem, dim>>();
1033 : //m_pThis->template register_func_for_elem_fvgeom< TElem, FV1CondensedGeometry<TElem, dim> >();
1034 : //m_pThis->template register_func_for_elem_fvgeom< TElem, HFV1Geometry<TElem, dim> >();
1035 0 : }
1036 : };
1037 :
1038 : /// registers the update function for an element type and a FV geometry
1039 : template <typename TElem, typename TFVGeom>
1040 : void register_func_for_elem_fvgeom()
1041 : {
1042 : typedef bool (this_type::*TFunc)(const TFVGeom *, const MathVector<dim> *, const MathMatrix<dim, dim> *, bool);
1043 0 : base_type::template register_update_func<TFVGeom, TFunc>(&this_type::template update<TFVGeom>);
1044 : }
1045 :
1046 : /// functor for registering the shapes for the reference-dimension-templated FV geometries
1047 : struct RegisterRefDimFunc
1048 : {
1049 : this_type *m_pThis;
1050 : RegisterRefDimFunc(this_type *pThis) : m_pThis(pThis) {}
1051 : template <typename TRefDim>
1052 : void operator()(TRefDim &) { m_pThis->register_func_for_refDim<TRefDim::value>(); }
1053 : };
1054 :
1055 : /// registers the update function for a reference dimension
1056 : template <int refDim>
1057 : void register_func_for_refDim()
1058 : {
1059 : //typedef DimFV1Geometry<refDim, dim> TGeom;
1060 : //typedef bool (this_type::*TFunc) (const TGeom*, const MathVector<dim>*, const MathMatrix<dim, dim>*, bool);
1061 : //base_type::template register_update_func<TGeom, TFunc>(&this_type::template update<TGeom>);
1062 : }
1063 : };
1064 :
1065 : template <int TDim>
1066 : template <typename TFVGeom>
1067 0 : bool ConvectionShapesLinearProfileSkewedUpwind<TDim>::
1068 : update(const TFVGeom *geo,
1069 : const MathVector<dim> *Velocity,
1070 : const MathMatrix<dim, dim> *DiffDisp,
1071 : bool computeDeriv)
1072 : {
1073 : UG_ASSERT(geo != nullptr, "Null pointer");
1074 : UG_ASSERT(Velocity != nullptr, "Null pointer");
1075 :
1076 : // \todo: think about: this should be something like scvf.num_sh()
1077 : const size_t numSH = geo->num_sh();
1078 :
1079 : // loop subcontrol volume faces
1080 0 : for (size_t ip = 0; ip < geo->num_scvf(); ++ip)
1081 : {
1082 : // get subcontrol volume face
1083 : const typename TFVGeom::SCVF &scvf = geo->scvf(ip);
1084 :
1085 : // Compute flux
1086 0 : const number flux = VecDot(scvf.normal(), Velocity[ip]);
1087 :
1088 : // Switch to no upwind in case of small velocity
1089 0 : if (VecTwoNorm(Velocity[ip]) < 1e-14)
1090 : {
1091 : // No upwind!
1092 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1093 0 : conv_shape(ip, sh) = flux * scvf.shape(sh);
1094 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1095 : conv_shape(ip, sh) = 0.0;
1096 0 : if (computeDeriv)
1097 : {
1098 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1099 : VecScale(D_vel(ip, sh), scvf.normal(), scvf.shape(sh));
1100 : // temporary, see above
1101 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1102 : VecSet(D_vel(ip, sh), 0.0);
1103 : }
1104 0 : continue;
1105 0 : }
1106 :
1107 : // initialize shapes with zeroes
1108 0 : for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1109 0 : conv_shape(ip, sh) = 0.0;
1110 :
1111 : // Hanging nodes
1112 : // this is introduced here, hopefully temporarily: The problem is, that
1113 : // for hanging nodes the number of shape function is not the number of
1114 : // corners, but scvf.num_sh() currently returns the number of corners.
1115 : // this is actually enough to interpolate the function, but still we
1116 : // should reset the interpolation adding for hanging dofs to zero
1117 : for (size_t sh = scvf.num_sh(); sh < numSH; sh++)
1118 : conv_shape(ip, sh) = 0.0;
1119 :
1120 : const MathVector<dim> *vCornerCoords = geo->corners();
1121 : // Reference element type
1122 : typedef typename TFVGeom::ref_elem_type TRefElem;
1123 0 : size_t side = 0;
1124 : MathVector<dim> globalIntersection;
1125 : MathVector<TRefElem::dim> localIntersection;
1126 :
1127 : try
1128 : {
1129 : ElementSideRayIntersection<TRefElem, dim>(side, globalIntersection, localIntersection,
1130 : scvf.global_ip(), Velocity[ip], false /* search upwind */, vCornerCoords);
1131 : }
1132 0 : UG_CATCH_THROW("GetLinearProfileSkewedUpwindShapes: Cannot find cut side.");
1133 :
1134 : // get linear trial space
1135 : static const ReferenceObjectID roid = TRefElem::REFERENCE_OBJECT_ID;
1136 : const LocalShapeFunctionSet<TRefElem::dim> &TrialSpace =
1137 0 : LocalFiniteElementProvider::get<TRefElem::dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
1138 :
1139 : // get Reference Element
1140 0 : static const TRefElem &rRefElem = Provider<TRefElem>::get();
1141 :
1142 : // Shape function values
1143 : std::vector<number> vShape;
1144 : // Values of shape function gradients
1145 : std::vector<MathVector<TRefElem::dim>> vShapeGrad;
1146 : // Get values at the intersection point
1147 0 : TrialSpace.shapes(vShape, localIntersection);
1148 0 : TrialSpace.grads(vShapeGrad, localIntersection);
1149 :
1150 0 : size_t num_corners_inters_side = rRefElem.num(dim - 1, side, 0);
1151 :
1152 : // loop corners of side
1153 0 : for (size_t j = 0; j < num_corners_inters_side; ++j)
1154 : {
1155 : // get corner
1156 0 : const size_t co = rRefElem.id(dim-1, side, 0, j);
1157 :
1158 : // write shape
1159 0 : conv_shape(ip, co) = flux * vShape[co];
1160 : }
1161 :
1162 : // Write Derivatives if wanted
1163 : if (computeDeriv)
1164 : {
1165 : // No derivatives yet!
1166 : }
1167 :
1168 : // The shapes do not depend of the diffusion tensor
1169 : }
1170 :
1171 : // we're done
1172 0 : return true;
1173 : }
1174 :
1175 :
1176 : } // end namespace ug
1177 :
1178 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE__ */
|
