Line data Source code
1 : /*
2 : * Copyright (c) 2010-2021: G-CSC, Goethe University Frankfurt
3 : * Authors: Andreas Vogel, Dmitrij Logashenko, Martin Stepniewski
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__DISC_HELPER__FINITE_VOLUME_UTIL__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_UTIL__
35 :
36 : // extern libraries
37 : #include <cmath>
38 : #include <vector>
39 :
40 : // other ug4 modules
41 : #include "common/common.h"
42 :
43 : #include "lib_disc/common/geometry_util.h"
44 : #include "lib_disc/local_finite_element/lagrange/lagrange.h"
45 : #include "lib_disc/spatial_disc/disc_util/fv_geom_base.h"
46 :
47 : namespace ug{
48 :
49 : /// averages positions by arithmetic mean
50 : /**
51 : * Arithmetic Mean of Positions
52 : * returns the arithmetic mean of positions
53 : *
54 : * \param[in] vCornerCoords positions
55 : * \param[in] num number of positions
56 : * \param[out] vOut arithmetic mean of positions
57 : */
58 : template <typename TPosition>
59 0 : void AveragePositions(TPosition& vOut, const TPosition* vCornerCoords, size_t num)
60 : {
61 : vOut = vCornerCoords[0];
62 0 : for(size_t j = 1; j < num; ++j)
63 : {
64 0 : vOut += vCornerCoords[j];
65 : }
66 0 : vOut *= 1./(number)num;
67 0 : }
68 :
69 :
70 :
71 : ////////////////////////////////////////////////////////////////////////////////
72 : // Finite Volume Traits
73 : ////////////////////////////////////////////////////////////////////////////////
74 :
75 : /// Base class, some fields are redefined in the instantiations for particular elements
76 : template <typename TRefElem> struct fv1_traits_most_common
77 : {
78 : /// type of reference element
79 : typedef TRefElem ref_elem_type;
80 :
81 : /// number of SubControlVolumes (for most of the elements - overridden for e.g. ROID_PYRAMID and ROID_OCTAHEDRON)
82 : static const size_t numSCV = ref_elem_type::numCorners;
83 :
84 : /// number of SubControlVolumeFaces (for most of the elements - overridden for e.g. ROID_PYRAMID and ROID_OCTAHEDRON)
85 : static const size_t numSCVF = ref_elem_type::numEdges;
86 :
87 : /// returns the 'from' and 'to' corner indices for a scvf
88 : static size_t scvf_from_to
89 : (
90 : const ref_elem_type& refElem, ///< the reference element object
91 : size_t i, ///< index of the scvf
92 : size_t ft ///< 0 = from, 1 = to
93 : )
94 : {
95 0 : return refElem.id(1, i, 0, ft);
96 : }
97 :
98 : /// returns the node id for a scv
99 : static size_t scv_node_id
100 : (
101 : const ref_elem_type& refElem, ///< the reference element object
102 : size_t i ///< index of the scv
103 : )
104 : {
105 : return i;
106 : }
107 :
108 : /// Calculates array vMidID for scvf i
109 : static void scvf_mid_id
110 : (
111 : const ref_elem_type& refElem, ///< the reference element object
112 : MidID *vMidID, ///< array vMidID for scvf i
113 : size_t i ///< index of the scvf
114 : )
115 : {
116 : static const int dim = TRefElem::dim;
117 :
118 : // set mid ids:
119 :
120 : // start at edge midpoint
121 : vMidID[0] = MidID(1,i);
122 :
123 : // loop up dimension
124 : if(dim == 2)
125 : {
126 : vMidID[1] = MidID(dim, 0); // center of element
127 : }
128 : else if (dim == 3)
129 : {
130 : vMidID[1] = MidID(2, refElem.id(1, i, 2, 0)); // side 0
131 : vMidID[2] = MidID(dim, 0); // center of element
132 : vMidID[3] = MidID(2, refElem.id(1, i, 2, 1)); // side 1
133 : }
134 : }
135 :
136 : /// Calculates array vMidID for scv i
137 : static void scv_mid_id
138 : (
139 : const ref_elem_type& refElem, ///< the reference element object
140 : MidID *vMidID, ///< array vMidID for scvf i
141 : size_t i ///< index of the scvf
142 : )
143 : {
144 : static const int dim = TRefElem::dim;
145 :
146 : // set mid ids:
147 :
148 : if(dim == 1)
149 : {
150 : vMidID[0] = MidID(0, i); // set node as corner of scv
151 : vMidID[1] = MidID(dim, 0); // center of element
152 : }
153 : else if(dim == 2)
154 : {
155 : vMidID[0] = MidID(0, i); // set node as corner of scv
156 : vMidID[1] = MidID(1, refElem.id(0, i, 1, 0)); // edge 1
157 : vMidID[2] = MidID(dim, 0); // center of element
158 : vMidID[3] = MidID(1, refElem.id(0, i, 1, 1)); // edge 2
159 : }
160 : else if(dim == 3)
161 : {
162 : vMidID[0] = MidID(0, i); // set node as corner of scv
163 : vMidID[1] = MidID(1, refElem.id(0, i, 1, 1)); // edge 1
164 : vMidID[2] = MidID(2, refElem.id(0, i, 2, 0)); // face 0
165 : vMidID[3] = MidID(1, refElem.id(0, i, 1, 0)); // edge 0
166 : vMidID[4] = MidID(1, refElem.id(0, i, 1, 2)); // edge 2
167 : vMidID[5] = MidID(2, refElem.id(0, i, 2, 2)); // face 2
168 : vMidID[6] = MidID(dim, 0); // center of element
169 : vMidID[7] = MidID(2, refElem.id(0, i, 2, 1)); // face 1
170 : }
171 : else {UG_THROW("Dimension higher than 3 not implemented.");}
172 : }
173 : };
174 :
175 : /// Traits for Finite Volumes (dummy implementation, s. the instantiations below)
176 : template <typename TRefElem, int TWorldDim> struct fv1_traits
177 : : public fv1_traits_most_common<TRefElem>
178 : {
179 : // maximum for dimension
180 : static const size_t maxNumSCVF;
181 : static const size_t maxNumSCV;
182 : static const size_t maxNSH;
183 :
184 : // number of corners of scvf
185 : const static size_t NumCornersOfSCVF;
186 :
187 : // maximum of corners of scv
188 : const static size_t NumCornersOfSCV;
189 :
190 : // maximum of corners of bf
191 : const static size_t NumCornersOfBF;
192 :
193 : // computes the normal to a scvf
194 : static void NormalOnSCVF(MathVector<TWorldDim>& outNormal,
195 : const MathVector<TWorldDim>* vSCVFCorner,
196 : const MathVector<TWorldDim>* vElemCorner);
197 :
198 : // computes the normal to a bf
199 : static void NormalOnBF(MathVector<TWorldDim>& outNormal,
200 : const MathVector<TWorldDim>* vSCVFCorner,
201 : const MathVector<TWorldDim>* vElemCorner);
202 :
203 : // types of scv and scvf and bf
204 : typedef void scv_type;
205 : typedef void scvf_type;
206 : typedef void bf_type;
207 : };
208 :
209 : /////////////////////////
210 : // 1D Reference Element
211 : /////////////////////////
212 :
213 : struct fv1_traits_ReferenceEdge
214 : {
215 : static const size_t maxNumSCVF = 1;
216 : static const size_t maxNumSCV = 2;
217 : static const size_t maxNSH = maxNumSCV;
218 :
219 : const static size_t NumCornersOfSCVF = 1;
220 : const static size_t NumCornersOfSCV = 2;
221 : const static size_t NumCornersOfBF = 1;
222 :
223 : typedef ReferenceEdge scv_type;
224 : typedef ReferenceVertex scvf_type;
225 : typedef ReferenceVertex bf_type;
226 : };
227 :
228 : template <> struct fv1_traits<ReferenceEdge, 1>
229 : : public fv1_traits_ReferenceEdge,
230 : public fv1_traits_most_common<ReferenceEdge>
231 : {
232 : static void NormalOnSCVF(MathVector<1>& outNormal,
233 : const MathVector<1>* vSCVFCorner,
234 : const MathVector<1>* vElemCorner)
235 : {
236 : ElementNormal<ReferenceVertex, 1>(outNormal, vSCVFCorner);
237 0 : VecNormalize(outNormal, outNormal);
238 : }
239 : static void NormalOnBF(MathVector<1>& outNormal,
240 : const MathVector<1>* vSCVFCorner,
241 : const MathVector<1>* vElemCorner)
242 : {
243 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
244 : }
245 : };
246 :
247 : template <> struct fv1_traits<ReferenceEdge, 2>
248 : : public fv1_traits_ReferenceEdge,
249 : public fv1_traits_most_common<ReferenceEdge>
250 : {
251 0 : static void NormalOnSCVF(MathVector<2>& outNormal,
252 : const MathVector<2>* vSCVFCorner,
253 : const MathVector<2>* vElemCorner)
254 : {
255 : VecSubtract(outNormal, vElemCorner[1], vElemCorner[0]);
256 0 : VecNormalize(outNormal, outNormal);
257 0 : }
258 : static void NormalOnBF(MathVector<2>& outNormal,
259 : const MathVector<2>* vSCVFCorner,
260 : const MathVector<2>* vElemCorner)
261 : {
262 0 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
263 : }
264 : };
265 :
266 : template <> struct fv1_traits<ReferenceEdge, 3>
267 : : public fv1_traits_ReferenceEdge,
268 : public fv1_traits_most_common<ReferenceEdge>
269 : {
270 0 : static void NormalOnSCVF(MathVector<3>& outNormal,
271 : const MathVector<3>* vSCVFCorner,
272 : const MathVector<3>* vElemCorner)
273 : {
274 : VecSubtract(outNormal, vElemCorner[1], vElemCorner[0]);
275 0 : VecNormalize(outNormal, outNormal);
276 0 : }
277 : static void NormalOnBF(MathVector<3>& outNormal,
278 : const MathVector<3>* vSCVFCorner,
279 : const MathVector<3>* vElemCorner)
280 : {
281 0 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
282 : }
283 : };
284 :
285 : /////////////////////////
286 : // 2D Reference Element
287 : /////////////////////////
288 :
289 : struct fv1_traits_ReferenceFace
290 : {
291 : static const size_t maxNumSCVF = 4;
292 : static const size_t maxNumSCV = 4;
293 : static const size_t maxNSH = maxNumSCV;
294 :
295 : const static size_t NumCornersOfSCVF = 2;
296 : const static size_t NumCornersOfSCV = 4;
297 : const static size_t NumCornersOfBF = 2;
298 :
299 : typedef ReferenceQuadrilateral scv_type;
300 : typedef ReferenceEdge scvf_type;
301 : typedef ReferenceEdge bf_type;
302 : };
303 :
304 : struct fv1_traits_ReferenceFace2d : public fv1_traits_ReferenceFace
305 : {
306 : static void NormalOnSCVF(MathVector<2>& outNormal,
307 : const MathVector<2>* vSCVFCorner,
308 : const MathVector<2>* vElemCorner)
309 : {ElementNormal<ReferenceEdge, 2>(outNormal, vSCVFCorner);}
310 : static void NormalOnBF(MathVector<2>& outNormal,
311 : const MathVector<2>* vSCVFCorner,
312 : const MathVector<2>* vElemCorner)
313 : {
314 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
315 : }
316 : };
317 :
318 : struct fv1_traits_ReferenceFace3d : public fv1_traits_ReferenceFace
319 : {
320 : template <typename TElem>
321 0 : static void NormalOnSCVF_Face(MathVector<3>& outNormal,
322 : const MathVector<3>* vSCVFCorner,
323 : const MathVector<3>* vElemCorner)
324 : {
325 : // compute Normal to Face (right-handed)
326 : MathVector<3> ElemNormal;
327 0 : ElementNormal<TElem,3>(ElemNormal, vElemCorner);
328 :
329 : // compute a Point such that the triangle given by vSCVFCorner and p
330 : // has the following property:
331 : // a) The new Triangle is normal to the ElementFace
332 : // b) The new Triangle contains the scvf
333 : MathVector<3> p;
334 : VecAdd(p, vSCVFCorner[0], ElemNormal);
335 :
336 0 : MathVector<3> vNewTriangleCorner[3];
337 : vNewTriangleCorner[0] = vSCVFCorner[0];
338 : vNewTriangleCorner[1] = vSCVFCorner[1];
339 : vNewTriangleCorner[2] = p;
340 :
341 : // now compute Normal to new triangle. This Normal will be normal to
342 : // the scvf and within the original element
343 0 : ElementNormal<ReferenceTriangle, 3>(outNormal, vNewTriangleCorner);
344 :
345 : // scale to normal to the size of the scvf
346 0 : const number size = VecDistance(vSCVFCorner[0], vSCVFCorner[1]);
347 0 : VecNormalize(outNormal, outNormal);
348 : VecScale(outNormal, outNormal, size);
349 0 : }
350 : static void NormalOnSCVF(MathVector<3>& outNormal,
351 : const MathVector<3>* vSCVFCorner,
352 : const MathVector<3>* vElemCorner)
353 : {
354 : UG_THROW("not implemented.")
355 : }
356 0 : static void NormalOnBF(MathVector<3>& outNormal,
357 : const MathVector<3>* vSCVFCorner,
358 : const MathVector<3>* vElemCorner)
359 : {
360 0 : UG_THROW("not implemented.")
361 : }
362 : };
363 :
364 : template <> struct fv1_traits<ReferenceTriangle, 2>
365 : : public fv1_traits_ReferenceFace2d,
366 : public fv1_traits_most_common<ReferenceTriangle>
367 : {};
368 : template <> struct fv1_traits<ReferenceTriangle, 3>
369 : : public fv1_traits_ReferenceFace,
370 : public fv1_traits_most_common<ReferenceTriangle>
371 : {
372 : static void NormalOnSCVF(MathVector<3>& outNormal,
373 : const MathVector<3>* vSCVFCorner,
374 : const MathVector<3>* vElemCorner)
375 0 : {fv1_traits_ReferenceFace3d::NormalOnSCVF_Face<ReferenceTriangle>(outNormal, vSCVFCorner, vElemCorner);}
376 : static void NormalOnBF(MathVector<3>& outNormal,
377 : const MathVector<3>* vSCVFCorner,
378 : const MathVector<3>* vElemCorner)
379 : {
380 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
381 : }
382 : };
383 :
384 : template <> struct fv1_traits<ReferenceQuadrilateral, 2>
385 : : public fv1_traits_ReferenceFace2d,
386 : public fv1_traits_most_common<ReferenceQuadrilateral>
387 : {};
388 : template <> struct fv1_traits<ReferenceQuadrilateral, 3>
389 : : public fv1_traits_ReferenceFace,
390 : public fv1_traits_most_common<ReferenceQuadrilateral>
391 : {
392 : static void NormalOnSCVF(MathVector<3>& outNormal,
393 : const MathVector<3>* vSCVFCorner,
394 : const MathVector<3>* vElemCorner)
395 0 : {fv1_traits_ReferenceFace3d::NormalOnSCVF_Face<ReferenceQuadrilateral>(outNormal, vSCVFCorner, vElemCorner);}
396 : static void NormalOnBF(MathVector<3>& outNormal,
397 : const MathVector<3>* vSCVFCorner,
398 : const MathVector<3>* vElemCorner)
399 : {
400 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
401 : }
402 : };
403 :
404 : /////////////////////////
405 : // 3D Reference Element
406 : /////////////////////////
407 :
408 : struct fv1_traits_ReferenceVolume
409 : {
410 : static const size_t maxNumSCVF = 24;
411 : static const size_t maxNumSCV = 32;
412 : static const size_t maxNSH = 8;
413 :
414 : const static size_t NumCornersOfSCVF = 4;
415 : const static size_t NumCornersOfSCV = 8;
416 : const static size_t NumCornersOfBF = 4;
417 :
418 : typedef ReferenceHexahedron scv_type;
419 : typedef ReferenceQuadrilateral scvf_type;
420 : typedef ReferenceQuadrilateral bf_type;
421 :
422 : static void NormalOnSCVF(MathVector<3>& outNormal,
423 : const MathVector<3>* vSCVFCorner,
424 : const MathVector<3>* vElemCorner)
425 0 : {ElementNormal<ReferenceQuadrilateral, 3>(outNormal, vSCVFCorner);}
426 : static void NormalOnBF(MathVector<3>& outNormal,
427 : const MathVector<3>* vSCVFCorner,
428 : const MathVector<3>* vElemCorner)
429 : {
430 : NormalOnSCVF(outNormal, vSCVFCorner, vElemCorner);
431 : }
432 : };
433 :
434 : template <> struct fv1_traits<ReferenceTetrahedron, 3>
435 : : public fv1_traits_ReferenceVolume,
436 : public fv1_traits_most_common<ReferenceTetrahedron>
437 : {};
438 : template <> struct fv1_traits<ReferencePrism, 3>
439 : : public fv1_traits_ReferenceVolume,
440 : public fv1_traits_most_common<ReferencePrism>
441 : {};
442 : template <> struct fv1_traits<ReferenceHexahedron, 3>
443 : : public fv1_traits_ReferenceVolume,
444 : public fv1_traits_most_common<ReferenceHexahedron>
445 : {};
446 :
447 : /// Pyramids: dimension-independent part of the FV1 traits
448 : /**
449 : * For Pyramids we use triangular scvf, since the quadrilateral scvf would not be
450 : * flat in general by the positions where its corners are placed
451 : */
452 : struct fv1_traits_ReferencePyramid
453 : : public fv1_traits_ReferenceVolume
454 : // Remark: Pyramid is a special case, fv1_traits_most_common<ReferencePyramid> is not inherited here!
455 : {
456 : static const size_t numSCV = 4 * ReferencePyramid::numEdges; ///< overridden field from fv1_traits_most_common
457 :
458 : static const size_t numSCVF = 2 * ReferencePyramid::numEdges; ///< overridden field from fv1_traits_most_common
459 :
460 : const static size_t NumCornersOfSCVF = 3; // triangles
461 : const static size_t NumCornersOfSCV = 4; // tetrahedrons
462 : const static size_t NumCornersOfBF = 4; // quadrilaterals
463 :
464 : typedef ReferenceTetrahedron scv_type;
465 : typedef ReferenceTriangle scvf_type;
466 : typedef ReferenceQuadrilateral bf_type;
467 :
468 : /// returns the 'from' and 'to' corner indices for a scvf (overridden function from fv1_traits_most_common)
469 : static size_t scvf_from_to
470 : (
471 : const ReferencePyramid& refElem, ///< the reference element object
472 : size_t i, ///< index of the scvf
473 : size_t ft ///< 0 = from, 1 = to
474 : )
475 : { // map according to the order defined in ComputeSCVFMidID
476 0 : return refElem.id(1, i/2, 0, ft);
477 : }
478 :
479 : /// returns the node id for a scv (overridden function from fv1_traits_most_common)
480 : static size_t scv_node_id
481 : (
482 : const ReferencePyramid& refElem, ///< the reference element object
483 : size_t i ///< index of the scv
484 : )
485 : { // map according to order defined in ComputeSCVMidID
486 0 : if(i%2 == 0){
487 0 : return refElem.id(1, i/4, 0, 0); // from
488 : } else {
489 0 : return refElem.id(1, i/4, 0, 1); // to
490 : }
491 : }
492 :
493 : /// Calculates array vMidID for scvf i
494 : static void scvf_mid_id
495 : (
496 : const ReferencePyramid& refElem, ///< the reference element object
497 : MidID *vMidID, ///< array vMidID for scvf i
498 : size_t i ///< index of the scvf
499 : )
500 : {
501 : static const int dim = 3;
502 :
503 : // set mid ids:
504 :
505 : // start at edge midpoint
506 0 : vMidID[0] = MidID(1,i/2);
507 :
508 : // there are 2 scvf per edge
509 0 : if(i%2 == 0){
510 0 : vMidID[1] = MidID(2, refElem.id(1, i/2, 2, 0)); // side 0
511 0 : vMidID[2] = MidID(dim, 0); // center of element
512 : } else {
513 0 : vMidID[1] = MidID(dim, 0); // center of element
514 0 : vMidID[2] = MidID(2, refElem.id(1, i/2, 2, 1)); // side 1
515 : }
516 : }
517 :
518 : /// Calculates array vMidID for scv i
519 0 : static void scv_mid_id
520 : (
521 : const ReferencePyramid& refElem, ///< the reference element object
522 : MidID *vMidID, ///< array vMidID for scvf i
523 : size_t i ///< index of the scvf
524 : )
525 : {
526 : static const int dim = 3;
527 :
528 : // set mid ids:
529 :
530 : // start at edge midpoint
531 0 : vMidID[3] = MidID(1,i/4);
532 :
533 : // there are 2 scvf per edge
534 0 : if(i%4 == 0 || i%4 == 1){
535 0 : vMidID[1] = MidID(2, refElem.id(1, i/4, 2, 0)); // side 0
536 0 : vMidID[2] = MidID(dim, 0); // center of element
537 : } else {
538 0 : vMidID[1] = MidID(dim, 0); // center of element
539 0 : vMidID[2] = MidID(2, refElem.id(1, i/4, 2, 1)); // side 1
540 : }
541 :
542 : // connect to from / to corners of edge
543 0 : if(i%2 == 0){
544 0 : vMidID[0] = MidID(0, refElem.id(1, i/4, 0, 0)); // from
545 : } else {
546 0 : vMidID[0] = MidID(0, refElem.id(1, i/4, 0, 1)); // to
547 : }
548 0 : }
549 : };
550 : /// Pyramids: the FV1 traits
551 : /**
552 : * For Pyramids we use triangular scvf, since the quadrilateral scvf would not be
553 : * flat in general by the positions where its corners are placed
554 : */
555 : template <> struct fv1_traits<ReferencePyramid, 3> : public fv1_traits_ReferencePyramid
556 : {
557 : static void NormalOnSCVF(MathVector<3>& outNormal,
558 : const MathVector<3>* vSCVFCorner,
559 : const MathVector<3>* vElemCorner)
560 0 : {ElementNormal<ReferenceTriangle, 3>(outNormal, vSCVFCorner);}
561 : static void NormalOnBF(MathVector<3>& outNormal,
562 : const MathVector<3>* vSCVFCorner,
563 : const MathVector<3>* vElemCorner)
564 : {
565 0 : ElementNormal<ReferenceQuadrilateral, 3>(outNormal, vSCVFCorner);
566 : }
567 : };
568 :
569 : /// Octahedra: dimension-independent part of the FV1 traits
570 : struct fv1_traits_ReferenceOctahedron
571 : : public fv1_traits_ReferenceVolume
572 : // Remark: Octahedron is a special case, fv1_traits_most_common<ReferenceOctahedron> is not inherited here!
573 : {
574 : /**
575 : * The octahedral reference element contains an implicit interior
576 : * substructure that is constructed by several geometric objects,
577 : * i.e. imaginary subfaces (8 triangles), subedges (2) and subvolumes (4 tetrahedra)
578 : * resulting from the division into 4 subtetrahedra alongside inner edge 3->1.
579 : * The dual fv1 geometry consists of the original hexahedral SCVs in each of the
580 : * 4 subtetrahedra.
581 : */
582 : static const size_t numSCV = 16; ///< overridden field from fv1_traits_most_common
583 :
584 : static const size_t numSCVF = 24; ///< overridden field from fv1_traits_most_common
585 : // Remark: Special case for octahedron, scvf not mappable by edges.
586 :
587 : // maximum dimension of substructure objects
588 : enum{MAXDIM = 3};
589 :
590 : // maximum number of substructure objects
591 : enum{MAXSUBSTRUCTOBJECTS = 8};
592 :
593 : // maximum number of substructure corners
594 : enum{MAXSUBSTRUCTCORNERS = 4};
595 :
596 : /// returns the id of corner j of obj i in dimension dim_i of the implicit interior substructure
597 : /**
598 : * The octahedral reference element contains an implicit interior
599 : * substructure that is constructed by several geometric objects, that
600 : * are mapped by a reference element by themselves. This method returns the
601 : * id (w.r.t. this reference element) of a sub-geometric object that is
602 : * part of a sub-geometric object of the implicit interior substructure of
603 : * this reference element.
604 : *
605 : * \param[in] dim_i dimension of sub geometric object
606 : * \param[in] i id of sub geometric object
607 : * \param[in] j number of obj contained in the sub-object
608 : * \returns id of the j'th corner that is
609 : * contained in the i*th (sub-)geom object of dimension dim_i
610 : */
611 0 : static int substruct_coID(int dim_i, size_t i, size_t j)
612 : {
613 : // Corner indices of implicit interior substructure Geometric Objects
614 : static int substruct_coID[MAXDIM+1][MAXSUBSTRUCTOBJECTS][MAXSUBSTRUCTCORNERS];
615 :
616 : // subedge 0 = (3,1)
617 0 : substruct_coID[EDGE][0][0] = 3;
618 0 : substruct_coID[EDGE][0][1] = 1;
619 : // subedge 1 = (1,3)
620 0 : substruct_coID[EDGE][1][0] = 1;
621 0 : substruct_coID[EDGE][1][1] = 3;
622 :
623 : // subface 0 = (1,2,3)
624 0 : substruct_coID[FACE][0][0] = 1;
625 0 : substruct_coID[FACE][0][1] = 2;
626 0 : substruct_coID[FACE][0][2] = 3;
627 : // subface 1 = (1,3,2)
628 0 : substruct_coID[FACE][1][0] = 1;
629 0 : substruct_coID[FACE][1][1] = 3;
630 0 : substruct_coID[FACE][1][2] = 2;
631 : // subface 2 = (1,3,4)
632 0 : substruct_coID[FACE][2][0] = 1;
633 0 : substruct_coID[FACE][2][1] = 3;
634 0 : substruct_coID[FACE][2][2] = 4;
635 : // subface 3 = (1,4,3)
636 0 : substruct_coID[FACE][3][0] = 1;
637 0 : substruct_coID[FACE][3][1] = 4;
638 0 : substruct_coID[FACE][3][2] = 3;
639 : // subface 4 = (1,3,5)
640 0 : substruct_coID[FACE][4][0] = 1;
641 0 : substruct_coID[FACE][4][1] = 3;
642 0 : substruct_coID[FACE][4][2] = 5;
643 : // subface 5 = (1,5,3)
644 0 : substruct_coID[FACE][5][0] = 1;
645 0 : substruct_coID[FACE][5][1] = 5;
646 0 : substruct_coID[FACE][5][2] = 3;
647 : // subface 6 = (1,0,3)
648 0 : substruct_coID[FACE][6][0] = 1;
649 0 : substruct_coID[FACE][6][1] = 0;
650 0 : substruct_coID[FACE][6][2] = 3;
651 : // subface 7 = (1,0,3)
652 0 : substruct_coID[FACE][7][0] = 1;
653 0 : substruct_coID[FACE][7][1] = 3;
654 0 : substruct_coID[FACE][7][2] = 0;
655 :
656 : // subvolume 0 = (1,2,3,5)
657 0 : substruct_coID[VOLUME][0][0] = 1;
658 0 : substruct_coID[VOLUME][0][1] = 2;
659 0 : substruct_coID[VOLUME][0][2] = 3;
660 0 : substruct_coID[VOLUME][0][3] = 5;
661 : // subvolume 1 = (1,3,4,5)
662 0 : substruct_coID[VOLUME][1][0] = 1;
663 0 : substruct_coID[VOLUME][1][1] = 3;
664 0 : substruct_coID[VOLUME][1][2] = 4;
665 0 : substruct_coID[VOLUME][1][3] = 5;
666 : // subvolume 2 = (1,2,3,0)
667 0 : substruct_coID[VOLUME][2][0] = 1;
668 0 : substruct_coID[VOLUME][2][1] = 2;
669 0 : substruct_coID[VOLUME][2][2] = 3;
670 0 : substruct_coID[VOLUME][2][3] = 0;
671 : // subvolume 3 = (1,3,4,0)
672 0 : substruct_coID[VOLUME][3][0] = 1;
673 0 : substruct_coID[VOLUME][3][1] = 3;
674 0 : substruct_coID[VOLUME][3][2] = 4;
675 0 : substruct_coID[VOLUME][3][3] = 0;
676 :
677 0 : return substruct_coID[dim_i][i][j];
678 : }
679 :
680 : /// returns the number of implicit interior substructure geometric objects of dim
681 : /**
682 : * The octahedral reference element contains an implicit interior
683 : * substructure that is constructed by several geometric objects, that
684 : * are mapped by a reference element by themselves. This method returns how
685 : * many (sub-)geometric objects of a given dimension are contained in the
686 : * implicit interior substructure of this reference element.
687 : *
688 : * \param[in] dim dimension
689 : * \returns number of objects of the dimension contained in the ref elem
690 : */
691 : static size_t substruct_num(int dim)
692 : {
693 : // number of interior substructure Geometric Objects
694 : size_t vSubStructNum[MAXDIM+1];
695 :
696 : vSubStructNum[VERTEX] = 0; // no additional vertices in the substructure
697 0 : vSubStructNum[EDGE] = 2;
698 0 : vSubStructNum[FACE] = 8;
699 0 : vSubStructNum[VOLUME] = 4;
700 :
701 0 : return vSubStructNum[dim];
702 : }
703 :
704 : /// returns the number of objects of dim for a sub-geometric object of the implicit interior substructure
705 : /**
706 : * The octahedral reference element contains an implicit interior
707 : * substructure that is constructed by several geometric objects, that
708 : * are mapped by a reference element by themselves. This method returns how
709 : * many (sub-)geometric objects of a given dimension are contained in the
710 : * (sub-)geometric object of the implicit interior substructure of
711 : * this reference element.
712 : *
713 : * \param[in] dim_i dimension of sub geometric object
714 : * \param[in] i number of sub geometric object
715 : * \param[in] dim_j dimension for elems contained in the sub-object
716 : * \returns number of objects of the dimension dim_j that are
717 : * contained in the i*th (sub-)geom object of dimension dim_i
718 : */
719 0 : static size_t substruct_num(int dim_i, size_t i, int dim_j)
720 : {
721 : // number of interior substructure Geometric Objects
722 : size_t vSubStructSubNum[MAXDIM+1][MAXSUBSTRUCTOBJECTS][MAXDIM+1];
723 :
724 0 : for(size_t i = 0; i < substruct_num(VOLUME); ++i)
725 : {
726 0 : vSubStructSubNum[VOLUME][i][VERTEX] = 4;
727 0 : vSubStructSubNum[VOLUME][i][EDGE] = 6;
728 0 : vSubStructSubNum[VOLUME][i][FACE] = 4;
729 0 : vSubStructSubNum[VOLUME][i][VOLUME] = 1;
730 : }
731 :
732 0 : for(size_t i = 0; i < substruct_num(FACE); ++i)
733 : {
734 0 : vSubStructSubNum[FACE][i][VERTEX] = 3;
735 0 : vSubStructSubNum[FACE][i][EDGE] = 3;
736 0 : vSubStructSubNum[FACE][i][FACE] = 1;
737 0 : vSubStructSubNum[FACE][i][VOLUME] = 1;
738 : }
739 :
740 0 : for(size_t i = 0; i < substruct_num(EDGE); ++i)
741 : {
742 0 : vSubStructSubNum[EDGE][i][VERTEX] = 2;
743 0 : vSubStructSubNum[EDGE][i][EDGE] = 1;
744 0 : vSubStructSubNum[EDGE][i][FACE] = 2;
745 0 : vSubStructSubNum[EDGE][i][VOLUME] = 1;
746 : }
747 :
748 : for(size_t i = 0; i < substruct_num(VERTEX); ++i)
749 : {
750 : vSubStructSubNum[VERTEX][i][VERTEX] = 1;
751 : vSubStructSubNum[VERTEX][i][EDGE] = 3;
752 : vSubStructSubNum[VERTEX][i][FACE] = 3;
753 : vSubStructSubNum[VERTEX][i][VOLUME] = 1;
754 : }
755 :
756 0 : return vSubStructSubNum[dim_i][i][dim_j];
757 : }
758 :
759 : /// returns the 'from' and 'to' corner indices for a scvf (overridden function from fv1_traits_most_common)
760 : static size_t scvf_from_to
761 : (
762 : const ReferenceOctahedron& refElem, ///< the reference element object
763 : size_t i, ///< index of the scvf
764 : size_t ft ///< 0 = from, 1 = to
765 : )
766 : { // map according to the order defined in ComputeSCVFMidID
767 : static size_t from_to_ind [24][2] =
768 : {
769 : {1, 2}, // face 0
770 : {2, 1}, // face 1
771 :
772 : {2, 3}, // face 2
773 : {3, 2}, // face 3
774 :
775 : {3, 1}, // face 4
776 : {1, 3}, // face 5
777 :
778 : {1, 5}, // face 6
779 : {1, 0}, // face 7
780 :
781 : {2, 5}, // face 8
782 : {2, 0}, // face 9
783 :
784 : {3, 5}, // face 10
785 : {3, 0}, // face 11
786 :
787 : {1, 3}, // face 12
788 : {3, 1}, // face 13
789 :
790 : {3, 4}, // face 14
791 : {4, 3}, // face 15
792 :
793 : {4, 1}, // face 16
794 : {1, 4}, // face 17
795 :
796 : {1, 5}, // face 18
797 : {1, 0}, // face 19
798 :
799 : {3, 5}, // face 20
800 : {3, 0}, // face 21
801 :
802 : {4, 5}, // face 22
803 : {4, 0} // face 23
804 : };
805 0 : return from_to_ind [i] [ft];
806 : }
807 :
808 : /// returns the node id for a scv (overridden function from fv1_traits_most_common)
809 : static size_t scv_node_id
810 : (
811 : const ReferenceOctahedron& refElem, ///< the reference element object
812 : size_t i ///< index of the scv
813 : )
814 : { // map according to order defined in ComputeSCVMidID
815 : static size_t node_id [16] =
816 : {
817 : 1, // scv 0
818 : 1, // scv 1
819 :
820 : 2, // scv 2
821 : 2, // scv 3
822 :
823 : 3, // scv 4
824 : 3, // scv 5
825 :
826 : 5, // scv 6
827 : 0, // scv 7
828 :
829 : 1, // scv 8
830 : 1, // scv 9
831 :
832 : 3, // scv 10
833 : 3, // scv 11
834 :
835 : 4, // scv 12
836 : 4, // scv 13
837 :
838 : 5, // scv 14
839 : 0 // scv 15
840 : };
841 0 : return node_id [i];
842 : }
843 :
844 : /// Calculates array vMidID for scvf i
845 0 : static void scvf_mid_id
846 : (
847 : const ReferenceOctahedron& refElem, ///< the reference element object
848 : MidID *vMidID, ///< array vMidID for scvf i
849 : size_t i ///< index of the scvf
850 : )
851 : {
852 0 : switch (i)
853 : {
854 : // scvf of edge 4 (top)
855 0 : case 0: vMidID[0] = MidID(1,4); // edge 4
856 0 : vMidID[1] = MidID(2,8); // subface 0 = 1,2,3
857 0 : vMidID[2] = MidID(3,1); // subvolume 0
858 0 : vMidID[3] = MidID(2,4); // face 4
859 0 : break;
860 : // scvf of edge 4 (bottom)
861 0 : case 1: vMidID[0] = MidID(1,4); // edge 4
862 0 : vMidID[1] = MidID(2,9); // subface 1 = 1,3,2
863 0 : vMidID[2] = MidID(3,3); // subvolume 2
864 0 : vMidID[3] = MidID(2,0); // face 0
865 0 : break;
866 :
867 :
868 : // scvf of edge 5 (top)
869 0 : case 2: vMidID[0] = MidID(1,5); // edge 5
870 0 : vMidID[1] = MidID(2,8); // subface 0 = 1,2,3
871 0 : vMidID[2] = MidID(3,1); // subvolume 0
872 0 : vMidID[3] = MidID(2,5); // face 5
873 0 : break;
874 : // scvf of edge 5 (bottom)
875 0 : case 3: vMidID[0] = MidID(1,5); // edge 5
876 0 : vMidID[1] = MidID(2,9); // subface 1 = 1,3,2
877 0 : vMidID[2] = MidID(3,3); // subvolume 2
878 0 : vMidID[3] = MidID(2,1); // face 1
879 0 : break;
880 :
881 :
882 : // scvf of diagonal 3->1 (top) in subvolume 0
883 0 : case 4: vMidID[0] = MidID(1,12);// diagonal 3->1
884 0 : vMidID[1] = MidID(2,8); // subface 0 = 1,2,3
885 0 : vMidID[2] = MidID(3,1); // subvolume 0
886 0 : vMidID[3] = MidID(2,13);// face 1,5,3
887 0 : break;
888 : // scvf of diagonal 1->3 (bottom) in subvolume 2
889 0 : case 5: vMidID[0] = MidID(1,13);// diagonal 1->3
890 0 : vMidID[1] = MidID(2,9); // subface 1 = 1,3,2
891 0 : vMidID[2] = MidID(3,3); // subvolume 2
892 0 : vMidID[3] = MidID(2,14);// face 1,0,3
893 0 : break;
894 :
895 :
896 : // scvf of edge 8 in subvolume 0
897 0 : case 6: vMidID[0] = MidID(1,8); // edge 8
898 0 : vMidID[1] = MidID(2,4); // face 4
899 0 : vMidID[2] = MidID(3,1); // subvolume 0
900 0 : vMidID[3] = MidID(2,13);// face 1,5,3
901 0 : break;
902 : // scvf of edge 0 in subvolume 2
903 0 : case 7: vMidID[0] = MidID(1,0); // edge 0
904 0 : vMidID[1] = MidID(2,15);// face 1,3,0
905 0 : vMidID[2] = MidID(3,3); // subvolume 2
906 0 : vMidID[3] = MidID(2,0); // face 0
907 0 : break;
908 :
909 :
910 : // scvf of edge 9 in subvolume 0
911 0 : case 8: vMidID[0] = MidID(1,9); // edge 9
912 0 : vMidID[1] = MidID(2,5); // face 5
913 0 : vMidID[2] = MidID(3,1); // subvolume 0
914 0 : vMidID[3] = MidID(2,4); // face 4
915 0 : break;
916 : // scvf of edge 1 in subvolume 2
917 0 : case 9: vMidID[0] = MidID(1,1); // edge 1
918 0 : vMidID[1] = MidID(2,0); // face 0
919 0 : vMidID[2] = MidID(3,3); // subvolume 2
920 0 : vMidID[3] = MidID(2,1); // face 1
921 0 : break;
922 :
923 :
924 : // scvf of edge 10 in subvolume 0
925 0 : case 10:vMidID[0] = MidID(1,10);// edge 10
926 0 : vMidID[1] = MidID(2,12);// face 1,3,5
927 0 : vMidID[2] = MidID(3,1); // subvolume 0
928 0 : vMidID[3] = MidID(2,5); // face 5
929 0 : break;
930 : // scvf of edge 2 in subvolume 2
931 0 : case 11:vMidID[0] = MidID(1,2); // edge 2
932 0 : vMidID[1] = MidID(2,1); // face 1
933 0 : vMidID[2] = MidID(3,3); // subvolume 2
934 0 : vMidID[3] = MidID(2,14);// face 1,0,3
935 0 : break;
936 :
937 :
938 : // scvf of diagonal 1->3 (top) in subvolume 1
939 0 : case 12:vMidID[0] = MidID(1,13);// diagonal 1->3
940 0 : vMidID[1] = MidID(2,10);// subface 2 = 1,3,4
941 0 : vMidID[2] = MidID(3,2); // subvolume 1
942 0 : vMidID[3] = MidID(2,12);// face 1,3,5
943 0 : break;
944 : // scvf of diagonal 3->1 (bottom) in subvolume 3
945 0 : case 13:vMidID[0] = MidID(1,12);// diagonal 3->1
946 0 : vMidID[1] = MidID(2,11);// subface 3 = 1,4,3
947 0 : vMidID[2] = MidID(3,4); // subvolume 3
948 0 : vMidID[3] = MidID(2,15);// face 1,3,0
949 0 : break;
950 :
951 :
952 : // scvf of edge 6 (top)
953 0 : case 14:vMidID[0] = MidID(1,6); // edge 6
954 0 : vMidID[1] = MidID(2,10);// subface 2 = 1,3,4
955 0 : vMidID[2] = MidID(3,2); // subvolume 1
956 0 : vMidID[3] = MidID(2,6); // face 6
957 0 : break;
958 : // scvf of edge 6 (bottom)
959 0 : case 15:vMidID[0] = MidID(1,6); // edge 6
960 0 : vMidID[1] = MidID(2,11);// subface 3 = 1,4,3
961 0 : vMidID[2] = MidID(3,4); // subvolume 3
962 0 : vMidID[3] = MidID(2,2); // face 2
963 0 : break;
964 :
965 :
966 : // scvf of edge 7 (top)
967 0 : case 16:vMidID[0] = MidID(1,7); // edge 7
968 0 : vMidID[1] = MidID(2,10);// subface 2 = 1,3,4
969 0 : vMidID[2] = MidID(3,2); // subvolume 1
970 0 : vMidID[3] = MidID(2,7); // face 7
971 0 : break;
972 : // scvf of edge 7 (bottom)
973 0 : case 17:vMidID[0] = MidID(1,7); // edge 7
974 0 : vMidID[1] = MidID(2,11);// subface 3 = 1,4,3
975 0 : vMidID[2] = MidID(3,4); // subvolume 3
976 0 : vMidID[3] = MidID(2,3); // face 3
977 0 : break;
978 :
979 :
980 : // scvf of edge 8 in subvolume 1
981 0 : case 18:vMidID[0] = MidID(1,8); // edge 8
982 0 : vMidID[1] = MidID(2,13);// face 1,5,3
983 0 : vMidID[2] = MidID(3,2); // subvolume 1
984 0 : vMidID[3] = MidID(2,7); // face 7
985 0 : break;
986 : // scvf of edge 0 in subvolume 3
987 0 : case 19:vMidID[0] = MidID(1,0); // edge 0
988 0 : vMidID[1] = MidID(2,3); // face 3
989 0 : vMidID[2] = MidID(3,4); // subvolume 3
990 0 : vMidID[3] = MidID(2,15);// face 1,3,0
991 0 : break;
992 :
993 :
994 : // scvf of edge 10 in subvolume 1
995 0 : case 20:vMidID[0] = MidID(1,10);// edge 10
996 0 : vMidID[1] = MidID(2,6); // face 6
997 0 : vMidID[2] = MidID(3,2); // subvolume 1
998 0 : vMidID[3] = MidID(2,12);// face 1,3,5
999 0 : break;
1000 : // scvf of edge 2 in subvolume 3
1001 0 : case 21:vMidID[0] = MidID(1,2); // edge 2
1002 0 : vMidID[1] = MidID(2,14);// face 1,0,3
1003 0 : vMidID[2] = MidID(3,4); // subvolume 3
1004 0 : vMidID[3] = MidID(2,2); // face 2
1005 0 : break;
1006 :
1007 :
1008 : // scvf of edge 11 in subvolume 1
1009 0 : case 22:vMidID[0] = MidID(1,11);// edge 11
1010 0 : vMidID[1] = MidID(2,7); // face 7
1011 0 : vMidID[2] = MidID(3,2); // subvolume 1
1012 0 : vMidID[3] = MidID(2,6); // face 6
1013 0 : break;
1014 : // scvf of edge 3 in subvolume 3
1015 0 : case 23:vMidID[0] = MidID(1,3); // edge 3
1016 0 : vMidID[1] = MidID(2,2); // face 2
1017 0 : vMidID[2] = MidID(3,4); // subvolume 3
1018 0 : vMidID[3] = MidID(2,3); // face 3
1019 0 : break;
1020 :
1021 :
1022 0 : default:UG_THROW("Octahedron only has 24 SCVFs (no. 0-23), but requested no. " << i << ".");
1023 : break;
1024 : }
1025 0 : }
1026 :
1027 : /// Calculates array vMidID for scv i
1028 0 : static void scv_mid_id
1029 : (
1030 : const ReferenceOctahedron& refElem, ///< the reference element object
1031 : MidID *vMidID, ///< array vMidID for scvf i
1032 : size_t i ///< index of the scvf
1033 : )
1034 : {
1035 0 : switch (i)
1036 : {
1037 : // scv of corner 1 in subvolume 0 (top)
1038 0 : case 0: vMidID[0] = MidID(0,1); // corner 1
1039 0 : vMidID[1] = MidID(1,4); // edge 4
1040 0 : vMidID[2] = MidID(2,8); // subface 0 = 1,2,3
1041 0 : vMidID[3] = MidID(1,12);// edge 3->1
1042 0 : vMidID[4] = MidID(1,8); // edge 8
1043 0 : vMidID[5] = MidID(2,4); // face 4
1044 0 : vMidID[6] = MidID(3,1); // subvolume 0
1045 0 : vMidID[7] = MidID(2,13);// face 1,5,3
1046 0 : break;
1047 : // scv of corner 1 in subvolume 2 (bottom)
1048 0 : case 1: vMidID[0] = MidID(0,1); // corner 1
1049 0 : vMidID[1] = MidID(1,13);// edge 1->3
1050 0 : vMidID[2] = MidID(2,9); // subface 1 = 1,3,2
1051 0 : vMidID[3] = MidID(1,4); // edge 4
1052 0 : vMidID[4] = MidID(1,0); // edge 0
1053 0 : vMidID[5] = MidID(2,15);// face 1,3,0
1054 0 : vMidID[6] = MidID(3,3); // subvolume 2
1055 0 : vMidID[7] = MidID(2,0); // face 0
1056 0 : break;
1057 :
1058 :
1059 : // scv of corner 2 in subvolume 0 (top)
1060 0 : case 2: vMidID[0] = MidID(0,2); // corner 2
1061 0 : vMidID[1] = MidID(1,5); // edge 5
1062 0 : vMidID[2] = MidID(2,8); // subface 0 = 1,2,3
1063 0 : vMidID[3] = MidID(1,4); // edge 4
1064 0 : vMidID[4] = MidID(1,9); // edge 9
1065 0 : vMidID[5] = MidID(2,5); // face 5
1066 0 : vMidID[6] = MidID(3,1); // subvolume 0
1067 0 : vMidID[7] = MidID(2,4); // face 4
1068 0 : break;
1069 : // scv of corner 2 in subvolume 2 (bottom)
1070 0 : case 3: vMidID[0] = MidID(0,2); // corner 2
1071 0 : vMidID[1] = MidID(1,4); // edge 4
1072 0 : vMidID[2] = MidID(2,9); // subface 1 = 1,3,2
1073 0 : vMidID[3] = MidID(1,5); // edge 5
1074 0 : vMidID[4] = MidID(1,1); // edge 1
1075 0 : vMidID[5] = MidID(2,0); // face 0
1076 0 : vMidID[6] = MidID(3,3); // subvolume 2
1077 0 : vMidID[7] = MidID(2,1); // face 1
1078 0 : break;
1079 :
1080 :
1081 : // scv of corner 3 in subvolume 0 (top)
1082 0 : case 4: vMidID[0] = MidID(0,3); // corner 3
1083 0 : vMidID[1] = MidID(1,12);// edge 3->1
1084 0 : vMidID[2] = MidID(2,8); // subface 0 = 1,2,3
1085 0 : vMidID[3] = MidID(1,5); // edge 5
1086 0 : vMidID[4] = MidID(1,10);// edge 10
1087 0 : vMidID[5] = MidID(2,13);// face 1,5,3
1088 0 : vMidID[6] = MidID(3,1); // subvolume 0
1089 0 : vMidID[7] = MidID(2,5); // face 5
1090 0 : break;
1091 : // scv of corner 3 in subvolume 2 (bottom)
1092 0 : case 5: vMidID[0] = MidID(0,3); // corner 3
1093 0 : vMidID[1] = MidID(1,5); // edge 5
1094 0 : vMidID[2] = MidID(2,9); // subface 0 = 1,3,2
1095 0 : vMidID[3] = MidID(1,13);// edge 1->3
1096 0 : vMidID[4] = MidID(1,2); // edge 2
1097 0 : vMidID[5] = MidID(2,1); // face 1
1098 0 : vMidID[6] = MidID(3,3); // subvolume 2
1099 0 : vMidID[7] = MidID(2,15);// face 1,3,0
1100 0 : break;
1101 :
1102 :
1103 : // scv of corner 5 in subvolume 0 (top)
1104 0 : case 6: vMidID[0] = MidID(0,5); // corner 5
1105 0 : vMidID[1] = MidID(1,9); // edge 9
1106 0 : vMidID[2] = MidID(2,4); // face 4
1107 0 : vMidID[3] = MidID(1,8); // edge 8
1108 0 : vMidID[4] = MidID(1,10);// edge 10
1109 0 : vMidID[5] = MidID(2,5); // face 5
1110 0 : vMidID[6] = MidID(3,1); // subvolume 0
1111 0 : vMidID[7] = MidID(2,13);// subface 1,5,3
1112 0 : break;
1113 : // scv of corner 0 in subvolume 2 (bottom)
1114 0 : case 7: vMidID[0] = MidID(0,0); // corner 0
1115 0 : vMidID[1] = MidID(1,0); // edge 0
1116 0 : vMidID[2] = MidID(2,0); // face 0
1117 0 : vMidID[3] = MidID(1,1); // edge 1
1118 0 : vMidID[4] = MidID(1,2); // edge 2
1119 0 : vMidID[5] = MidID(2,14);// subface 1,0,3
1120 0 : vMidID[6] = MidID(3,3); // subvolume 2
1121 0 : vMidID[7] = MidID(2,1); // face 1
1122 0 : break;
1123 :
1124 :
1125 : // scv of corner 1 in subvolume 1 (top)
1126 0 : case 8: vMidID[0] = MidID(0,1); // corner 1
1127 0 : vMidID[1] = MidID(1,13);// edge 1->3
1128 0 : vMidID[2] = MidID(2,10);// subface 2 = 1,3,4
1129 0 : vMidID[3] = MidID(1,7); // edge 7
1130 0 : vMidID[4] = MidID(1,8); // edge 8
1131 0 : vMidID[5] = MidID(2,12);// face 1,3,5
1132 0 : vMidID[6] = MidID(3,2); // subvolume 1
1133 0 : vMidID[7] = MidID(2,7); // face 7
1134 0 : break;
1135 : // scv of corner 1 in subvolume 3 (bottom)
1136 0 : case 9: vMidID[0] = MidID(0,1); // corner 1
1137 0 : vMidID[1] = MidID(1,7); // edge 7
1138 0 : vMidID[2] = MidID(2,11);// subface 3 = 1,4,3
1139 0 : vMidID[3] = MidID(1,12);// edge 3->1
1140 0 : vMidID[4] = MidID(1,0); // edge 0
1141 0 : vMidID[5] = MidID(2,3); // face 3
1142 0 : vMidID[6] = MidID(3,4); // subvolume 3
1143 0 : vMidID[7] = MidID(2,14);// face 1,0,3
1144 0 : break;
1145 :
1146 :
1147 : // scv of corner 3 in subvolume 1 (top)
1148 0 : case 10:vMidID[0] = MidID(0,3); // corner 3
1149 0 : vMidID[1] = MidID(1,6); // edge 6
1150 0 : vMidID[2] = MidID(2,10);// subface 2 = 1,3,4
1151 0 : vMidID[3] = MidID(1,13);// edge 1->3
1152 0 : vMidID[4] = MidID(1,10);// edge 10
1153 0 : vMidID[5] = MidID(2,6); // face 6
1154 0 : vMidID[6] = MidID(3,2); // subvolume 1
1155 0 : vMidID[7] = MidID(2,12);// face 1,3,5
1156 0 : break;
1157 : // scv of corner 3 in subvolume 3 (bottom)
1158 0 : case 11:vMidID[0] = MidID(0,3); // corner 3
1159 0 : vMidID[1] = MidID(1,12);// edge 3->1
1160 0 : vMidID[2] = MidID(2,11);// subface 3 = 1,4,3
1161 0 : vMidID[3] = MidID(1,6); // edge 6
1162 0 : vMidID[4] = MidID(1,2); // edge 2
1163 0 : vMidID[5] = MidID(2,14);// face 1,0,3
1164 0 : vMidID[6] = MidID(3,4); // subvolume 3
1165 0 : vMidID[7] = MidID(2,2); // face 2
1166 0 : break;
1167 :
1168 :
1169 : // scv of corner 4 in subvolume 1 (top)
1170 0 : case 12:vMidID[0] = MidID(0,4); // corner 4
1171 0 : vMidID[1] = MidID(1,7); // edge 7
1172 0 : vMidID[2] = MidID(2,10);// subface 2 = 1,3,4
1173 0 : vMidID[3] = MidID(1,6); // edge 6
1174 0 : vMidID[4] = MidID(1,11);// edge 11
1175 0 : vMidID[5] = MidID(2,7); // face 7
1176 0 : vMidID[6] = MidID(3,2); // subvolume 1
1177 0 : vMidID[7] = MidID(2,6); // face 6
1178 0 : break;
1179 : // scv of corner 4 in subvolume 3 (bottom)
1180 0 : case 13:vMidID[0] = MidID(0,4); // corner 4
1181 0 : vMidID[1] = MidID(1,6); // edge 6
1182 0 : vMidID[2] = MidID(2,11);// subface 3 = 1,4,3
1183 0 : vMidID[3] = MidID(1,7); // edge 7
1184 0 : vMidID[4] = MidID(1,3); // edge 3
1185 0 : vMidID[5] = MidID(2,2); // face 2
1186 0 : vMidID[6] = MidID(3,4); // subvolume 3
1187 0 : vMidID[7] = MidID(2,3); // face 3
1188 0 : break;
1189 :
1190 :
1191 : // scv of corner 5 in subvolume 1 (top)
1192 0 : case 14:vMidID[0] = MidID(0,5); // corner 5
1193 0 : vMidID[1] = MidID(1,10);// edge 10
1194 0 : vMidID[2] = MidID(2,12);// subface 1,3,5
1195 0 : vMidID[3] = MidID(1,8); // edge 8
1196 0 : vMidID[4] = MidID(1,11);// edge 11
1197 0 : vMidID[5] = MidID(2,6); // face 6
1198 0 : vMidID[6] = MidID(3,2); // subvolume 1
1199 0 : vMidID[7] = MidID(2,7); // face 7
1200 0 : break;
1201 : // scv of corner 0 in subvolume 3 (bottom)
1202 0 : case 15:vMidID[0] = MidID(0,0); // corner 0
1203 0 : vMidID[1] = MidID(1,0); // edge 0
1204 0 : vMidID[2] = MidID(2,15);// subface 1,3,0
1205 0 : vMidID[3] = MidID(1,2); // edge 2
1206 0 : vMidID[4] = MidID(1,3); // edge 3
1207 0 : vMidID[5] = MidID(2,3); // face 3
1208 0 : vMidID[6] = MidID(3,4); // subvolume 3
1209 0 : vMidID[7] = MidID(2,2); // face 2
1210 0 : break;
1211 :
1212 :
1213 0 : default:UG_THROW("Octahedron only has 16 SCVs (no. 0-15), but requested no. " << i << ".");
1214 : break;
1215 : }
1216 0 : }
1217 : };
1218 : /// Octahedra: the FV1 traits
1219 : template <> struct fv1_traits<ReferenceOctahedron, 3> : public fv1_traits_ReferenceOctahedron {};
1220 :
1221 : ////////////////////////////////////////////////////////////////////////////////
1222 : // Dimension dependent traits DIM FV1
1223 : ////////////////////////////////////////////////////////////////////////////////
1224 :
1225 : /// Base for the Traits for Finite Volumes for a generic element of the fixed dimensionalities
1226 : template <int TDim, int TWorldDim> struct fv1_dim_traits_base
1227 : {
1228 : /// dimension of reference element
1229 : static const int dim = TDim;
1230 :
1231 : /// generic reference element type
1232 : typedef DimReferenceElement<dim> ref_elem_type;
1233 :
1234 : /// returns the number of the SCV
1235 0 : static void dim_get_num_SCV_and_SCVF
1236 : (
1237 : const ref_elem_type& refElem, ///< reference element object
1238 : ReferenceObjectID roid, ///< reference element object id
1239 : size_t& numSCV, ///< to write the number of the SCVs
1240 : size_t& numSCVF ///< to write the number of the SCVFs
1241 : )
1242 : {
1243 0 : if(roid != ROID_PYRAMID && roid != ROID_OCTAHEDRON)
1244 : {
1245 0 : numSCV = refElem.num(0);
1246 0 : numSCVF = refElem.num(1);
1247 : }
1248 0 : else if(dim == 3 && roid == ROID_PYRAMID)
1249 : {
1250 : UG_WARNING("Pyramid Finite Volume Geometry for 1st order currently "
1251 : "implemented in DimFV1Geom EXPERIMENTATLLY. Please contact "
1252 : "Martin Stepniewski or Andreas Vogel if you see this message.")
1253 :
1254 0 : numSCV = 4*refElem.num(1);
1255 0 : numSCVF = 2*refElem.num(1);
1256 : }
1257 : else if(dim == 3 && roid == ROID_OCTAHEDRON)
1258 : {
1259 : // Case octahedron
1260 0 : numSCV = 16;
1261 0 : numSCVF = 24;
1262 : }
1263 : else
1264 0 : UG_THROW ("fv1_dim_traits_base: Unsupported combination of dimension and reference element.");
1265 0 : }
1266 :
1267 : /// returns the 'from' and 'to' corner indices for a scvf
1268 0 : static void get_dim_scvf_from_to
1269 : (
1270 : const ref_elem_type& refElem, ///< reference element object
1271 : ReferenceObjectID roid, ///< reference element object id
1272 : size_t i, ///< index of the scvf
1273 : size_t& From, ///< to write the from-index
1274 : size_t& To ///< to write the to-index
1275 : )
1276 : {
1277 0 : if (roid != ROID_PYRAMID && roid != ROID_OCTAHEDRON)
1278 : {
1279 0 : From = refElem.id(1, i, 0, 0);
1280 0 : To = refElem.id(1, i, 0, 1);
1281 : }
1282 : // special case pyramid (scvf not mappable by edges)
1283 0 : else if (dim == 3 && roid == ROID_PYRAMID)
1284 : {
1285 : // map according to order defined in ComputeSCVFMidID
1286 0 : From = fv1_traits_ReferencePyramid::scvf_from_to ((const ReferencePyramid&) refElem, i, 0);
1287 0 : To = fv1_traits_ReferencePyramid::scvf_from_to ((const ReferencePyramid&) refElem, i, 1);
1288 : }
1289 : // special case octahedron (scvf not mappable by edges)
1290 : else if(dim == 3 && roid == ROID_OCTAHEDRON)
1291 : {
1292 : // map according to order defined in ComputeSCVFMidID
1293 0 : From = fv1_traits_ReferenceOctahedron::scvf_from_to ((const ReferenceOctahedron&) refElem, i, 0);
1294 0 : To = fv1_traits_ReferenceOctahedron::scvf_from_to ((const ReferenceOctahedron&) refElem, i, 1);
1295 : }
1296 : else
1297 0 : UG_THROW ("fv1_dim_traits_base: Unsupported combination of dimension and reference element.");
1298 0 : }
1299 :
1300 : /// returns the node id for a scv
1301 0 : static size_t dim_scv_node_id
1302 : (
1303 : const ref_elem_type& refElem, ///< reference element object
1304 : ReferenceObjectID roid, ///< reference element object id
1305 : size_t i ///< index of the scv
1306 : )
1307 : {
1308 : // "classical" elements
1309 0 : if (roid != ROID_PYRAMID && roid != ROID_OCTAHEDRON)
1310 : return i;
1311 :
1312 : // special case pyramid (scv not mappable by corners)
1313 0 : if(dim == 3 && roid == ROID_PYRAMID)
1314 : return fv1_traits_ReferencePyramid::scv_node_id ((const ReferencePyramid&) refElem, i);
1315 :
1316 : // special case octahedron (scvf not mappable by edges)
1317 : if(dim == 3 && roid == ROID_OCTAHEDRON)
1318 0 : return fv1_traits_ReferenceOctahedron::scv_node_id ((const ReferenceOctahedron&) refElem, i);
1319 :
1320 0 : UG_THROW ("fv1_dim_traits_base: Unsupported combination of dimension and reference element.");
1321 : }
1322 : };
1323 :
1324 : /// Traits for Finite Volumes for a generic element of the fixed dimensionalities
1325 : template <int TDim, int TWorldDim> struct fv1_dim_traits;
1326 :
1327 : template <> struct fv1_dim_traits<1, 1> : public fv1_traits<ReferenceEdge, 1>, public fv1_dim_traits_base<1, 1> {};
1328 : template <> struct fv1_dim_traits<1, 2> : public fv1_traits<ReferenceEdge, 2>, public fv1_dim_traits_base<1, 2> {};
1329 : template <> struct fv1_dim_traits<1, 3> : public fv1_traits<ReferenceEdge, 3>, public fv1_dim_traits_base<1, 3> {};
1330 :
1331 : template <> struct fv1_dim_traits<2, 2> : public fv1_traits_ReferenceFace2d, public fv1_dim_traits_base<2, 2> {};
1332 : template <> struct fv1_dim_traits<2, 3> : public fv1_traits_ReferenceFace3d, public fv1_dim_traits_base<2, 3>
1333 : {
1334 : static void NormalOnSCVF(MathVector<3>& outNormal,
1335 : const MathVector<3>* vSCVFCorner,
1336 : const MathVector<3>* vElemCorner)
1337 : {
1338 : // Little bit dirty, but should be correct:
1339 : // Even if the true element has more than three vertices (quadrilateral),
1340 : // we only need three to compute the direction of the normal ElemNormal in NormalOnSCVF_Face,
1341 : // the norm is not needed!
1342 0 : fv1_traits_ReferenceFace3d::NormalOnSCVF_Face<ReferenceTriangle>(outNormal, vSCVFCorner, vElemCorner);
1343 : //UG_THROW("Not implemented.")
1344 : }
1345 :
1346 : };
1347 :
1348 : template <> struct fv1_dim_traits<3, 3> : public fv1_traits_ReferenceVolume, public fv1_dim_traits_base<3, 3> {};
1349 :
1350 : ////////////////////////////////////////////////////////////////////////////////
1351 : // Hanging Finite Volume Traits
1352 : ////////////////////////////////////////////////////////////////////////////////
1353 :
1354 : /// Traits for hanging finite volume (dummy implementation)
1355 : template <typename TRefElem, int TWorldDim> struct hfv1_traits
1356 : {
1357 : const static size_t NumCornersOfSCVF;
1358 : const static size_t MaxNumCornersOfSCV;
1359 :
1360 : static void NormalOnSCVF(MathVector<TWorldDim>& outNormal, const MathVector<TWorldDim>* vCornerCoords);
1361 :
1362 : typedef void scv_type;
1363 : };
1364 :
1365 : /////////////////////////
1366 : // 1D Reference Element
1367 : /////////////////////////
1368 :
1369 : struct hfv1_traits_ReferenceEdge
1370 : {
1371 : const static size_t NumCornersOfSCVF = 1;
1372 : const static size_t MaxNumCornersOfSCV = 2;
1373 : typedef ReferenceEdge scv_type;
1374 : };
1375 :
1376 : template <> struct hfv1_traits<ReferenceEdge, 1> : public hfv1_traits_ReferenceEdge
1377 : {
1378 : static void NormalOnSCVF(MathVector<1>& outNormal, const MathVector<1>* vCornerCoords)
1379 : {ElementNormal<ReferenceVertex, 1>(outNormal, vCornerCoords);}
1380 : };
1381 :
1382 : template <> struct hfv1_traits<ReferenceEdge, 2> : public hfv1_traits_ReferenceEdge
1383 : {
1384 0 : static void NormalOnSCVF(MathVector<2>& outNormal, const MathVector<2>* vCornerCoords)
1385 0 : {UG_THROW("Not implemented");}
1386 : };
1387 :
1388 : template <> struct hfv1_traits<ReferenceEdge, 3> : public hfv1_traits_ReferenceEdge
1389 : {
1390 0 : static void NormalOnSCVF(MathVector<3>& outNormal, const MathVector<3>* vCornerCoords)
1391 0 : {UG_THROW("Not implemented");}
1392 : };
1393 :
1394 : /////////////////////////
1395 : // 2D Reference Element
1396 : /////////////////////////
1397 :
1398 : struct hfv1_traits_ReferenceFace
1399 : {
1400 : const static size_t NumCornersOfSCVF = 2;
1401 : const static size_t MaxNumCornersOfSCV = 4;
1402 : typedef ReferenceQuadrilateral scv_type;
1403 : };
1404 :
1405 : template <> struct hfv1_traits<ReferenceTriangle, 2> : public hfv1_traits_ReferenceFace
1406 : {
1407 : static void NormalOnSCVF(MathVector<2>& outNormal, const MathVector<2>* vCornerCoords)
1408 : {ElementNormal<ReferenceEdge, 2>(outNormal, vCornerCoords);}
1409 : };
1410 :
1411 : template <> struct hfv1_traits<ReferenceTriangle, 3> : public hfv1_traits_ReferenceFace
1412 : {
1413 0 : static void NormalOnSCVF(MathVector<3>& outNormal, const MathVector<3>* vCornerCoords)
1414 0 : {UG_THROW("Not implemented");}
1415 : };
1416 :
1417 : template <> struct hfv1_traits<ReferenceQuadrilateral, 2> : public hfv1_traits_ReferenceFace
1418 : {
1419 : static void NormalOnSCVF(MathVector<2>& outNormal, const MathVector<2>* vCornerCoords)
1420 : {ElementNormal<ReferenceEdge, 2>(outNormal, vCornerCoords);}
1421 : };
1422 :
1423 : template <> struct hfv1_traits<ReferenceQuadrilateral, 3> : public hfv1_traits_ReferenceFace
1424 : {
1425 0 : static void NormalOnSCVF(MathVector<3>& outNormal, const MathVector<3>* vCornerCoords)
1426 0 : {UG_THROW("Not implemented");}
1427 : };
1428 :
1429 : /////////////////////////
1430 : // 3D Reference Element
1431 : /////////////////////////
1432 :
1433 : struct hfv1_traits_ReferenceVolume
1434 : {
1435 : static void NormalOnSCVF(MathVector<3>& outNormal, const MathVector<3>* vCornerCoords)
1436 0 : {ElementNormal<ReferenceTriangle, 3>(outNormal, vCornerCoords);}
1437 :
1438 : typedef ReferenceTetrahedron scv_type;
1439 : };
1440 :
1441 : template <> struct hfv1_traits<ReferenceTetrahedron, 3> : public hfv1_traits_ReferenceVolume
1442 : {
1443 : const static size_t NumCornersOfSCVF = 3;
1444 : const static size_t MaxNumCornersOfSCV = 8;
1445 : };
1446 :
1447 : template <> struct hfv1_traits<ReferencePrism, 3> : public hfv1_traits_ReferenceVolume
1448 : {
1449 : const static size_t NumCornersOfSCVF = 3;
1450 : const static size_t MaxNumCornersOfSCV = 8;
1451 : };
1452 :
1453 : template <> struct hfv1_traits<ReferencePyramid, 3> : public hfv1_traits_ReferenceVolume
1454 : {
1455 : const static size_t NumCornersOfSCVF = 3;
1456 : const static size_t MaxNumCornersOfSCV = 10;
1457 : };
1458 :
1459 : template <> struct hfv1_traits<ReferenceHexahedron, 3> : public hfv1_traits_ReferenceVolume
1460 : {
1461 : const static size_t NumCornersOfSCVF = 3;
1462 : const static size_t MaxNumCornersOfSCV = 8;
1463 : };
1464 :
1465 : template <> struct hfv1_traits<ReferenceOctahedron, 3> : public hfv1_traits_ReferenceVolume
1466 : {
1467 : const static size_t NumCornersOfSCVF = 3;
1468 : const static size_t MaxNumCornersOfSCV = 8;
1469 : };
1470 :
1471 : template <int TDim> struct hdimfv1_traits
1472 : {
1473 : typedef void scv_type;
1474 : typedef void elem_type_0;
1475 : typedef void elem_type_1;
1476 : typedef void elem_type_2;
1477 : typedef void elem_type_3;
1478 : typedef void elem_type_4;
1479 : const static size_t NumCornersOfSCVF;
1480 : const static size_t MaxNumCornersOfSCV;
1481 : };
1482 :
1483 : template <> struct hdimfv1_traits<1>
1484 : {
1485 : typedef ReferenceEdge scv_type;
1486 : typedef ReferenceEdge elem_type_0;
1487 : typedef ReferenceEdge elem_type_1;
1488 : typedef ReferenceEdge elem_type_2;
1489 : typedef ReferenceEdge elem_type_3;
1490 : typedef ReferenceEdge elem_type_4;
1491 : const static size_t NumCornersOfSCVF = 1;
1492 : const static size_t MaxNumCornersOfSCV = 2;
1493 : };
1494 :
1495 : template <> struct hdimfv1_traits<2>
1496 : {
1497 : typedef ReferenceQuadrilateral scv_type;
1498 : typedef ReferenceTriangle elem_type_0;
1499 : typedef ReferenceQuadrilateral elem_type_1;
1500 : typedef ReferenceTriangle elem_type_2;
1501 : typedef ReferenceQuadrilateral elem_type_3;
1502 : typedef ReferenceQuadrilateral elem_type_4;
1503 : const static size_t NumCornersOfSCVF = 2;
1504 : const static size_t MaxNumCornersOfSCV = 4;
1505 : };
1506 :
1507 : template <> struct hdimfv1_traits<3>
1508 : {
1509 : typedef ReferenceTetrahedron scv_type;
1510 : typedef ReferenceTetrahedron elem_type_0;
1511 : typedef ReferencePyramid elem_type_1;
1512 : typedef ReferencePrism elem_type_2;
1513 : typedef ReferenceHexahedron elem_type_3;
1514 : typedef ReferenceOctahedron elem_type_4;
1515 : const static size_t NumCornersOfSCVF = 3;
1516 : const static size_t MaxNumCornersOfSCV = 10;
1517 : };
1518 :
1519 :
1520 : ////////////////////////////////////////////////////////////////////////////////
1521 : // Functions
1522 : ////////////////////////////////////////////////////////////////////////////////
1523 :
1524 : template <typename TRefElem, int TWorldDim>
1525 : void HangingNormalOnSCVF(MathVector<TWorldDim>& outNormal, const MathVector<TWorldDim>* vCornerCoords)
1526 : {
1527 0 : hfv1_traits<TRefElem, TWorldDim>::NormalOnSCVF(outNormal, vCornerCoords);
1528 0 : }
1529 :
1530 : ////////////////////////////////////////////////////////////////////////////////
1531 : // FVHO: Finite Volume for Higher Order traits
1532 : ////////////////////////////////////////////////////////////////////////////////
1533 :
1534 : /// Traits for Finite Volumes of higher order
1535 : template <int TOrder, typename TRefElem, int TWorldDim> struct fvho_traits
1536 : {
1537 : // can be inherited from fv1_traits (since the same)
1538 : const static size_t NumCornersOfSCVF;
1539 : const static size_t MaxNumCornersOfSCV;
1540 : static void NormalOnSCVF(MathVector<TWorldDim>& outNormal, const MathVector<TWorldDim>* vCornerCoords);
1541 : typedef void scv_type;
1542 :
1543 : // own data
1544 : const static size_t NumSubElem;
1545 : };
1546 :
1547 : } // end namespace ug
1548 :
1549 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_UTIL__ */
|
