Initial implementation of templated frontal marking algorithm

git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@27 fdcce57e-7e00-11e2-b579-49867b4cea03

meshLibrary/utilities/helperFunctions/helperFunctions.H

@@ -41,6 +41,7 @@ SourceFiles
 #include "helperFunctionsTopologyManipulation.H"
 #include "helperFunctionsStringConversion.H"
 #include "helperFunctionsPar.H"
+#include "helperFunctionsFrontalMarking.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

meshLibrary/utilities/helperFunctions/helperFunctionsFrontalMarking.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2005-2007 Franjo Juretic
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "helperFunctionsFrontalMarking.H"
+#include "DynList.H"
+#include "labelPair.H"
+#include "HashSet.H"
+
+# ifdef USE_OMP
+#include <omp.h>
+# endif
+
+#define DEBUGFrontalMarking
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+namespace help
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+template<class labelListType, class neiOp, class filterOp>
+void frontalMarking
+(
+    labelListType& result,
+    const label startingIndex,
+    const neiOp& neighbourCalculator,
+    const filterOp& selector
+)
+{
+    //- add the starting element
+    result.clear();
+    result.append(startingIndex);
+
+    //- add the starting element to the front
+    labelListPMG front;
+    front.append(startingIndex);
+
+    //- store information which element were already visited
+    boolList alreadySelected(neighbourCalculator.size(), false);
+
+    //- start with frontal marking
+    while( front.size() )
+    {
+        const label eLabel = front.removeLastElement();
+
+        //- find neighbours of the current element
+        DynList<label> neighbours;
+        neighbourCalculator(eLabel, neighbours);
+
+        forAll(neighbours, neiI)
+        {
+            const label nei = neighbours[neiI];
+
+            if( nei < 0 )
+                continue;
+            if( alreadySelected[nei] )
+                continue;
+
+            if( selector(nei) )
+            {
+                alreadySelected[nei] = true;
+                front.append(nei);
+                result.append(nei);
+            }
+        }
+    }
+}
+
+template<class labelListType, class neiOp, class filterOp>
+label groupMarking
+(
+    labelListType& elementInGroup,
+    const neiOp& neighbourCalculator,
+    const filterOp& selector
+)
+{
+    label nGroups(0);
+
+    elementInGroup.setSize(neighbourCalculator.size());
+    elementInGroup = -1;
+
+    forAll(neighbourCalculator, elI)
+    {
+        if( elementInGroup[elI] != -1 )
+            continue;
+        if( !selector(elI) )
+            continue;
+
+        elementInGroup[elI] = nGroups;
+        labelListPMG front;
+        front.append(elI);
+
+        while( front.size() )
+        {
+            const label eLabel = front.removeLastElement();
+
+            DynList<label> neighbours;
+            neighbourCalculator(eLabel, neighbours);
+
+            forAll(neighbours, neiI)
+            {
+                const label nei = neighbours[neiI];
+
+                if( (nei < 0) || (elementInGroup[nei] != -1) )
+                    continue;
+
+                if( selector(nei) )
+                {
+                    elementInGroup[nei] = nGroups;
+                    front.append(nei);
+                }
+            }
+        }
+
+        //- increment the number of groups
+        ++nGroups;
+    }
+
+    return nGroups;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+} // End namespace help
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

meshLibrary/utilities/helperFunctions/helperFunctionsFrontalMarking.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    helperFunctionsPar
+
+Description
+    Functions used for exchanging data between processors
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef helperFunctionsFrontalMarking_H
+#define helperFunctionsFrontalMarking_H
+
+#include "labelList.H"
+#include "LongList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace help
+{
+
+//- template implementation of the fron-marking algorithm
+//- labelListType contains indices of elements satisfying the filtering criteria
+//- startingIndex is the first element from where the search is started
+//- neiOp determines which neighbouring elements are available. It requires
+//- an operator neiOp(const label, DynList<label>&) to be defined
+//- filterOp determines which neighbouring elements shall be stored
+//- in the result and used in the front. It requires an operator
+//- bool filterOp(const label) which return true if the element shall be part
+//- of the front
+template<class labelListType, class neiOp, class filterOp>
+void frontalMarking
+(
+    labelListType& result,
+    const label startingIndex,
+    const neiOp& neighbourCalculator,
+    const filterOp& selector
+);
+
+//- templated implementation of functionality for finding discrete groups
+//- of elements depending on the user-defined criteria
+//- the first argument is the result of the operation. Each element is assigned
+//- to a group amd the ones failing the critearia are set to -1
+//- neiOp is the functionality for determining neighbours of the currently
+//- processed front element. This class requires an operator
+//- neiOp(const label, DynList<label>&)
+//- filterOp represents the user-defined criteria and shall have an operator
+//- bool filterOp(const label)
+template<class labelListType, class neiOp, class filterOp>
+label groupMarking
+(
+    labelListType& elementInGroup,
+    const neiOp& neighbourCalculator,
+    const filterOp& selector
+);
+
+} // End namespace help
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+#ifdef NoRepository
+#   include "helperFunctionsFrontalMarking.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

meshLibrary/utilities/helperFunctions/helperFunctionsGeometryQueries.H

@@ -39,6 +39,7 @@ SourceFiles
 #include "plane.H"
 #include "face.H"
 
+#include "triSurf.H"
 #include "triangle.H"
 #include "tetrahedron.H"
 
@@ -47,7 +48,6 @@ SourceFiles
 class boolList;
 class labelHashSet;
 class boundBox;
-class triSurf;
 
 namespace Foam
 {

meshLibrary/utilities/helperFunctions/helperFunctionsPar.C

@@ -32,9 +32,9 @@ Description
 #include "labelPair.H"
 #include "HashSet.H"
 
+# ifdef USE_OMP
 #include <omp.h>
-
-#define DEBUGExchangeMap
+# endif
 
 namespace Foam
 {

meshLibrary/utilities/octrees/meshOctree/meshOctreeAddressing/meshOctreeAddressing.C

@@ -356,7 +356,7 @@ void meshOctreeAddressing::edgeIntersections
             //- check for geometric intersection
             point intersection;
 
-            if(
+            const bool intersectionExists =
                 help::triLineIntersection
                 (
                     octree_.surface(),
@@ -364,8 +364,9 @@ void meshOctreeAddressing::edgeIntersections
                     points[edges[eI].start()],
                     points[edges[eI].end()],
                     intersection
-                )
-            )
+                );
+
+            if( intersectionExists )
             {
                 bool store(true);
                 forAll(intersections, i)

meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctreeCube/surfaceIntersectionsOctreeCubeIntersections.C

@@ -198,13 +198,15 @@ void surfaceIntersectionsOctreeCube::intersectionCandidates
         Info << "Cube " << cubeBox_ << " contains elements "
             << containedElements_ << endl;
         # endif
-        point intersection;
+
         //- add elements within the cube
         for(SLList<label>::const_iterator eIter = containedElements_.begin();
             eIter != containedElements_.end();
             ++eIter
         )
         {
+            point intersection;
+
             if(
                 help::triLineIntersection
                 (
@@ -215,13 +217,15 @@ void surfaceIntersectionsOctreeCube::intersectionCandidates
                     intersection
                 )
             )
-                cp.append(pointIndexHit(true,intersection, eIter()));
+                cp.append(pointIndexHit(true, intersection, eIter()));
         }
 
         //- check subCubes if there are any
         if( subCubesPtr_ )
         {
-            const FixedList<surfaceIntersectionsOctreeCube*, 8>& subCubes_ = *subCubesPtr_;
+            const FixedList<surfaceIntersectionsOctreeCube*, 8>& subCubes_ =
+                *subCubesPtr_;
+
             forAll(subCubes_, cubeI)
                 subCubes_[cubeI]->intersectionCandidates(cp, s, e);
         }

testingInterfaces/testFrontalMarking/Make/files

+
+testFrontalMarking.C
+
+EXE = $(FOAM_USER_APPBIN)/testFrontalMarking

testingInterfaces/testFrontalMarking/Make/options

+EXE_INC = -I$(LIB_SRC)/triSurface/lnInclude -I$(LIB_SRC)/meshTools/lnInclude -I../../meshLibrary/lnInclude
+EXE_LIBS = -lMeshLibrary -ltriSurface -lmeshTools

testingInterfaces/testFrontalMarking/testFrontalMarking.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2005-2007 Franjo Juretic
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    Test of the cartesian mesh
+
+Description
+    - creates an octree and creates cartesian mesh
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "meshOctreeCreator.H"
+#include "meshOctreeAutomaticRefinement.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "polyMeshGen.H"
+#include "cartesianMeshExtractor.H"
+#include "triSurf.H"
+#include "helperFunctions.H"
+#include "findCellsIntersectingSurface.H"
+
+#include <sstream>
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- find octree leaves connected over faces
+class octreeNeighbours
+{
+    const meshOctree& octree_;
+
+public:
+
+    octreeNeighbours(const meshOctree& octree)
+    :
+      octree_(octree)
+    {}
+
+    label size() const
+    {
+        return octree_.numberOfLeaves();
+    }
+
+    void operator()(const label leafI, DynList<label>& neighbours) const
+    {
+        octree_.findNeighboursForLeaf(leafI, neighbours);
+    }
+};
+
+class octreeSelectOperator
+{
+    const meshOctree& octree_;
+
+public:
+
+    octreeSelectOperator(const meshOctree& octree)
+    :
+      octree_(octree)
+    {}
+
+    bool operator()(const label leafI) const
+    {
+        if( octree_.returnLeaf(leafI).cubeType() & meshOctreeCube::DATA )
+            return false;
+
+        return true;
+    }
+};
+
+class meshNeighbourOperator
+{
+    const polyMeshGen& mesh_;
+
+public:
+
+    meshNeighbourOperator(const polyMeshGen& mesh)
+    :
+        mesh_(mesh)
+    {}
+
+    label size() const
+    {
+        return mesh_.cells().size();
+    }
+
+    void operator()(const label cellI, DynList<label>& neighbourCells) const
+    {
+        neighbourCells.clear();
+
+        const labelList& owner = mesh_.owner();
+        const labelList& neighbour = mesh_.neighbour();
+
+        const cell& c = mesh_.cells()[cellI];
+
+        forAll(c, fI)
+        {
+            label nei = owner[c[fI]];
+
+            if( nei == cellI )
+                nei = neighbour[c[fI]];
+
+            if( nei >= 0 )
+                neighbourCells.append(nei);
+        }
+    }
+};
+
+class meshSelectorOperator
+{
+    const VRWGraph& patchesIntersectingCell_;
+
+public:
+
+    meshSelectorOperator(const VRWGraph& patchesIntersectingCell)
+    :
+        patchesIntersectingCell_(patchesIntersectingCell)
+    {}
+
+    bool operator()(const label cellI) const
+    {
+        if( patchesIntersectingCell_[cellI].size() )
+            return false;
+
+        return true;
+    }
+};
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    IOdictionary meshDict
+    (
+        IOobject
+        (
+            "meshDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    fileName surfaceFile = meshDict.lookup("surfaceFile");
+    if( Pstream::parRun() )
+        surfaceFile = ".."/surfaceFile;
+
+    triSurf surf;
+    surf.readSurface(runTime.path()/surfaceFile);
+
+    // construct the octree
+    meshOctree mo(surf);
+    meshOctreeCreator(mo, meshDict).createOctreeBoxes();
+
+    meshOctreeAutomaticRefinement(mo, meshDict, false).automaticRefinement();
+
+    octreeNeighbours on(mo);
+    octreeSelectOperator oso(mo);
+    labelListPMG result;
+    label start(0);
+
+    help::frontalMarking(result, start, on, oso);
+
+    Info << "Number of octree elements meeting the criteria " << result << endl;
+
+    Info<< "Execution time for octree creation = "
+        << runTime.elapsedCpuTime()
+        << " s\n" << endl << endl;
+
+    polyMeshGen pmg(runTime);
+    cartesianMeshExtractor cmg(mo, meshDict, pmg);
+
+    //cmg.decomposeSplitHexes();
+    cmg.createMesh();
+
+    findCellsIntersectingSurface fis(pmg, mo);
+
+    meshNeighbourOperator mnop(pmg);
+    meshSelectorOperator msop(fis.facetsIntersectingCells());
+
+    help::frontalMarking(result, 0, mnop, msop);
+    Info << "Cells in the group " << result.size() << endl;
+
+    const label nGroups = help::groupMarking(result, mnop, msop);
+
+    Info << "Number of groups " << nGroups << endl;
+    labelList groupId(nGroups);
+    forAll(groupId, i)
+        groupId[i] = pmg.addCellSubset("group_"+help::scalarToText(i));
+
+    forAll(result, cellI)
+    {
+        if( result[cellI] < 0 )
+            continue;
+
+        pmg.addCellToSubset(groupId[result[cellI]], cellI);
+    }
+
+    pmg.write();
+
+    Info << "End\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //