View Issue Details

IDProjectCategoryView StatusLast Update
0003611OpenFOAMPatchpublic2021-01-08 16:50
Reporterhandrake0724Assigned Towill 
PrioritynormalSeverityminorReproducibilityhave not tried
Status closedResolutionno change required 
Product Versiondev 
Fixed in Version 
Summary0003611: support centre position information in triCut and tetCut
DescriptionI think triCut and tetCut function is useful to calculate the part of scalar field in a cell or a face.
One thing I am missing in triCut and tetCut is that it does not provide centre information of the partial shape of a cell or a face.
I am working on irregular wave generation with triCut and tetCut function such that wave velocity is evaluated on the centre of wet side of a cell or a face.
for that, I added several function like triCutAndCentre, tetCutand Centre and dependent functions in cut.H and found the functions are working as I expected.
It would be nice if these functions are supported in triCut and tetCut functions.
please review the attached file and let me know your opinion.
TagsNo tags attached.

Activities

handrake0724

2021-01-08 14:53

reporter  

patch.diff (16,020 bytes)
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
index 1e412a7cd..c215db2d3 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
@@ -576,6 +576,15 @@ typename cut::opAddResult<AboveOp, BelowOp>::type triCut
     const BelowOp& belowOp
 );
 
+template<class Point, class AboveOp, class BelowOp>
+Tuple2<typename cut::opAddResult<AboveOp, BelowOp>::type, Point> triCutAndCentre
+(
+    const FixedList<vector, 3>& tri,
+    const FixedList<scalar, 3>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+);
+
 //- As above, but with a plane specifying the location of the cut
 template<class AboveOp, class BelowOp>
 typename cut::opAddResult<AboveOp, BelowOp>::type triCut
@@ -596,6 +605,15 @@ typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
     const BelowOp& belowOp
 );
 
+template<class Point, class AboveOp, class BelowOp>
+Tuple2<typename cut::opAddResult<AboveOp, BelowOp>::type, Point> tetCutAndCentre
+(
+    const FixedList<vector, 4>& tet,
+    const FixedList<scalar, 4>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+);
+
 //- As above, but with a plane specifying the location of the cut
 template<class AboveOp, class BelowOp>
 typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
index 9b7ca4794..263b2a352 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
@@ -286,6 +286,27 @@ inline typename Op::result triCutTri
     return Op(triCutTri(op, f))(triCutTri(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> triCutTriAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 3>& p,
+    const Pair<scalar>& f
+)
+{
+    FixedList<Point, 3> cutFace = triCutTri(p, f);
+    Point centre = Zero;
+    forAll(cutFace, pi)
+    {
+        centre += cutFace[pi];
+    }
+    centre /= 3;
+    return Tuple2<typename Op::result, Point>
+        (
+            Op(triCutTri(op, f))(cutFace),
+            centre
+        );
+}
 
 //- Apply an operation to a quad. Splits the quad into two tris.
 template<class Op, class OpData, class Point>
@@ -317,6 +338,27 @@ inline typename Op::result triCutQuad
     return quadOp<Op>(triCutQuad(op, f), triCutQuad(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> triCutQuadAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 3>& p,
+    const FixedList<scalar, 2>& f
+)
+{
+    FixedList<Point, 4> cutFace = triCutQuad(p, f);
+    Point centre = Zero;
+    forAll(cutFace, pi)
+    {
+        centre += cutFace[pi];
+    }
+    centre /= 4;
+    return Tuple2<typename Op::result, Point>
+        (
+            quadOp<Op>(triCutQuad(op, f), cutFace),
+            centre
+        );
+}
 
 //- Cut a tet from a tet and apply an operation to the result. The cut is made
 //  along the three edges connected to vertex 0, and the cut locations are given
@@ -332,6 +374,28 @@ inline typename Op::result tetCutTet
     return Op(tetCutTet(op, f))(tetCutTet(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutTetAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 3>& f
+)
+{
+    FixedList<Point, 4> cutCell = tetCutTet(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 4;
+
+    return Tuple2<typename Op::result, Point>
+        (
+            Op(tetCutTet(op, f))(cutCell),
+            centre
+        );
+}
 
 //- Apply an operation to a prism. Splits the prism into three tets.
 template<class Op, class OpData, class Point>
@@ -364,6 +428,27 @@ inline typename Op::result tetCutPrism0
     return prismOp<Op>(tetCutPrism0(op, f), tetCutPrism0(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism0AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 3>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism0(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism0(op, f), cutCell),
+            centre
+        );
+}
 
 //- Cut a prism from a tet and apply an operation to the result. The cut is made
 //  along four edges, not edges 01 or 23, and the cut locations are given as
@@ -379,6 +464,27 @@ inline typename Op::result tetCutPrism01
     return prismOp<Op>(tetCutPrism01(op, f), tetCutPrism01(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism01AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 4>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism01(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism01(op, f), cutCell),
+            centre
+        );
+}
 
 //- Cut a prism from a tet and apply an operation to the result. The cuts are
 //  the same as for tetCutPrism01. The result is the side connected to edge 23.
@@ -393,6 +499,28 @@ inline typename Op::result tetCutPrism23
     return prismOp<Op>(tetCutPrism23(op, f), tetCutPrism23(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism23AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 4>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism23(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism23(op, f), cutCell),
+            centre
+        );
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
index 11bff4eba..25d1a1106 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
@@ -97,6 +97,96 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
     }
 }
 
+template<class Point, class AboveOp, class BelowOp>
+Foam::Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+Foam::triCutAndCentre
+(
+    const FixedList<Point, 3>& tri,
+    const FixedList<scalar, 3>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+)
+{
+    // Quick return if all levels are zero or have the same sign
+    if (level[0] == 0 && level[1] == 0 && level[2] == 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(), Zero);
+    }
+    if (level[0] >= 0 && level[1] >= 0 && level[2] >= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp(tri) + belowOp(), Zero);
+    }
+    if (level[0] <= 0 && level[1] <= 0 && level[2] <= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(tri), triPointRef(tri[0], tri[1], tri[2]).centre());
+    }
+
+    // There will be just one edge without a sign change. Find it, and put it
+    // opposite the first vertex. This may change the sign of the tri.
+    FixedList<label, 3> indices({0, 1, 2});
+    label i;
+    for (i = 0; i < 3; ++ i)
+    {
+        if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
+        {
+            Swap(indices[0], indices[i]);
+            break;
+        }
+    }
+    if (i == 3)
+    {
+        FatalErrorInFunction
+            << "The number of tri vertices above the level set should always "
+            << "be 1" << exit(FatalError);
+    }
+
+    // Correct the sign
+    if (indices[0] != 0)
+    {
+        Swap(indices[1], indices[2]);
+    }
+
+    // Permute the data
+    const FixedList<Point, 3> p = triReorder(tri, indices);
+    const FixedList<scalar, 3> l = triReorder(level, indices);
+    AboveOp a = triReorder(aboveOp, indices);
+    BelowOp b = triReorder(belowOp, indices);
+
+    // Slice off one corner to form a tri and a quad
+    Pair<scalar> f;
+    for (label i = 0; i < 2; ++ i)
+    {
+        f[i] = l[0]/(l[0] - l[i+1]);
+    }
+    if (l[0] > 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (
+            triCutTriAndCentre(a, p, f).first() + triCutQuadAndCentre(b, p, f).first(),
+            (
+                triCutTriAndCentre(a, p, f).first()*triCutTriAndCentre(a, p, f).second()
+                + triCutQuadAndCentre(b, p, f).first()*triCutQuadAndCentre(b, p, f).second()
+            )
+          / (triCutTriAndCentre(a, p, f).first() + triCutQuadAndCentre(b, p, f).first())
+        );
+    }
+    else
+    {
+
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (
+            triCutQuadAndCentre(a, p, f).first() + triCutTriAndCentre(b, p, f).first(),
+            (
+                triCutQuadAndCentre(a, p, f).first()*triCutQuadAndCentre(a, p, f).second()
+                + triCutTriAndCentre(b, p, f).first()*triCutTriAndCentre(b, p, f).second()
+            )
+          / (triCutQuadAndCentre(a, p, f).first() + triCutTriAndCentre(b, p, f).first())
+        );
+    }
+}
 
 template<class AboveOp, class BelowOp>
 typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
@@ -246,6 +336,169 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
     return aboveOp() + belowOp();
 }
 
+template<class Point, class AboveOp, class BelowOp>
+Foam::Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+Foam::tetCutAndCentre
+(
+    const FixedList<Point, 4>& tet,
+    const FixedList<scalar, 4>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+)
+{
+    // Get the min and max over all four vertices and quick return if
+    // all levels are zero or have the same sign
+    scalar levelMin = vGreat, levelMax = - vGreat;
+    for (label i = 0; i < 4; ++ i)
+    {
+        levelMin = min(levelMin, level[i]);
+        levelMax = max(levelMax, level[i]);
+    }
+    if (levelMin == 0 && levelMax == 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(), Zero);
+    }
+    if (levelMin >= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp(tet) + belowOp(), Zero);
+    }
+    if (levelMax <= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(tet), tetPointRef(tet[0], tet[1], tet[2], tet[3]).centre());
+    }
+
+    // Partition the level so that positive values are at the start. This is
+    // like a single iteration of quick-sort, except that the pivot is a hard-
+    // coded zero, rather than an element of the array. This can change the sign
+    // of the tet.
+    FixedList<label, 4> indices({0, 1, 2, 3});
+    bool signChange = false;
+    label i = 0, j = 3;
+    while (true)
+    {
+        while (i < j && level[indices[i]] > 0)
+        {
+            i ++;
+        }
+        while (j > i && level[indices[j]] <= 0)
+        {
+            j --;
+        }
+        if (i == j)
+        {
+            break;
+        }
+        Swap(indices[i], indices[j]);
+        signChange = !signChange;
+    }
+
+    // The number of vertices above the slice
+    label n = i;
+
+    // If there are more positives than negatives then reverse the order so that
+    // the negatives are at the start
+    if (n > 2)
+    {
+        n = 4 - n;
+        for (label i = 0; i < 2; ++ i)
+        {
+            Swap(indices[i], indices[3-i]);
+        }
+    }
+
+    // Correct the sign
+    if (signChange)
+    {
+        Swap(indices[2], indices[3]);
+    }
+
+    // Permute the data
+    const FixedList<Point, 4> p = tetReorder(tet, indices);
+    const FixedList<scalar, 4> l = tetReorder(level, indices);
+    AboveOp a = tetReorder(aboveOp, indices);
+    BelowOp b = tetReorder(belowOp, indices);
+
+    // Calculate the integrals above and below the level set
+    if (n == 1)
+    {
+        // Slice off one corner to form a tet and a prism
+        FixedList<scalar, 3> f;
+        for (label i = 0; i < 3; ++ i)
+        {
+            f[i] = l[0]/(l[0] - l[i+1]);
+        }
+        if (l[0] > 0)
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutTetAndCentre(a, p, f).first() + tetCutPrism0AndCentre(b, p, f).first(),
+                (
+                    tetCutTetAndCentre(a, p, f).first()*tetCutTetAndCentre(a, p, f).second()
+                    + tetCutPrism0AndCentre(b, p, f).first()*tetCutPrism0AndCentre(b, p, f).second()
+                )
+              / (tetCutTetAndCentre(a, p, f).first() + tetCutPrism0AndCentre(b, p, f).first())
+            );
+        }
+        else
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism0AndCentre(a, p, f).first() + tetCutTetAndCentre(b, p, f).first(),
+                (
+                    tetCutPrism0AndCentre(a, p, f).first()*tetCutPrism0AndCentre(a, p, f).second()
+                    + tetCutTetAndCentre(b, p, f).first()*tetCutTetAndCentre(b, p, f).second()
+                )
+              / (tetCutPrism0AndCentre(a, p, f).first() + tetCutTetAndCentre(b, p, f).first())
+            );
+        }
+    }
+    else if (n == 2)
+    {
+        // Slice off two corners to form two prisms
+        FixedList<scalar, 4> f;
+        for (label i = 0; i < 2; ++ i)
+        {
+            for (label j = 0; j < 2; ++ j)
+            {
+                f[2*i+j] = l[i]/(l[i] - l[j+2]);
+            }
+        }
+        if (l[0] > 0)
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism01AndCentre(a, p, f).first() + tetCutPrism23AndCentre(b, p, f).first(),
+                (
+                    tetCutPrism01AndCentre(a, p, f).first()*tetCutPrism01AndCentre(a, p, f).second()
+                    + tetCutPrism23AndCentre(b, p, f).first()*tetCutPrism23AndCentre(b, p, f).second()
+                )
+              / (tetCutPrism01AndCentre(a, p, f).first() + tetCutPrism23AndCentre(b, p, f).first())
+            );
+        }
+        else
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism23AndCentre(a, p, f).first() + tetCutPrism01AndCentre(b, p, f).first(),
+                (
+                    tetCutPrism23AndCentre(a, p, f).first()*tetCutPrism23AndCentre(a, p, f).second()
+                    + tetCutPrism01AndCentre(b, p, f).first()*tetCutPrism01AndCentre(b, p, f).second()
+                )
+              / (tetCutPrism23AndCentre(a, p, f).first() + tetCutPrism01AndCentre(b, p, f).first())
+            );
+        }
+    }
+
+    FatalErrorInFunction
+        << "The number of tet vertices above the level set should always be "
+        << "either 1 or 2" << exit(FatalError);
+
+    return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (aboveOp() + belowOp(), Zero);
+}
 
 template<class AboveOp, class BelowOp>
 typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
patch.diff (16,020 bytes)

henry

2021-01-08 15:03

manager   ~0011809

Are you now in a position to fund maintenance to cover the cost of checking, testing, releasing and maintaining your code? See

https://openfoam.org/news/funding-2021/

If not do you know of any companies that are interested in supporting this work?

handrake0724

2021-01-08 16:19

reporter   ~0011810

KRISO (Korean Research Institute of Ship and Ocean Engineering) has a plan to develop numerical wave tank (NWT) with OpenFOAM.
There is possibility they start the project based on my code, which is in discussion with KRISO and me.
In that sense, KRISO might be interested in supporting the code but not sure.
Anyway, it is fine to close this issue if you want.

will

2021-01-08 16:28

manager   ~0011811

Last edited: 2021-01-08 16:28

View 2 revisions

It's already supported. The whole point of the templated cut operations is that you don't have to replicate all the cutting logic every time you want to calculate a new property. You can calculate the area or volume of the cut, and you can calculate the area or volume weighted integral of the vertex positions, and if you divide those values then that will give you the centroid. E.g.,

FixedList<point, 4> tet = ...;
FixedList<scalar, 4> level = ...;
const scalar volumeAbove = tetCut(tet, level, cut::volumeOp(), cut::noOp());
const point volumeCentreAbove = tetCut(tet, level, cut::volumeIntegrateOp<point>(tet), cut::noOp());
const point centreAbove = volumeCentreAbove/volumeAbove;

handrake0724

2021-01-08 16:32

reporter   ~0011812

Thank you for enlighten me.

Issue History

Date Modified Username Field Change
2021-01-08 14:53 handrake0724 New Issue
2021-01-08 14:53 handrake0724 File Added: patch.diff
2021-01-08 15:03 henry Note Added: 0011809
2021-01-08 16:19 handrake0724 Note Added: 0011810
2021-01-08 16:28 will Note Added: 0011811
2021-01-08 16:28 will Note Edited: 0011811 View Revisions
2021-01-08 16:32 handrake0724 Note Added: 0011812
2021-01-08 16:50 will Assigned To => will
2021-01-08 16:50 will Status new => closed
2021-01-08 16:50 will Resolution open => no change required