View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0003611 | OpenFOAM | Patch | public | 2021-01-08 14:53 | 2021-01-08 16:50 |
Reporter | handrake0724 | Assigned To | will | ||
Priority | normal | Severity | minor | Reproducibility | have not tried |
Status | closed | Resolution | no change required | ||
Product Version | dev | ||||
Summary | 0003611: support centre position information in triCut and tetCut | ||||
Description | I 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. | ||||
Tags | No tags attached. | ||||
|
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 |
|
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? |
|
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. |
|
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; |
|
Thank you for enlighten me. |
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 | |
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 |