3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dinteractionvolumecontaineradaptive.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
25#define DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
26
27// dumux environment
32
33namespace Dumux {
34
45template<class TypeTag>
47{
48 friend class FvMpfaL3dInteractionVolumeContainer<TypeTag>;
50 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
51
52 enum
53 {
54 dim = GridView::dimension, dimWorld = GridView::dimensionworld
55 };
56
57 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
58 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
59
60 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
61
62 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
63
64 using Element = typename GridView::Traits::template Codim<0>::Entity;
65 using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
66 using ElementGeometry = typename Element::Geometry;
67
68 using GlobalPosition = typename ElementGeometry::GlobalCoordinate;
69 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
70
71 using DimVector = Dune::FieldVector<Scalar, dim>;
72 using FaceVerticesVector = std::vector<Dune::FieldVector<std::set<int>, 2*dim> >;
73
74 enum
75 {
76 pressureEqIdx = Indices::pressureEqIdx,
77 };
78
79 using IndexTranslator = IndexTranslatorAdaptive;
80
81public:
85
86private:
87
88 void storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex);
89 void storeInnerInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex);
90
91protected:
93
94public:
95
104 std::set<int>& faceVerticeIndices(int eIdxGlobal, int fIdx)
105 {
106 return faceVertices_[eIdxGlobal][fIdx];
107 }
108
114 ParentType(problem), problem_(problem)
115 {
116 if (dim != 3)
117 {
118 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
119 }
120 }
121
122private:
123 using Implementation = typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolumeContainer);
124
126 Implementation &asImp_()
127 { return *static_cast<Implementation *>(this);}
128
130 const Implementation &asImp_() const
131 { return *static_cast<const Implementation *>(this);}
132
133 Problem& problem_;
134 FaceVerticesVector faceVertices_;
135};
136
160template<class TypeTag>
161void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume,
162 const Vertex& vertex)
163{
164 // if cells of different grid level appear, the geometric information is constructed going from the
165 // coarsest level to the finest level. This ensures that coarser information ins always overwritten by
166 // finer information.
167 if (!interactionVolume.sameLevel())
168 {
169 // sort the local element indices according to the grid level
170 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
171 for (int i = 0; i < 8; i++)
172 {
173 levelIdx[i][0] = i;
174 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
175 }
176
177 std::sort(levelIdx.begin(), levelIdx.end(), [](const auto& a, const auto& b) { return (a[1]<b[1]); });
178
179 // Generate and store the geometric information going from the coarsest to the finest level.
180 // For the calculation we take advantage from the fact that the ordering inside the interaction volume
181 // with respect to the DUNE reference element is known due to the storage process of the elements in
182 // FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements
183 // (const Element& element, std::vector < std::vector<int> >& elemVertMap)
184 for (int i = 0; i < 8; i++)
185 {
186 int idx = levelIdx[i][0];
187 auto element = interactionVolume.getSubVolumeElement(idx);
188
189 const ElementGeometry& geometry = element.geometry();
190
191 const auto referenceElement = ReferenceElements::general(geometry.type());
192
193 switch (idx)
194 {
195 case 0:
196 {
197 DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
198 interactionVolume.setEdgePosition(edgeCoord, 2);
199 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
200 interactionVolume.setEdgePosition(edgeCoord, 0);
201 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
202 interactionVolume.setEdgePosition(edgeCoord, 5);
203
204 break;
205 }
206 case 1:
207 {
208 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
209 interactionVolume.setEdgePosition(edgeCoord, 0);
210 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
211 interactionVolume.setEdgePosition(edgeCoord, 2);
212 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
213 interactionVolume.setEdgePosition(edgeCoord, 3);
214
215 break;
216 }
217 case 2:
218 {
219 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
220 interactionVolume.setEdgePosition(edgeCoord, 0);
221 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
222 interactionVolume.setEdgePosition(edgeCoord, 4);
223 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
224 interactionVolume.setEdgePosition(edgeCoord, 5);
225
226 break;
227 }
228 case 3:
229 {
230 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
231 interactionVolume.setEdgePosition(edgeCoord, 0);
232 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
233 interactionVolume.setEdgePosition(edgeCoord, 4);
234 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
235 interactionVolume.setEdgePosition(edgeCoord, 3);
236
237 break;
238 }
239 case 4:
240 {
241 DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
242 interactionVolume.setEdgePosition(edgeCoord, 1);
243 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
244 interactionVolume.setEdgePosition(edgeCoord, 2);
245 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
246 interactionVolume.setEdgePosition(edgeCoord, 5);
247
248 break;
249 }
250 case 5:
251 {
252 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
253 interactionVolume.setEdgePosition(edgeCoord, 1);
254 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
255 interactionVolume.setEdgePosition(edgeCoord, 2);
256 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
257 interactionVolume.setEdgePosition(edgeCoord, 3);
258
259 break;
260 }
261 case 6:
262 {
263 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
264 interactionVolume.setEdgePosition(edgeCoord, 1);
265 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
266 interactionVolume.setEdgePosition(edgeCoord, 4);
267 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
268 interactionVolume.setEdgePosition(edgeCoord, 5);
269
270 break;
271 }
272 case 7:
273 {
274 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
275 interactionVolume.setEdgePosition(edgeCoord, 1);
276 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
277 interactionVolume.setEdgePosition(edgeCoord, 4);
278 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
279 interactionVolume.setEdgePosition(edgeCoord, 3);
280
281 break;
282 }
283 }
284 }
285 ParentType::storeInnerInteractionVolume(interactionVolume, vertex, false);
286 }
287 else
288 {
289 ParentType::storeInnerInteractionVolume(interactionVolume, vertex, true);
290 }
291}
292
315template<class TypeTag>
316void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume,
317 const Vertex& vertex)
318{
319 const DimVector& centerPos = vertex.geometry().center();
320
321 interactionVolume.setCenterPosition(centerPos);
322
323 // sort the local element indices according to the grid level
324 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
325 for (int i = 0; i < 8; i++)
326 {
327 levelIdx[i][0] = i;
328 if (interactionVolume.hasSubVolumeElement(i))
329 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
330 else
331 levelIdx[i][1] = -1;
332 }
333
334 std::sort(levelIdx.begin(), levelIdx.end(), [](const auto& a, const auto& b) { return (a[1]<b[1]); });
335
336 // Generate and store the geometric information going from the coarsest to the finest level.
337 // For the calculation we take advantage from the fact that the ordering inside the interaction volume
338 // with respect to the DUNE reference element is known due to the storage process of the elements in
339 // FvMpfaL3dInteractionVolumeContainer<TypeTag>::
340 // storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap)
341 for (int i = 0; i < 8; i++)
342 {
343 if (levelIdx[i][1] < 0)
344 continue;
345
346 int idx = levelIdx[i][0];
347
348 auto element = interactionVolume.getSubVolumeElement(idx);
349
350 const ElementGeometry& geometry = element.geometry();
351
352 const auto referenceElement = ReferenceElements::general(geometry.type());
353
354 switch (idx)
355 {
356 case 0:
357 {
358 DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
359 interactionVolume.setEdgePosition(edgeCoord, 2);
360 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
361 interactionVolume.setEdgePosition(edgeCoord, 0);
362 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
363 interactionVolume.setEdgePosition(edgeCoord, 5);
364
365 break;
366 }
367 case 1:
368 {
369 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
370 interactionVolume.setEdgePosition(edgeCoord, 0);
371 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
372 interactionVolume.setEdgePosition(edgeCoord, 2);
373 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
374 interactionVolume.setEdgePosition(edgeCoord, 3);
375
376 break;
377 }
378 case 2:
379 {
380 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
381 interactionVolume.setEdgePosition(edgeCoord, 0);
382 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
383 interactionVolume.setEdgePosition(edgeCoord, 4);
384 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
385 interactionVolume.setEdgePosition(edgeCoord, 5);
386
387 break;
388 }
389 case 3:
390 {
391 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
392 interactionVolume.setEdgePosition(edgeCoord, 0);
393 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
394 interactionVolume.setEdgePosition(edgeCoord, 4);
395 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
396 interactionVolume.setEdgePosition(edgeCoord, 3);
397
398 break;
399 }
400 case 4:
401 {
402 DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
403 interactionVolume.setEdgePosition(edgeCoord, 1);
404 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
405 interactionVolume.setEdgePosition(edgeCoord, 2);
406 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
407 interactionVolume.setEdgePosition(edgeCoord, 5);
408
409 break;
410 }
411 case 5:
412 {
413 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
414 interactionVolume.setEdgePosition(edgeCoord, 1);
415 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
416 interactionVolume.setEdgePosition(edgeCoord, 2);
417 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
418 interactionVolume.setEdgePosition(edgeCoord, 3);
419
420 break;
421 }
422 case 6:
423 {
424 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
425 interactionVolume.setEdgePosition(edgeCoord, 1);
426 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
427 interactionVolume.setEdgePosition(edgeCoord, 4);
428 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
429 interactionVolume.setEdgePosition(edgeCoord, 5);
430
431 break;
432 }
433 case 7:
434 {
435 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
436 interactionVolume.setEdgePosition(edgeCoord, 1);
437 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
438 interactionVolume.setEdgePosition(edgeCoord, 4);
439 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
440 interactionVolume.setEdgePosition(edgeCoord, 3);
441
442 break;
443 }
444 }
445 }
446
447 // Choose the type of the hanging-node-interaction-volume depending on the number of stored fine
448 // elements (see dissertation M. Wolff, http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/)
449 // and add missing elements and geometric information
450 switch (interactionVolume.getElementNumber())
451 {
452 // hanging-node interaction volume of type 5 or 7
453 case 2:
454 {
455 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
456
457 std::vector<int> elemIdxOld;
458 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
459 {
460 if (interactionVolume.hasSubVolumeElement(i))
461 elemIdxOld.push_back(i);
462 }
463
464 int zeroFaceIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[0], elemIdxOld[1]);
465
466 std::vector<int> elemIdxNew(2);
467 elemIdxNew[0] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[0]);
468 elemIdxNew[1] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[1]);
469
470
471 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[0]),
472 elemIdxNew[0]);
473 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[1]),
474 elemIdxNew[1]);
475
476 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
477 {
478 for (int i = 0; i < 3; i++)
479 {
480 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
481 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
482
483 for (int j = 0; j < 3; j++)
484 {
485 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
486 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
487
488 if (faceIdxNew == faceIdxNewTest)
489 {
490 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
491 elemNew, j);
492 hangingNodeVolume.setIndexOnElement(
493 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
494 }
495 }
496 }
497 }
498
499 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
500 {
501 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
502 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
503 }
504
505 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
506 {
507 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
508 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
509 }
510
511 interactionVolume = hangingNodeVolume;
512
513 interactionVolume.setHangingNodeType(InteractionVolume::twoSmallCells);
514
515 interactionVolume.setCenterPosition(centerPos);
516
517 auto element1 = interactionVolume.getSubVolumeElement(0);
518
519 for (const auto& intersection : intersections(problem_.gridView(), element1))
520 {
521 int idxInInside = intersection.indexInInside();
522
523 if (idxInInside == interactionVolume.getIndexOnElement(0, 2))
524 {
525 if (intersection.neighbor())
526 {
527 auto outside = intersection.outside();
528 if (element1.level() > outside.level())
529 {
530 interactionVolume.setSubVolumeElement(outside, 4);
531 interactionVolume.setSubVolumeElement(outside, 5);
532 }
533 }
534 }
535 else if (idxInInside == interactionVolume.getIndexOnElement(0, 1))
536 {
537 if (intersection.neighbor())
538 {
539 auto outside = intersection.outside();
540 if (element1.level() > outside.level())
541 {
542 interactionVolume.setSubVolumeElement(outside, 2);
543 interactionVolume.setSubVolumeElement(outside, 3);
544 }
545 }
546 }
547 }
548 auto element2 = interactionVolume.getSubVolumeElement(1);
549 auto element45 = interactionVolume.getSubVolumeElement(4);
550 auto element23 = interactionVolume.getSubVolumeElement(2);
551
552 for (const auto& intersection1 : intersections(problem_.gridView(), element45))
553 {
554 if (intersection1.neighbor())
555 {
556 auto element45Outside = intersection1.outside();
557
558 for (const auto& intersection2 : intersections(problem_.gridView(), element23))
559 {
560 if (intersection2.neighbor())
561 {
562 auto element23Outside = intersection2.outside();
563
564 if (element45Outside == element23Outside && element45Outside != element1
565 && element45Outside != element2)
566 {
567 interactionVolume.setSubVolumeElement(element45Outside, 6);
568 interactionVolume.setSubVolumeElement(element45Outside, 7);
569 DimVector normal = intersection2.centerUnitOuterNormal();
570 interactionVolume.setNormal(normal, 2, 2);
571 interactionVolume.setNormal(normal, 3, 2);
572 normal *= -1;
573 interactionVolume.setNormal(normal, 6, 0);
574 interactionVolume.setNormal(normal, 7, 0);
575
576 GlobalPosition globalPosFace = intersection2.geometry().center();
577 interactionVolume.setFacePosition(globalPosFace, 10);
578 interactionVolume.setFacePosition(globalPosFace, 11);
579 interactionVolume.setEdgePosition(centerPos, 4);
580
581 Scalar faceArea = intersection2.geometry().volume()/4.0;
582
583 interactionVolume.setFaceArea(faceArea, 10);
584 interactionVolume.setFaceArea(faceArea, 11);
585
586 normal = intersection1.centerUnitOuterNormal();
587 interactionVolume.setNormal(normal, 4, 2);
588 interactionVolume.setNormal(normal, 5, 1);
589 normal *= -1;
590 interactionVolume.setNormal(normal, 6, 1);
591 interactionVolume.setNormal(normal, 7, 2);
592
593 globalPosFace = intersection1.geometry().center();
594 interactionVolume.setFacePosition(globalPosFace, 5);
595 interactionVolume.setFacePosition(globalPosFace, 7);
596 interactionVolume.setEdgePosition(centerPos, 1);
597
598 faceArea = intersection1.geometry().volume()/4.0;
599
600 interactionVolume.setFaceArea(faceArea, 5);
601 interactionVolume.setFaceArea(faceArea, 7);
602 }
603 }
604 }
605 }
606 }
607
608 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
609 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
610 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
611 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
612 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
613 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
614
615 DimVector crossProductVector1(0);
616 DimVector crossProductVector2(0);
617
618 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
619 crossProductVector1 = centerPos - globalPosFace;
620 crossProductVector2 = edgeCoord1 - edgeCoord3;
621 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
622 interactionVolume.setFaceArea(faceArea, 0);
623
624 globalPosFace = interactionVolume.getFacePosition(1);
625 crossProductVector1 = centerPos - globalPosFace;
626 crossProductVector2 = edgeCoord1 - edgeCoord4;
627 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
628 interactionVolume.setFaceArea(faceArea, 1);
629
630 globalPosFace = interactionVolume.getFacePosition(3);
631 crossProductVector1 = centerPos - globalPosFace;
632 crossProductVector2 = edgeCoord1 - edgeCoord6;
633 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
634 interactionVolume.setFaceArea(faceArea, 3);
635
636 globalPosFace = interactionVolume.getFacePosition(8);
637 crossProductVector1 = centerPos - globalPosFace;
638 crossProductVector2 = edgeCoord3 - edgeCoord6;
639 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
640 interactionVolume.setFaceArea(faceArea, 8);
641
642 globalPosFace = interactionVolume.getFacePosition(9);
643 crossProductVector1 = centerPos - globalPosFace;
644 crossProductVector2 = edgeCoord3 - edgeCoord4;
645 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
646 interactionVolume.setFaceArea(faceArea, 9);
647
648 break;
649 }
650 //hanging-node interaction volume of type 1, 3 or 4
651 case 4:
652 {
653 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
654
655 std::vector<int> elemIdxOld;
656 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
657 {
658 if (interactionVolume.hasSubVolumeElement(i))
659 elemIdxOld.push_back(i);
660 }
661
662 std::set<int> zeroFaceIdxVec;
663 for (int i = 0; i < 4; i++)
664 {
665 for (int j = 0; j < 4; j++)
666 {
667 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
668 if (fIdx >= 0)
669 zeroFaceIdxVec.insert(fIdx);
670 }
671 }
672
673 std::vector<int> elemIdxNew(4);
674 int zeroFaceIdx = 0;
675
676 if (zeroFaceIdxVec.size() == 2)
677 {
678 hangingNodeVolume.setHangingNodeType(InteractionVolume::fourSmallCellsDiag);
679
680 if (zeroFaceIdxVec.find(0) == zeroFaceIdxVec.end())
681 zeroFaceIdx = *zeroFaceIdxVec.begin();
682
683 for (int i = 0; i < 4; i++)
684 {
685 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
686 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
687 elemIdxNew[i]);
688 }
689 }
690 else if (zeroFaceIdxVec.size() == 4)
691 {
692 std::set<int>::iterator it = zeroFaceIdxVec.begin();
693 for (; it != zeroFaceIdxVec.end(); ++it)
694 {
695 zeroFaceIdx = *it;
696
697 bool isFace = true;
698 for (int i = 0; i < 4; i++)
699 {
700 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx,
701 elemIdxOld[i]);
702 if (elemIdxNew[i] == 4 || elemIdxNew[i] == 5 || elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
703 {
704 isFace = false;
705 break;
706 }
707 }
708 if (isFace)
709 {
710 for (int i = 0; i < 4; i++)
711 {
712 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]), elemIdxNew[i]);
713 }
714
715 break;
716 }
717 }
718 }
719 else
720 {
721 std::set<int>::iterator it = zeroFaceIdxVec.begin();
722 std::cout<<"zeroFaceIdxVec = ";
723 for(;it != zeroFaceIdxVec.end(); ++it)
724 {
725 std::cout<<" "<<*it;
726 }
727 std::cout<<"\n";
728 }
729
730 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
731 {
732 for (int i = 0; i < 3; i++)
733 {
734 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
735 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
736
737 for (int j = 0; j < 3; j++)
738 {
739 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
740 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
741
742 if (faceIdxNew == faceIdxNewTest)
743 {
744 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
745 elemNew, j);
746 hangingNodeVolume.setIndexOnElement(interactionVolume.getIndexOnElement(elem, i), elemNew, j);
747 }
748 }
749 }
750 }
751
752 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
753 {
754 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
755 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
756 }
757
758 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
759 {
760 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
761 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
762 }
763
764 interactionVolume = hangingNodeVolume;
765
766 interactionVolume.setCenterPosition(centerPos);
767
768 if (zeroFaceIdxVec.size() == 4)
769 {
770 auto element1 = interactionVolume.getSubVolumeElement(0);
771 auto element4 = interactionVolume.getSubVolumeElement(3);
772
773 auto outside1 = element1;
774 auto outside4 = element4;
775
776 for (const auto& intersection1 : intersections(problem_.gridView(), element1))
777 {
778 if (intersection1.neighbor())
779 {
780 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
781 {
782 outside1 = intersection1.outside();
783 break;
784 }
785 }
786 }
787 for (const auto& intersection4 : intersections(problem_.gridView(), element4))
788 {
789 if (intersection4.neighbor())
790 {
791 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
792 {
793 outside4 = intersection4.outside();
794 break;
795 }
796 }
797 }
798 if (outside1 != outside4)
799 {
800 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsEdge);
801 }
802 else if (outside1 == outside4)
803 {
804 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsFace);
805 }
806 }
807
808 if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsFace)
809 {
810 auto element = interactionVolume.getSubVolumeElement(0);
811
812 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
813 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
814 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
815 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
816 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
817 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
818
819 DimVector crossProductVector1(0);
820 DimVector crossProductVector2(0);
821
822 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
823 crossProductVector1 = centerPos - globalPosFace;
824 crossProductVector2 = edgeCoord1 - edgeCoord3;
825 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
826 interactionVolume.setFaceArea(faceArea, 0);
827
828 globalPosFace = interactionVolume.getFacePosition(1);
829 crossProductVector1 = centerPos - globalPosFace;
830 crossProductVector2 = edgeCoord1 - edgeCoord4;
831 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
832 interactionVolume.setFaceArea(faceArea, 1);
833
834 globalPosFace = interactionVolume.getFacePosition(2);
835 crossProductVector1 = centerPos - globalPosFace;
836 crossProductVector2 = edgeCoord1 - edgeCoord5;
837 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
838 interactionVolume.setFaceArea(faceArea, 2);
839
840 globalPosFace = interactionVolume.getFacePosition(3);
841 crossProductVector1 = centerPos - globalPosFace;
842 crossProductVector2 = edgeCoord1 - edgeCoord6;
843 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
844 interactionVolume.setFaceArea(faceArea, 3);
845
846 globalPosFace = interactionVolume.getFacePosition(8);
847 crossProductVector1 = centerPos - globalPosFace;
848 crossProductVector2 = edgeCoord3 - edgeCoord6;
849 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
850 interactionVolume.setFaceArea(faceArea, 8);
851
852 globalPosFace = interactionVolume.getFacePosition(9);
853 crossProductVector1 = centerPos - globalPosFace;
854 crossProductVector2 = edgeCoord3 - edgeCoord4;
855 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
856 interactionVolume.setFaceArea(faceArea, 9);
857
858 globalPosFace = interactionVolume.getFacePosition(10);
859 crossProductVector1 = centerPos - globalPosFace;
860 crossProductVector2 = edgeCoord4 - edgeCoord5;
861 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
862 interactionVolume.setFaceArea(faceArea, 10);
863
864 globalPosFace = interactionVolume.getFacePosition(11);
865 crossProductVector1 = centerPos - globalPosFace;
866 crossProductVector2 = edgeCoord5 - edgeCoord6;
867 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
868 interactionVolume.setFaceArea(faceArea, 11);
869
870 for (const auto& intersection : intersections(problem_.gridView(), element))
871 {
872 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
873 {
874 if (intersection.neighbor())
875 {
876 auto outside = intersection.outside();
877 if (element.level() > outside.level())
878 {
879 interactionVolume.setSubVolumeElement(outside, 4);
880 interactionVolume.setSubVolumeElement(outside, 5);
881 interactionVolume.setSubVolumeElement(outside, 6);
882 interactionVolume.setSubVolumeElement(outside, 7);
883
884 break;
885 }
886 }
887 }
888
889 }
890 }
891 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsEdge)
892 {
893 auto element1 = interactionVolume.getSubVolumeElement(0);
894 auto element2 = interactionVolume.getSubVolumeElement(1);
895 auto element3 = interactionVolume.getSubVolumeElement(2);
896 auto element4 = interactionVolume.getSubVolumeElement(3);
897
898 for (const auto& intersection1 : intersections(problem_.gridView(), element1))
899 {
900 if (intersection1.neighbor())
901 {
902 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
903 {
904 auto outside = intersection1.outside();
905 if (element1.level() > outside.level())
906 {
907 interactionVolume.setSubVolumeElement(outside, 4);
908 }
909 }
910 }
911 }
912 for (const auto& intersection2 : intersections(problem_.gridView(), element2))
913 {
914 if (intersection2.neighbor())
915 {
916 if (intersection2.indexInInside() == interactionVolume.getIndexOnElement(1, 2))
917 {
918 auto outside = intersection2.outside();
919 if (element2.level() > outside.level())
920 {
921 interactionVolume.setSubVolumeElement(outside, 5);
922
923 break;
924 }
925 }
926 }
927 }
928 for (const auto& intersection3 : intersections(problem_.gridView(), element3))
929 {
930 if (intersection3.neighbor())
931 {
932 if (intersection3.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
933 {
934 auto outside = intersection3.outside();
935 if (element3.level() > outside.level())
936 {
937 interactionVolume.setSubVolumeElement(outside, 6);
938 break;
939 }
940 }
941 }
942 }
943 for (const auto& intersection4 : intersections(problem_.gridView(), element4))
944 {
945 if (intersection4.neighbor())
946 {
947 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
948 {
949 auto outside = intersection4.outside();
950 if (element4.level() > outside.level())
951 {
952 interactionVolume.setSubVolumeElement(outside, 7);
953
954 break;
955 }
956 }
957 }
958 }
959
960 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
961 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
962 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
963 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
964 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
965
966 DimVector crossProductVector1(0);
967 DimVector crossProductVector2(0);
968
969 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
970 crossProductVector1 = centerPos - globalPosFace;
971 crossProductVector2 = edgeCoord1 - edgeCoord3;
972 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
973 interactionVolume.setFaceArea(faceArea, 0);
974
975 globalPosFace = interactionVolume.getFacePosition(1);
976 crossProductVector1 = centerPos - globalPosFace;
977 crossProductVector2 = edgeCoord1 - edgeCoord4;
978 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
979 interactionVolume.setFaceArea(faceArea, 1);
980
981 globalPosFace = interactionVolume.getFacePosition(2);
982 crossProductVector1 = centerPos - globalPosFace;
983 crossProductVector2 = edgeCoord1 - edgeCoord5;
984 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
985 interactionVolume.setFaceArea(faceArea, 2);
986
987 globalPosFace = interactionVolume.getFacePosition(3);
988 crossProductVector1 = centerPos - globalPosFace;
989 crossProductVector2 = edgeCoord1 - edgeCoord6;
990 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
991 interactionVolume.setFaceArea(faceArea, 3);
992
993 globalPosFace = interactionVolume.getFacePosition(8);
994 crossProductVector1 = centerPos - globalPosFace;
995 crossProductVector2 = edgeCoord3 - edgeCoord6;
996 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
997 interactionVolume.setFaceArea(faceArea, 8);
998
999 globalPosFace = interactionVolume.getFacePosition(9);
1000 crossProductVector1 = centerPos - globalPosFace;
1001 crossProductVector2 = edgeCoord3 - edgeCoord4;
1002 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1003 interactionVolume.setFaceArea(faceArea, 9);
1004
1005 globalPosFace = interactionVolume.getFacePosition(10);
1006 crossProductVector1 = centerPos - globalPosFace;
1007 crossProductVector2 = edgeCoord4 - edgeCoord5;
1008 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1009 interactionVolume.setFaceArea(faceArea, 10);
1010
1011 globalPosFace = interactionVolume.getFacePosition(11);
1012 crossProductVector1 = centerPos - globalPosFace;
1013 crossProductVector2 = edgeCoord5 - edgeCoord6;
1014 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1015 interactionVolume.setFaceArea(faceArea, 11);
1016
1017 auto element5 = interactionVolume.getSubVolumeElement(4);
1018 auto element6 = interactionVolume.getSubVolumeElement(5);
1019 auto element7 = interactionVolume.getSubVolumeElement(6);
1020 auto element8 = interactionVolume.getSubVolumeElement(7);
1021
1022 if (element5 == element6)
1023 {
1024 interactionVolume.setFacePosition(element5.geometry().center(), 4);
1025 interactionVolume.setFacePosition(element7.geometry().center(), 6);
1026
1027 for (const auto& intersection : intersections(problem_.gridView(), element5))
1028 {
1029 if (intersection.neighbor())
1030 {
1031 auto outside = intersection.outside();
1032
1033 if (outside == element7 || outside == element8)
1034 {
1035 int indexInInside = intersection.indexInInside();
1036 interactionVolume.setIndexOnElement(indexInInside, 4, 2);
1037 interactionVolume.setIndexOnElement(indexInInside, 5, 1);
1038 DimVector normal = intersection.centerUnitOuterNormal();
1039 interactionVolume.setNormal(normal, 4, 2);
1040 interactionVolume.setNormal(normal, 5, 1);
1041 globalPosFace = intersection.geometry().center();
1042 interactionVolume.setFacePosition(globalPosFace, 5);
1043 interactionVolume.setFacePosition(globalPosFace, 7);
1044 int indexInOutside = intersection.indexInOutside();
1045 interactionVolume.setIndexOnElement(indexInOutside, 6, 1);
1046 interactionVolume.setIndexOnElement(indexInOutside, 7, 2);
1047 normal *= -1;
1048 interactionVolume.setNormal(normal, 6, 1);
1049 interactionVolume.setNormal(normal, 7, 2);
1050
1051 break;
1052 }
1053 }
1054 }
1055
1056 DimVector edgeCoord2(interactionVolume.getFacePosition(7));
1057
1058 globalPosFace = interactionVolume.getFacePosition(5);
1059 crossProductVector1 = centerPos - globalPosFace;
1060 crossProductVector2 = edgeCoord2 - edgeCoord4;
1061 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1062 interactionVolume.setFaceArea(faceArea, 5);
1063 // interactionVolume.setFaceArea(0.0, 5);
1064 globalPosFace = interactionVolume.getFacePosition(7);
1065 crossProductVector1 = centerPos - globalPosFace;
1066 crossProductVector2 = edgeCoord2 - edgeCoord6;
1067 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1068 interactionVolume.setFaceArea(faceArea, 7);
1069 // interactionVolume.setFaceArea(0.0, 7);
1070 }
1071 else if (element5 == element7)
1072 {
1073 interactionVolume.setFacePosition(element6.geometry().center(), 5);
1074 interactionVolume.setFacePosition(element5.geometry().center(), 7);
1075 interactionVolume.setFacePosition(globalPosFace, 6);
1076
1077 for (const auto& intersection : intersections(problem_.gridView(), element5))
1078 {
1079 if (intersection.neighbor())
1080 {
1081 auto outside = intersection.outside();
1082
1083 if (outside == element6 || outside == element8)
1084 {
1085 int indexInInside = intersection.indexInInside();
1086 interactionVolume.setIndexOnElement(indexInInside, 4, 1);
1087 interactionVolume.setIndexOnElement(indexInInside, 6, 2);
1088 DimVector normal = intersection.centerUnitOuterNormal();
1089 interactionVolume.setNormal(normal, 4, 1);
1090 interactionVolume.setNormal(normal, 6, 2);
1091 globalPosFace = intersection.geometry().center();
1092 interactionVolume.setFacePosition(globalPosFace, 4);
1093 interactionVolume.setFacePosition(globalPosFace, 6);
1094 int indexInOutside = intersection.indexInOutside();
1095 interactionVolume.setIndexOnElement(indexInOutside, 5, 2);
1096 interactionVolume.setIndexOnElement(indexInOutside, 7, 1);
1097 normal *= -1;
1098 interactionVolume.setNormal(normal, 5, 2);
1099 interactionVolume.setNormal(normal, 7, 1);
1100
1101 break;
1102 }
1103 }
1104 }
1105 DimVector edgeCoord2(interactionVolume.getFacePosition(4));
1106
1107 globalPosFace = interactionVolume.getFacePosition(4);
1108 crossProductVector1 = centerPos - globalPosFace;
1109 crossProductVector2 = edgeCoord2 - edgeCoord3;
1110 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1111 interactionVolume.setFaceArea(faceArea, 4);
1112
1113 globalPosFace = interactionVolume.getFacePosition(6);
1114 crossProductVector1 = centerPos - globalPosFace;
1115 crossProductVector2 = edgeCoord2 - edgeCoord5;
1116 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1117 interactionVolume.setFaceArea(faceArea, 6);
1118 // interactionVolume.setFaceArea(0.0, 4);
1119 // interactionVolume.setFaceArea(0.0, 6);
1120 }
1121
1122 const ElementGeometry& geometry = element5.geometry();
1123
1124 const auto referenceElement = ReferenceElements::general(geometry.type());
1125
1126 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
1127 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
1128
1129 DimVector edgeCoord(0.0);
1130 switch (oldSubVolumElemIdx)
1131 {
1132 case 0:
1133 {
1134 switch (oldEdgeIdx)
1135 {
1136 case 2:
1137 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1138 break;
1139 case 0:
1140 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1141 break;
1142 case 5:
1143 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1144 break;
1145 }
1146
1147 break;
1148 }
1149 case 1:
1150 {
1151 switch (oldEdgeIdx)
1152 {
1153 case 0:
1154 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1155 break;
1156 case 2:
1157 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1158 break;
1159 case 3:
1160 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1161 break;
1162 }
1163
1164 break;
1165 }
1166 case 2:
1167 {
1168 switch (oldEdgeIdx)
1169 {
1170 case 0:
1171 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1172 break;
1173 case 4:
1174 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1175 break;
1176 case 5:
1177 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1178 break;
1179 }
1180
1181 break;
1182 }
1183 case 3:
1184 {
1185 switch (oldEdgeIdx)
1186 {
1187 case 0:
1188 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1189 break;
1190 case 4:
1191 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1192 break;
1193 case 3:
1194 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1195 break;
1196 }
1197
1198 break;
1199 }
1200 case 4:
1201 {
1202 switch (oldEdgeIdx)
1203 {
1204 case 1:
1205 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1206 break;
1207 case 2:
1208 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1209 break;
1210 case 5:
1211 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1212 break;
1213 }
1214
1215 break;
1216 }
1217 case 5:
1218 {
1219 switch (oldEdgeIdx)
1220 {
1221 case 1:
1222 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1223 break;
1224 case 2:
1225 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1226 break;
1227 case 3:
1228 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1229 break;
1230 }
1231
1232 break;
1233 }
1234 case 6:
1235 {
1236 switch (oldEdgeIdx)
1237 {
1238 case 1:
1239 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1240 break;
1241 case 4:
1242 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1243 break;
1244 case 5:
1245 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1246 break;
1247 }
1248
1249 break;
1250 }
1251 case 7:
1252 {
1253 switch (oldEdgeIdx)
1254 {
1255 case 1:
1256 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1257 break;
1258 case 4:
1259 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1260 break;
1261 case 3:
1262 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1263 break;
1264 }
1265
1266 break;
1267 }
1268 }
1269
1270 interactionVolume.setEdgePosition(edgeCoord, 1);
1271 }
1272 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsDiag)
1273 {
1274 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1275 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1276 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1277 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1278 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1279 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1280
1281 DimVector crossProductVector1(0);
1282 DimVector crossProductVector2(0);
1283
1284 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1285 crossProductVector1 = centerPos - globalPosFace;
1286 crossProductVector2 = edgeCoord1 - edgeCoord3;
1287 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1288 interactionVolume.setFaceArea(faceArea, 0);
1289
1290 globalPosFace = interactionVolume.getFacePosition(1);
1291 crossProductVector1 = centerPos - globalPosFace;
1292 crossProductVector2 = edgeCoord1 - edgeCoord4;
1293 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1294 interactionVolume.setFaceArea(faceArea, 1);
1295
1296 globalPosFace = interactionVolume.getFacePosition(3);
1297 crossProductVector1 = centerPos - globalPosFace;
1298 crossProductVector2 = edgeCoord1 - edgeCoord6;
1299 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1300 interactionVolume.setFaceArea(faceArea, 3);
1301
1302 globalPosFace = interactionVolume.getFacePosition(5);
1303 crossProductVector1 = centerPos - globalPosFace;
1304 crossProductVector2 = edgeCoord2 - edgeCoord4;
1305 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1306 interactionVolume.setFaceArea(faceArea, 5);
1307
1308 globalPosFace = interactionVolume.getFacePosition(6);
1309 crossProductVector1 = centerPos - globalPosFace;
1310 crossProductVector2 = edgeCoord2 - edgeCoord5;
1311 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1312 interactionVolume.setFaceArea(faceArea, 6);
1313
1314 globalPosFace = interactionVolume.getFacePosition(7);
1315 crossProductVector1 = centerPos - globalPosFace;
1316 crossProductVector2 = edgeCoord2 - edgeCoord6;
1317 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1318 interactionVolume.setFaceArea(faceArea, 7);
1319
1320 globalPosFace = interactionVolume.getFacePosition(8);
1321 crossProductVector1 = centerPos - globalPosFace;
1322 crossProductVector2 = edgeCoord3 - edgeCoord6;
1323 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1324 interactionVolume.setFaceArea(faceArea, 8);
1325
1326 globalPosFace = interactionVolume.getFacePosition(9);
1327 crossProductVector1 = centerPos - globalPosFace;
1328 crossProductVector2 = edgeCoord3 - edgeCoord4;
1329 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1330 interactionVolume.setFaceArea(faceArea, 9);
1331
1332 globalPosFace = interactionVolume.getFacePosition(10);
1333 crossProductVector1 = centerPos - globalPosFace;
1334 crossProductVector2 = edgeCoord4 - edgeCoord5;
1335 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1336 interactionVolume.setFaceArea(faceArea, 10);
1337
1338 globalPosFace = interactionVolume.getFacePosition(11);
1339 crossProductVector1 = centerPos - globalPosFace;
1340 crossProductVector2 = edgeCoord5 - edgeCoord6;
1341 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1342 interactionVolume.setFaceArea(faceArea, 11);
1343
1344 auto element1 = interactionVolume.getSubVolumeElement(0);
1345
1346 bool hasFaceOne = false;
1347 bool hasFaceTwo = false;
1348 for (const auto& intersection : intersections(problem_.gridView(), element1))
1349 {
1350 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 1))
1351 {
1352 if (intersection.neighbor())
1353 {
1354 auto outside = intersection.outside();
1355 if (element1.level() > outside.level())
1356 {
1357 interactionVolume.setSubVolumeElement(outside, 2);
1358 interactionVolume.setSubVolumeElement(outside, 3);
1359
1360 hasFaceOne = true;
1361 if (hasFaceTwo)
1362 break;
1363 }
1364 }
1365 }
1366 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
1367 {
1368 if (intersection.neighbor())
1369 {
1370 auto outside = intersection.outside();
1371 if (element1.level() > outside.level())
1372 {
1373 interactionVolume.setSubVolumeElement(outside, 4);
1374 interactionVolume.setSubVolumeElement(outside, 5);
1375
1376 hasFaceTwo = true;
1377 if (hasFaceOne)
1378 break;
1379 }
1380 }
1381 }
1382 }
1383 }
1384
1385 break;
1386 }
1387 //hanging-node interaction volume of type 2 or 6
1388 case 6:
1389 {
1390 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
1391
1392 std::vector<int> elemIdxOld;
1393 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
1394 {
1395 if (interactionVolume.hasSubVolumeElement(i))
1396 elemIdxOld.push_back(i);
1397 }
1398
1399 std::set<int> zeroFaceIdxVec;
1400 for (int i = 0; i < 6; i++)
1401 {
1402 for (int j = 0; j < 6; j++)
1403 {
1404 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
1405 if (fIdx >= 0)
1406 zeroFaceIdxVec.insert(fIdx);
1407 }
1408 }
1409
1410 std::vector<int> elemIdxNew(6);
1411 int zeroFaceIdx = 0;
1412
1413 bool isFace = true;
1414 std::set<int>::iterator it = zeroFaceIdxVec.begin();
1415 for (; it != zeroFaceIdxVec.end(); ++it)
1416 {
1417 zeroFaceIdx = *it;
1418 isFace = true;
1419 // bool isFace = true;
1420 for (int i = 0; i < 6; i++)
1421 {
1422 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
1423 if (elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
1424 {
1425 isFace = false;
1426 break;
1427 }
1428 }
1429 if (isFace)
1430 {
1431 for (int i = 0; i < 6; i++)
1432 {
1433 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
1434 elemIdxNew[i]);
1435 }
1436
1437 break;
1438 }
1439 }
1440
1441 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
1442 {
1443 for (int i = 0; i < 3; i++)
1444 {
1445 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
1446 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
1447
1448 for (int j = 0; j < 3; j++)
1449 {
1450 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
1451 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
1452
1453 if (faceIdxNew == faceIdxNewTest)
1454 {
1455 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
1456 elemNew, j);
1457 hangingNodeVolume.setIndexOnElement(
1458 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
1459 }
1460 }
1461 }
1462 }
1463
1464 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
1465 {
1466 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
1467 hangingNodeVolume.setFaceArea(interactionVolume.getFaceArea(i), idxNew);
1468 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
1469 }
1470
1471 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
1472 {
1473 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
1474 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
1475 }
1476
1477 interactionVolume = hangingNodeVolume;
1478
1479 interactionVolume.setCenterPosition(centerPos);
1480
1481 interactionVolume.setHangingNodeType(InteractionVolume::sixSmallCells);
1482
1483 auto element3 = interactionVolume.getSubVolumeElement(2);
1484
1485 for (const auto& intersection : intersections(problem_.gridView(), element3))
1486 {
1487 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
1488 {
1489 if (intersection.neighbor())
1490 {
1491 auto outside = intersection.outside();
1492 if (element3.level() > outside.level())
1493 {
1494 interactionVolume.setSubVolumeElement(outside, 6);
1495 interactionVolume.setSubVolumeElement(outside, 7);
1496
1497 break;
1498 }
1499 }
1500 }
1501 }
1502
1503 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1504 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1505 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1506 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1507 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1508 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1509
1510 DimVector crossProductVector1(0);
1511 DimVector crossProductVector2(0);
1512
1513 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1514 crossProductVector1 = centerPos - globalPosFace;
1515 crossProductVector2 = edgeCoord1 - edgeCoord3;
1516 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1517 interactionVolume.setFaceArea(faceArea, 0);
1518
1519 globalPosFace = interactionVolume.getFacePosition(1);
1520 crossProductVector1 = centerPos - globalPosFace;
1521 crossProductVector2 = edgeCoord1 - edgeCoord4;
1522 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1523 interactionVolume.setFaceArea(faceArea, 1);
1524
1525 globalPosFace = interactionVolume.getFacePosition(2);
1526 crossProductVector1 = centerPos - globalPosFace;
1527 crossProductVector2 = edgeCoord1 - edgeCoord5;
1528 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1529 interactionVolume.setFaceArea(faceArea, 2);
1530
1531 globalPosFace = interactionVolume.getFacePosition(3);
1532 crossProductVector1 = centerPos - globalPosFace;
1533 crossProductVector2 = edgeCoord1 - edgeCoord6;
1534 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1535 interactionVolume.setFaceArea(faceArea, 3);
1536
1537 globalPosFace = interactionVolume.getFacePosition(4);
1538 crossProductVector1 = centerPos - globalPosFace;
1539 crossProductVector2 = edgeCoord2 - edgeCoord3;
1540 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1541 interactionVolume.setFaceArea(faceArea, 4);
1542
1543 globalPosFace = interactionVolume.getFacePosition(5);
1544 crossProductVector1 = centerPos - globalPosFace;
1545 crossProductVector2 = edgeCoord2 - edgeCoord4;
1546 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1547 interactionVolume.setFaceArea(faceArea, 5);
1548
1549 globalPosFace = interactionVolume.getFacePosition(7);
1550 crossProductVector1 = centerPos - globalPosFace;
1551 crossProductVector2 = edgeCoord2 - edgeCoord6;
1552 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1553 interactionVolume.setFaceArea(faceArea, 7);
1554
1555 globalPosFace = interactionVolume.getFacePosition(8);
1556 crossProductVector1 = centerPos - globalPosFace;
1557 crossProductVector2 = edgeCoord3 - edgeCoord6;
1558 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1559 interactionVolume.setFaceArea(faceArea, 8);
1560
1561 globalPosFace = interactionVolume.getFacePosition(9);
1562 crossProductVector1 = centerPos - globalPosFace;
1563 crossProductVector2 = edgeCoord3 - edgeCoord4;
1564 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1565 interactionVolume.setFaceArea(faceArea, 9);
1566
1567 globalPosFace = interactionVolume.getFacePosition(10);
1568 crossProductVector1 = centerPos - globalPosFace;
1569 crossProductVector2 = edgeCoord4 - edgeCoord5;
1570 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1571 interactionVolume.setFaceArea(faceArea, 10);
1572
1573 globalPosFace = interactionVolume.getFacePosition(11);
1574 crossProductVector1 = centerPos - globalPosFace;
1575 crossProductVector2 = edgeCoord5 - edgeCoord6;
1576 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1577 interactionVolume.setFaceArea(faceArea, 11);
1578 break;
1579 }
1580 default:
1581 {
1582 interactionVolume.printInteractionVolumeInfo();
1583 DUNE_THROW(Dune::NotImplemented, "Hanging node shape not implemented");
1584 break;
1585 }
1586 }
1587}
1588
1590template<class TypeTag>
1592{
1593 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
1594
1595 //Add elements to the interaction volumes and store element-vertex map
1596 for (const auto& element : elements(problem_.gridView()))
1597 asImp_().storeSubVolumeElements(element, elemVertMap);
1598
1599 for (unsigned int i = 0; i < asImp_().interactionVolumes_.size(); i++)
1600 if (asImp_().interactionVolumes_[i].getElementNumber() == 0)
1601 asImp_().interactionVolumes_[i].printInteractionVolumeInfo();
1602
1603 // Store information related to DUNE intersections for all interaction volumes
1604 for (const auto& element : elements(problem_.gridView()))
1605 asImp_().storeIntersectionInfo(element, elemVertMap);
1606
1607 faceVertices_.clear();
1608 faceVertices_.resize(problem_.gridView().size(0));
1609
1610 // Complete storage of the interaction volumes using the previously stored information
1611 // about the orientation and relationship of the DUNE elements in the interaction volumes (see doc/docextra/3dmpfa)
1612 for (const auto& vertex : vertices(problem_.gridView()))
1613 {
1614 int vIdxGlobal = problem_.variables().index(vertex);
1615
1616 InteractionVolume& interactionVolume = asImp_().interactionVolumes_[vIdxGlobal];
1617
1618 if (interactionVolume.getElementNumber() == 8)
1619 {
1620 asImp_().storeInnerInteractionVolume(interactionVolume, vertex);
1621 }
1622 else if (interactionVolume.isBoundaryInteractionVolume())
1623 {
1624 asImp_().storeBoundaryInteractionVolume(interactionVolume, vertex);
1625 }
1626 //hanging node!
1627 else
1628 {
1629 storeHangingNodeInteractionVolume(interactionVolume, vertex);
1630 }
1631
1632 if (!interactionVolume.isBoundaryInteractionVolume())
1633 {
1634 auto element1 = interactionVolume.getSubVolumeElement(0);
1635 auto element2 = interactionVolume.getSubVolumeElement(1);
1636 auto element3 = interactionVolume.getSubVolumeElement(2);
1637 auto element4 = interactionVolume.getSubVolumeElement(3);
1638 auto element5 = interactionVolume.getSubVolumeElement(4);
1639 auto element6 = interactionVolume.getSubVolumeElement(5);
1640 auto element7 = interactionVolume.getSubVolumeElement(6);
1641 auto element8 = interactionVolume.getSubVolumeElement(7);
1642
1643 int globalIdx1 = problem_.variables().index(element1);
1644 int globalIdx2 = problem_.variables().index(element2);
1645 int globalIdx3 = problem_.variables().index(element3);
1646 int globalIdx4 = problem_.variables().index(element4);
1647 int globalIdx5 = problem_.variables().index(element5);
1648 int globalIdx6 = problem_.variables().index(element6);
1649 int globalIdx7 = problem_.variables().index(element7);
1650 int globalIdx8 = problem_.variables().index(element8);
1651
1652 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 0), globalIdx1, interactionVolume.getIndexOnElement(0, 0));
1653 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 1), globalIdx1, interactionVolume.getIndexOnElement(0, 1));
1654 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 2), globalIdx1, interactionVolume.getIndexOnElement(0, 2));
1655 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 0), globalIdx2, interactionVolume.getIndexOnElement(1, 0));
1656 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 1), globalIdx2, interactionVolume.getIndexOnElement(1, 1));
1657 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 2), globalIdx2, interactionVolume.getIndexOnElement(1, 2));
1658 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 0), globalIdx3, interactionVolume.getIndexOnElement(2, 0));
1659 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 1), globalIdx4, interactionVolume.getIndexOnElement(3, 1));
1660 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 0), globalIdx5, interactionVolume.getIndexOnElement(4, 0));
1661 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 0), globalIdx6, interactionVolume.getIndexOnElement(5, 0));
1662
1663 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 0)].insert(vIdxGlobal);
1664 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 1)].insert(vIdxGlobal);
1665 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 2)].insert(vIdxGlobal);
1666 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 0)].insert(vIdxGlobal);
1667 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 1)].insert(vIdxGlobal);
1668 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 2)].insert(vIdxGlobal);
1669 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 0)].insert(vIdxGlobal);
1670 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 1)].insert(vIdxGlobal);
1671 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 0)].insert(vIdxGlobal);
1672 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 0)].insert(vIdxGlobal);
1673
1674 if (interactionVolume.getHangingNodeType() != InteractionVolume::twoSmallCells)
1675 {
1676 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 1), globalIdx3, interactionVolume.getIndexOnElement(2, 1));
1677 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 2), globalIdx3, interactionVolume.getIndexOnElement(2, 2));
1678 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 0), globalIdx4, interactionVolume.getIndexOnElement(3, 0));
1679 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 2), globalIdx4, interactionVolume.getIndexOnElement(3, 2));
1680 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 1), globalIdx5, interactionVolume.getIndexOnElement(4, 1));
1681 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 2), globalIdx5, interactionVolume.getIndexOnElement(4, 2));
1682 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 1), globalIdx6, interactionVolume.getIndexOnElement(5, 1));
1683 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 2), globalIdx6, interactionVolume.getIndexOnElement(5, 2));
1684 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 0), globalIdx7, interactionVolume.getIndexOnElement(6, 0));
1685 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 1), globalIdx7, interactionVolume.getIndexOnElement(6, 1));
1686 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 2), globalIdx7, interactionVolume.getIndexOnElement(6, 2));
1687 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 0), globalIdx8, interactionVolume.getIndexOnElement(7, 0));
1688 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 1), globalIdx8, interactionVolume.getIndexOnElement(7, 1));
1689 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 2), globalIdx8, interactionVolume.getIndexOnElement(7, 2));
1690
1691 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 1)].insert(vIdxGlobal);
1692 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 2)].insert(vIdxGlobal);
1693 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 0)].insert(vIdxGlobal);
1694 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 2)].insert(vIdxGlobal);
1695 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 1)].insert(vIdxGlobal);
1696 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 2)].insert(vIdxGlobal);
1697 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 1)].insert(vIdxGlobal);
1698 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 2)].insert(vIdxGlobal);
1699 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 0)].insert(vIdxGlobal);
1700 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 1)].insert(vIdxGlobal);
1701 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 2)].insert(vIdxGlobal);
1702 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 0)].insert(vIdxGlobal);
1703 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 1)].insert(vIdxGlobal);
1704 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 2)].insert(vIdxGlobal);
1705 }
1706 }
1707
1708 }
1709}
1710} // end namespace Dumux
1711#endif
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Interactionvolume container for 3-d MPFA L-method.
Class including the information of an interaction volume of a MPFA 3D method that does not change wit...
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:631
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag Indices
Definition: porousmediumflow/sequential/properties.hh:59
Property tag MPFAInteractionVolumeContainer
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:97
Property tag MPFAInteractionVolume
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:96
Interactionvolume container for 3-d MPFA L-method.
Definition: 3dinteractionvolumecontainer.hh:46
typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume
Type for storing an MPFA-interaction-volume. (Usually of type FvMpfaL3dInteractionVolume or FvMpfaL3d...
Definition: 3dinteractionvolumecontainer.hh:93
InteractionVolume & interactionVolume(int vertexIdx)
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:142
Interactionvolume container for 3-d MPFA L-method on an h-adaptive grid.
Definition: 3dinteractionvolumecontaineradaptive.hh:47
void storeInteractionVolumeInfo()
Stores interaction volumes for each grid vertex.
Definition: 3dinteractionvolumecontaineradaptive.hh:1591
std::set< int > & faceVerticeIndices(int eIdxGlobal, int fIdx)
Returns the set of vertices on an element face.
Definition: 3dinteractionvolumecontaineradaptive.hh:104
FvMpfaL3dInteractionVolumeContainerAdaptive(Problem &problem)
Constructs a FvMpfaL3dInteractionVolumeContainerAdaptive object.
Definition: 3dinteractionvolumecontaineradaptive.hh:113
Properties for a MPFA method.
Base file for properties related to sequential IMPET algorithms.