From c4e19274696a9ba6ebc201daf31bc87a52dd9ca8 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 31 Oct 2017 10:28:05 +0000 Subject: [PATCH] AMI: Added sweptFaceAreaWeightAMI method This method projects the source patch to the target using the point normals. The projection fills space, which results in target weights that correctly sum to unity. A source patch face can still project onto an area larger or smaller than the face, so the source weights do not (in general) sum to unity as a result of this method. This has not been made the default AMI method. Further investigation is needed to asses the benefits of this sort of projection. --- .../AMIInterpolation/AMIInterpolation.C | 12 +- .../AMIInterpolation/AMIInterpolation.H | 3 +- .../sweptFaceAreaWeightAMI.C | 602 ++++++++++++++++++ .../sweptFaceAreaWeightAMI.H | 273 ++++++++ .../AMIPatchToPatchInterpolation.C | 5 +- 5 files changed, 892 insertions(+), 3 deletions(-) create mode 100644 src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C create mode 100644 src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.H diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index 825de3f335..9d27328d0f 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -62,6 +62,11 @@ Foam::AMIInterpolation::interpolationMethodToWord method = "partialFaceAreaWeightAMI"; break; } + case imSweptFaceAreaWeight: + { + method = "sweptFaceAreaWeightAMI"; + break; + } default: { FatalErrorInFunction @@ -91,7 +96,8 @@ Foam::AMIInterpolation::wordTointerpolationMethod "directAMI " "mapNearestAMI " "faceAreaWeightAMI " - "partialFaceAreaWeightAMI" + "partialFaceAreaWeightAMI " + "sweptFaceAreaWeightAMI" ")" )() ); @@ -112,6 +118,10 @@ Foam::AMIInterpolation::wordTointerpolationMethod { method = imPartialFaceAreaWeight; } + else if (im == "sweptFaceAreaWeightAMI") + { + method = imSweptFaceAreaWeight; + } else { FatalErrorInFunction diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index 4840f82fb5..1237a1d690 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -89,7 +89,8 @@ public: imDirect, imMapNearest, imFaceAreaWeight, - imPartialFaceAreaWeight + imPartialFaceAreaWeight, + imSweptFaceAreaWeight }; //- Convert interpolationMethod to word representation diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C new file mode 100644 index 0000000000..288ebce6ef --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C @@ -0,0 +1,602 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 . + +\*---------------------------------------------------------------------------*/ + +#include "sweptFaceAreaWeightAMI.H" +#include "cut.H" +#include "linearEqn.H" +#include "quadraticEqn.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +const Foam::scalar Foam::sweptFaceAreaWeightAMI:: +minCutRatio_ = 10*SMALL; + +template +const Foam::scalar Foam::sweptFaceAreaWeightAMI:: +maxDot_ = - cos(degToRad(89.9)); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +template +void Foam::sweptFaceAreaWeightAMI::writeCutTrisVTK +( + const cutTriList& tris, + const word& name +) const +{ + OFstream obj(typeName + "_" + name + ".vtk"); + + obj << "# vtk DataFile Version 2.0" << endl + << "triSurface" << endl + << "ASCII" << endl + << "DATASET POLYDATA" << endl + << "POINTS " << 3*tris.size() << " float" << endl; + + for (label i = 0; i < tris.size(); ++ i) + { + for (label j = 0; j < 3; ++ j) + { + const vector& x = tris[i][j]; + obj << x.x() << ' ' << x.y() << ' ' << x.z() << endl; + } + } + + obj << "POLYGONS " << tris.size() << ' ' << 4*tris.size() << endl; + + for (label i = 0; i < tris.size(); ++ i) + { + obj << 3 << ' ' << 3*i << ' ' << 3*i + 1 << ' ' << 3*i + 2 << endl; + } +} + + +template +void Foam::sweptFaceAreaWeightAMI::writeFaceOBJ +( + const face& f, + const pointField& ps, + const string& name +) const +{ + OFstream obj(typeName + "_" + name + ".obj"); + + for (label i = 0; i < f.size(); ++ i) + { + const vector& x = ps[f[i]]; + obj << "v " << x.x() << ' ' << x.y() << ' ' << x.z() << endl; + } + + obj << 'f'; + for (label i = 0; i < f.size(); ++ i) + { + obj << ' ' << i + 1; + } + obj << endl; +} + + +template +void Foam::sweptFaceAreaWeightAMI::writeProjectionOBJ +( + const label srcN, + const FixedList& srcTri, + const FixedList& srcNrm +) const +{ + scalar l = 0; + for (label i = 0; i < srcN - 1; ++ i) + { + l = max(l, mag(srcTri[i] - srcTri[i + 1])); + } + + const label nu = 20, nv = 20; + const scalar u0 = 0, u1 = 1, v0 = -l, v1 = l; + + OFstream obj(typeName + "_projection.obj"); + + for (label i = 0; i < srcN; ++ i) + { + const point& p0 = srcTri[i], p1 = srcTri[(i + 1) % srcN]; + const vector& n0 = srcNrm[i], n1 = srcNrm[(i + 1) % srcN]; + + for (label iu = 0; iu <= nu; ++ iu) + { + const scalar u = u0 + (u1 - u0)*scalar(iu)/nu; + for (label iv = 0; iv <= nv; ++ iv) + { + const scalar v = v0 + (v1 - v0)*scalar(iv)/nv; + const vector x = p0 + (p1 - p0)*u + (n0 + (n1 - n0)*u)*v; + obj << "v " << x.x() << ' ' << x.y() << ' ' << x.z() << endl; + } + } + } + + for (label i = 0; i < srcN; ++ i) + { + for (label iu = 0; iu < nu; ++ iu) + { + for (label iv = 0; iv < nv; ++ iv) + { + obj << "f " + << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 1 << ' ' + << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 2 << ' ' + << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 2 << ' ' + << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 1<< endl; + } + } + } +} + + +template +Foam::label +Foam::sweptFaceAreaWeightAMI::getSourceProjection +( + FixedList& srcTri, + FixedList& srcNrm, + const FixedList& tgtTri +) const +{ + // The target normal + const vector np = plane(tgtTri[0], tgtTri[1], tgtTri[2]).normal(); + + // Dot products between the projection normals and the target plane + FixedList dots; + for (label i = 0; i < 3; ++ i) + { + dots[i] = srcNrm[i] & np; + } + dots[3] = 0; + + // Search though for the number of reversed directions and an index for + // both a reverse and a forward direction + label nR = 0, iR = -1, iF = -1; + for (label i = 0; i < 3; ++ i) + { + if (dots[i] > maxDot_) + { + ++ nR; + iR = i; + } + else + { + iF = i; + } + } + + // If all the normals hit in the forward or reverse direction then return + // the existing projection or no projection respectively + if (nR == 0 || nR == 3) + { + return 3 - nR; + } + + // One normal hits in the reverse direction + if (nR == 1) + { + /* + o j1 + / \ + / \ + / \ + / o jR + / / . + / / . + / / . + o - > - o . . . . (old iR) + i0 iR + */ + + // Make a gap by duplicating the reverse the point + for (label j = 2; j >= iR; -- j) + { + srcTri[j + 1] = srcTri[j]; + srcNrm[j + 1] = srcNrm[j]; + dots[j + 1] = dots[j]; + } + + const label jR = iR + 1; + const label i0 = (iR + 3) % 4, j1 = (jR + 1) % 4; + + const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]); + const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[j1]); + + srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi; + srcTri[jR] = srcTri[jR] + (srcTri[j1] - srcTri[jR])*wj; + + srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi; + srcNrm[jR] = srcNrm[jR] + (srcNrm[j1] - srcNrm[jR])*wj; + + return 4; + } + + // Two normals hit in the reverse direction + if (nR == 2) + { + /* + . (old jR) + . . + jR o . + / \ . + / \ . + / \ . + / \ . + / \ . + o - - > - - o . . (old iR) + iF iR + */ + + const label iR = (iF + 1) % 3, jR = (iF + 2) % 3; + const label i0 = (iR + 2) % 3, j1 = (jR + 1) % 3; + + const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]); + const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[j1]); + + srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi; + srcTri[jR] = srcTri[jR] + (srcTri[j1] - srcTri[jR])*wj; + + srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi; + srcNrm[jR] = srcNrm[jR] + (srcNrm[j1] - srcNrm[jR])*wj; + + return 3; + } + + // This cannot happen + return -1; +} + + +template +Foam::plane Foam::sweptFaceAreaWeightAMI::getCutPlane +( + const point& p0, + const point& p1, + const vector& n0, + const vector& n1, + const FixedList& tgtTri +) const +{ + // The target plane + const plane tgtPln(tgtTri[0], tgtTri[1], tgtTri[2]); + const point& pp = tgtPln.refPoint(); + const vector& np = tgtPln.normal(); + + // Calculate the bounding intersection points. These are the locations at + // which the bounding lines of the projected surface intersect with the + // target plane. + const scalar v0 = ((pp - p0) & np)/(n0 & np); + const scalar v1 = ((pp - p1) & np)/(n1 & np); + Pair cutPnts(p0 + v0*n0, p1 + v1*n1); + Pair cutNrms((p1 - p0 + n1*v0) ^ n0, (p1 - p0 - n0*v1) ^ n1); + + // Calculate edge intersections with the surface + for (label i = 0; i < 3; ++ i) + { + // The current target edge + const point& l0 = tgtTri[i], l1 = tgtTri[(i + 1)%3]; + + // Form the quadratic in the surface parameter u + const vector k = l0 - p0; + const vector a = l0 - l1, b = p1 - p0, c = n0, d = n1 - n0; + const vector ka = k ^ a, ba = b ^ a, ca = c ^ a, da = d ^ a; + const scalar A = d & ba; + const scalar B = (c & ba) - (d & ka); + const scalar C = - c & ka; + + // Solve and extract the other parameters + const Roots<2> us = quadraticEqn(A, B, C).roots(); + Pair vs, ws; + Pair ns; + Pair cuts(false, false); + for (label j = 0; j < 2; ++ j) + { + if (us.type(j) == roots::real) + { + const vector den = ca + da*us[j]; + if (magSqr(den) > VSMALL) + { + const vector vNum = ka - ba*us[j]; + const vector wNum = (- k + b*us[j]) ^ (c + d*us[j]); + vs[j] = (vNum & den)/magSqr(den); + ws[j] = (wNum & den)/magSqr(den); + ns[j] = (b ^ (c + d*us[j])) + (n1 ^ n0)*vs[j]; + const bool cutu = 0 < us[j] && us[j] < 1; + const bool cutw = 0 < ws[j] && ws[j] < 1; + cuts[j] = cutu && cutw; + } + } + } + + // If we have just one intersection in bounds then store it in the + // result list based on its direction + if (cuts[0] != cuts[1]) + { + const label j = cuts[0] ? 0 : 1; + const scalar na = ns[j] & a; + cutPnts[na < 0] = l0 + (l1 - l0)*ws[j]; + cutNrms[na < 0] = ns[j]; + } + } + + // Generate and return the cut plane. If the cut points are not coincident + // then form a plane normal by crossing the displacement between the points + // by the target plane normal. If the points are coincident then use the + // projected surface normal evaluated at the first cut point. + const vector cutDelta = cutPnts[1] - cutPnts[0]; + const bool coincident = mag(cutDelta) < minCutRatio_*mag(p1 - p0); + return plane(cutPnts[0], coincident ? cutNrms[0] : np ^ cutDelta); +}; + + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +template +Foam::scalar Foam::sweptFaceAreaWeightAMI::interArea +( + const label srcFacei, + const label tgtFacei +) const +{ + const label debugTgtFace = + (- 1 - debug) % this->tgtPatch_.size() == tgtFacei + ? (- 1 - debug) / this->tgtPatch_.size() + 1 : 0; + const bool debugPrint = debug > 0 || debugTgtFace > 0; + const bool debugWrite = debug > 1 || debugTgtFace > 1; + + if (debugPrint) + { + Info<< "Inter area between source face #" << srcFacei + << " and target face #" << tgtFacei << (debugWrite ? "\n" : ": "); + } + + // Patch data + const pointField& srcPoints = this->srcPatch_.localPoints(); + const pointField& tgtPoints = this->tgtPatch_.localPoints(); + const vectorField& srcPointNormals = this->srcPatch_.pointNormals(); + + // Faces + const face& srcFace = this->srcPatch_.localFaces()[srcFacei]; + const face& tgtFace = this->tgtPatch_.localFaces()[tgtFacei]; + + // Write out the faces + if (debugWrite) + { + writeFaceOBJ(srcFace, srcPoints, "source"); + writeFaceOBJ(tgtFace, tgtPoints, "target"); + } + if (debugTgtFace) + { + writeFaceOBJ(tgtFace, tgtPoints, "target" + std::to_string(tgtFacei)); + } + + // Triangulate the faces + const faceAreaIntersect::triangulationMode triMode = this->triMode_; + faceList srcFaceTris, tgtFaceTris; + faceAreaIntersect::triangulate(srcFace, srcPoints, triMode, srcFaceTris); + faceAreaIntersect::triangulate(tgtFace, tgtPoints, triMode, tgtFaceTris); + + // Area sum + scalar areaMag = Zero; + + // Loop the target triangles + forAllConstIter(faceList, tgtFaceTris, tgtIter) + { + const FixedList + tgtTri = + { + tgtPoints[(*tgtIter)[0]], + tgtPoints[(*tgtIter)[this->reverseTarget_ ? 2 : 1]], + tgtPoints[(*tgtIter)[this->reverseTarget_ ? 1 : 2]] + }; + + // Loop the source triangles + forAllConstIter(faceList, srcFaceTris, srcIter) + { + FixedList + srcTri = + { + srcPoints[(*srcIter)[0]], + srcPoints[(*srcIter)[1]], + srcPoints[(*srcIter)[2]], + vector::zero + }; + FixedList + srcNrm = + { + srcPointNormals[(*srcIter)[0]], + srcPointNormals[(*srcIter)[1]], + srcPointNormals[(*srcIter)[2]], + vector::zero + }; + + // Get the source projection modified for any reverse intersections + const label srcN = getSourceProjection(srcTri, srcNrm, tgtTri); + if (srcN <= 0) + { + continue; + } + + // Write the source projection + if (debugWrite) + { + writeProjectionOBJ(srcN, srcTri, srcNrm); + } + + // Create the initial cut triangle list + cutTriList<8> cutTris; + cutTris.append(tgtTri); + + // Write the initial triangle + if (debugWrite) + { + writeCutTrisVTK(cutTris, "tris0"); + } + + // Do all but one of the cuts + for + ( + label i = 0; + i < srcN - 1 && (debugWrite || cutTris.size()); + ++ i + ) + { + // Do the cut + const plane cutPlane = + getCutPlane + ( + srcTri[i], + srcTri[i+1], + srcNrm[i], + srcNrm[i+1], + tgtTri + ); + cutTriList<8> cutTrisTmp; + for (label i = 0; i < cutTris.size(); ++ i) + { + triCut + ( + cutTris[i], + cutPlane, + cut::noOp(), + cut::appendOp>(cutTrisTmp) + ); + } + Swap(cutTris, cutTrisTmp); + + // Write the triangles resulting from the cut + if (debugWrite) + { + writeCutTrisVTK(cutTris, "tris" + std::to_string(i + 1)); + } + } + + // Do the last cut + const plane cutPlane = + getCutPlane + ( + srcTri[srcN - 1], + srcTri[0], + srcNrm[srcN - 1], + srcNrm[0], + tgtTri + ); + cutTriList<8> cutTrisTmp; + for (label i = 0; i < cutTris.size(); ++ i) + { + // Sum the area of the cut triangle + areaMag += + mag + ( + triCut + ( + cutTris[i], + cutPlane, + cut::noOp(), + cut::areaOp() + ) + ); + // Store the cut triangle if it is needed for output + if (debugWrite || debugTgtFace) + { + triCut + ( + cutTris[i], + cutPlane, + cut::noOp(), + cut::appendOp>(cutTrisTmp) + ); + } + } + Swap(cutTris, cutTrisTmp); + + // Write the triangles resulting from the cuts + if (debugTgtFace && cutTris.size()) + { + static label writeCutTrisVTKIndex = 0; + writeCutTrisVTK + ( + cutTris, + "target" + std::to_string(tgtFacei) + "_" + + "tris" + std::to_string(writeCutTrisVTKIndex ++) + ); + } + if (debugWrite) + { + writeCutTrisVTK(cutTris, "tris" + std::to_string(srcN)); + + Info << "view triangles then press ENTER to continue ..."; + getchar(); + } + } + } + + // Print the difference between this inter-area and that obtained by the + // standard algorithm built around faceAreaIntersect + if (debugPrint) + { + const scalar standardAreaMag = + faceAreaWeightAMI::interArea + ( + srcFacei, + tgtFacei + ); + + Info<<"standard=" << standardAreaMag << ", swept=" << areaMag + << endl; + } + + return areaMag; +} + + +template +Foam::scalar +Foam::sweptFaceAreaWeightAMI::minWeight() const +{ + return SMALL; +} + + +template +Foam::scalar +Foam::sweptFaceAreaWeightAMI::maxWalkAngle() const +{ + return degToRad(180); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::sweptFaceAreaWeightAMI::~sweptFaceAreaWeightAMI +( +) +{} + + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.H new file mode 100644 index 0000000000..519e68752d --- /dev/null +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.H @@ -0,0 +1,273 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 . + +Class + Foam::sweptFaceAreaWeightAMI + +Description + Swept face area weighted Arbitrary Mesh Interface (AMI) method + + This method uses the point normals of the source patch to sweep the source + faces over the target patches. This creates a projection which fills space, + and which therefore generates overlap areas which are consistent between + neighbouring faces. The projection of a source edge is shown below: + \verbatim + / + / + n_1 ^ + / + / + p_1 o + \\ + \\ ^ u + \\ \ + \\ + - > v + \\ + o - - - > - - - + p_0 n_0 + \endverbatim + + The surface is, in general, not flat. Any deviation between the two end + normals generates expansion, contraction, and twist in the surface. The + surface connects with surfaces eminating from connected edges along the end + normals. This is what makes the projection fill space and generate + consistent overlaps between neighbouring faces. + + The projected surface is parameterised by the local coordinates, \f$u\f$ + and \f$v\f$, and a position on this plane is calculated from \f$u\f$ and + \f$v\f$ as follows: + \f[ + x(u, v) = p_0 + (p_1 - p_0) u + [ q_0 + (q_1 - q_0) u ] v + \f] + + To calculate an intersection with a line between points \f$l_0\f$ to + \f$l_1\f$, we define a local coordinate, \f$w\f$, along the line, and + subtract it's equation from that of the surface: + \f[ + 0 = (p_0 - l_0) - (l_1 - l_0) w + (p_1 - p_0) u + + [ q_0 + (q_1 - q_0) u ] v + \f] + + This is a system of three equations in three unknowns. It is non-linear, + courtesy of the \f$u v\f$ term at the end. It can be reduced to a single + quadratic in any of the three variables. We choose to solve for \f$u\f$ by + taking the dot product of the above equation with the following vector: + \f[ + (l_1 - l_0) \times [ q_0 + (q_1 - q_0) u ] + \f] + + The sign of the intersection (i.e., whether the line crosses from below the + surface to above or vice versa) can be determined by taking the dot product + of the line vector with the surface normal at the intersection. The surface + normal is as follows: + \f[ + n(u, v) = (p_1 - p_0) \times [q_0 + (q_1 - q_0) u] + (q_1 \times q_0) v + \f] + +SourceFiles + sweptFaceAreaWeightAMI.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sweptFaceAreaWeightAMI_H +#define sweptFaceAreaWeightAMI_H + +#include "faceAreaWeightAMI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class sweptFaceAreaWeightAMI Declaration +\*---------------------------------------------------------------------------*/ + +template +class sweptFaceAreaWeightAMI +: + public faceAreaWeightAMI +{ +private: + + // Private static data + + //- Minimum length of a cut as a ratio of the overall projected edge + // length + static const scalar minCutRatio_; + + //- Maximum allowable dot product between a source point normal and a + // target triangle + static const scalar maxDot_; + + + // Private classes + + //- A fixed list of tris which simulates a dynamic list by incrementing + // a counter whenever its append method is called. This is used as an + // optimisation so that the tri cutting does not allocate memory. + template + class cutTriList + : + public FixedList, Size> + { + private: + + //- The number of stored elements + label n_; + + + public: + + //- Construct null + cutTriList() + : + n_(0) + {} + + //- Clear the array + void clear() + { + n_ = 0; + } + + //- Get the current size + label size() const + { + return n_; + } + + //- Add a new tet to the end of the array + void append(const FixedList& t) + { + this->operator[](n_) = t; + ++ n_; + } + }; + + + // Private member functions + + // Debugging + + //- Write a VTK file of cut triangles + template + void writeCutTrisVTK + ( + const cutTriList& tris, + const word& name + ) const; + + //- Write an OBJ file of a face + void writeFaceOBJ + ( + const face& f, + const pointField& ps, + const string& name + ) const; + + //- Write an OBJ file of the source projection + void writeProjectionOBJ + ( + const label srcN, + const FixedList& srcTri, + const FixedList& srcPrj + ) const; + + + //- Convert the source tris and normals to a projection. Most of the + // time this does nothing, but if some of the normals point in the + // reverse direction the projection will be reduced to span only the + // region in which the projection points forward throught the target + // plane. Returns the number of edges in the projection (0, 3 or 4). + label getSourceProjection + ( + FixedList& srcTri, + FixedList& srcNrm, + const FixedList& tgtTri + ) const; + + //- Get the cutting plane, for an edge of the source projection. + plane getCutPlane + ( + const point& p0, + const point& p1, + const vector& n0, + const vector& n1, + const FixedList& tgtTri + ) const; + + //- The minimum weight below which connections are discarded + virtual scalar minWeight() const; + + //- The maximum edge angle that the walk will cross + virtual scalar maxWalkAngle() const; + + + +protected: + + // Protected Member Functions + + // Evaluation + + //- Area of intersection between source and target faces + virtual scalar interArea + ( + const label srcFacei, + const label tgtFacei + ) const; + + +public: + + //- Runtime type information + TypeName("sweptFaceAreaWeightAMI"); + + + // Constructors + + //- Use parent constructors + using faceAreaWeightAMI::faceAreaWeightAMI; + + + //- Destructor + virtual ~sweptFaceAreaWeightAMI(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "sweptFaceAreaWeightAMI.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C index 319f837593..2b93aa78fc 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIPatchToPatchInterpolation.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,6 +29,8 @@ License #include "mapNearestAMI.H" #include "faceAreaWeightAMI.H" #include "partialFaceAreaWeightAMI.H" +#include "sweptFaceAreaWeightAMI.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -40,6 +42,7 @@ namespace Foam makeAMIMethodType(AMIPatchToPatchInterpolation, mapNearestAMI); makeAMIMethodType(AMIPatchToPatchInterpolation, faceAreaWeightAMI); makeAMIMethodType(AMIPatchToPatchInterpolation, partialFaceAreaWeightAMI); + makeAMIMethodType(AMIPatchToPatchInterpolation, sweptFaceAreaWeightAMI); }