Merge branch 'feature-vf-ext' into 'develop'

Update of view factor generation using 2AI and 2LI methods plus CGAL for ray tracing

See merge request Development/openfoam!551

applications/utilities/preProcessing/viewFactorsGen/Allwmake

@@ -0,0 +1,19 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments # (error catching)
+. ${WM_PROJECT_DIR:?}/wmake/scripts/have_cgal
+
+#------------------------------------------------------------------------------
+unset COMP_FLAGS LINK_FLAGS
+
+if have_cgal
+then
+    echo "    found CGAL -- enabling CGAL support."
+else
+    echo "    did not find CGAL -- disabling CGAL functionality"
+    export COMP_FLAGS="-DNO_CGAL"
+fi
+
+wmake $targetType
+
+#------------------------------------------------------------------------------

applications/utilities/preProcessing/viewFactorsGen/Make/options

@@ -1,10 +1,18 @@
+include $(GENERAL_RULES)/cgal-header-only
+
 EXE_INC = \
+    -Wno-old-style-cast \
+    $(COMP_FLAGS) \
+    ${CGAL_INC} \
+    -DCGAL_HEADER_ONLY \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/distributed/lnInclude \
+    -I$(LIB_SRC)/parallel/distributed/lnInclude
+
 
 EXE_LIBS = \
+    /* ${CGAL_LIBS} */ \
     -lfiniteVolume \
     -lsurfMesh \
     -lmeshTools \

applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H

@@ -0,0 +1,84 @@
+Random rndGen(653213);
+
+// Determine mesh bounding boxes:
+List<treeBoundBox> meshBb
+(
+    1,
+    treeBoundBox
+    (
+        boundBox(coarseMesh.points(), false)
+    ).extend(rndGen, 1e-3)
+);
+
+// Dummy bounds dictionary
+dictionary dict;
+dict.add("bounds", meshBb);
+dict.add
+(
+    "distributionType",
+    distributedTriSurfaceMesh::distributionTypeNames_
+    [
+        distributedTriSurfaceMesh::FROZEN
+    ]
+);
+dict.add("mergeDistance", SMALL);
+
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
+// CGAL surface
+
+distributedTriSurfaceMesh surfacesMesh
+(
+    IOobject
+    (
+        "wallSurface.stl",
+        runTime.constant(),     // directory
+        "triSurface",           // instance
+        runTime,                // registry
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    localSurface,
+    dict
+);
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
+surfacesMesh.setField(triSurfaceToAgglom);
+
+const pointField& pts = surfacesMesh.localPoints();
+
+std::list<Triangle> triangles;
+
+for (const auto& triLocal : surfacesMesh.localFaces())
+{
+    point p1l = pts[triLocal[0]];
+    point p2l = pts[triLocal[1]];
+    point p3l = pts[triLocal[2]];
+
+    Point p1(p1l[0], p1l[1], p1l[2]);
+    Point p2(p2l[0], p2l[1], p2l[2]);
+    Point p3(p3l[0], p3l[1], p3l[2]);
+
+    Triangle tri(p1, p2, p3);
+
+    if (tri.is_degenerate())
+    {
+        std::cout << tri << std::endl;
+    }
+    triangles.push_back(tri);
+}
+
+// constructs AABB tree
+Tree tree(triangles.begin(), triangles.end());

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -1,6 +1,15 @@
 // All rays expressed as start face (local) index end end face (global)
 // Pre-size by assuming a certain percentage is visible.
 
+
+if (Pstream::master())
+{
+    Info<< "\nShooting rays using distributedTriSurfaceMesh (no CGAL support)"
+        << endl;
+}
+
+
+
 // Maximum length for dynamicList
 const label maxDynListLength
 (

applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H

@@ -0,0 +1,131 @@
+// Maximum length for dynamicList
+const label maxDynListLength
+(
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+);
+
+for (const int proci : Pstream::allProcs())
+{
+    std::vector<Point> start;
+    start.reserve(coarseMesh.nFaces());
+
+    std::vector<Point> end(start.size());
+    end.reserve(start.size());
+
+    DynamicList<label> startIndex(start.size());
+    DynamicList<label> endIndex(start.size());
+
+    DynamicList<label> startAgg(start.size());
+    DynamicList<label> endAgg(start.size());
+
+    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
+    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
+
+    const pointField& remoteArea = remoteCoarseSf[proci];
+    const pointField& remoteFc = remoteCoarseCf[proci];
+    const labelField& remoteAgg = remoteCoarseAgg[proci];
+
+    label i = 0;
+    label j = 0;
+    do
+    {
+        for (; i < myFc.size(); i++)
+        {
+            const point& fc = myFc[i];
+            const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
+
+            for (; j < remoteFc.size(); j++)
+            {
+                if (proci != Pstream::myProcNo() || i != j)
+                {
+                    const point& remFc = remoteFc[j];
+                    const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
+                    const vector d(remFc - fc);
+
+
+                    const vector nd = d/mag(d);
+                    const vector nfA = fA/mag(fA);
+                    const vector nremA = remA/mag(remA);
+
+                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
+                    {
+                        vector direction(d[0], d[1], d[2]);
+                        vector s(fc[0], fc[1], fc[2]);
+                        vector rayEnd(s + (1-intTol)*direction);
+                        end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
+
+                        s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
+
+                        start.push_back(Point(s[0], s[1], s[2]));
+                        startIndex.append(i);
+                        if (useAgglomeration)
+                        {
+                            startAgg.append
+                            (
+                                globalNumbering.toGlobal(proci, fAgg)
+                            );
+                            endAgg.append
+                            (
+                                globalNumbering.toGlobal(proci, remAgg)
+                            );
+                        }
+
+                        label globalI = globalNumbering.toGlobal(proci, j);
+                        endIndex.append(globalI);
+
+                        if (startIndex.size() > maxDynListLength)
+                        {
+                            FatalErrorInFunction
+                                << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
+                        }
+                    }
+                }
+            }
+
+            if (j == remoteFc.size())
+            {
+                j = 0;
+            }
+        }
+
+    } while (i < myFc.size());
+
+    label totalRays(startIndex.size());
+    reduce(totalRays, sumOp<scalar>());
+
+    if (debug)
+    {
+        Pout<< "Number of rays :" << totalRays << endl;
+    }
+
+    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+    {
+        Segment ray(start[rayI], end[rayI]);
+        Segment_intersection intersects = tree.any_intersection(ray);
+
+        if (!intersects)
+        {
+            rayStartFace.append(startIndex[rayI]);
+            rayEndFace.append(endIndex[rayI]);
+        }
+    }
+    start.clear();
+
+    if (debug)
+    {
+        Pout << "hits : " <<  rayStartFace.size()<< endl;
+    }
+
+    startIndex.clear();
+    end.clear();
+    endIndex.clear();
+    startAgg.clear();
+    endAgg.clear();
+}

applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

@@ -31,14 +31,34 @@ Group
     grpPreProcessingUtilities
 
 Description
-    View factors are calculated based on a face agglomeration array
+    This view factors generation application uses a combined approach of
-    (finalAgglom generated by faceAgglomerate utility).
+    double area integral (2AI) and double linear integral (2LI). 2AI is used
+    when the two surfaces are 'far' apart and 2LI when they are 'close'.
+    2LI is integrated along edges using Gaussian quadrature.
+    The distance between faces is calculating a ratio between averaged areas
+    and the distance between face centres.
 
-    Each view factor between the agglomerated faces i and j (Fij) is calculated
+    The input from viewFactorsDict are:
-    using a double integral of the sub-areas composing the agglomeration.
 
-    The patches involved in the view factor calculation are taken from the
+        GaussQuadTol              0.1;      // GaussQuad  error
-    boundary file and should be part on the group viewFactorWall. ie.:
+        distTol                   8;        // R/Average(rm)
+        alpha                     0.22;     // Use for common edges for 2LI
+
+
+    For debugging purposes, the following entries can be set in viewFactorsDict:
+
+        writeViewFactorMatrix     true;
+        writeFacesAgglomeration   false;
+        dumpRays		          false;
+
+
+        writeViewFactorMatrix   writes the sum of the VF on each face.
+        writeFacesAgglomeration writes the agglomeration
+        dumpRays                dumps rays
+
+
+    The participating patches in the VF calculation have to be in the
+    'viewFactorWall' patch group (in the polyMesh/boundary file), e.g.
 
     floor
     {
@@ -48,6 +68,9 @@ Description
         startFace       3100;
     }
 
+    Compile with -DNO_CGAL only if no CGAL present - CGAL AABB tree performs
+    better than the built-in octree.
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
@@ -56,21 +79,52 @@ Description
 #include "surfaceFields.H"
 #include "distributedTriSurfaceMesh.H"
 #include "meshTools.H"
+#include "constants.H"
 
 #include "indirectPrimitivePatch.H"
 #include "DynamicField.H"
-#include "unitConversion.H"
 
 #include "scalarMatrices.H"
 #include "labelListIOList.H"
 #include "scalarListIOList.H"
 
 #include "singleCellFvMesh.H"
-
 #include "IOmapDistribute.H"
 
-using namespace Foam;
+#ifndef NO_CGAL
 
+// Silence boost bind deprecation warnings (before CGAL-5.2.1)
+#include "CGAL/version.h"
+#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
+#define BOOST_BIND_GLOBAL_PLACEHOLDERS
+#endif
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/Surface_mesh.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point;
+typedef K::Direction_3 Vector3C;
+typedef K::Triangle_3 Triangle;
+typedef K::Segment_3 Segment;
+
+typedef std::list<Triangle>::iterator Iterator;
+typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+typedef boost::optional
+<
+    Tree::Intersection_and_primitive_id<Segment>::Type
+> Segment_intersection;
+
+#endif // NO_CGAL
+
+using namespace Foam;
+using namespace Foam::constant;
+using namespace Foam::constant::mathematical;
 
 triSurface triangulate
 (
@@ -127,20 +181,56 @@ triSurface triangulate
         newPatchI++;
     }
 
-    //striSurfaceToAgglom.resize(localTriFaceI-1);
+    //triSurfaceToAgglom.resize(localTriFaceI-1);
 
     triangles.shrink();
+    triSurface surface(triangles, mesh.points());
+    surface.compactPoints();
 
-    // Create globally numbered tri surface
-    triSurface rawSurface(triangles, mesh.points());
 
-    // Create locally numbered tri surface
+#ifndef NO_CGAL
-    triSurface surface
+
+    // CGAL : every processor has whole surface
+
+    globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
+    globalIndex globalPointIdx
     (
-        rawSurface.localFaces(),
+        surface.points().size(), globalIndex::gatherOnly()
-        rawSurface.localPoints()
     );
 
+    List<labelledTri> globalSurfaceTris(globalFaceIdx.gather(surface));
+    pointField globalSurfacePoints(globalPointIdx.gather(surface.points()));
+
+    //label offset = 0;
+    for (const label proci : globalPointIdx.allProcs())
+    {
+        const label offset = globalPointIdx.localStart(proci);
+
+        if (offset)
+        {
+            for
+            (
+                labelledTri& tri
+             :  globalSurfaceTris.slice(globalFaceIdx.range(proci))
+            )
+            {
+                tri[0] += offset;
+                tri[1] += offset;
+                tri[2] += offset;
+            }
+        }
+    }
+
+    surface =
+        triSurface
+        (
+            std::move(globalSurfaceTris),
+            std::move(globalSurfacePoints)
+        );
+
+    Pstream::broadcast(surface);
+#endif
+
     // Add patch names to surface
     surface.patches().setSize(newPatchI);
 
@@ -193,7 +283,7 @@ void writeRays
 }
 
 
-scalar calculateViewFactorFij
+scalar calculateViewFactorFij2AI
 (
     const vector& i,
     const vector& j,
@@ -250,6 +340,74 @@ void insertMatrixElements
 }
 
 
+scalar GaussQuad
+(
+
+    const scalarList& w,
+    const scalarList& p,
+    const scalar& magSi,
+    const scalar& magSj,
+    const vector& di,
+    const vector& dj,
+    const vector& ci,
+    const vector& cj,
+    const scalar cosij,
+    const scalar alpha,
+    label gi
+)
+{
+    scalar dIntFij = 0;
+    if (gi == 0)
+    {
+        vector r(ci - cj);
+        if (mag(r) < SMALL)
+        {
+            r = (alpha*magSi)*di;
+        }
+        dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
+    }
+    else
+    {
+        List<vector> pi(w.size());
+        forAll (pi, i)
+        {
+            pi[i] = ci + p[i]*(magSi/2)*di;
+        }
+
+        List<vector> pj(w.size());
+        forAll (pj, i)
+        {
+            pj[i] = cj + p[i]*(magSj/2)*dj;
+        }
+
+        forAll (w, i)
+        {
+            forAll (w, j)
+            {
+                vector r(pi[i] - pj[j]);
+
+                if (mag(r) < SMALL)
+                {
+                    r = (alpha*magSi)*di;
+                     dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+                else
+                {
+                    dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+
+            }
+        }
+
+        dIntFij *= (magSi/2) * (magSj/2);
+
+    }
+    return dIntFij;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -273,7 +431,7 @@ int main(int argc, char *argv[])
             "viewFactorsDict",
             runTime.constant(),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE
        )
     );
@@ -288,6 +446,23 @@ int main(int argc, char *argv[])
 
     const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
 
+    const scalar GaussQuadTol =
+        viewFactorDict.getOrDefault<scalar>("GaussQuadTol", 0.01);
+
+    const scalar distTol =
+         viewFactorDict.getOrDefault<scalar>("distTol", 8);
+
+    const scalar alpha =
+         viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
+
+    const scalar intTol =
+         viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
+
+    bool useAgglomeration(true);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const labelList viewFactorsPatches(patches.indices(viewFactorWall));
+
     // Read agglomeration map
     labelListIOList finalAgglom
     (
@@ -296,12 +471,22 @@ int main(int argc, char *argv[])
             "finalAgglom",
             mesh.facesInstance(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     );
 
+    if (!finalAgglom.typeHeaderOk<labelListIOList>())
+    {
+        finalAgglom.setSize(patches.size());
+        for (label patchi=0;  patchi < patches.size(); patchi++)
+        {
+            finalAgglom[patchi] = identity(patches[patchi].size());
+        }
+        useAgglomeration = false;
+    }
+
     // Create the coarse mesh  using agglomeration
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -336,10 +521,8 @@ int main(int argc, char *argv[])
     label nCoarseFaces = 0;      //total number of coarse faces
     label nFineFaces = 0;        //total number of fine faces
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
     const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
 
-    labelList viewFactorsPatches(patches.indices(viewFactorWall));
     for (const label patchi : viewFactorsPatches)
     {
         nCoarseFaces += coarsePatches[patchi].size();
@@ -465,7 +648,13 @@ int main(int argc, char *argv[])
 
     // Set up searching engine for obstacles
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifdef NO_CGAL
+    // Using octree
     #include "searchingEngine.H"
+#else
+    // Using CGAL aabbtree (faster, more robust)
+    #include "searchingEngine_CGAL.H"
+#endif
 
     // Determine rays between coarse face centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -473,11 +662,15 @@ int main(int argc, char *argv[])
 
     DynamicList<label> rayEndFace(rayStartFace.size());
 
-
     // Return rayStartFace in local index and rayEndFace in global index
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+#ifdef NO_CGAL
+    // Using octree, distributedTriSurfaceMesh
     #include "shootRays.H"
+#else
+    // Using CGAL aabbtree (faster, more robust)
+    #include "shootRays_CGAL.H"
+#endif
 
     // Calculate number of visible faces from local index
     labelList nVisibleFaceFaces(nCoarseFaces, Zero);
@@ -518,9 +711,9 @@ int main(int argc, char *argv[])
         visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
     }
 
-
     // Construct data in compact addressing
-    // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2LI) need edges (li)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     pointField compactCoarseCf(map.constructSize(), Zero);
@@ -528,12 +721,16 @@ int main(int argc, char *argv[])
     List<List<point>> compactFineSf(map.constructSize());
     List<List<point>> compactFineCf(map.constructSize());
 
+    DynamicList<List<point>> compactPoints(map.constructSize());
+
     DynamicList<label> compactPatchId(map.constructSize());
 
     // Insert my coarse local values
     SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
     SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
 
+    const faceList& faces = mesh.faces();
+
     // Insert my fine local values
     label compactI = 0;
     forAll(viewFactorsPatches, i)
@@ -548,11 +745,21 @@ int main(int argc, char *argv[])
             const labelList& coarsePatchFace =
                 coarseMesh.patchFaceMap()[patchID];
 
+            const polyPatch& pp = patches[patchID];
+
             forAll(coarseToFine, coarseI)
             {
                 compactPatchId.append(patchID);
                 List<point>& fineCf = compactFineCf[compactI];
-                List<point>& fineSf = compactFineSf[compactI++];
+                List<point>& fineSf = compactFineSf[compactI];
+
+                label startFace = pp.start();
+
+                const vectorField locPoints
+                (
+                    mesh.points(),
+                    faces[coarseI + startFace]
+                );
 
                 const label coarseFaceI = coarsePatchFace[coarseI];
                 const labelList& fineFaces = coarseToFine[coarseFaceI];
@@ -560,6 +767,8 @@ int main(int argc, char *argv[])
                 fineCf.setSize(fineFaces.size());
                 fineSf.setSize(fineFaces.size());
 
+                compactPoints.append(locPoints);
+
                 fineCf = UIndirectList<point>
                 (
                     mesh.Cf().boundaryField()[patchID],
@@ -570,19 +779,25 @@ int main(int argc, char *argv[])
                     mesh.Sf().boundaryField()[patchID],
                     coarseToFine[coarseFaceI]
                 );
+
+                compactI++;
             }
         }
     }
 
+    if (Pstream::master() && debug)
+    {
+        Info<< "map distribute..."  << endl;
+    }
+
     // Do all swapping
     map.distribute(compactCoarseSf);
     map.distribute(compactCoarseCf);
     map.distribute(compactFineCf);
     map.distribute(compactFineSf);
-
+    map.distribute(compactPoints);
     map.distribute(compactPatchId);
 
-
     // Plot all rays between visible faces.
     if (dumpRays)
     {
@@ -598,8 +813,7 @@ int main(int argc, char *argv[])
 
     // Fill local view factor matrix
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+    scalarListIOList F2LI
-    scalarListIOList F
     (
         IOobject
         (
@@ -630,6 +844,61 @@ int main(int argc, char *argv[])
         Info<< "\nCalculating view factors..." << endl;
     }
 
+    FixedList<scalarList, 5> GaussPoints;
+    GaussPoints[0].setSize(1);
+    GaussPoints[0] = 0;
+
+    GaussPoints[1].setSize(2);
+    GaussPoints[1][0] =  1/std::sqrt(3);
+    GaussPoints[1][1] = -1/std::sqrt(3);
+
+    GaussPoints[2].setSize(3);
+    GaussPoints[2][0] =  0;
+    GaussPoints[2][1] =  std::sqrt(3.0/5.0);
+    GaussPoints[2][2] = -std::sqrt(3.0/5.0);
+
+    GaussPoints[3].setSize(4);
+    GaussPoints[3][0] = std::sqrt(3.0/7.0 - (2.0/7.0)*std::sqrt(6.0/5.0));
+    GaussPoints[3][1] = -GaussPoints[3][0];
+    GaussPoints[3][2] = std::sqrt(3.0/7.0 + (2.0/7.0)*std::sqrt(6.0/5.0));
+    GaussPoints[3][3] = -GaussPoints[3][2];
+
+    GaussPoints[4].setSize(5);
+    GaussPoints[4][0] =  0;
+    GaussPoints[4][1] = (1.0/3.0)*std::sqrt(5.0 - 2.0*std::sqrt(10.0/7.0));
+    GaussPoints[4][2] = -GaussPoints[4][1];
+    GaussPoints[4][3] = (1.0/3.0)*std::sqrt(5.0 + 2.0*std::sqrt(10.0/7.0));
+    GaussPoints[4][4] = -GaussPoints[4][3];
+
+
+    List<scalarList> GaussWeights(5);
+    GaussWeights[0].setSize(1);
+    GaussWeights[0] = 2;
+
+    GaussWeights[1].setSize(2);
+    GaussWeights[1][0] =  1;
+    GaussWeights[1][1] =  1;
+
+    GaussWeights[2].setSize(3);
+    GaussWeights[2][0] =  8.0/9.0;
+    GaussWeights[2][1] =  5.0/9.0;
+    GaussWeights[2][2] =  5.0/9.0;
+
+    GaussWeights[3].setSize(4);
+    GaussWeights[3][0] =  (18.0 + std::sqrt(30))/36.0;
+    GaussWeights[3][1] =  (18.0 + std::sqrt(30))/36.0;
+    GaussWeights[3][2] =  (18.0 - std::sqrt(30))/36.0;
+    GaussWeights[3][3] =  (18.0 - std::sqrt(30))/36.0;
+
+    GaussWeights[4].setSize(5);
+    GaussWeights[4][0] =  128.0/225.0;
+    GaussWeights[4][1] =  (322.0 + 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][2] =  (322.0 + 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][3] =  (322.0 - 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][4] =  (322.0 - 13.0*std::sqrt(70))/900.0;
+
+    const label maxQuadOrder = 5;
+
     if (mesh.nSolutionD() == 3)
     {
         forAll(localCoarseSf, coarseFaceI)
@@ -638,119 +907,207 @@ int main(int argc, char *argv[])
             const vector Ai = sum(localFineSf);
             const List<point>& localFineCf = compactFineCf[coarseFaceI];
             const label fromPatchId = compactPatchId[coarseFaceI];
+
+            const List<point>& lPoints = compactPoints[coarseFaceI];
+
             patchArea[fromPatchId] += mag(Ai);
 
             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 
             forAll(visCoarseFaces, visCoarseFaceI)
             {
-                F[coarseFaceI].setSize(visCoarseFaces.size());
+                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
+                F2LI[coarseFaceI].setSize(visCoarseFaces.size());
                 label compactJ = visCoarseFaces[visCoarseFaceI];
                 const List<point>& remoteFineSj = compactFineSf[compactJ];
                 const List<point>& remoteFineCj = compactFineCf[compactJ];
 
+                const List<point>& rPointsCj = compactPoints[compactJ];
+
                 const label toPatchId = compactPatchId[compactJ];
 
-                scalar Fij = 0;
+                bool far(false);
+                // Relative distance
                 forAll(localFineSf, i)
                 {
-                    const vector& dAi = localFineSf[i];
+                    const scalar dAi =
+                        Foam::sqrt
+                        (
+                            mag(localFineSf[i])/constant::mathematical::pi
+                        );
                     const vector& dCi = localFineCf[i];
 
                     forAll(remoteFineSj, j)
                     {
-                        const vector& dAj = remoteFineSj[j];
+                        const scalar dAj =
+                            Foam::sqrt
+                            (
+                                mag(remoteFineSj[j])/constant::mathematical::pi
+                            );
                         const vector& dCj = remoteFineCj[j];
 
-                        scalar dIntFij = calculateViewFactorFij
+                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
-                        (
-                            dCi,
-                            dCj,
-                            dAi,
-                            dAj
-                        );
 
-                        Fij += dIntFij;
+                        if (dist > distTol)
+                        {
+                            far = true;
+                        }
                     }
                 }
-                F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
+
-                sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
+                if (far)
+                {
+                    // 2AI method
+                    scalar F2AIij = 0;
+
+                    forAll(localFineSf, i)
+                    {
+                        const vector& dAi = localFineSf[i];
+                        const vector& dCi = localFineCf[i];
+
+                        forAll(remoteFineSj, j)
+                        {
+                            const vector& dAj = remoteFineSj[j];
+                            const vector& dCj = remoteFineCj[j];
+
+                            scalar dIntFij = calculateViewFactorFij2AI
+                            (
+                                dCi,
+                                dCj,
+                                dAi,
+                                dAj
+                            );
+
+                            F2AIij += dIntFij;
+                        }
+                    }
+                    F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
+                }
+                else
+                {
+                    // 2LI method
+                    label nLocal = lPoints.size();
+                    label nRemote = rPointsCj.size();
+
+                    // Using sub-divisions (quadrature)
+                    scalar oldEToeInt = 0;
+                    for (label gi=0; gi < maxQuadOrder; gi++)
+                    {
+                        scalar F2LIij = 0;
+                        for(label i=0; i<nLocal; i++)
+                        {
+                            vector si;
+                            vector ci;
+
+                            vector sj;
+                            vector cj;
+
+                            if (i == 0)
+                            {
+                                si = lPoints[i] - lPoints[nLocal-1];
+                                ci = (lPoints[i] + lPoints[nLocal-1])/2;
+                            }
+                            else
+                            {
+                                si = lPoints[i] - lPoints[i-1];
+                                ci = (lPoints[i] + lPoints[i-1])/2;
+                            }
+
+                            for(label j=0; j<nRemote; j++)
+                            {
+                                if (j == 0)
+                                {
+                                    sj = rPointsCj[j]-rPointsCj[nRemote-1];
+                                    cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
+                                }
+                                else
+                                {
+                                    sj = rPointsCj[j] - rPointsCj[j-1];
+                                    cj = (rPointsCj[j] + rPointsCj[j-1])/2;
+                                }
+
+
+                                scalar magSi = mag(si);
+                                scalar magSj = mag(sj);
+                                scalar cosij = (si & sj)/(magSi * magSj);
+
+                                vector di = si/magSi;
+                                vector dj = sj/magSj;
+
+                                label quadOrder = gi;
+                                const vector r(ci - cj);
+                                // Common edges use n = 0
+                                if (mag(r) < SMALL)
+                                {
+                                    quadOrder = 0;
+                                }
+
+                                scalar dIntFij =
+                                    GaussQuad
+                                    (
+                                        GaussWeights[gi],
+                                        GaussPoints[gi],
+                                        magSi,
+                                        magSj,
+                                        di,
+                                        dj,
+                                        ci,
+                                        cj,
+                                        cosij,
+                                        alpha,
+                                        quadOrder
+                                    );
+
+                                F2LIij += dIntFij;
+                            }
+                        }
+
+                        scalar err = (F2LIij-oldEToeInt)/F2LIij;
+
+                        if
+                        (
+                            (mag(err) < GaussQuadTol && gi > 0)
+                         || gi == maxQuadOrder-1
+                        )
+                        {
+                            F2LI[coarseFaceI][visCoarseFaceI] =
+                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
+                            break;
+                        }
+                        else
+                        {
+                            oldEToeInt = F2LIij;
+                        }
+                    }
+                }
+
+                sumViewFactorPatch[fromPatchId][toPatchId] +=
+                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
             }
         }
     }
-    else if (mesh.nSolutionD() == 2)
+    else
     {
-        const boundBox& box = mesh.bounds();
+         FatalErrorInFunction
-        const Vector<label>& dirs = mesh.geometricD();
+            << " View factors are not available in 2D "
-        vector emptyDir = Zero;
+            << exit(FatalError);
-        forAll(dirs, i)
-        {
-            if (dirs[i] == -1)
-            {
-                emptyDir[i] = 1.0;
-            }
-        }
-
-        scalar wideBy2 = (box.span() & emptyDir)*2.0;
-
-        forAll(localCoarseSf, coarseFaceI)
-        {
-            const vector& Ai = localCoarseSf[coarseFaceI];
-            const vector& Ci = localCoarseCf[coarseFaceI];
-            vector Ain = Ai/mag(Ai);
-            vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
-            vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
-
-            const label fromPatchId = compactPatchId[coarseFaceI];
-            patchArea[fromPatchId] += mag(Ai);
-
-            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
-            forAll(visCoarseFaces, visCoarseFaceI)
-            {
-                F[coarseFaceI].setSize(visCoarseFaces.size());
-                label compactJ = visCoarseFaces[visCoarseFaceI];
-                const vector& Aj = compactCoarseSf[compactJ];
-                const vector& Cj = compactCoarseCf[compactJ];
-
-                const label toPatchId = compactPatchId[compactJ];
-
-                vector Ajn = Aj/mag(Aj);
-                vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
-                vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
-
-                scalar d1 = mag(R1i - R2j);
-                scalar d2 = mag(R2i - R1j);
-                scalar s1 = mag(R1i - R1j);
-                scalar s2 = mag(R2i - R2j);
-
-                scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
-
-                F[coarseFaceI][visCoarseFaceI] = Fij;
-                sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
-            }
-        }
-    }
-
-    if (Pstream::master())
-    {
-        Info << "Writing view factor matrix..." << endl;
     }
 
     // Write view factors matrix in listlist form
-    F.write();
+    F2LI.write();
 
     reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
     reduce(patchArea, sumOp<scalarList>());
 
-
     if (Pstream::master() && debug)
     {
         forAll(viewFactorsPatches, i)
         {
             label patchI =  viewFactorsPatches[i];
-            forAll(viewFactorsPatches, i)
+            for (label j=i; j<viewFactorsPatches.size(); j++)
             {
-                label patchJ =  viewFactorsPatches[i];
+                label patchJ =  viewFactorsPatches[j];
+
                 Info << "F" << patchI << patchJ << ": "
                      << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
                      << endl;
@@ -758,9 +1115,13 @@ int main(int argc, char *argv[])
         }
     }
 
-
     if (writeViewFactors)
     {
+        if (Pstream::master())
+        {
+            Info << "Writing view factor matrix..." << endl;
+        }
+
         volScalarField viewFactorField
         (
             IOobject
@@ -791,13 +1152,14 @@ int main(int argc, char *argv[])
 
                 forAll(coarseToFine, coarseI)
                 {
-                    const scalar Fij = sum(F[compactI]);
+                    const scalar FiSum = sum(F2LI[compactI]);
+
                     const label coarseFaceID = coarsePatchFace[coarseI];
                     const labelList& fineFaces = coarseToFine[coarseFaceID];
                     forAll(fineFaces, fineId)
                     {
                         const label faceID = fineFaces[fineId];
-                        vfbf[patchID][faceID] = Fij;
+                        vfbf[patchID][faceID] = FiSum;
                     }
                     compactI++;
                 }

src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C

@@ -35,6 +35,8 @@ License
 #include "boundaryRadiationProperties.H"
 #include "lduCalculatedProcessorField.H"
 
+
+
 using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -55,7 +57,36 @@ const Foam::word Foam::radiation::viewFactor::viewFactorWalls
 
 void Foam::radiation::viewFactor::initialise()
 {
-    const polyBoundaryMesh& coarsePatches = coarseMesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    if (!finalAgglom_.typeHeaderOk<labelListIOList>())
+    {
+        finalAgglom_.setSize(patches.size());
+        for (label patchi=0;  patchi < patches.size(); patchi++)
+        {
+            finalAgglom_[patchi] = identity(patches[patchi].size());
+        }
+    }
+
+    coarseMesh_.reset
+    (
+        new singleCellFvMesh
+        (
+            IOobject
+            (
+                "coarse:" + mesh_.name(),
+                mesh_.polyMesh::instance(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            finalAgglom_
+        )
+    );
+
+    const polyBoundaryMesh& coarsePatches = coarseMesh_->boundaryMesh();
 
     selectedPatches_ = mesh_.boundaryMesh().indices(viewFactorWalls);
 
@@ -80,6 +111,8 @@ void Foam::radiation::viewFactor::initialise()
 
     useDirect_ = coeffs_.getOrDefault<bool>("useDirectSolver", true);
 
+
+
     map_.reset
     (
         new IOmapDistribute
@@ -186,7 +219,7 @@ void Foam::radiation::viewFactor::initialise()
         labelList upper(rays_.size(), -1);
         labelList lower(rays_.size(), -1);
 
-        const edgeList& raysLst = rays_.sortedToc();
+        const edgeList raysLst(rays_.sortedToc());
         label rayI = 0;
         for (const auto& e : raysLst)
         {
@@ -408,8 +441,8 @@ void Foam::radiation::viewFactor::initialise()
             totalDelta /= myF.size();
             reduce(totalDelta, sumOp<scalar>());
             reduce(maxDelta, maxOp<scalar>());
-            Info << "Smoothng average delta : " << totalDelta << endl;
+            Info << "Smoothing average delta : " << totalDelta << endl;
-            Info << "Smoothng maximum delta : " << maxDelta << nl << endl;
+            Info << "Smoothing maximum delta : " << maxDelta << nl << endl;
         }
     }
 
@@ -510,26 +543,26 @@ Foam::radiation::viewFactor::viewFactor(const volScalarField& T)
             "finalAgglom",
             mesh_.facesInstance(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     ),
     map_(),
-    coarseMesh_
+    coarseMesh_(),
-    (
+//     (
-        IOobject
+//         IOobject
-        (
+//         (
-            "coarse:" + mesh_.name(),
+//             "coarse:" + mesh_.name(),
-            mesh_.polyMesh::instance(),
+//             mesh_.polyMesh::instance(),
-            mesh_.time(),
+//             mesh_.time(),
-            IOobject::NO_READ,
+//             IOobject::NO_READ,
-            IOobject::NO_WRITE,
+//             IOobject::NO_WRITE,
-            false
+//             false
-        ),
+//         ),
-        mesh_,
+//         mesh_,
-        finalAgglom_
+//         finalAgglom_
-    ),
+//     ),
     qr_
     (
         IOobject
@@ -573,26 +606,26 @@ Foam::radiation::viewFactor::viewFactor
             "finalAgglom",
             mesh_.facesInstance(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     ),
     map_(),
-    coarseMesh_
+    coarseMesh_(),
-    (
+//     (
-        IOobject
+//         IOobject
-        (
+//         (
-            "coarse:" + mesh_.name(),
+//             "coarse:" + mesh_.name(),
-            mesh_.polyMesh::instance(),
+//             mesh_.polyMesh::instance(),
-            mesh_.time(),
+//             mesh_.time(),
-            IOobject::NO_READ,
+//             IOobject::NO_READ,
-            IOobject::NO_WRITE,
+//             IOobject::NO_WRITE,
-            false
+//             false
-        ),
+//         ),
-        mesh_,
+//         mesh_,
-        finalAgglom_
+//         finalAgglom_
-    ),
+//     ),
     qr_
     (
         IOobject
@@ -719,10 +752,10 @@ void Foam::radiation::viewFactor::calculate()
             const tmp<scalarField> tHoi = qrp.qro(bandI);
             const scalarField& Hoi = tHoi();
 
-            const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
+            const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
 
             const labelList& coarsePatchFace =
-                coarseMesh_.patchFaceMap()[patchID];
+                coarseMesh_->patchFaceMap()[patchID];
 
             scalarList T4ave(pp.size(), 0.0);
             scalarList Eave(pp.size(), 0.0);
@@ -820,7 +853,7 @@ void Foam::radiation::viewFactor::calculate()
             // Local matrix coefficients
             if (!constEmissivity_ || iterCounter_ == 0)
             {
-                const edgeList& raysLst = rays_.sortedToc();
+                const edgeList raysLst(rays_.sortedToc());
 
                 label rayI = 0;
                 for (const auto& e : raysLst)
@@ -1044,7 +1077,7 @@ void Foam::radiation::viewFactor::calculate()
     label globCoarseId = 0;
     for (const label patchID : selectedPatches_)
     {
-        const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
+        const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
 
         if (pp.size() > 0)
         {
@@ -1056,7 +1089,7 @@ void Foam::radiation::viewFactor::calculate()
             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 
             const labelList& coarsePatchFace =
-                coarseMesh_.patchFaceMap()[patchID];
+                coarseMesh_->patchFaceMap()[patchID];
 
             //scalar heatFlux = 0.0;
             forAll(coarseToFine, coarseI)

src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -107,7 +107,7 @@ protected:
         autoPtr<IOmapDistribute> map_;
 
         //- Coarse mesh
-        singleCellFvMesh coarseMesh_;
+        autoPtr<singleCellFvMesh> coarseMesh_;
 
         //- Net radiative heat flux [W/m2]
         volScalarField qr_;

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun

@@ -9,11 +9,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 #-- Run on single processor
 
 # Agglomerate patch faces
-for region in air
+#for region in air
-do
+#do
-    runApplication -s $region \
+#    runApplication -s $region \
-        faceAgglomerate -region $region -dict constant/viewFactorsDict
+#        faceAgglomerate -region $region -dict constant/viewFactorsDict
-done
+#done
 
 # Generate view factors
 for region in air

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun-parallel

@@ -12,11 +12,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 runApplication decomposePar -allRegions -constant
 
 # Agglomerate patch faces
-for region in air
+#for region in air
-do
+#do
-    runParallel -s $region \
+#    runParallel -s $region \
-        faceAgglomerate -region $region -dict constant/viewFactorsDict
+#        faceAgglomerate -region $region -dict constant/viewFactorsDict
-done
+#done
 
 # Generate view factors
 for region in air

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/viewFactorsDict

@@ -15,11 +15,10 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
 writePatchViewFactors     false;
 // dumpRays    true;
 
-maxDynListLength     200000;
+maxDynListLength          200000;
 
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun

@@ -9,11 +9,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 #-- Run on single processor
 
 # Agglomerate patch faces
-for region in bottomAir topAir
+#for region in bottomAir topAir
-do
+#do
-    runApplication -s "$region" \
+#    runApplication -s "$region" \
-        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
+#        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
-done
+#done
 
 # Generate view factors
 for region in bottomAir topAir

tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun-parallel

@@ -12,11 +12,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 runApplication decomposePar -allRegions
 
 # Agglomerate patch faces
-for region in bottomAir topAir
+#for region in bottomAir topAir
-do
+#do
-    runParallel -s "$region" \
+#    runParallel -s "$region" \
-        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
+#        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
-done
+#done
 
 # Generate view factors
 for region in bottomAir topAir

tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/bottomAir/viewFactorsDict

@@ -15,56 +15,6 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
-writePatchViewFactors     false;
-
-bottomAir_to_heater
-{
-    nFacesInCoarsestLevel     30;
-    featureAngle              10;
-}
-
-bottomAir_to_leftSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-bottomAir_to_rightSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-minX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-minY
-{
-    nFacesInCoarsestLevel     30;
-    featureAngle              10;
-}
-
-minZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-maxX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
 
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/topAir/viewFactorsDict

@@ -15,56 +15,5 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
-writePatchViewFactors     false;
-
-topAir_to_heater
-{
-    nFacesInCoarsestLevel     24;
-    featureAngle              10;
-}
-
-topAir_to_leftSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-topAir_to_rightSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-minX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxY
-{
-    nFacesInCoarsestLevel     40;
-    featureAngle              10;
-}
-
-minZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-maxX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
 
 // ************************************************************************* //