From 5edba1d109d377a5d675e836211890550aca268f Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Tue, 10 Jan 2012 15:38:07 +0000
Subject: [PATCH] ENH: Adding mapRegion option to meshToMesh and volume
weighted interpolation method
---
src/sampling/Make/files | 4 +
.../meshToMesh/calculateMeshToMeshWeights.C | 123 ++
.../meshToMesh/meshToMesh.C | 10 +-
.../meshToMesh/meshToMesh.H | 29 +-
.../meshToMesh/meshToMeshInterpolate.C | 68 +-
.../tetOverlapVolume/tetOverlapVolume.C | 795 +++++++++++++
.../tetOverlapVolume/tetOverlapVolume.H | 218 ++++
.../tetOverlapVolume/tetPoints.H | 114 ++
.../tetOverlapVolume/tetrahedron.C | 341 ++++++
.../tetOverlapVolume/tetrahedron.H | 315 +++++
.../tetOverlapVolume/tetrahedronI.H | 1012 +++++++++++++++++
11 files changed, 3020 insertions(+), 9 deletions(-)
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H
create mode 100644 src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 0c89693a8ea..17099071a57 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -69,4 +69,8 @@ $(meshToMesh)/meshToMesh.C
$(meshToMesh)/calculateMeshToMeshAddressing.C
$(meshToMesh)/calculateMeshToMeshWeights.C
+tetOverlapVolume = meshToMeshInterpolation/tetOverlapVolume
+$(tetOverlapVolume)/tetOverlapVolume.C
+
+
LIB = $(FOAM_LIBBIN)/libsampling
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C
index 532fee3f85c..c3db6586ad5 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C
@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "meshToMesh.H"
+#include "tetOverlapVolume.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -99,6 +100,108 @@ void Foam::meshToMesh::calculateInverseDistanceWeights() const
}
+void Foam::meshToMesh::calculateInverseVolumeWeights() const
+{
+ if (debug)
+ {
+ Info<< "meshToMesh::calculateInverseVolumeWeights() : "
+ << "calculating inverse volume weighting factors" << endl;
+ }
+
+ if (inverseVolumeWeightsPtr_)
+ {
+ FatalErrorIn("meshToMesh::calculateInverseVolumeWeights()")
+ << "weighting factors already calculated"
+ << exit(FatalError);
+ }
+
+ inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
+ scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
+
+ labelListList& cellToCell = *cellToCellAddressingPtr_;
+
+ tetOverlapVolume overlapEngine;
+
+ forAll(cellToCell, celli)
+ {
+ const labelList& overlapCells = cellToCell[celli];
+
+ if (overlapCells.size() > 0)
+ {
+ invVolCoeffs[celli].setSize(overlapCells.size());
+ scalar v(0);
+ forAll (overlapCells, j)
+ {
+ label cellFrom = overlapCells[j];
+ treeBoundBox bbFromMesh
+ (
+ pointField
+ (
+ fromMesh_.points(),
+ fromMesh_.cellPoints()[cellFrom]
+ )
+ );
+
+ v = overlapEngine.cellCellOverlapVolumeMinDecomp
+ (
+ toMesh_,
+ celli,
+
+ fromMesh_,
+ cellFrom,
+ bbFromMesh
+ );
+ invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
+ }
+ if (celli == 2)
+ {
+ Info << "cellToCell :" << cellToCell[celli] << endl;
+ Info << "invVolCoeffs :" << invVolCoeffs[celli] << endl;
+ }
+ }
+ }
+}
+
+
+void Foam::meshToMesh::calculateCellToCellAddressing() const
+{
+ if (debug)
+ {
+ Info<< "meshToMesh::calculateCellToCellAddressing() : "
+ << "calculating cell to cell addressing" << endl;
+ }
+
+ if (cellToCellAddressingPtr_)
+ {
+ FatalErrorIn("meshToMesh::calculateCellToCellAddressing()")
+ << "addressing already calculated"
+ << exit(FatalError);
+ }
+
+ tetOverlapVolume overlapEngine;
+
+ cellToCellAddressingPtr_ = new labelListList(toMesh_.nCells());
+ labelListList& cellToCell = *cellToCellAddressingPtr_;
+
+
+ forAll(cellToCell, iTo)
+ {
+ const labelList overLapCells =
+ overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo);
+ if (overLapCells.size() > 0)
+ {
+ //Info << "To " << iTo << endl;
+ //Info << "cellToCell " << overLapCells << endl;
+
+ cellToCell[iTo].setSize(overLapCells.size());
+ forAll(overLapCells, j)
+ {
+ cellToCell[iTo][j] = overLapCells[j];
+ }
+ }
+ }
+}
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
const Foam::scalarListList& Foam::meshToMesh::inverseDistanceWeights() const
@@ -112,4 +215,24 @@ const Foam::scalarListList& Foam::meshToMesh::inverseDistanceWeights() const
}
+const Foam::scalarListList& Foam::meshToMesh::inverseVolumeWeights() const
+{
+ if (!inverseVolumeWeightsPtr_)
+ {
+ calculateInverseVolumeWeights();
+ }
+
+ return *inverseVolumeWeightsPtr_;
+}
+
+const Foam::labelListList& Foam::meshToMesh::cellToCellAddressing() const
+{
+ if (!cellToCellAddressingPtr_)
+ {
+ calculateCellToCellAddressing();
+ }
+
+ return *cellToCellAddressingPtr_;
+}
+
// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
index c5f95879501..fa495c40ce5 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
@@ -48,7 +48,9 @@ Foam::meshToMesh::meshToMesh
patchMap_(patchMap),
cellAddressing_(toMesh_.nCells()),
boundaryAddressing_(toMesh_.boundaryMesh().size()),
- inverseDistanceWeightsPtr_(NULL)
+ inverseDistanceWeightsPtr_(NULL),
+ inverseVolumeWeightsPtr_(NULL),
+ cellToCellAddressingPtr_(NULL)
{
forAll(fromMesh_.boundaryMesh(), patchi)
{
@@ -118,7 +120,9 @@ Foam::meshToMesh::meshToMesh
toMesh_(meshTo),
cellAddressing_(toMesh_.nCells()),
boundaryAddressing_(toMesh_.boundaryMesh().size()),
- inverseDistanceWeightsPtr_(NULL)
+ inverseDistanceWeightsPtr_(NULL),
+ inverseVolumeWeightsPtr_(NULL),
+ cellToCellAddressingPtr_(NULL)
{
// check whether both meshes have got the same number
// of boundary patches
@@ -198,6 +202,8 @@ Foam::meshToMesh::meshToMesh
Foam::meshToMesh::~meshToMesh()
{
deleteDemandDrivenData(inverseDistanceWeightsPtr_);
+ deleteDemandDrivenData(inverseVolumeWeightsPtr_);
+ deleteDemandDrivenData(cellToCellAddressingPtr_);
}
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 83d1ee902ac..e60871be259 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -88,6 +88,12 @@ class meshToMesh
//- Inverse-distance interpolation weights
mutable scalarListList* inverseDistanceWeightsPtr_;
+ //- Inverse-volume interpolation weights
+ mutable scalarListList* inverseVolumeWeightsPtr_;
+
+ //- Cell to cell overlap addressing
+ mutable labelListList* cellToCellAddressingPtr_;
+
// Private Member Functions
@@ -104,8 +110,16 @@ class meshToMesh
void calculateInverseDistanceWeights() const;
+ void calculateInverseVolumeWeights() const;
+
+ void calculateCellToCellAddressing() const;
+
const scalarListList& inverseDistanceWeights() const;
+ const scalarListList& inverseVolumeWeights() const;
+
+ const labelListList& cellToCellAddressing() const;
+
// Private static data members
@@ -124,7 +138,8 @@ public:
{
MAP,
INTERPOLATE,
- CELL_POINT_INTERPOLATE
+ CELL_POINT_INTERPOLATE,
+ CELL_VOLUME_WEIGHT
};
@@ -239,6 +254,18 @@ public:
const CombineOp& cop
) const;
+ //- Interpolate field using inverse-volume weights
+ template<class Type, class CombineOp>
+ void interpolateField
+ (
+ Field<Type>&,
+ const GeometricField<Type, fvPatchField, volMesh>&,
+ const labelListList& adr,
+ const scalarListList& weights,
+ const CombineOp& cop
+ ) const;
+
+
//- Interpolate field using cell-point interpolation
template<class Type, class CombineOp>
void interpolateField
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
index 7ed543ed3cd..372194e3225 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
@@ -54,6 +54,33 @@ void Foam::meshToMesh::mapField
}
+template<class Type, class CombineOp>
+void Foam::meshToMesh::interpolateField
+(
+ Field<Type>& toF,
+ const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+ const labelListList& adr,
+ const scalarListList& weights,
+ const CombineOp& cop
+) const
+{
+ // Inverse volume weighted interpolation
+ forAll(toF, celli)
+ {
+ const labelList& overlapCells = adr[celli];
+ const scalarList& w = weights[celli];
+
+ Type f = pTraits<Type>::zero;
+ forAll(overlapCells, i)
+ {
+ label fromCelli = overlapCells[i];
+ f += fromVf[fromCelli]*w[i];
+ cop(toF[celli], f);
+ }
+ }
+}
+
+
template<class Type, class CombineOp>
void Foam::meshToMesh::interpolateField
(
@@ -162,6 +189,7 @@ void Foam::meshToMesh::interpolateInternalField
break;
case INTERPOLATE:
+ {
interpolateField
(
toF,
@@ -170,9 +198,10 @@ void Foam::meshToMesh::interpolateInternalField
inverseDistanceWeights(),
cop
);
- break;
-
+ break;
+ }
case CELL_POINT_INTERPOLATE:
+ {
interpolateField
(
toF,
@@ -181,8 +210,24 @@ void Foam::meshToMesh::interpolateInternalField
toMesh_.cellCentres(),
cop
);
- break;
+ break;
+ }
+ case CELL_VOLUME_WEIGHT:
+ {
+ const labelListList& cellToCell = cellToCellAddressing();
+ const scalarListList& invVolWeights = inverseVolumeWeights();
+
+ interpolateField
+ (
+ toF,
+ fromVf,
+ cellToCell,
+ invVolWeights,
+ cop
+ );
+ break;
+ }
default:
FatalErrorIn
(
@@ -229,6 +274,7 @@ void Foam::meshToMesh::interpolate
switch(ord)
{
case MAP:
+ {
mapField
(
toVf.boundaryField()[patchi],
@@ -236,9 +282,11 @@ void Foam::meshToMesh::interpolate
boundaryAddressing_[patchi],
cop
);
- break;
+ break;
+ }
case INTERPOLATE:
+ {
interpolateField
(
toVf.boundaryField()[patchi],
@@ -247,9 +295,11 @@ void Foam::meshToMesh::interpolate
toPatch.Cf(),
cop
);
- break;
+ break;
+ }
case CELL_POINT_INTERPOLATE:
+ {
interpolateField
(
toVf.boundaryField()[patchi],
@@ -258,7 +308,13 @@ void Foam::meshToMesh::interpolate
toPatch.Cf(),
cop
);
- break;
+ break;
+ }
+ case CELL_VOLUME_WEIGHT:
+ {
+ // Do nothing
+ break;
+ }
default:
FatalErrorIn
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
new file mode 100644
index 00000000000..d36730a0eea
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
@@ -0,0 +1,795 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetOverlapVolume.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::tetOverlapVolume, 0);
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::tetOverlapVolume::tetOverlapVolume()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::tetOverlapVolume::~tetOverlapVolume()
+{}
+
+
+// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+Foam::point Foam::tetOverlapVolume::planeIntersection
+(
+ const FixedList<scalar, 4>& d,
+ const tetPoints& t,
+ const label negI,
+ const label posI
+)
+{
+ return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI]+d[posI]);
+}
+
+
+template <class TetOp>
+inline void Foam::tetOverlapVolume::decomposePrism
+(
+ const FixedList<point, 6>& points,
+ TetOp& op
+)
+{
+ op(tetPoints(points[1], points[3], points[2], points[0]));
+ op(tetPoints(points[1], points[2], points[3], points[4]));
+ op(tetPoints(points[4], points[2], points[3], points[5]));
+}
+
+
+template <class AboveTetOp, class BelowTetOp>
+inline void Foam::tetOverlapVolume::tetSliceWithPlane
+(
+ const tetPoints& tet,
+ const plane& pl,
+
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+)
+{
+ // Distance to plane
+ FixedList<scalar, 4> d;
+ label nPos = 0;
+ forAll(tet, i)
+ {
+ d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
+ if (d[i] > 0)
+ {
+ nPos++;
+ }
+ }
+
+ if (nPos == 4)
+ {
+ aboveOp(tet);
+ }
+ else if (nPos == 3)
+ {
+ // Sliced into below tet and above prism. Prism gets split into
+ // two tets.
+
+ // Find the below tet
+ label i0 = -1;
+ forAll(d, i)
+ {
+ if (d[i] <= 0)
+ {
+ i0 = i;
+ break;
+ }
+ }
+
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+ label i3 = d.fcIndex(i2);
+
+ point p01 = planeIntersection(d, tet, i0, i1);
+ point p02 = planeIntersection(d, tet, i0, i2);
+ point p03 = planeIntersection(d, tet, i0, i3);
+
+ // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
+ // ,, 1 : ,, inwards pointing triad
+ // ,, 2 : ,, outwards pointing triad
+ // ,, 3 : ,, inwards pointing triad
+
+ //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
+
+ if (i0 == 0 || i0 == 2)
+ {
+ tetPoints t(tet[i0], p01, p02, p03);
+ //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 3, belowTet i0==0 or 2");
+ belowOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i1];
+ p[1] = tet[i3];
+ p[2] = tet[i2];
+ p[3] = p01;
+ p[4] = p03;
+ p[5] = p02;
+ //Pout<< " aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ else
+ {
+ tetPoints t(p01, p02, p03, tet[i0]);
+ //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 3, belowTet i0==1 or 3");
+ belowOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i3];
+ p[1] = tet[i1];
+ p[2] = tet[i2];
+ p[3] = p03;
+ p[4] = p01;
+ p[5] = p02;
+ //Pout<< " aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ }
+ else if (nPos == 2)
+ {
+ // Tet cut into two prisms. Determine the positive one.
+ label pos0 = -1;
+ label pos1 = -1;
+ label neg0 = -1;
+ label neg1 = -1;
+ forAll(d, i)
+ {
+ if (d[i] > 0)
+ {
+ if (pos0 == -1)
+ {
+ pos0 = i;
+ }
+ else
+ {
+ pos1 = i;
+ }
+ }
+ else
+ {
+ if (neg0 == -1)
+ {
+ neg0 = i;
+ }
+ else
+ {
+ neg1 = i;
+ }
+ }
+ }
+
+ //Pout<< "Split 2pos tet " << tet << " d:" << d
+ // << " around pos0:" << pos0 << " pos1:" << pos1
+ // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
+
+ const edge posEdge(pos0, pos1);
+
+ if (posEdge == edge(0, 1))
+ {
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[0];
+ p[1] = p02;
+ p[2] = p03;
+ p[3] = tet[1];
+ p[4] = p12;
+ p[5] = p13;
+ //Pout<< " 01 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p02;
+ p[2] = p12;
+ p[3] = tet[3];
+ p[4] = p03;
+ p[5] = p13;
+ //Pout<< " 01 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(1, 2))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p01;
+ p[2] = p13;
+ p[3] = tet[2];
+ p[4] = p02;
+ p[5] = p23;
+ //Pout<< " 12 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[3];
+ p[1] = p23;
+ p[2] = p13;
+ p[3] = tet[0];
+ p[4] = p02;
+ p[5] = p01;
+ //Pout<< " 12 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(2, 0))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p12;
+ p[2] = p23;
+ p[3] = tet[0];
+ p[4] = p01;
+ p[5] = p03;
+ //Pout<< " 20 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p12;
+ p[2] = p01;
+ p[3] = tet[3];
+ p[4] = p23;
+ p[5] = p03;
+ //Pout<< " 20 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(0, 3))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[3];
+ p[1] = p23;
+ p[2] = p13;
+ p[3] = tet[0];
+ p[4] = p02;
+ p[5] = p01;
+ //Pout<< " 03 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p23;
+ p[2] = p02;
+ p[3] = tet[1];
+ p[4] = p13;
+ p[5] = p01;
+ //Pout<< " 03 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(1, 3))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p12;
+ p[2] = p01;
+ p[3] = tet[3];
+ p[4] = p23;
+ p[5] = p03;
+ //Pout<< " 13 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p12;
+ p[2] = p23;
+ p[3] = tet[0];
+ p[4] = p01;
+ p[5] = p03;
+ //Pout<< " 13 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(2, 3))
+ {
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p02;
+ p[2] = p12;
+ p[3] = tet[3];
+ p[4] = p03;
+ p[5] = p13;
+ //Pout<< " 23 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[0];
+ p[1] = p02;
+ p[2] = p03;
+ p[3] = tet[1];
+ p[4] = p12;
+ p[5] = p13;
+ //Pout<< " 23 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else
+ {
+ FatalErrorIn("tetSliceWithPlane(..)") << "Missed edge:" << posEdge
+ << abort(FatalError);
+ }
+ }
+ else if (nPos == 1)
+ {
+ // Find the positive tet
+ label i0 = -1;
+ forAll(d, i)
+ {
+ if (d[i] > 0)
+ {
+ i0 = i;
+ break;
+ }
+ }
+
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+ label i3 = d.fcIndex(i2);
+
+ point p01 = planeIntersection(d, tet, i0, i1);
+ point p02 = planeIntersection(d, tet, i0, i2);
+ point p03 = planeIntersection(d, tet, i0, i3);
+
+ //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
+
+ if (i0 == 0 || i0 == 2)
+ {
+ tetPoints t(tet[i0], p01, p02, p03);
+ //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
+ aboveOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i1];
+ p[1] = tet[i3];
+ p[2] = tet[i2];
+ p[3] = p01;
+ p[4] = p03;
+ p[5] = p02;
+ //Pout<< " belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ else
+ {
+ tetPoints t(p01, p02, p03, tet[i0]);
+ //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
+ aboveOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i3];
+ p[1] = tet[i1];
+ p[2] = tet[i2];
+ p[3] = p03;
+ p[4] = p01;
+ p[5] = p02;
+ //Pout<< " belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else // nPos == 0
+ {
+ belowOp(tet);
+ }
+}
+
+
+void Foam::tetOverlapVolume::tetTetOverlap
+(
+ const tetPoints& tetA,
+ const tetPoints& tetB,
+ FixedList<tetPoints, 200>& insideTets,
+ label& nInside,
+ FixedList<tetPoints, 200>& outsideTets,
+ label& nOutside
+)
+{
+ // Work storage
+ FixedList<tetPoints, 200> cutInsideTets;
+ label nCutInside = 0;
+
+ storeTetOp inside(insideTets, nInside);
+ storeTetOp cutInside(cutInsideTets, nCutInside);
+ dummyTetOp outside;
+
+
+
+ // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
+ // in aboveTets are inside tetB.
+
+ {
+ // face0
+ plane pl0(tetB[1], tetB[3], tetB[2]);
+
+ // Cut and insert subtets into cutInsideTets (either by getting
+ // an index from freeSlots or by appending to insideTets) or
+ // insert into outsideTets
+ tetSliceWithPlane(tetA, pl0, cutInside, outside);
+ }
+
+ if (nCutInside == 0)
+ {
+ nInside = nCutInside;
+ return;
+ }
+
+ {
+ // face1
+ plane pl1(tetB[0], tetB[2], tetB[3]);
+
+ nInside = 0;
+
+ for (label i = 0; i < nCutInside; i++)
+ {
+ tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside);
+ }
+
+ if (nInside == 0)
+ {
+ return;
+ }
+ }
+
+ {
+ // face2
+ plane pl2(tetB[0], tetB[3], tetB[1]);
+
+ nCutInside = 0;
+
+ for (label i = 0; i < nInside; i++)
+ {
+ tetSliceWithPlane(insideTets[i], pl2, cutInside, outside);
+ }
+
+ if (nCutInside == 0)
+ {
+ nInside = nCutInside;
+ return;
+ }
+ }
+
+ {
+ // face3
+ plane pl3(tetB[0], tetB[1], tetB[2]);
+
+ nInside = 0;
+
+ for (label i = 0; i < nCutInside; i++)
+ {
+ tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside);
+ }
+ }
+}
+
+
+inline Foam::scalar
+Foam::tetOverlapVolume::tetTetOverlapVol
+(
+ const tetPoints& tetA,
+ const tetPoints& tetB
+)
+{
+ FixedList<tetPoints, 200> insideTets;
+ label nInside = 0;
+ FixedList<tetPoints, 200> cutInsideTets;
+ label nCutInside = 0;
+
+ storeTetOp inside(insideTets, nInside);
+ storeTetOp cutInside(cutInsideTets, nCutInside);
+ sumTetVolOp volInside;
+ dummyTetOp outside;
+
+ // face0
+ plane pl0(tetB[1], tetB[3], tetB[2]);
+ tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+ if (nCutInside == 0)
+ {
+ return 0.0;
+ }
+
+ // face1
+ plane pl1(tetB[0], tetB[2], tetB[3]);
+ nInside = 0;
+ for (label i = 0; i < nCutInside; i++)
+ {
+ const tetPointRef t = cutInsideTets[i].tet();
+ t.sliceWithPlane(pl1, inside, outside);
+ }
+ if (nInside == 0)
+ {
+ return 0.0;
+ }
+
+ // face2
+ plane pl2(tetB[0], tetB[3], tetB[1]);
+ nCutInside = 0;
+ for (label i = 0; i < nInside; i++)
+ {
+ const tetPointRef t = insideTets[i].tet();
+ t.sliceWithPlane(pl2, cutInside, outside);
+ }
+ if (nCutInside == 0)
+ {
+ return 0.0;
+ }
+
+ // face3
+ plane pl3(tetB[0], tetB[1], tetB[2]);
+ for (label i = 0; i < nCutInside; i++)
+ {
+ const tetPointRef t = cutInsideTets[i].tet();
+ t.sliceWithPlane(pl3, volInside, outside);
+ }
+
+ return volInside.vol_;
+}
+
+
+inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
+(
+ const pointField& points,
+ const face& f,
+ const point& fc
+)
+{
+ treeBoundBox bb(fc, fc);
+ forAll(f, fp)
+ {
+ const point& pt = points[f[fp]];
+ bb.min() = min(bb.min(), pt);
+ bb.max() = max(bb.max(), pt);
+ }
+ return bb;
+}
+
+
+// * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
+(
+ const primitiveMesh& meshA,
+ const label cellAI,
+
+ const primitiveMesh& meshB,
+ const label cellBI,
+ const treeBoundBox& cellBbB
+)
+{
+ const cell& cFacesA = meshA.cells()[cellAI];
+ const point& ccA = meshA.cellCentres()[cellAI];
+
+ const cell& cFacesB = meshB.cells()[cellBI];
+ const point& ccB = meshB.cellCentres()[cellBI];
+
+ scalar vol = 0.0;
+
+ forAll(cFacesA, cFA)
+ {
+ label faceAI = cFacesA[cFA];
+
+ const face& fA = meshA.faces()[faceAI];
+ const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+ if (!pyrA.overlaps(cellBbB))
+ {
+ continue;
+ }
+
+ bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+ label tetBasePtAI = 0;
+
+ const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+ for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+ {
+ label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+ label otherFacePtAI = fA.fcIndex(facePtAI);
+
+ label pt0I = -1;
+ label pt1I = -1;
+
+ if (ownA)
+ {
+ pt0I = fA[facePtAI];
+ pt1I = fA[otherFacePtAI];
+ }
+ else
+ {
+ pt0I = fA[otherFacePtAI];
+ pt1I = fA[facePtAI];
+ }
+
+ const tetPoints tetA
+ (
+ ccA,
+ tetBasePtA,
+ meshA.points()[pt0I],
+ meshA.points()[pt1I]
+ );
+ const treeBoundBox tetABb(tetA.bounds());
+
+
+ // Loop over tets of cellB
+ forAll(cFacesB, cFB)
+ {
+ label faceBI = cFacesB[cFB];
+
+ const face& fB = meshB.faces()[faceBI];
+ const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+ if (!pyrB.overlaps(pyrA))
+ {
+ continue;
+ }
+
+ bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+ label tetBasePtBI = 0;
+
+ const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+ for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+ {
+ label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+ label otherFacePtBI = fB.fcIndex(facePtBI);
+
+ label pt0I = -1;
+ label pt1I = -1;
+
+ if (ownB)
+ {
+ pt0I = fB[facePtBI];
+ pt1I = fB[otherFacePtBI];
+ }
+ else
+ {
+ pt0I = fB[otherFacePtBI];
+ pt1I = fB[facePtBI];
+ }
+
+ const tetPoints tetB
+ (
+ ccB,
+ tetBasePtB,
+ meshB.points()[pt0I],
+ meshB.points()[pt1I]
+ );
+ if (!tetB.bounds().overlaps(tetABb))
+ {
+ continue;
+ }
+
+ vol += tetTetOverlapVol(tetA, tetB);
+ }
+ }
+ }
+ }
+
+ return vol;
+}
+
+
+Foam::labelList Foam::tetOverlapVolume::overlappingCells
+(
+ const fvMesh& fromMesh,
+ const fvMesh& toMesh,
+ const label iTo
+) const
+{
+ const indexedOctree<treeDataCell>& treeA = fromMesh.cellTree();
+
+ treeBoundBox bbB
+ (
+ pointField(toMesh.points(), toMesh.cellPoints()[iTo])
+ );
+
+ return treeA.findBox(bbB);
+}
+
+
+/*
+
+ forAll(cellsA, i)
+ {
+ label cellAI = cellsA[i];
+ treeBoundBox bbA
+ (
+ pointField(meshA.points(), meshA.cellPoints()[cellAI])
+ );
+
+ scalar v = cellCellOverlapVolumeMinDecomp
+ (
+ meshA,
+ cellAI,
+ bbA,
+ meshB,
+ cellBI,
+ bbB
+ );
+
+ overlapVol += v;
+ nOverlapTests++;
+ }
+
+*/
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
new file mode 100644
index 00000000000..d0891c14d14
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+
+Class
+ Foam::tetOverlapVolume
+
+Description
+ Calculates overlap volume of two tets.
+
+SourceFiles
+ tetOverlapVolume.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetOverlapVolume_H
+#define tetOverlapVolume_H
+
+#include "tetrahedron.H"
+#include "fvMesh.H"
+#include "plane.H"
+#include "tetPointRef.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "tetPoints.H"
+#include "tetCell.H"
+#include "EdgeMap.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class tetOverlapVolume Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetOverlapVolume
+{
+ // Private member functions
+
+ //- Plane intersection
+ inline point planeIntersection
+ (
+ const FixedList<scalar, 4>& d,
+ const tetPoints& t,
+ const label negI,
+ const label posI
+ );
+
+
+ //- Decompose prism
+ template <class TetOp> inline void decomposePrism
+ (
+ const FixedList<point, 6>& points,
+ TetOp& op
+ );
+
+
+ //- Helping cľasses
+ class dummyTetOp
+ {
+ public:
+
+ inline void operator()(const tetPoints&){}
+ };
+
+
+ class sumTetVolOp
+ {
+ public:
+ scalar vol_;
+
+ inline sumTetVolOp()
+ :
+ vol_(0.0)
+ {}
+
+ inline void operator()(const tetPoints& tet)
+ {
+ vol_ += tet.tet().mag();
+ }
+ };
+
+
+ class storeTetOp
+ {
+ FixedList<tetPoints, 200>& tets_;
+ label& nTets_;
+
+ public:
+
+ inline storeTetOp(FixedList<tetPoints, 200>& tets, label& nTets)
+ :
+ tets_(tets),
+ nTets_(nTets)
+ {}
+
+ inline void operator()(const tetPoints& tet)
+ {
+ tets_[nTets_++] = tet;
+ }
+ };
+
+
+ //- Slice. Split tet into subtets above and below plane
+ template <class AboveTetOp, class BelowTetOp>
+ inline void tetSliceWithPlane
+ (
+ const tetPoints& tet,
+ const plane& pl,
+
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+ );
+
+
+ //- Tet overlap
+ void tetTetOverlap
+ (
+ const tetPoints& tetA,
+ const tetPoints& tetB,
+ FixedList<tetPoints, 200>& insideTets,
+ label& nInside,
+ FixedList<tetPoints, 200>& outsideTets,
+ label& nOutside
+ );
+
+
+ //- Tet Overlap Vol
+ inline scalar tetTetOverlapVol
+ (
+ const tetPoints& tetA,
+ const tetPoints& tetB
+ );
+
+
+ //- Return a const treeBoundBox
+ inline const treeBoundBox pyrBb
+ (
+ const pointField& points,
+ const face& f,
+ const point& fc
+ );
+
+
+public:
+
+
+ //- Runtime type information
+ TypeName("tetOverlapVolume");
+
+
+ // Constructors
+
+ //- Null constructor
+ tetOverlapVolume();
+
+
+ //- Destructor
+ virtual ~tetOverlapVolume();
+
+
+ // Public members
+
+ //- Return a list of cells in meshA which overlaps with cellBI in
+ // meshB
+ labelList overlappingCells
+ (
+ const fvMesh& meshA,
+ const fvMesh& meshB,
+ const label cellBI
+ ) const;
+
+
+ //- Calculates the overlap volume
+ scalar cellCellOverlapVolumeMinDecomp
+ (
+ const primitiveMesh& meshA,
+ const label cellAI,
+
+ const primitiveMesh& meshB,
+ const label cellBI,
+ const treeBoundBox& cellBbB
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H
new file mode 100644
index 00000000000..81499832688
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::tetPoints
+
+Description
+ Tet storage. Null constructable (unfortunately tetrahedron<point, point>
+ is not)
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetPoints_H
+#define tetPoints_H
+
+#include "FixedList.H"
+#include "treeBoundBox.H"
+#include "tetPointRef.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class tetPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetPoints
+:
+ public FixedList<point, 4>
+{
+public:
+
+ // Constructors
+
+ //- Construct null
+ inline tetPoints()
+ {}
+
+ //- Construct from four points
+ inline tetPoints
+ (
+ const point& a,
+ const point& b,
+ const point& c,
+ const point& d
+ )
+ {
+ operator[](0) = a;
+ operator[](1) = b;
+ operator[](2) = c;
+ operator[](3) = d;
+ }
+
+ // Member Functions
+
+ //- Return the tetrahedron
+ inline tetPointRef tet() const
+ {
+ return tetPointRef
+ (
+ operator[](0),
+ operator[](1),
+ operator[](2),
+ operator[](3)
+ );
+ }
+
+ //- Calculate the bounding box
+ inline treeBoundBox bounds() const
+ {
+ treeBoundBox bb(operator[](0), operator[](0));
+ for (label i = 1; i < size(); i++)
+ {
+ bb.min() = min(bb.min(), operator[](i));
+ bb.max() = max(bb.max(), operator[](i));
+ }
+ return bb;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C
new file mode 100644
index 00000000000..14c59d5eca8
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C
@@ -0,0 +1,341 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Description
+ Calculation of shape function product for a tetrahedron
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetrahedron.H"
+#include "triPointRef.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+// (Probably very inefficient) minimum containment sphere calculation.
+// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
+// Sphere ctr is smallest one of
+// - tet circumcentre
+// - triangle circumcentre
+// - edge mids
+template<class Point, class PointRef>
+Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
+(
+ const scalar tol
+) const
+{
+ const scalar fac = 1 + tol;
+
+ // Halve order of tolerance for comparisons of sqr.
+ const scalar facSqr = Foam::sqrt(fac);
+
+
+ // 1. Circumcentre itself.
+
+ pointHit pHit(circumCentre());
+ pHit.setHit();
+ scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
+
+
+ // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
+ // check if 4th point is inside.
+
+ {
+ point ctr = triPointRef(a_, b_, c_).circumCentre();
+ scalar radiusSqr = magSqr(ctr - a_);
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+ {
+ point ctr = triPointRef(a_, b_, d_).circumCentre();
+ scalar radiusSqr = magSqr(ctr - a_);
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+ {
+ point ctr = triPointRef(a_, c_, d_).circumCentre();
+ scalar radiusSqr = magSqr(ctr - a_);
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+ {
+ point ctr = triPointRef(b_, c_, d_).circumCentre();
+ scalar radiusSqr = magSqr(ctr - b_);
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+
+ // 3. Try midpoints of edges
+
+ // mid of edge A-B
+ {
+ point ctr = 0.5*(a_ + b_);
+ scalar radiusSqr = magSqr(a_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(c_-ctr) <= testRadSrq
+ && magSqr(d_-ctr) <= testRadSrq)
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+ // mid of edge A-C
+ {
+ point ctr = 0.5*(a_ + c_);
+ scalar radiusSqr = magSqr(a_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(b_-ctr) <= testRadSrq
+ && magSqr(d_-ctr) <= testRadSrq
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+ // mid of edge A-D
+ {
+ point ctr = 0.5*(a_ + d_);
+ scalar radiusSqr = magSqr(a_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(b_-ctr) <= testRadSrq
+ && magSqr(c_-ctr) <= testRadSrq
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+ // mid of edge B-C
+ {
+ point ctr = 0.5*(b_ + c_);
+ scalar radiusSqr = magSqr(b_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(a_-ctr) <= testRadSrq
+ && magSqr(d_-ctr) <= testRadSrq
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+ // mid of edge B-D
+ {
+ point ctr = 0.5*(b_ + d_);
+ scalar radiusSqr = magSqr(b_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(a_-ctr) <= testRadSrq
+ && magSqr(c_-ctr) <= testRadSrq)
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+ // mid of edge C-D
+ {
+ point ctr = 0.5*(c_ + d_);
+ scalar radiusSqr = magSqr(c_ - ctr);
+ scalar testRadSrq = facSqr*radiusSqr;
+
+ if
+ (
+ radiusSqr < minRadiusSqr
+ && magSqr(a_-ctr) <= testRadSrq
+ && magSqr(b_-ctr) <= testRadSrq
+ )
+ {
+ pHit.setMiss(false);
+ pHit.setPoint(ctr);
+ minRadiusSqr = radiusSqr;
+ }
+ }
+
+
+ pHit.setDistance(sqrt(minRadiusSqr));
+
+ return pHit;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiSquared
+(
+ scalarField& buffer
+) const
+{
+ // Change of sign between face area vector and gradient
+ // does not matter because of square
+
+ // Warning: Added a mag to produce positive coefficients even if
+ // the tetrahedron is twisted inside out. This is pretty
+ // dangerous, but essential for mesh motion.
+ scalar magVol = Foam::mag(mag());
+
+ buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
+ buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
+ buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
+ buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiDotGradNj
+(
+ scalarField& buffer
+) const
+{
+ // Warning. Ordering of edges needs to be the same for a tetrahedron
+ // class, a tetrahedron cell shape model and a tetCell
+
+ // Warning: Added a mag to produce positive coefficients even if
+ // the tetrahedron is twisted inside out. This is pretty
+ // dangerous, but essential for mesh motion.
+
+ // Double change of sign between face area vector and gradient
+
+ scalar magVol = Foam::mag(mag());
+ vector sa = Sa();
+ vector sb = Sb();
+ vector sc = Sc();
+ vector sd = Sd();
+
+ buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
+ buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
+ buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
+ buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
+ buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
+ buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiGradNi
+(
+ tensorField& buffer
+) const
+{
+ // Change of sign between face area vector and gradient
+ // does not matter because of square
+
+ scalar magVol = Foam::mag(mag());
+
+ buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
+ buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
+ buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
+ buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
+(
+ tensorField& buffer
+) const
+{
+ // Warning. Ordering of edges needs to be the same for a tetrahedron
+ // class, a tetrahedron cell shape model and a tetCell
+
+ // Double change of sign between face area vector and gradient
+
+ scalar magVol = Foam::mag(mag());
+ vector sa = Sa();
+ vector sb = Sb();
+ vector sc = Sc();
+ vector sd = Sd();
+
+ buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
+ buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
+ buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
+ buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
+ buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
+ buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H
new file mode 100644
index 00000000000..0424f91b118
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H
@@ -0,0 +1,315 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::tetrahedron
+
+Description
+ A tetrahedron primitive.
+
+ Ordering of edges needs to be the same for a tetrahedron
+ class, a tetrahedron cell shape model and a tetCell.
+
+SourceFiles
+ tetrahedronI.H
+ tetrahedron.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetrahedron_H
+#define tetrahedron_H
+
+#include "point.H"
+#include "primitiveFieldsFwd.H"
+#include "pointHit.H"
+#include "cachedRandom.H"
+#include "Random.H"
+#include "FixedList.H"
+#include "UList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class Istream;
+class Ostream;
+class tetPoints;
+class plane;
+
+// Forward declaration of friend functions and operators
+
+template<class Point, class PointRef> class tetrahedron;
+
+template<class Point, class PointRef>
+inline Istream& operator>>
+(
+ Istream&,
+ tetrahedron<Point, PointRef>&
+);
+
+template<class Point, class PointRef>
+inline Ostream& operator<<
+(
+ Ostream&,
+ const tetrahedron<Point, PointRef>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+ class tetrahedron Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Point, class PointRef>
+class tetrahedron
+{
+public:
+
+ // Classes for use in sliceWithPlane. What to do with decomposition
+ // of tet.
+
+ //- Dummy
+ class dummyOp
+ {
+ public:
+ inline void operator()(const tetPoints&);
+ };
+
+ //- Sum resulting volumes
+ class sumVolOp
+ {
+ public:
+ scalar vol_;
+
+ inline sumVolOp();
+
+ inline void operator()(const tetPoints&);
+ };
+
+ //- Store resulting tets
+ class storeOp
+ {
+ FixedList<tetPoints, 200>& tets_;
+ label& nTets_;
+
+ public:
+ inline storeOp(FixedList<tetPoints, 200>&, label&);
+
+ inline void operator()(const tetPoints&);
+ };
+
+private:
+
+ // Private data
+
+ PointRef a_, b_, c_, d_;
+
+ inline static point planeIntersection
+ (
+ const FixedList<scalar, 4>&,
+ const tetPoints&,
+ const label,
+ const label
+ );
+
+ template<class TetOp>
+ inline static void decomposePrism
+ (
+ const FixedList<point, 6>& points,
+ TetOp& op
+ );
+
+ template<class AboveTetOp, class BelowTetOp>
+ inline static void tetSliceWithPlane
+ (
+ const plane& pl,
+ const tetPoints& tet,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+ );
+
+
+public:
+
+ // Member constants
+
+ enum
+ {
+ nVertices = 4, // Number of vertices in tetrahedron
+ nEdges = 6 // Number of edges in tetrahedron
+ };
+
+
+ // Constructors
+
+ //- Construct from points
+ inline tetrahedron
+ (
+ const Point& a,
+ const Point& b,
+ const Point& c,
+ const Point& d
+ );
+
+ //- Construct from four points in the list of points
+ inline tetrahedron
+ (
+ const UList<Point>&,
+ const FixedList<label, 4>& indices
+ );
+
+ //- Construct from Istream
+ inline tetrahedron(Istream&);
+
+
+ // Member Functions
+
+ // Access
+
+ //- Return vertices
+ inline const Point& a() const;
+
+ inline const Point& b() const;
+
+ inline const Point& c() const;
+
+ inline const Point& d() const;
+
+
+ // Properties
+
+ //- Return face normal
+ inline vector Sa() const;
+
+ inline vector Sb() const;
+
+ inline vector Sc() const;
+
+ inline vector Sd() const;
+
+ //- Return centre (centroid)
+ inline Point centre() const;
+
+ //- Return volume
+ inline scalar mag() const;
+
+ //- Return circum-centre
+ inline Point circumCentre() const;
+
+ //- Return circum-radius
+ inline scalar circumRadius() const;
+
+ //- Return quality: Ratio of tetrahedron and circum-sphere
+ // volume, scaled so that a regular tetrahedron has a
+ // quality of 1
+ inline scalar quality() const;
+
+ //- Return a random point in the tetrahedron from a
+ // uniform distribution
+ inline Point randomPoint(Random& rndGen) const;
+
+ //- Return a random point in the tetrahedron from a
+ // uniform distribution
+ inline Point randomPoint(cachedRandom& rndGen) const;
+
+ //- Calculate the barycentric coordinates of the given
+ // point, in the same order as a, b, c, d. Returns the
+ // determinant of the solution.
+ inline scalar barycentric
+ (
+ const point& pt,
+ List<scalar>& bary
+ ) const;
+
+ //- Return nearest point to p on tetrahedron
+ inline pointHit nearestPoint(const point& p) const;
+
+ //- Return true if point is inside tetrahedron
+ inline bool inside(const point& pt) const;
+
+ //- Decompose tet into tets above and below plane
+ template<class AboveTetOp, class BelowTetOp>
+ inline void sliceWithPlane
+ (
+ const plane& pl,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+ ) const;
+
+
+ //- Return (min)containment sphere, i.e. the smallest sphere with
+ // all points inside. Returns pointHit with:
+ // - hit : if sphere is equal to circumsphere
+ // (biggest sphere)
+ // - point : centre of sphere
+ // - distance : radius of sphere
+ // - eligiblemiss: false
+ // Tol (small compared to 1, e.g. 1E-9) is used to determine
+ // whether point is inside: mag(pt - ctr) < (1+tol)*radius.
+ pointHit containmentSphere(const scalar tol) const;
+
+ //- Fill buffer with shape function products
+ void gradNiSquared(scalarField& buffer) const;
+
+ void gradNiDotGradNj(scalarField& buffer) const;
+
+ void gradNiGradNi(tensorField& buffer) const;
+
+ void gradNiGradNj(tensorField& buffer) const;
+
+
+ // IOstream operators
+
+ friend Istream& operator>> <Point, PointRef>
+ (
+ Istream&,
+ tetrahedron&
+ );
+
+ friend Ostream& operator<< <Point, PointRef>
+ (
+ Ostream&,
+ const tetrahedron&
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "tetrahedronI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "tetrahedron.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H
new file mode 100644
index 00000000000..6bd780aa550
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H
@@ -0,0 +1,1012 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triangle.H"
+#include "IOstreams.H"
+#include "triPointRef.H"
+#include "tetPoints.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::tetrahedron
+(
+ const Point& a,
+ const Point& b,
+ const Point& c,
+ const Point& d
+)
+:
+ a_(a),
+ b_(b),
+ c_(c),
+ d_(d)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::tetrahedron
+(
+ const UList<Point>& points,
+ const FixedList<label, 4>& indices
+)
+:
+ a_(points[indices[0]]),
+ b_(points[indices[1]]),
+ c_(points[indices[2]]),
+ d_(points[indices[3]])
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::tetrahedron(Istream& is)
+{
+ is >> *this;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
+{
+ return a_;
+}
+
+
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
+{
+ return b_;
+}
+
+
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
+{
+ return c_;
+}
+
+
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
+{
+ return d_;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sa() const
+{
+ return triangle<Point, PointRef>(b_, c_, d_).normal();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sb() const
+{
+ return triangle<Point, PointRef>(a_, d_, c_).normal();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sc() const
+{
+ return triangle<Point, PointRef>(a_, b_, d_).normal();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::vector Foam::tetrahedron<Point, PointRef>::Sd() const
+{
+ return triangle<Point, PointRef>(a_, c_, b_).normal();
+}
+
+
+template<class Point, class PointRef>
+inline Point Foam::tetrahedron<Point, PointRef>::centre() const
+{
+ return 0.25*(a_ + b_ + c_ + d_);
+}
+
+
+template<class Point, class PointRef>
+inline Foam::scalar Foam::tetrahedron<Point, PointRef>::mag() const
+{
+ return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_));
+}
+
+
+template<class Point, class PointRef>
+inline Point Foam::tetrahedron<Point, PointRef>::circumCentre() const
+{
+ vector a = b_ - a_;
+ vector b = c_ - a_;
+ vector c = d_ - a_;
+
+ scalar lambda = magSqr(c) - (a & c);
+ scalar mu = magSqr(b) - (a & b);
+
+ vector ba = b ^ a;
+ vector ca = c ^ a;
+
+ vector num = lambda*ba - mu*ca;
+ scalar denom = (c & ba);
+
+ if (Foam::mag(denom) < ROOTVSMALL)
+ {
+ // Degenerate tetrahedron, returning centre instead of circumCentre.
+
+ return centre();
+ }
+
+ return a_ + 0.5*(a + num/denom);
+}
+
+
+template<class Point, class PointRef>
+inline Foam::scalar Foam::tetrahedron<Point, PointRef>::circumRadius() const
+{
+ vector a = b_ - a_;
+ vector b = c_ - a_;
+ vector c = d_ - a_;
+
+ scalar lambda = magSqr(c) - (a & c);
+ scalar mu = magSqr(b) - (a & b);
+
+ vector ba = b ^ a;
+ vector ca = c ^ a;
+
+ vector num = lambda*ba - mu*ca;
+ scalar denom = (c & ba);
+
+ if (Foam::mag(denom) < ROOTVSMALL)
+ {
+ // Degenerate tetrahedron, returning GREAT for circumRadius.
+ return GREAT;
+ }
+
+ return Foam::mag(0.5*(a + num/denom));
+}
+
+
+template<class Point, class PointRef>
+inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
+{
+ return
+ mag()
+ /(
+ 8.0/(9.0*sqrt(3.0))
+ *pow3(min(circumRadius(), GREAT))
+ + ROOTVSMALL
+ );
+}
+
+
+template<class Point, class PointRef>
+inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
+(
+ Random& rndGen
+) const
+{
+ // Adapted from
+ // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
+
+ scalar s = rndGen.scalar01();
+ scalar t = rndGen.scalar01();
+ scalar u = rndGen.scalar01();
+
+ if (s + t > 1.0)
+ {
+ s = 1.0 - s;
+ t = 1.0 - t;
+ }
+
+ if (t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = 1.0 - s - t;
+ t = 1.0 - tmp;
+ }
+ else if (s + t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = s + t + u - 1.0;
+ s = 1.0 - t - tmp;
+ }
+
+ return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
+}
+
+
+template<class Point, class PointRef>
+inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
+(
+ cachedRandom& rndGen
+) const
+{
+ // Adapted from
+ // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
+
+ scalar s = rndGen.sample01<scalar>();
+ scalar t = rndGen.sample01<scalar>();
+ scalar u = rndGen.sample01<scalar>();
+
+ if (s + t > 1.0)
+ {
+ s = 1.0 - s;
+ t = 1.0 - t;
+ }
+
+ if (t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = 1.0 - s - t;
+ t = 1.0 - tmp;
+ }
+ else if (s + t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = s + t + u - 1.0;
+ s = 1.0 - t - tmp;
+ }
+
+ return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
+}
+
+
+template<class Point, class PointRef>
+Foam::scalar Foam::tetrahedron<Point, PointRef>::barycentric
+(
+ const point& pt,
+ List<scalar>& bary
+) const
+{
+ // From:
+ // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
+
+ vector e0(a_ - d_);
+ vector e1(b_ - d_);
+ vector e2(c_ - d_);
+
+ tensor t
+ (
+ e0.x(), e1.x(), e2.x(),
+ e0.y(), e1.y(), e2.y(),
+ e0.z(), e1.z(), e2.z()
+ );
+
+ scalar detT = det(t);
+
+ if (Foam::mag(detT) < SMALL)
+ {
+ // Degenerate tetrahedron, returning 1/4 barycentric coordinates.
+
+ bary = List<scalar>(4, 0.25);
+
+ return detT;
+ }
+
+ vector res = inv(t, detT) & (pt - d_);
+
+ bary.setSize(4);
+
+ bary[0] = res.x();
+ bary[1] = res.y();
+ bary[2] = res.z();
+ bary[3] = (1.0 - res.x() - res.y() - res.z());
+
+ return detT;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
+(
+ const point& p
+) const
+{
+ // Adapted from:
+ // Real-time collision detection, Christer Ericson, 2005, p142-144
+
+ // Assuming initially that the point is inside all of the
+ // halfspaces, so closest to itself.
+
+ point closestPt = p;
+
+ scalar minOutsideDistance = VGREAT;
+
+ bool inside = true;
+
+ if (((p - b_) & Sa()) >= 0)
+ {
+ // p is outside halfspace plane of tri
+ pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
+
+ inside = false;
+
+ if (info.distance() < minOutsideDistance)
+ {
+ closestPt = info.rawPoint();
+
+ minOutsideDistance = info.distance();
+ }
+ }
+
+ if (((p - a_) & Sb()) >= 0)
+ {
+ // p is outside halfspace plane of tri
+ pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
+
+ inside = false;
+
+ if (info.distance() < minOutsideDistance)
+ {
+ closestPt = info.rawPoint();
+
+ minOutsideDistance = info.distance();
+ }
+ }
+
+ if (((p - a_) & Sc()) >= 0)
+ {
+ // p is outside halfspace plane of tri
+ pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
+
+ inside = false;
+
+ if (info.distance() < minOutsideDistance)
+ {
+ closestPt = info.rawPoint();
+
+ minOutsideDistance = info.distance();
+ }
+ }
+
+ if (((p - a_) & Sd()) >= 0)
+ {
+ // p is outside halfspace plane of tri
+ pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
+
+ inside = false;
+
+ if (info.distance() < minOutsideDistance)
+ {
+ closestPt = info.rawPoint();
+
+ minOutsideDistance = info.distance();
+ }
+ }
+
+ // If the point is inside, then the distance to the closest point
+ // is zero
+ if (inside)
+ {
+ minOutsideDistance = 0;
+ }
+
+ return pointHit
+ (
+ inside,
+ closestPt,
+ minOutsideDistance,
+ !inside
+ );
+}
+
+
+template<class Point, class PointRef>
+bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
+{
+ // For robustness, assuming that the point is in the tet unless
+ // "definitively" shown otherwise by obtaining a positive dot
+ // product greater than a tolerance of SMALL.
+
+ // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
+ // vectors and base points for the half-space planes are:
+ // area[0] = Sa();
+ // area[1] = Sb();
+ // area[2] = Sc();
+ // area[3] = Sd();
+ // planeBase[0] = tetBasePt = b_
+ // planeBase[1] = ptA = c_
+ // planeBase[2] = tetBasePt = b_
+ // planeBase[3] = tetBasePt = b_
+
+ vector n = vector::zero;
+
+ {
+ // 0, a
+ const point& basePt = b_;
+
+ n = Sa();
+ n /= (Foam::mag(n) + VSMALL);
+
+ if (((pt - basePt) & n) > SMALL)
+ {
+ return false;
+ }
+ }
+
+ {
+ // 1, b
+ const point& basePt = c_;
+
+ n = Sb();
+ n /= (Foam::mag(n) + VSMALL);
+
+ if (((pt - basePt) & n) > SMALL)
+ {
+ return false;
+ }
+ }
+
+ {
+ // 2, c
+ const point& basePt = b_;
+
+ n = Sc();
+ n /= (Foam::mag(n) + VSMALL);
+
+ if (((pt - basePt) & n) > SMALL)
+ {
+ return false;
+ }
+ }
+
+ {
+ // 3, d
+ const point& basePt = b_;
+
+ n = Sd();
+ n /= (Foam::mag(n) + VSMALL);
+
+ if (((pt - basePt) & n) > SMALL)
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::dummyOp::operator()
+(
+ const tetPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::sumVolOp::sumVolOp()
+:
+ vol_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
+(
+ const tetPoints& tet
+)
+{
+ vol_ += tet.tet().mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
+(
+ FixedList<tetPoints, 200>& tets,
+ label& nTets
+)
+:
+ tets_(tets),
+ nTets_(nTets)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::storeOp::operator()
+(
+ const tetPoints& tet
+)
+{
+ tets_[nTets_++] = tet;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::tetrahedron<Point, PointRef>::planeIntersection
+(
+ const FixedList<scalar, 4>& d,
+ const tetPoints& t,
+ const label negI,
+ const label posI
+)
+{
+ return
+ (d[posI]*t[negI] - d[negI]*t[posI])
+ / (-d[negI]+d[posI]);
+}
+
+
+template<class Point, class PointRef>
+template<class TetOp>
+inline void Foam::tetrahedron<Point, PointRef>::decomposePrism
+(
+ const FixedList<point, 6>& points,
+ TetOp& op
+)
+{
+ op(tetPoints(points[1], points[3], points[2], points[0]));
+ op(tetPoints(points[1], points[2], points[3], points[4]));
+ op(tetPoints(points[4], points[2], points[3], points[5]));
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::
+tetSliceWithPlane
+(
+ const plane& pl,
+ const tetPoints& tet,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+)
+{
+ // Distance to plane
+ FixedList<scalar, 4> d;
+ label nPos = 0;
+ forAll(tet, i)
+ {
+ d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
+ if (d[i] > 0)
+ {
+ nPos++;
+ }
+ }
+
+ if (nPos == 4)
+ {
+ aboveOp(tet);
+ }
+ else if (nPos == 3)
+ {
+ // Sliced into below tet and above prism. Prism gets split into
+ // two tets.
+
+ // Find the below tet
+ label i0 = -1;
+ forAll(d, i)
+ {
+ if (d[i] <= 0)
+ {
+ i0 = i;
+ break;
+ }
+ }
+
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+ label i3 = d.fcIndex(i2);
+
+ point p01 = planeIntersection(d, tet, i0, i1);
+ point p02 = planeIntersection(d, tet, i0, i2);
+ point p03 = planeIntersection(d, tet, i0, i3);
+
+ // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
+ // ,, 1 : ,, inwards pointing triad
+ // ,, 2 : ,, outwards pointing triad
+ // ,, 3 : ,, inwards pointing triad
+
+ //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
+
+ if (i0 == 0 || i0 == 2)
+ {
+ tetPoints t(tet[i0], p01, p02, p03);
+ //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 3, belowTet i0==0 or 2");
+ belowOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i1];
+ p[1] = tet[i3];
+ p[2] = tet[i2];
+ p[3] = p01;
+ p[4] = p03;
+ p[5] = p02;
+ //Pout<< " aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ else
+ {
+ tetPoints t(p01, p02, p03, tet[i0]);
+ //Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 3, belowTet i0==1 or 3");
+ belowOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i3];
+ p[1] = tet[i1];
+ p[2] = tet[i2];
+ p[3] = p03;
+ p[4] = p01;
+ p[5] = p02;
+ //Pout<< " aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ }
+ else if (nPos == 2)
+ {
+ // Tet cut into two prisms. Determine the positive one.
+ label pos0 = -1;
+ label pos1 = -1;
+ label neg0 = -1;
+ label neg1 = -1;
+ forAll(d, i)
+ {
+ if (d[i] > 0)
+ {
+ if (pos0 == -1)
+ {
+ pos0 = i;
+ }
+ else
+ {
+ pos1 = i;
+ }
+ }
+ else
+ {
+ if (neg0 == -1)
+ {
+ neg0 = i;
+ }
+ else
+ {
+ neg1 = i;
+ }
+ }
+ }
+
+ //Pout<< "Split 2pos tet " << tet << " d:" << d
+ // << " around pos0:" << pos0 << " pos1:" << pos1
+ // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
+
+ const edge posEdge(pos0, pos1);
+
+ if (posEdge == edge(0, 1))
+ {
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[0];
+ p[1] = p02;
+ p[2] = p03;
+ p[3] = tet[1];
+ p[4] = p12;
+ p[5] = p13;
+ //Pout<< " 01 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p02;
+ p[2] = p12;
+ p[3] = tet[3];
+ p[4] = p03;
+ p[5] = p13;
+ //Pout<< " 01 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(1, 2))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p01;
+ p[2] = p13;
+ p[3] = tet[2];
+ p[4] = p02;
+ p[5] = p23;
+ //Pout<< " 12 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[3];
+ p[1] = p23;
+ p[2] = p13;
+ p[3] = tet[0];
+ p[4] = p02;
+ p[5] = p01;
+ //Pout<< " 12 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(2, 0))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p12;
+ p[2] = p23;
+ p[3] = tet[0];
+ p[4] = p01;
+ p[5] = p03;
+ //Pout<< " 20 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p12;
+ p[2] = p01;
+ p[3] = tet[3];
+ p[4] = p23;
+ p[5] = p03;
+ //Pout<< " 20 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(0, 3))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[3];
+ p[1] = p23;
+ p[2] = p13;
+ p[3] = tet[0];
+ p[4] = p02;
+ p[5] = p01;
+ //Pout<< " 03 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p23;
+ p[2] = p02;
+ p[3] = tet[1];
+ p[4] = p13;
+ p[5] = p01;
+ //Pout<< " 03 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(1, 3))
+ {
+ point p01 = planeIntersection(d, tet, 0, 1);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p23 = planeIntersection(d, tet, 2, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[1];
+ p[1] = p12;
+ p[2] = p01;
+ p[3] = tet[3];
+ p[4] = p23;
+ p[5] = p03;
+ //Pout<< " 13 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p12;
+ p[2] = p23;
+ p[3] = tet[0];
+ p[4] = p01;
+ p[5] = p03;
+ //Pout<< " 13 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else if (posEdge == edge(2, 3))
+ {
+ point p02 = planeIntersection(d, tet, 0, 2);
+ point p12 = planeIntersection(d, tet, 1, 2);
+ point p03 = planeIntersection(d, tet, 0, 3);
+ point p13 = planeIntersection(d, tet, 1, 3);
+ // Split the resulting prism
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[2];
+ p[1] = p02;
+ p[2] = p12;
+ p[3] = tet[3];
+ p[4] = p03;
+ p[5] = p13;
+ //Pout<< " 23 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList<point, 6> p;
+ p[0] = tet[0];
+ p[1] = p02;
+ p[2] = p03;
+ p[3] = tet[1];
+ p[4] = p12;
+ p[5] = p13;
+ //Pout<< " 23 belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else
+ {
+ FatalErrorIn("tetSliceWithPlane(..)")
+ << "Missed edge:" << posEdge
+ << abort(FatalError);
+ }
+ }
+ else if (nPos == 1)
+ {
+ // Find the positive tet
+ label i0 = -1;
+ forAll(d, i)
+ {
+ if (d[i] > 0)
+ {
+ i0 = i;
+ break;
+ }
+ }
+
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+ label i3 = d.fcIndex(i2);
+
+ point p01 = planeIntersection(d, tet, i0, i1);
+ point p02 = planeIntersection(d, tet, i0, i2);
+ point p03 = planeIntersection(d, tet, i0, i3);
+
+ //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
+
+ if (i0 == 0 || i0 == 2)
+ {
+ tetPoints t(tet[i0], p01, p02, p03);
+ //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
+ aboveOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i1];
+ p[1] = tet[i3];
+ p[2] = tet[i2];
+ p[3] = p01;
+ p[4] = p03;
+ p[5] = p02;
+ //Pout<< " belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ else
+ {
+ tetPoints t(p01, p02, p03, tet[i0]);
+ //Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
+ //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
+ aboveOp(t);
+
+ // Prism
+ FixedList<point, 6> p;
+ p[0] = tet[i3];
+ p[1] = tet[i1];
+ p[2] = tet[i2];
+ p[3] = p03;
+ p[4] = p01;
+ p[5] = p02;
+ //Pout<< " belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else // nPos == 0
+ {
+ belowOp(tet);
+ }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::sliceWithPlane
+(
+ const plane& pl,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+) const
+{
+ tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+inline Foam::Istream& Foam::operator>>
+(
+ Istream& is,
+ tetrahedron<Point, PointRef>& t
+)
+{
+ is.readBegin("tetrahedron");
+ is >> t.a_ >> t.b_ >> t.c_ >> t.d_;
+ is.readEnd("tetrahedron");
+
+ is.check("Istream& operator>>(Istream&, tetrahedron&)");
+
+ return is;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const tetrahedron<Point, PointRef>& t
+)
+{
+ os << nl
+ << token::BEGIN_LIST
+ << t.a_ << token::SPACE
+ << t.b_ << token::SPACE
+ << t.c_ << token::SPACE
+ << t.d_
+ << token::END_LIST;
+
+ return os;
+}
+
+
+// ************************************************************************* //
--
GitLab