stitchMesh: Replacement utility based on the new patchIntersection algorithm

The mergePatchPairs functionality in blockMesh also now uses patchIntersection.

The new mergePatchPairs and patchIntersection replaces the old, fragile and
practically unusable polyTopoChanger::slidingInterface functionality the removal
of which has allowed the deletion of a lot of other ancient and otherwise unused
clutter including polyTopoChanger, polyMeshModifier, polyTopoChange::setAction
and associated addObject/*, modifyObject/* and removeObject/*.  This
rationalisation paves the way for the completion of the update of zone handling
allowing mesh points, faces and cells to exist in multiple zones which is
currently not supported with mesh topology change.

Application
    stitchMesh

Description
    Utility to stitch or conform pairs of patches,
    converting the patch faces either into internal faces
    or conformal faces or another patch.

Usage
    \b stitchMesh (\<list of patch pairs\>)

    E.g. to stitch patches \c top1 to \c top2 and \c bottom1 to \c bottom2
        stitchMesh "((top1 top2) (bottom1 bottom2))"

    Options:
      - \par -overwrite \n
        Replace the old mesh with the new one, rather than writing the new one
        into a separate time directory

      - \par -region \<name\>
        Specify an alternative mesh region.

      - \par -fields
        Update vol and point fields

      - \par -tol
        Merge tolerance relative to local edge length (default 1e-4)

See also
    Foam::mergePatchPairs
This commit is contained in:
Henry Weller
2023-12-25 13:32:39 +00:00
parent 602347f6c1
commit 69da8f3d7b
57 changed files with 7729 additions and 9852 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,24 +48,16 @@ Usage
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "IOdictionary.H"
#include "IOPtrList.H"
#include "systemDict.H"
#include "blockMesh.H"
#include "attachPolyTopoChanger.H"
#include "mergePatchPairs.H"
#include "polyTopoChange.H"
#include "emptyPolyPatch.H"
#include "cyclicPolyPatch.H"
#include "argList.H"
#include "OSspecific.H"
#include "OFstream.H"
#include "Pair.H"
#include "slidingInterface.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -236,22 +228,23 @@ int main(int argc, char *argv[])
// Read in a list of dictionaries for the merge patch pairs
if (meshDict.found("mergePatchPairs"))
{
List<Pair<word>> mergePatchPairs
List<Pair<word>> patchPairNames
(
meshDict.lookup("mergePatchPairs")
);
#include "mergePatchPairs.H"
if (patchPairNames.size())
{
const word oldInstance = mesh.pointsInstance();
mergePatchPairs(mesh, patchPairNames, 1e-4);
mesh.setInstance(oldInstance);
}
}
else
{
Info<< nl << "There are no merge patch pairs edges" << endl;
Info<< nl << "No patch pairs to merge" << endl;
}
// Set any cellZones (note: cell labelling unaffected by above
// mergePatchPairs)
label nZones = blocks.numZonedBlocks();
if (nZones > 0)

View File

@ -1,119 +0,0 @@
if (mergePatchPairs.size())
{
Info<< "Creating merge patch pairs" << nl << endl;
// Create and add point and face zones and mesh modifiers
List<pointZone*> pz(mergePatchPairs.size());
List<faceZone*> fz(3*mergePatchPairs.size());
List<cellZone*> cz(0);
forAll(mergePatchPairs, pairI)
{
const word mergeName
(
mergePatchPairs[pairI].first()
+ mergePatchPairs[pairI].second()
+ name(pairI)
);
pz[pairI] = new pointZone
(
mergeName + "CutPointZone",
labelList(0),
0,
mesh.pointZones()
);
// Master patch
const word masterPatchName(mergePatchPairs[pairI].first());
const polyPatch& masterPatch =
mesh.boundaryMesh()[masterPatchName];
labelList isf(masterPatch.size());
forAll(isf, i)
{
isf[i] = masterPatch.start() + i;
}
fz[3*pairI] = new faceZone
(
mergeName + "MasterZone",
isf,
boolList(masterPatch.size(), false),
0,
mesh.faceZones()
);
// Slave patch
const word slavePatchName(mergePatchPairs[pairI].second());
const polyPatch& slavePatch =
mesh.boundaryMesh()[slavePatchName];
labelList osf(slavePatch.size());
forAll(osf, i)
{
osf[i] = slavePatch.start() + i;
}
fz[3*pairI + 1] = new faceZone
(
mergeName + "SlaveZone",
osf,
boolList(slavePatch.size(), false),
1,
mesh.faceZones()
);
// Add empty zone for cut faces
fz[3*pairI + 2] = new faceZone
(
mergeName + "CutFaceZone",
labelList(0),
boolList(0, false),
2,
mesh.faceZones()
);
} // end of all merge pairs
Info<< "Adding point and face zones" << endl;
mesh.addZones(pz, fz, cz);
Info<< "Creating attachPolyTopoChanger" << endl;
attachPolyTopoChanger polyMeshAttacher(mesh);
polyMeshAttacher.setSize(mergePatchPairs.size());
forAll(mergePatchPairs, pairI)
{
const word mergeName
(
mergePatchPairs[pairI].first()
+ mergePatchPairs[pairI].second()
+ name(pairI)
);
// Add the sliding interface mesh modifier
polyMeshAttacher.set
(
pairI,
new slidingInterface
(
"couple" + name(pairI),
mesh,
mergeName + "MasterZone",
mergeName + "SlaveZone",
mergeName + "CutPointZone",
mergeName + "CutFaceZone",
mergePatchPairs[pairI].first(),
mergePatchPairs[pairI].second(),
slidingInterface::INTEGRAL, // always integral
false,
intersection::algorithm::visible
)
);
}
polyMeshAttacher.attach(true);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,185 +25,55 @@ Application
stitchMesh
Description
'Stitches' a mesh.
Utility to stitch or conform pairs of patches,
converting the patch faces either into internal faces
or conformal faces or another patch.
Takes a mesh and two patches and merges the faces on the two patches
(if geometrically possible) so the faces become internal.
Usage
\b stitchMesh (\<list of patch pairs\>)
Can do
- 'perfect' match: faces and points on patches align exactly. Order might
be different though.
- 'integral' match: where the surfaces on both patches exactly
match but the individual faces not
- 'partial' match: where the non-overlapping part of the surface remains
in the respective patch.
E.g. to stitch patches \c top1 to \c top2 and \c bottom1 to \c bottom2
stitchMesh "((top1 top2) (bottom1 bottom2))"
Note : Is just a front-end to perfectInterface/slidingInterface.
Options:
- \par -overwrite \n
Replace the old mesh with the new one, rather than writing the new one
into a separate time directory
Comparable to running a meshModifier of the form
(if masterPatch is called "M" and slavePatch "S"):
- \par -region \<name\>
Specify an alternative mesh region.
\verbatim
couple
{
type slidingInterface;
masterFaceZoneName MSMasterZone
slaveFaceZoneName MSSlaveZone
cutPointZoneName MSCutPointZone
cutFaceZoneName MSCutFaceZone
masterPatchName M;
slavePatchName S;
typeOfMatch partial or integral
}
\endverbatim
- \par -fields
Update vol and point fields
- \par -tol
Merge tolerance relative to local edge length (default 1e-4)
See also
Foam::mergePatchPairs
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "polyTopoChanger.H"
#include "polyTopoChangeMap.H"
#include "slidingInterface.H"
#include "perfectInterface.H"
#include "fvMesh.H"
#include "timeSelector.H"
#include "ReadFields.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "pointFields.H"
#include "mergePatchPairs.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
label addPointZone(const polyMesh& mesh, const word& name)
{
label zoneID = mesh.pointZones().findIndex(name);
if (zoneID != -1)
{
Info<< "Reusing existing pointZone "
<< mesh.pointZones()[zoneID].name()
<< " at index " << zoneID << endl;
}
else
{
meshPointZones& pointZones = const_cast<polyMesh&>(mesh).pointZones();
zoneID = pointZones.size();
Info<< "Adding pointZone " << name << " at index " << zoneID << endl;
pointZones.setSize(zoneID+1);
pointZones.set
(
zoneID,
new pointZone
(
name,
labelList(0),
zoneID,
pointZones
)
);
}
return zoneID;
}
label addFaceZone(const polyMesh& mesh, const word& name)
{
label zoneID = mesh.faceZones().findIndex(name);
if (zoneID != -1)
{
Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
<< " at index " << zoneID << endl;
}
else
{
meshFaceZones& faceZones = const_cast<polyMesh&>(mesh).faceZones();
zoneID = faceZones.size();
Info<< "Adding faceZone " << name << " at index " << zoneID << endl;
faceZones.setSize(zoneID+1);
faceZones.set
(
zoneID,
new faceZone
(
name,
labelList(0),
boolList(),
zoneID,
faceZones
)
);
}
return zoneID;
}
label addCellZone(const polyMesh& mesh, const word& name)
{
label zoneID = mesh.cellZones().findIndex(name);
if (zoneID != -1)
{
Info<< "Reusing existing cellZone " << mesh.cellZones()[zoneID].name()
<< " at index " << zoneID << endl;
}
else
{
meshCellZones& cellZones = const_cast<polyMesh&>(mesh).cellZones();
zoneID = cellZones.size();
Info<< "Adding cellZone " << name << " at index " << zoneID << endl;
cellZones.setSize(zoneID+1);
cellZones.set
(
zoneID,
new cellZone
(
name,
labelList(0),
zoneID,
cellZones
)
);
}
return zoneID;
}
// Checks whether patch present
void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
{
const label patchi = bMesh.findIndex(name);
if (patchi == -1)
{
FatalErrorInFunction
<< "Cannot find patch " << name << endl
<< "It should be present and of non-zero size" << endl
<< "Valid patches are " << bMesh.names()
<< exit(FatalError);
}
if (bMesh[patchi].empty())
{
FatalErrorInFunction
<< "Patch " << name << " is present but zero size"
<< exit(FatalError);
}
}
int main(int argc, char *argv[])
{
timeSelector::addOptions(true, false);
argList::addNote
(
"Merge the faces on the specified patches (if geometrically possible)\n"
"so the faces become internal.\n"
"Integral matching is used when the options -partial and -perfect are "
"omitted.\n"
"Stitch the faces on the specified patch pairs\n"
"converting the overlapping patch faces into internal faces.\n"
);
argList::noParallel();
@ -211,276 +81,61 @@ int main(int argc, char *argv[])
#include "addRegionOption.H"
argList::addBoolOption
(
"noFields",
"do not update fields"
);
argList::validArgs.append("masterPatch");
argList::validArgs.append("slavePatch");
argList::addBoolOption
(
"partial",
"couple partially overlapping patches (optional)"
);
argList::addBoolOption
(
"perfect",
"couple perfectly aligned patches (optional)"
"fields",
"update fields"
);
argList::addOption
(
"toleranceDict",
"file",
"dictionary file with tolerances"
"tol",
"scalar",
"merge tolerance relative to local edge length (default 1e-4)"
);
argList::validArgs.append("patchPairs");
#include "setRootCase.H"
#include "createTimeNoFunctionObjects.H"
// Select time if specified
timeSelector::selectIfPresent(runTime, args);
#include "createNamedMesh.H"
const scalar snapTol = args.optionLookupOrDefault("tol", 1e-4);
const bool overwrite = args.optionFound("overwrite");
const bool fields = args.optionFound("fields");
const List<Pair<word>> patchPairNames((IStringStream(args[1])()));
const word oldInstance = mesh.pointsInstance();
const word masterPatchName = args[1];
const word slavePatchName = args[2];
const bool partialCover = args.optionFound("partial");
const bool perfectCover = args.optionFound("perfect");
const bool overwrite = args.optionFound("overwrite");
const bool fields = !args.optionFound("noFields");
if (partialCover && perfectCover)
{
FatalErrorInFunction
<< "Cannot supply both partial and perfect." << endl
<< "Use perfect match option if the patches perfectly align"
<< " (both vertex positions and face centres)" << endl
<< exit(FatalError);
}
const word mergePatchName(masterPatchName + slavePatchName);
const word cutZoneName(mergePatchName + "CutFaceZone");
slidingInterface::typeOfMatch tom = slidingInterface::INTEGRAL;
if (partialCover)
{
Info<< "Coupling partially overlapping patches "
<< masterPatchName << " and " << slavePatchName << nl
<< "Resulting internal faces will be in faceZone " << cutZoneName
<< nl
<< "Any uncovered faces will remain in their patch"
<< endl;
tom = slidingInterface::PARTIAL;
}
else if (perfectCover)
{
Info<< "Coupling perfectly aligned patches "
<< masterPatchName << " and " << slavePatchName << nl
<< "Resulting (internal) faces will be in faceZone " << cutZoneName
<< nl << nl
<< "Note: both patches need to align perfectly." << nl
<< "Both the vertex"
<< " positions and the face centres need to align to within" << nl
<< "a tolerance given by the minimum edge length on the patch"
<< endl;
}
else
{
Info<< "Coupling patches " << masterPatchName << " and "
<< slavePatchName << nl
<< "Resulting (internal) faces will be in faceZone " << cutZoneName
<< nl << nl
<< "Note: the overall area covered by both patches should be"
<< " identical (\"integral\" interface)." << endl
<< "If this is not the case use the -partial option" << nl << endl;
}
// set up the tolerances for the sliding mesh
dictionary slidingTolerances;
if (args.options().found("toleranceDict"))
{
IOdictionary toleranceFile
(
IOobject
(
args.options()["toleranceDict"],
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
slidingTolerances += toleranceFile;
}
// Check for non-empty master and slave patches
checkPatch(mesh.boundaryMesh(), masterPatchName);
checkPatch(mesh.boundaryMesh(), slavePatchName);
// Create and add face zones and mesh modifiers
// Master patch
const polyPatch& masterPatch = mesh.boundaryMesh()[masterPatchName];
// Make list of masterPatch faces
labelList isf(masterPatch.size());
forAll(isf, i)
{
isf[i] = masterPatch.start() + i;
}
polyTopoChanger stitcher(mesh);
stitcher.setSize(1);
mesh.pointZones().clearAddressing();
mesh.faceZones().clearAddressing();
mesh.cellZones().clearAddressing();
if (perfectCover)
{
// Add empty zone for resulting internal faces
label cutZoneID = addFaceZone(mesh, cutZoneName);
mesh.faceZones()[cutZoneID].resetAddressing
(
isf,
boolList(masterPatch.size(), false)
);
// Add the perfect interface mesh modifier
stitcher.set
(
0,
new perfectInterface
(
"couple",
mesh,
cutZoneName,
masterPatchName,
slavePatchName
)
);
}
else
{
label pointZoneID = addPointZone(mesh, mergePatchName + "CutPointZone");
mesh.pointZones()[pointZoneID] = labelList(0);
label masterZoneID = addFaceZone(mesh, mergePatchName + "MasterZone");
mesh.faceZones()[masterZoneID].resetAddressing
(
isf,
boolList(masterPatch.size(), false)
);
// Slave patch
const polyPatch& slavePatch = mesh.boundaryMesh()[slavePatchName];
labelList osf(slavePatch.size());
forAll(osf, i)
{
osf[i] = slavePatch.start() + i;
}
label slaveZoneID = addFaceZone(mesh, mergePatchName + "SlaveZone");
mesh.faceZones()[slaveZoneID].resetAddressing
(
osf,
boolList(slavePatch.size(), false)
);
// Add empty zone for cut faces
label cutZoneID = addFaceZone(mesh, cutZoneName);
mesh.faceZones()[cutZoneID].resetAddressing
(
labelList(0),
boolList(0, false)
);
// Add the sliding interface mesh modifier
stitcher.set
(
0,
new slidingInterface
(
"couple",
mesh,
mergePatchName + "MasterZone",
mergePatchName + "SlaveZone",
mergePatchName + "CutPointZone",
cutZoneName,
masterPatchName,
slavePatchName,
tom, // integral or partial
true // couple/decouple mode
)
);
static_cast<slidingInterface&>(stitcher[0]).setTolerances
(
slidingTolerances,
true
);
}
// Search for list of objects for this time
IOobjectList objects(mesh, runTime.name());
IOobjectList objects(mesh, runTime.timeName());
if (fields) Info<< "Reading geometric fields" << nl << endl;
#include "readVolFields.H"
// #include "readSurfaceFields.H"
#include "readPointFields.H"
Info<< endl;
if (!overwrite)
{
runTime++;
}
// Execute all polyMeshModifiers
autoPtr<polyTopoChangeMap> map = stitcher.changeMesh(true);
mesh.setPoints(map->preMotionPoints());
// Stitch the patch-pairs
mergePatchPairs(mesh, patchPairNames, snapTol);
// Write mesh
if (overwrite)
{
mesh.setInstance(oldInstance);
}
Info<< nl << "Writing polyMesh to time " << runTime.name() << endl;
IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
Info<< "Writing mesh to " << mesh.facesInstance() << endl;
// Bypass runTime write (since only writes at writeTime)
if
(
!runTime.objectRegistry::writeObject
(
runTime.writeFormat(),
IOstream::currentVersion,
runTime.writeCompression(),
true
)
)
{
FatalErrorInFunction
<< "Failed writing polyMesh."
<< exit(FatalError);
}
mesh.faceZones().write();
mesh.pointZones().write();
mesh.cellZones().write();
// Write fields
runTime.write();
mesh.write();
Info<< "End\n" << endl;

View File

@ -36,7 +36,6 @@ Description
#include "polyMesh.H"
#include "faceSet.H"
#include "polyTopoChange.H"
#include "polyModifyFace.H"
#include "globalMeshData.H"
using namespace Foam;

View File

@ -1,44 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object toleranceDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// For complete details on what this dictionary file should provide, see
// $FOAM_SRC/src/polyTopoChange/slidingInterface/slidingInterface.C
// method: Foam::slidingInterface::setTolerances
//- Point merge tolerance
pointMergeTol 0.05;
//- Edge merge tolerance
edgeMergeTol 0.01;
//- Estimated number of faces an edge goes through
nFacesPerSlaveEdge 5;
//- Edge-face interaction escape limit
edgeFaceEscapeLimit 10;
//- Integral match point adjustment tolerance
integralAdjTol 0.05;
//- Edge intersection master catch fraction
edgeMasterCatchFraction 0.4;
//- Edge intersection co-planar tolerance
edgeCoPlanarTol 0.8;
//- Edge end cut-off tolerance
edgeEndCutoffTol 0.0001;
// ************************************************************************* //

View File

@ -355,7 +355,7 @@ _checkMesh_ ()
local line=${COMP_LINE}
local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
opts="-allGeometry -allTopology -case -constant -doc -fileHandler -help -hostRoots -latestTime -libs -meshQuality -noFunctionObjects -noTopology -noZero -parallel -region -roots -srcDoc -time -writeSets"
opts="-allGeometry -allTopology -case -constant -doc -fileHandler -help -hostRoots -latestTime -libs -meshQuality -noFunctionObjects -nonOrthThreshold -noTopology -noZero -parallel -region -roots -setFormat -skewThreshold -srcDoc -surfaceFormat -time -writeSets -writeSurfaces"
for o in $used ; do opts="${opts/$o/}" ; done
extra=""
@ -367,7 +367,7 @@ _checkMesh_ ()
opts="uncollated collated masterUncollated" ; extra="" ;;
-time)
opts="$(foamListTimes -withZero 2> /dev/null)" ; extra="" ;;
-hostRoots|-libs|-region|-roots|-writeSets)
-hostRoots|-libs|-nonOrthThreshold|-region|-roots|-setFormat|-skewThreshold|-surfaceFormat)
opts="" ; extra="" ;;
*) ;;
esac
@ -1989,31 +1989,6 @@ _mixtureAdiabaticFlameT_ ()
}
complete -o filenames -o nospace -F _mixtureAdiabaticFlameT_ mixtureAdiabaticFlameT
_modifyMesh_ ()
{
local cur="${COMP_WORDS[COMP_CWORD]}"
local prev="${COMP_WORDS[COMP_CWORD-1]}"
local line=${COMP_LINE}
local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -overwrite -parallel -roots -srcDoc"
for o in $used ; do opts="${opts/$o/}" ; done
extra=""
[ "$COMP_CWORD" = 1 ] || \
case "$prev" in
-case)
opts="" ; extra="-d" ;;
-fileHandler)
opts="uncollated collated masterUncollated" ; extra="" ;;
-hostRoots|-libs|-roots)
opts="" ; extra="" ;;
*) ;;
esac
COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
}
complete -o filenames -o nospace -F _modifyMesh_ modifyMesh
_mshToFoam_ ()
{
local cur="${COMP_WORDS[COMP_CWORD]}"
@ -2994,31 +2969,6 @@ _splitCells_ ()
}
complete -o filenames -o nospace -F _splitCells_ splitCells
_splitMesh_ ()
{
local cur="${COMP_WORDS[COMP_CWORD]}"
local prev="${COMP_WORDS[COMP_CWORD-1]}"
local line=${COMP_LINE}
local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
opts="-case -doc -fileHandler -help -libs -noFunctionObjects -overwrite -srcDoc"
for o in $used ; do opts="${opts/$o/}" ; done
extra=""
[ "$COMP_CWORD" = 1 ] || \
case "$prev" in
-case)
opts="" ; extra="-d" ;;
-fileHandler)
opts="uncollated collated masterUncollated" ; extra="" ;;
-libs)
opts="" ; extra="" ;;
*) ;;
esac
COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
}
complete -o filenames -o nospace -F _splitMesh_ splitMesh
_splitMeshRegions_ ()
{
local cur="${COMP_WORDS[COMP_CWORD]}"
@ -3132,7 +3082,7 @@ _stitchMesh_ ()
local line=${COMP_LINE}
local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
opts="-case -doc -fileHandler -help -libs -noFields -noFunctionObjects -overwrite -partial -perfect -region -srcDoc -toleranceDict"
opts="-case -constant -doc -fields -fileHandler -help -latestTime -libs -noFunctionObjects -noZero -overwrite -region -srcDoc -time -tol"
for o in $used ; do opts="${opts/$o/}" ; done
extra=""
@ -3140,11 +3090,11 @@ _stitchMesh_ ()
case "$prev" in
-case)
opts="" ; extra="-d" ;;
-toleranceDict)
opts="" ; extra="-d -f" ;;
-fileHandler)
opts="uncollated collated masterUncollated" ; extra="" ;;
-libs|-region)
-time)
opts="$(foamListTimes -withZero 2> /dev/null)" ; extra="" ;;
-libs|-region|-tol)
opts="" ; extra="" ;;
*) ;;
esac

View File

@ -275,4 +275,9 @@ nonConformal/boundary/nonConformalBoundary.C
pointDist/pointEdgeDist.C
pointDist/pointDist.C
patchIntersection/star/star.C
patchIntersection/patchIntersection.C
patchIntersection/primitivePatchIntersection.C
patchIntersection/polyPatchIntersection.C
LIB = $(FOAM_LIBBIN)/libmeshTools

View File

@ -0,0 +1,510 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "FacePatchIntersection.H"
#include "cpuTime.H"
#include "primitiveTriPatch.H"
#include "polygonTriangulate.H"
#include "TriPatchIntersection.H"
#include "vtkWritePolyData.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class T, class ListType>
void inplaceRemove(ListType& lst, const T& value)
{
if (lst.size())
{
inplaceRemove(lst, value, lst[0]);
}
}
template<class T, class NotT, class ListType>
void inplaceRemove(ListType& lst, const T& value, const NotT&)
{
forAll(lst, i)
{
if (lst[i].size())
{
inplaceRemove(lst[i], value, lst[i][0]);
}
}
}
template<class T, class ListType>
void inplaceRemove(ListType& lst, const T& value, const T&)
{
label iNew = 0;
forAll(lst, iOld)
{
if (lst[iOld] != value)
{
lst[iNew] = lst[iOld];
++ iNew;
}
}
lst.resize(iNew);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SrcPatchType, class TgtPatchType>
Foam::FacePatchIntersection<SrcPatchType, TgtPatchType>::FacePatchIntersection
(
const SrcPatchType& srcPatch,
const TgtPatchType& tgtPatch,
const scalar snapTol
)
:
FacePatchIntersection(srcPatch, srcPatch.pointNormals(), tgtPatch, snapTol)
{}
template<class SrcPatchType, class TgtPatchType>
Foam::FacePatchIntersection<SrcPatchType, TgtPatchType>::FacePatchIntersection
(
const SrcPatchType& srcPatch,
const vectorField& srcPointNormals,
const TgtPatchType& tgtPatch,
const scalar snapTol
)
:
PatchIntersection<SrcPatchType, TgtPatchType>(srcPatch, tgtPatch)
{
cpuTime time;
Info<< type() << ": Intersecting "
<< this->srcPatch_.size() << " source faces and "
<< this->tgtPatch_.size() << " target faces" << incrIndent << endl;
if (this->debug)
{
Info<< indent << "Writing patches" << incrIndent << endl;
const fileName srcFileName = type() + "_srcPatch.vtk";
Info<< indent << "Writing patch to " << srcFileName << endl;
vtkWritePolyData::write
(
srcFileName,
"source",
false,
this->srcPatch_.localPoints(),
labelList(),
labelListList(),
this->srcPatch_.localFaces()
);
const fileName tgtFileName = type() + "_tgtPatch.vtk";
Info<< indent << "Writing patch to " << tgtFileName << endl;
vtkWritePolyData::write
(
tgtFileName,
"target",
false,
this->tgtPatch_.localPoints(),
labelList(),
labelListList(),
this->tgtPatch_.localFaces()
);
Info<< decrIndent;
}
// Step 1: Triangulate the source and target patches and create addressing
DynamicList<triFace> srcTris, tgtTris;
DynamicList<label> srcTriSrcFaces, tgtTriTgtFaces;
DynamicList<FixedList<label, 3>> srcTriSrcEdges, tgtTriTgtEdges;
polygonTriangulate triEngine;
forAll(this->srcPatch_, srcFacei)
{
const face& srcFacePoints = this->srcPatch_[srcFacei];
const labelList& srcFaceEdges = this->srcPatch_.faceEdges()[srcFacei];
triEngine.triangulate
(
UIndirectList<point>(this->srcPatch_.points(), srcFacePoints)
);
srcTris.append(triEngine.triPoints(srcFacePoints));
srcTriSrcFaces.resize(srcTris.size(), srcFacei);
srcTriSrcEdges.append(triEngine.triEdges(srcFaceEdges));
}
forAll(this->tgtPatch_, tgtFacei)
{
const face tgtFacePoints = this->tgtPatch_[tgtFacei].reverseFace();
const labelList tgtFaceEdges =
reverseList(this->tgtPatch_.faceEdges()[tgtFacei]);
triEngine.triangulate
(
UIndirectList<point>(this->tgtPatch_.points(), tgtFacePoints)
);
tgtTris.append(triEngine.triPoints(tgtFacePoints));
tgtTriTgtFaces.resize(tgtTris.size(), tgtFacei);
tgtTriTgtEdges.append(triEngine.triEdges(tgtFaceEdges));
const label n = triEngine.triPoints().size();
for (label i = 0; i < n; ++ i)
{
tgtTris[tgtTris.size() - n + i] =
tgtTris[tgtTris.size() - n + i].reverseFace();
tgtTriTgtEdges[tgtTris.size() - n + i] =
reverseList(tgtTriTgtEdges[tgtTris.size() - n + i]);
}
}
const primitiveTriPatch srcTriPatch
(
SubList<triFace>(srcTris, srcTris.size()),
this->srcPatch_.points()
);
const primitiveTriPatch tgtTriPatch
(
SubList<triFace>(tgtTris, tgtTris.size()),
this->tgtPatch_.points()
);
labelList srcTriPatchPointSrcPatchPoints(srcTriPatch.nPoints());
labelList tgtTriPatchPointTgtPatchPoints(tgtTriPatch.nPoints());
{
labelList map
(
max(max(srcPatch.meshPoints()), max(tgtPatch.meshPoints())) + 1,
-1
);
UIndirectList<label>(map, srcPatch.meshPoints()) =
identityMap(srcPatch.nPoints());
srcTriPatchPointSrcPatchPoints =
renumber(map, srcTriPatch.meshPoints());
UIndirectList<label>(map, srcPatch.meshPoints()) = -1;
UIndirectList<label>(map, tgtPatch.meshPoints()) =
identityMap(tgtPatch.nPoints());
tgtTriPatchPointTgtPatchPoints =
renumber(map, tgtTriPatch.meshPoints());
UIndirectList<label>(map, tgtPatch.meshPoints()) = -1;
}
labelList srcTriPatchEdgeSrcPatchEdges(srcTriPatch.nEdges());
labelList tgtTriPatchEdgeTgtPatchEdges(tgtTriPatch.nEdges());
forAll(srcTriPatch, srcTrii)
{
forAll(srcTriPatch[srcTrii], srcTriEdgei)
{
const label srcTriPatchEdgei =
srcTriPatch.faceEdges()[srcTrii][srcTriEdgei];
const label srcPatchEdgei = srcTriSrcEdges[srcTrii][srcTriEdgei];
srcTriPatchEdgeSrcPatchEdges[srcTriPatchEdgei] = srcPatchEdgei;
}
}
forAll(tgtTriPatch, tgtTrii)
{
forAll(tgtTriPatch[tgtTrii], tgtTriEdgei)
{
const label tgtTriPatchEdgei =
tgtTriPatch.faceEdges()[tgtTrii][tgtTriEdgei];
const label tgtPatchEdgei = tgtTriTgtEdges[tgtTrii][tgtTriEdgei];
tgtTriPatchEdgeTgtPatchEdges[tgtTriPatchEdgei] = tgtPatchEdgei;
}
}
// Step 2: Do the tri-patch intersection
const TriPatchIntersection<primitiveTriPatch, primitiveTriPatch> tpi
(
srcTriPatch,
vectorField(srcPointNormals, srcTriPatchPointSrcPatchPoints),
tgtTriPatch,
snapTol
);
// Step 3: Map points and point associations back into the class data
{
this->points_ = tpi.points();
// Point-point connectivity
this->srcPointPoints_ =
reorder(srcTriPatchPointSrcPatchPoints, tpi.srcPointPoints());
this->tgtPointPoints_ =
reorder(tgtTriPatchPointTgtPatchPoints, tpi.tgtPointPoints());
this->pointSrcPoints_ = tpi.pointSrcPoints();
inplaceRenumber(srcTriPatchPointSrcPatchPoints, this->pointSrcPoints_);
this->pointTgtPoints_ = tpi.pointTgtPoints();
inplaceRenumber(tgtTriPatchPointTgtPatchPoints, this->pointTgtPoints_);
// Point-edge connectivity
this->srcEdgePoints_ =
reorder(srcTriPatchEdgeSrcPatchEdges, tpi.srcEdgePoints());
this->srcEdgePoints_.resize(this->srcPatch_.nEdges());
this->tgtEdgePoints_ =
reorder(tgtTriPatchEdgeTgtPatchEdges, tpi.tgtEdgePoints());
this->tgtEdgePoints_.resize(this->tgtPatch_.nEdges());
this->pointSrcEdges_ = tpi.pointSrcEdges();
inplaceRenumber(srcTriPatchEdgeSrcPatchEdges, this->pointSrcEdges_);
this->pointTgtEdges_ = tpi.pointTgtEdges();
inplaceRenumber(tgtTriPatchEdgeTgtPatchEdges, this->pointTgtEdges_);
// Point-tri connectivity
this->pointSrcFaces_ = tpi.pointSrcFaces();
inplaceRenumber(srcTriSrcFaces, this->pointSrcFaces_);
this->pointTgtFaces_ = tpi.pointTgtFaces();
inplaceRenumber(tgtTriTgtFaces, this->pointTgtFaces_);
// Point-face-diagonal connectivity
auto setFaceDiagonalPoints = []
(
const List<DynamicList<label>>& edgePoints,
const labelListList& edgeTris,
const labelUList& triFaces,
DynamicList<label>& pointFaces
)
{
forAll(edgePoints, edgei)
{
const labelList& triis = edgeTris[edgei];
if (triis.size() != 2)
{
continue;
}
const label facei0 = triFaces[triis[0]];
const label facei1 = triFaces[triis[1]];
if (facei0 != facei1)
{
continue;
}
const label facei = facei0;
for
(
label edgePointi = 1;
edgePointi < edgePoints[edgei].size() - 1;
++ edgePointi
)
{
const label pointi = edgePoints[edgei][edgePointi];
pointFaces[pointi] = facei;
}
}
};
setFaceDiagonalPoints
(
tpi.srcEdgePoints(),
srcTriPatch.edgeFaces(),
srcTriSrcFaces,
this->pointSrcFaces_
);
setFaceDiagonalPoints
(
tpi.tgtEdgePoints(),
tgtTriPatch.edgeFaces(),
tgtTriTgtFaces,
this->pointTgtFaces_
);
}
// Step 4: Construct faces and face associations by merging the faces that
// result from the tri-patch intersection
{
// Loop the tri intersection faces, initialising a star at each one,
// and propagating out to include everything connected with the same
// source and target face. This is then part of the face intersection.
boolList tpiFaceCombined(tpi.faces().size(), false);
star tpiStar;
forAll(tpi.faces(), tpiFacei)
{
// If this face has been done then skip to the next
if (tpiFaceCombined[tpiFacei]) continue;
// Get the source and target faces
const label srcFacei =
tpi.faceSrcFaces()[tpiFacei] != -1
? srcTriSrcFaces[tpi.faceSrcFaces()[tpiFacei]]
: -1;
const label tgtFacei =
tpi.faceTgtFaces()[tpiFacei] != -1
? tgtTriTgtFaces[tpi.faceTgtFaces()[tpiFacei]]
: -1;
// Get the relevant set of face-edge and edge-face addressing
const UList<labelList>& tpiFaceEdges = tpi.faceEdges();
const UList<labelPair>& tpiEdgeFaces =
srcFacei != -1 && tgtFacei != -1
? tpi.intersectEdgeFaces()
: tpi.nonIntersectEdgeFaces();
// Fill the star with all connected faces with the same source and
// target face associations
star::context tpiStarContext = tpiStar.populate
(
tpiFacei,
true,
[&](const label tpiEdgei, const label tpiFacei)
{
const label tpiSrcFacei = tpi.faceSrcFaces()[tpiFacei];
const label tpiTgtFacei = tpi.faceTgtFaces()[tpiFacei];
return
(
tpiSrcFacei != -1
? srcTriSrcFaces[tpiSrcFacei] == srcFacei
: srcFacei == -1
)
&& (
tpiTgtFacei != -1
? tgtTriTgtFaces[tpiTgtFacei] == tgtFacei
: tgtFacei == -1
);
},
tpiFaceEdges,
tpiEdgeFaces
);
// Mark everything that has been visited
forAllStarFaces(tpiStar, starFacei, tpiFacei)
{
tpiFaceCombined[tpiFacei] = true;
}
// Create the new face
const label facei = this->faces_.size();
this->faces_.append(face(tpiStar.starEdgeEdges().size()));
if (srcFacei != -1)
{
this->srcFaceFaces_[srcFacei].append(facei);
}
if (tgtFacei != -1)
{
this->tgtFaceFaces_[tgtFacei].append(facei);
}
this->faceSrcFaces_.append(srcFacei);
this->faceTgtFaces_.append(tgtFacei);
forAllStarEdges(tpiStar, i, starEdgei, tpiEdgei)
{
const label tpiEdgeFacei =
tpiEdgeFaces[tpiEdgei][0] == -1
|| tpiStar.faceStarFaces()[tpiEdgeFaces[tpiEdgei][0]] == -1;
const label tpiFacei = tpiEdgeFaces[tpiEdgei][tpiEdgeFacei];
const label tpiFaceEdgei =
findIndex(tpiFaceEdges[tpiFacei], tpiEdgei);
this->faces_.last()[i] =
tpi.faces()[tpiFacei].faceEdge(tpiFaceEdgei).start();
}
}
}
// Step 5: Filter out points that aren't used
{
// Determine how many faces each point is a part of
labelList pointNFaces(this->points_.size(), 0);
forAll(this->faces_, facei)
{
const face& f = this->faces_[facei];
forAll(f, fPointi)
{
++ pointNFaces[f[fPointi]];
}
}
// Remove points that are referenced by two faces or fewer, and which
// do not correspond to original geometry or an intersection
labelList oldPointNewPoints(this->points_.size(), -1);
label pointi = 0;
forAll(this->points_, pointj)
{
if
(
pointNFaces[pointj] > 2
|| (
pointNFaces[pointj] > 0
&& (
this->pointSrcPoints_[pointj] != -1
|| this->pointTgtPoints_[pointj] != -1
|| (
this->pointSrcEdges_[pointj] != -1
&& this->pointTgtEdges_[pointj] != -1
)
)
)
)
{
oldPointNewPoints[pointj] = pointi;
this->points_[pointi] = this->points_[pointj];
this->pointSrcPoints_[pointi] = this->pointSrcPoints_[pointj];
this->pointTgtPoints_[pointi] = this->pointTgtPoints_[pointj];
this->pointSrcEdges_[pointi] = this->pointSrcEdges_[pointj];
this->pointTgtEdges_[pointi] = this->pointTgtEdges_[pointj];
this->pointSrcFaces_[pointi] = this->pointSrcFaces_[pointj];
this->pointTgtFaces_[pointi] = this->pointTgtFaces_[pointj];
++ pointi;
}
}
this->points_.resize(pointi);
this->pointSrcPoints_.resize(pointi);
this->pointTgtPoints_.resize(pointi);
this->pointSrcEdges_.resize(pointi);
this->pointTgtEdges_.resize(pointi);
this->pointSrcFaces_.resize(pointi);
this->pointTgtFaces_.resize(pointi);
// Map
inplaceRenumber(oldPointNewPoints, this->faces_);
inplaceRenumber(oldPointNewPoints, this->srcPointPoints_);
inplaceRenumber(oldPointNewPoints, this->tgtPointPoints_);
inplaceRenumber(oldPointNewPoints, this->srcEdgePoints_);
inplaceRenumber(oldPointNewPoints, this->tgtEdgePoints_);
// Remove deleted points
inplaceRemove(this->faces_, -1);
inplaceRemove(this->srcEdgePoints_, -1);
inplaceRemove(this->tgtEdgePoints_, -1);
}
// Step 6: Reporting
this->report();
Info<< indent << this->faces_.size() << " couplings generated in "
<< time.cpuTimeIncrement() << 's' << endl;
Info<< decrIndent;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SrcPatchType, class TgtPatchType>
Foam::FacePatchIntersection<SrcPatchType, TgtPatchType>::
~FacePatchIntersection()
{}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,103 +22,71 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyTopoChanger
Foam::FacePatchIntersection
Description
List of mesh modifiers
Patch intersection based on polygonal faces. This triangulates the supplied
patches and passes them to TriPatchIntersection. It then re-combines the
triangles that TriPatchIntersection generates.
SourceFiles
polyTopoChanger.C
FacePatchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef polyTopoChanger_H
#define polyTopoChanger_H
#ifndef FacePatchIntersection_H
#define FacePatchIntersection_H
#include "polyMeshModifier.H"
#include "PtrList.H"
#include "PatchIntersection.H"
#include "triFace.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class polyMesh;
class polyTopoChangeMap;
class polyBoundaryMesh;
/*---------------------------------------------------------------------------*\
Class polyTopoChanger Declaration
Class FacePatchIntersection Declaration
\*---------------------------------------------------------------------------*/
class polyTopoChanger
template<class SrcPatchType, class TgtPatchType>
class FacePatchIntersection
:
public PtrList<polyMeshModifier>
public PatchIntersection<SrcPatchType, TgtPatchType>
{
// Private member functions
//- Is topology change required
bool changeTopology() const;
//- Return topology change request
autoPtr<polyTopoChange> topoChangeRequest() const;
//- Force recalculation of locally stored data on topological change
void update(const polyTopoChangeMap& m);
protected:
// Protected member data
//- Reference to mesh
polyMesh& mesh_;
public:
//- Runtime type information
TypeName("polyTopoChanger");
// Runtime type information
virtual word type() const
{
return "face" + patchIntersection::typeName.capitalise();
}
// Constructors
//- Read constructor for given polyMesh
explicit polyTopoChanger(polyMesh&);
//- Disallow default bitwise copy construction
polyTopoChanger(const polyTopoChanger&) = delete;
//- Destructor
virtual ~polyTopoChanger()
{}
// Member Functions
//- Return the mesh reference
const polyMesh& mesh() const
{
return mesh_;
}
autoPtr<polyTopoChangeMap> changeMesh
//- Construct from a source and a target patch
FacePatchIntersection
(
const bool inflate,
const bool syncParallel = true,
const bool orderCells = false,
const bool orderPoints = false
const SrcPatchType& srcPatch,
const TgtPatchType& tgtPatch,
const scalar snapTol
);
//- Construct from a source and a target patch, and specified source
// point normals
FacePatchIntersection
(
const SrcPatchType& srcPatch,
const vectorField& srcPointNormals,
const TgtPatchType& tgtPatch,
const scalar snapTol
);
// Member Operators
//- Disallow default bitwise assignment
void operator=(const polyTopoChanger&) = delete;
//- Destructor
virtual ~FacePatchIntersection();
};
@ -128,6 +96,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "FacePatchIntersection.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,46 +21,41 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Virtual base class for mesh modifiers.
\*---------------------------------------------------------------------------*/
#include "polyMeshModifier.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polyMeshModifier, 0);
}
#include "PatchIntersection.H"
#include "primitivePatch.H"
#include "vtkWritePolyData.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyMeshModifier::polyMeshModifier
template<class SrcPatchType, class TgtPatchType>
Foam::PatchIntersection<SrcPatchType, TgtPatchType>::PatchIntersection
(
const word& name,
const polyMesh& mesh
const SrcPatchType& srcPatch,
const TgtPatchType& tgtPatch
)
:
name_(name),
mesh_(mesh)
patchIntersection
(
srcPatch.nPoints(),
tgtPatch.nPoints(),
srcPatch.nEdges(),
tgtPatch.nEdges(),
srcPatch.size(),
tgtPatch.size()
),
srcPatch_(srcPatch),
tgtPatch_(tgtPatch)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::polyMeshModifier::~polyMeshModifier()
template<class SrcPatchType, class TgtPatchType>
Foam::PatchIntersection<SrcPatchType, TgtPatchType>::~PatchIntersection()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::polyMesh& Foam::polyMeshModifier::mesh() const
{
return mesh_;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,18 +22,25 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyRemoveFace
Foam::PatchIntersection
Description
Class containing data for face removal.
Base class for patch intersections. Provides storage and access to the
intersection points and faces and their relationship between the source and
target patches.
SourceFiles
PatchIntersectionI.H
PatchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef polyRemoveFace_H
#define polyRemoveFace_H
#ifndef PatchIntersection_H
#define PatchIntersection_H
#include "label.H"
#include "topoAction.H"
#include "patchIntersection.H"
#include "pointField.H"
#include "faceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -41,68 +48,50 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyRemoveFace Declaration
Class PatchIntersection Declaration
\*---------------------------------------------------------------------------*/
class polyRemoveFace
template<class SrcPatchType, class TgtPatchType>
class PatchIntersection
:
public topoAction
public patchIntersection
{
// Private Data
protected:
//- Face ID
label faceIndex_;
// Protected Data
//- Reference to the source patch
const SrcPatchType& srcPatch_;
//- Reference to the target patch
const TgtPatchType& tgtPatch_;
//- Merge faceID or -1
label mergeFaceIndex_;
public:
// Static Data Members
//- Runtime type information
TypeName("removeFace");
// Constructors
//- Construct null. Used for constructing lists
polyRemoveFace()
:
faceIndex_(-1),
mergeFaceIndex_(-1)
{}
//- Construct from components
polyRemoveFace(const label faceID, const label mergeFaceID = -1)
:
faceIndex_(faceID),
mergeFaceIndex_(mergeFaceID)
{}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyRemoveFace(*this));
}
//- Construct from a source and a target patch
PatchIntersection
(
const SrcPatchType& srcPatch,
const TgtPatchType& tgtPatch
);
// Default Destructor
//- Destructor
virtual ~PatchIntersection();
// Member Functions
//- Return face ID
label faceIndex() const
{
return faceIndex_;
}
// Access
//- Return merge face ID
label mergeFaceIndex() const
{
return mergeFaceIndex_;
}
//- The source patch
inline const SrcPatchType& srcPatch() const;
//- The target patch
inline const TgtPatchType& tgtPatch() const;
};
@ -112,6 +101,14 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "PatchIntersectionI.H"
#ifdef NoRepository
#include "PatchIntersection.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,27 +23,23 @@ License
\*---------------------------------------------------------------------------*/
#include "topoAction.H"
#include "polyAddPoint.H"
#include "polyAddFace.H"
#include "polyModifyPoint.H"
#include "polyModifyFace.H"
#include "polyRemovePoint.H"
#include "polyRemoveFace.H"
#include "PatchIntersection.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
template<class SrcPatchType, class TgtPatchType>
inline const SrcPatchType&
Foam::PatchIntersection<SrcPatchType, TgtPatchType>::srcPatch() const
{
defineTypeNameAndDebug(topoAction, 0);
return srcPatch_;
}
defineTypeNameAndDebug(polyAddPoint, 0);
defineTypeNameAndDebug(polyModifyPoint, 0);
defineTypeNameAndDebug(polyRemovePoint, 0);
defineTypeNameAndDebug(polyAddFace, 0);
defineTypeNameAndDebug(polyModifyFace, 0);
defineTypeNameAndDebug(polyRemoveFace, 0);
template<class SrcPatchType, class TgtPatchType>
inline const TgtPatchType&
Foam::PatchIntersection<SrcPatchType, TgtPatchType>::tgtPatch() const
{
return tgtPatch_;
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,488 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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::TriPatchIntersection
Description
Patch intersection based on triangular faces. Intersects and combines two
triangulated patches incrementally. The intersected surface is valid at
every stage of the process. Failure to intersect does not produce a
catastrophic error. Rather, it results in regions of the surface remaining
associated with only one of the source or the target patch.
SourceFiles
TriPatchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef TriPatchIntersection_H
#define TriPatchIntersection_H
#include "PatchIntersection.H"
#include "star.H"
#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TriPatchIntersection Declaration
\*---------------------------------------------------------------------------*/
template<class SrcPatchType, class TgtPatchType>
class TriPatchIntersection
:
public PatchIntersection<SrcPatchType, TgtPatchType>
{
private:
// Private Data
// Points
//- The source points
DynamicField<point> srcPoints_;
//- The source point normals
DynamicField<vector> srcPointNormals_;
//- The target points. Reference to base class data.
DynamicField<point>& tgtPoints_;
//- Point-point addressing. Facilitates quick point removal.
// Instead of shuffling up, when points are combined this just
// acts as a redirect. Everything gets resolved on clean.
DynamicList<label> pointPoints_;
// Edges
//- The two triangles adjacent to each edge
DynamicList<labelPair> edgeTris_;
//- ...
DynamicList<labelPair> intersectEdgeFaces_;
//- ...
DynamicList<labelPair> nonIntersectEdgeFaces_;
// Tris
//- The triangles' points
DynamicList<triFace> triPoints_;
//- The triangles' edges
DynamicList<FixedList<label, 3>> triEdges_;
//- The source patch face associated with each triangle. -1 if not
// associated with a source face.
DynamicList<label> triSrcFace_;
//- The target patch face associated with each triangle. -1 if not
// associated with a target face.
DynamicList<label> triTgtFace_;
//- The triangles constructed from each source patch face
List<DynamicList<label>> srcFaceTris_;
//- The triangles constructed from each target patch face
List<DynamicList<label>> tgtFaceTris_;
// Faces
//- ...
DynamicList<labelList> faceEdges_;
// Removal
//- Indices of triangles that have been removed. Prevents shuffling
// up and re-allocation during intersection. Cleared out on clean.
DynamicList<label> removedTris_;
//- Indices of edges that have been removed. Prevents shuffling
// up and re-allocation during intersection. Cleared out on clean.
DynamicList<label> removedEdges_;
// Front propagation
//- The edges in the front along which the intersection propagates
DynamicList<label> frontEdgeEdges_;
//- Inverse of the above
DynamicList<label> edgeFrontEdges_;
// Insertion
//- Indices of triangles considered candidates for insertion
DynamicList<label> candidateTriTris_;
//- Inverse of the above
DynamicList<label> triCandidateTris_;
// Marking
//- Indices of triangles that have been "marked". Used for multiple
// purposes.
DynamicList<label> markedTriTris_;
//- Inverse of the above
DynamicList<label> triMarkedTris_;
// Triangulation
//- Triangulation engine
polygonTriangulate polygonTriangulate_;
// Polygon generation
//- Star propagation engine
star star_;
// Debugging
//- Iteration counter with which to number surface files
label writei_;
// Private Member Functions
//- Check consistency of the mesh
void checkPatchFace(const label patchFacei, const bool isSrc) const;
void checkPatchEdge(const label patchFacei, const bool isSrc) const;
void checkPatchFaces(const bool isSrc) const;
//- Remove an edge
void removeEdge(const label edgei);
//- Remove a tri
void removeTri(const label trii);
//- Create a new triangle and return it's index
label newTrii();
//- Create a new edge and return it's index
label newEdgei();
//- Create a new point and return it's index
label newPointi();
//- Return a point of a given triangle
label triPoint(const label trii, const label triPointi) const;
//- Return the points of a given triangle
triFace triPoints(const label trii) const;
//- Get the values on the points of a given triangle
template <class Type>
FixedList<Type, 3> triPointValues
(
const label trii,
const UList<Type> values
) const;
//- Get a list indicating whether a tri "owns" it's edges; i.e.,
// whether they are numbered in the same order as the points
FixedList<bool, 3> triOwns(const label trii) const;
//- Get the shared point indices of another tri in a tri, or -1 if
// the point is not shared
FixedList<label, 3> triOtherTriPoints
(
const label trii,
const label otherTrii
) const;
//- Get the shared edge indices of another tri in a tri, or -1 if
// the point is not shared
FixedList<label, 3> triOtherTriEdges
(
const label trii,
const label otherTrii
) const;
//- Return the tri-edge-points for a given tri-edge
edge triEdgePoints(const label trii, const label triEdgei) const;
//- Return the edge-points for a given edge
edge edgePoints(const label edgei) const;
//- Return a point for the given patch face
label patchFacePoint
(
const label patchFacei,
const label patchFacePointi,
const bool isSrc
) const;
//- Return the points for the given patch face
triFace patchFacePoints
(
const label patchFacei,
const bool isSrc
) const;
//- Get the values on the points of a given patch face
template <class Type>
FixedList<Type, 3> patchFacePointValues
(
const label patchFacei,
const bool isSrc,
const UList<Type>& values
) const;
//- Get a list indicating whether a patch face "owns" it's edges; i.e.,
// whether they are numbered in the same order as the points
FixedList<bool, 3> patchFaceOwns
(
const label patchFacei,
const bool isSrc
) const;
//- Get the shared point indices of another patch face in a patch face,
// or -1 if the point is not shared
FixedList<label, 3> patchFaceOtherPatchFacePoints
(
const label patchFacei,
const label otherPatchFacei,
const bool isSrc
) const;
//- Return a patch point for the given patch face
label patchFacePatchPoint
(
const label patchFacei,
const label patchFacePointi,
const bool isSrc
) const;
//- Return the patch points for the given patch face
triFace patchFacePatchPoints
(
const label patchFacei,
const bool isSrc
) const;
//- Return a patch edge for the given patch face
label patchFacePatchEdge
(
const label patchFacei,
const label patchFaceEdgei,
const bool isSrc
) const;
//- Return the patch edges for the given patch face
triFace patchFacePatchEdges
(
const label patchFacei,
const bool isSrc
) const;
//- Compute the patch edge that a given edge lies along
label edgePatchEdge(const label edgei, const bool isSrc) const;
//- Compute the patch edges that a given edge lies along
labelPair edgePatchEdges(const label edgei) const;
//- Add a tri. Return the new tri label.
label addTri
(
const triFace& pointis,
const FixedList<label, 3>& edgeis,
const label patchFacei,
const bool isSrc
);
//- Flip an edge
void flipEdge(const label edgei);
//- Return the signed distance squared of the given point outside of a
// triangle's circumcircle. A negative value means the point is inside
// the circle.
scalar circumDistSqr(const label trii, const label pointi) const;
//- Insert points into the given tri or edge. If an edge, points are to
// be given in order along the edge.
void insertPoints
(
const label triOrEdgei,
const bool isTri,
const UList<label>& pointis,
UList<label>& insertionEdgeis,
const UList<label>& fixedEdgeis
);
//- Return whether or not a point can be intersected with the other side
bool pointCanIntersect(const label pointi) const;
//- Return whether or not an edge can be intersected with the other side
bool edgeCanIntersect(const label pointi) const;
/*
//- Snap points to points and edges between two tris
// !!! Not a good idea. Gets tangled up. Need to snap entire faces.
bool snapTris(const label srcTrii, const label tgtTrii);
*/
//- Snap points to points and edges between two patch faces
void snapPatchFaceTris
(
const label srcFacei,
const label tgtFacei,
const scalar snapTol
);
//- Intersect two tris
bool intersectTris(const label srcTrii, const label tgtTrii);
//- Intersect the triangulations of the given patch faces
void intersectPatchFaceTris
(
const label srcFacei,
const label tgtFacei
);
//- Make the triangulations of the given patch faces conform
bool conformPatchFaceTris
(
const label patchFacei,
const label otherPatchFacei,
const bool isSrc
);
//- Make the triangulations of the given patch faces conform
bool conformPatchFaceTris
(
const label srcFacei,
const label tgtFacei
);
//- Combine conformal parts of the given patch faces into intersection
// faces
bool combinePatchFaceTris(const label srcFacei, const label tgtFacei);
//- Initialise the member data
void initialise(const vectorField& srcPointNormals);
//- Clean up removed tris and edges
void clean();
//- Finalise the data in the base class
void finalise();
//- Undo the finalise process so the intersection process can continue
void unFinalise();
//- Write the current state of the surfaces for debugging purposes
void write();
//- Write the current state of the given patch face triangulation for
// debugging purposes
void writePatchFace(const label patchFacei, const bool isSrc) const;
public:
// Runtime type information
virtual word type() const
{
return "tri" + patchIntersection::typeName.capitalise();
}
// Constructors
//- Construct from a source and a target patch
TriPatchIntersection
(
const SrcPatchType& srcPatch,
const TgtPatchType& tgtPatch,
const scalar snapTol
);
//- Construct from a source and a target patch, and specified source
// point normals
TriPatchIntersection
(
const SrcPatchType& srcPatch,
const vectorField& srcPointNormals,
const TgtPatchType& tgtPatch,
const scalar snapTol
);
// Member Functions
// Access
//- ...
inline const DynamicList<labelList>& faceEdges() const
{
return faceEdges_;
}
//- ...
inline const DynamicList<labelPair>& intersectEdgeFaces() const
{
return intersectEdgeFaces_;
}
//- ...
inline const DynamicList<labelPair>& nonIntersectEdgeFaces() const
{
return nonIntersectEdgeFaces_;
}
//- Destructor
virtual ~TriPatchIntersection();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "TriPatchIntersection.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,189 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "patchIntersection.H"
#include "primitivePatch.H"
#include "vtkWritePolyData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const bool Foam::patchIntersection::orientToSource_ = true;
namespace Foam
{
defineTypeNameAndDebug(patchIntersection, 0);
}
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::patchIntersection::report(const word& writeSuffix)
{
{
const primitivePatch patch
(
SubList<face>(faces_, faces_.size()),
points_
);
scalar area = 0, srcArea = 0, tgtArea = 0;
forAll(faces_, facei)
{
const scalar a = faces_[facei].mag(points_);
area += a;
srcArea += faceSrcFaces_[facei] != -1 ? a : 0;
tgtArea += faceTgtFaces_[facei] != -1 ? a : 0;
}
Info<< indent << "Source/target coverage = " << srcArea/area
<< "/" << tgtArea/area << endl;
DynamicList<label> nEdgesNFaces, nFacesNEdges;
forAll(faces_, facei)
{
const label n = faces_[facei].size();
nEdgesNFaces.resize(max(nEdgesNFaces.size(), n + 1), 0);
++ nEdgesNFaces[n];
}
forAll(patch.edgeFaces(), edgei)
{
const label n = patch.edgeFaces()[edgei].size();
nFacesNEdges.resize(max(nFacesNEdges.size(), n + 1), 0);
++ nFacesNEdges[n];
}
Info<< indent << "Faces by number of edges = (";
forAll(nEdgesNFaces, n)
{
Info<< (n ? " " : "") << nEdgesNFaces[n];
}
Info<< ")" << endl << indent << "Edges by number of faces = (";
forAll(nFacesNEdges, n)
{
Info<< (n ? " " : "") << nFacesNEdges[n];
}
Info<< ")" << endl;
}
if (debug)
{
Info<< indent << "Writing intersected patch" << incrIndent << endl;
const fileName patchFileName =
type() + "_patch" + (writeSuffix.empty() ? "" : "_")
+ writeSuffix + ".vtk";
Info<< indent << "Writing patch to " << patchFileName << endl;
vtkWritePolyData::write
(
patchFileName,
"intersectedPatch",
false,
points_,
labelList(),
labelListList(),
faces_,
"srcFace",
false,
labelField(faceSrcFaces_),
"tgtFace",
false,
labelField(faceTgtFaces_)
);
const fileName patchEdgesFileName =
type() + "_patchEdges" + (writeSuffix.empty() ? "" : "_")
+ writeSuffix + ".vtk";
Info<< indent << "Writing patch edges to " << patchEdgesFileName
<< endl;
const primitivePatch patch
(
SubList<face>(faces_, faces_.size()),
points_
);
labelField edgeNFaces(patch.nEdges());
forAll(patch.edgeFaces(), edgei)
{
edgeNFaces[edgei] = patch.edgeFaces()[edgei].size();
}
vtkWritePolyData::write
(
patchEdgesFileName,
"intersectedPatchEdges",
false,
patch.localPoints(),
labelList(),
patch.edges(),
faceList(),
"nFaces",
false,
edgeNFaces
);
Info<< decrIndent;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::patchIntersection::patchIntersection
(
const label srcNPoints,
const label tgtNPoints,
const label srcNEdges,
const label tgtNEdges,
const label srcNFaces,
const label tgtNFaces
)
:
points_(),
srcPointPoints_(srcNPoints),
tgtPointPoints_(tgtNPoints),
pointSrcPoints_(),
pointTgtPoints_(),
srcEdgePoints_(srcNEdges),
tgtEdgePoints_(tgtNEdges),
pointSrcEdges_(),
pointTgtEdges_(),
pointSrcFaces_(),
pointTgtFaces_(),
faces_(),
srcFaceFaces_(srcNFaces),
tgtFaceFaces_(tgtNFaces),
faceSrcFaces_(),
faceTgtFaces_()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::patchIntersection::~patchIntersection()
{}
// ************************************************************************* //

View File

@ -0,0 +1,242 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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::patchIntersection
Description
Base class for patch intersections. Provides type name and debugging. See
templated derivatives for actual functionality.
SourceFiles
patchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef patchIntersection_H
#define patchIntersection_H
#include "DynamicField.H"
#include "face.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class patchIntersection Declaration
\*---------------------------------------------------------------------------*/
class patchIntersection
{
protected:
// Protected static data
//- Flag to set whether the faces of the intersection are oriented the
// same way as the source patch (true) or the target patch (false).
static const bool orientToSource_;
// Protected Data
// Points
//- The intersection points
DynamicField<point> points_;
//- The source points' corresponding intersection points
labelList srcPointPoints_;
//- The target points' corresponding intersection points
labelList tgtPointPoints_;
//- The intersection points' corresponding source points, or -1
// if the point does not correspond to a source point
DynamicList<label> pointSrcPoints_;
//- The intersection points' corresponding target points, or -1
// if the point does not correspond to a target point
DynamicList<label> pointTgtPoints_;
//- The source edges' intersection points. Ordered along the edge.
List<DynamicList<label>> srcEdgePoints_;
//- The target edges' intersection points. Ordered along the edge.
List<DynamicList<label>> tgtEdgePoints_;
//- The intersection points' source edges, or -1 if the point
// is not on a source edge
DynamicList<label> pointSrcEdges_;
//- The intersection points' target edges, or -1 if the point
// is not on a target edge
DynamicList<label> pointTgtEdges_;
// !!! We don't store srcFacePoints and tgtFacePoints. These are
// not needed at present. They also don't have any ordering
// associated with them, so they don't need to be maintained by the
// intersection process. They could be generated on demand as a
// post-processing step from from pointSrcFaces and pointTgtFaces.
//- The intersection points' source faces, or -1 if the point
// did not project into a source face
DynamicList<label> pointSrcFaces_;
//- The intersection points' target faces, or -1 if the point
// did not project into a target face
DynamicList<label> pointTgtFaces_;
// Faces
//- The intersection faces
DynamicList<face> faces_;
//- The source faces' intersection faces
List<DynamicList<label>> srcFaceFaces_;
//- The target faces' intersection faces
List<DynamicList<label>> tgtFaceFaces_;
//- The intersection faces' corresponding source faces, or -1
// if the face does not correspond to source face
DynamicList<label> faceSrcFaces_;
//- The intersection faces' corresponding target faces, or -1
// if the face does not correspond to target face
DynamicList<label> faceTgtFaces_;
// Protected Member Functions
//- Report properties of the intersection process
void report(const word& writeSuffix=word::null);
public:
// Runtime type information
ClassName("patchIntersection");
virtual word type() const = 0;
// Constructors
//- Construct given sizes
patchIntersection
(
const label srcNPoints,
const label tgtNPoints,
const label srcNEdges,
const label tgtNEdges,
const label srcNFaces,
const label tgtNFaces
);
//- Destructor
virtual ~patchIntersection();
// Member Functions
// Points
//- The intersection points
inline const pointField& points() const;
//- The source points' corresponding intersection points
inline const labelList& srcPointPoints() const;
//- The target points' corresponding intersection points
inline const labelList& tgtPointPoints() const;
//- The intersection points' corresponding source points, or -1
// if the point does not correspond to a source point
inline const DynamicList<label>& pointSrcPoints() const;
//- The intersection points' corresponding target points, or -1
// if the point does not correspond to a target point
inline const DynamicList<label>& pointTgtPoints() const;
//- The source edges' intersection points. Ordered along the edge.
inline const List<DynamicList<label>>& srcEdgePoints() const;
//- The target edges' intersection points. Ordered along the edge.
inline const List<DynamicList<label>>& tgtEdgePoints() const;
//- The intersection points' source edges, or -1 if the point
// is not on a source edge
inline const DynamicList<label>& pointSrcEdges() const;
//- The intersection points' target edges, or -1 if the point
// is not on a target edge
inline const DynamicList<label>& pointTgtEdges() const;
//- The intersection points' source faces, or -1 if the point
// did not project into a source face
inline const DynamicList<label>& pointSrcFaces() const;
//- The intersection points' target faces, or -1 if the point
// did not project into a target face
inline const DynamicList<label>& pointTgtFaces() const;
// Faces
//- The intersection faces
inline const faceList& faces() const;
//- The source faces' intersection faces
inline const List<DynamicList<label>>& srcFaceFaces() const;
//- The target faces' intersection faces
inline const List<DynamicList<label>>& tgtFaceFaces() const;
//- The intersection faces' corresponding source faces, or -1
// if the face does not correspond to source face
inline const DynamicList<label>& faceSrcFaces() const;
//- The intersection faces' corresponding target faces, or -1
// if the face does not correspond to target face
inline const DynamicList<label>& faceTgtFaces() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "patchIntersectionI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "patchIntersection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline const Foam::pointField& Foam::patchIntersection::points() const
{
return points_;
}
inline const Foam::labelList&
Foam::patchIntersection::srcPointPoints() const
{
return srcPointPoints_;
}
inline const Foam::labelList&
Foam::patchIntersection::tgtPointPoints() const
{
return tgtPointPoints_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointSrcPoints() const
{
return pointSrcPoints_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointTgtPoints() const
{
return pointTgtPoints_;
}
inline const Foam::List<Foam::DynamicList<Foam::label>>&
Foam::patchIntersection::srcEdgePoints() const
{
return srcEdgePoints_;
}
inline const Foam::List<Foam::DynamicList<Foam::label>>&
Foam::patchIntersection::tgtEdgePoints() const
{
return tgtEdgePoints_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointSrcEdges() const
{
return pointSrcEdges_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointTgtEdges() const
{
return pointTgtEdges_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointSrcFaces() const
{
return pointSrcFaces_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::pointTgtFaces() const
{
return pointTgtFaces_;
}
inline const Foam::faceList& Foam::patchIntersection::faces() const
{
return faces_;
}
inline const Foam::List<Foam::DynamicList<Foam::label>>&
Foam::patchIntersection::srcFaceFaces() const
{
return srcFaceFaces_;
}
inline const Foam::List<Foam::DynamicList<Foam::label>>&
Foam::patchIntersection::tgtFaceFaces() const
{
return tgtFaceFaces_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::faceSrcFaces() const
{
return faceSrcFaces_;
}
inline const Foam::DynamicList<Foam::label>&
Foam::patchIntersection::faceTgtFaces() const
{
return faceTgtFaces_;
}
// ************************************************************************* //

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "polyPatchIntersection.H"
#include "FacePatchIntersection.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyPatchIntersection::polyPatchIntersection
(
const polyPatch& srcPatch,
const polyPatch& tgtPatch,
const scalar snapTol
)
:
PatchIntersection<polyPatch, polyPatch>
(
FacePatchIntersection<polyPatch, polyPatch>
(
srcPatch,
tgtPatch,
snapTol
)
)
{}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,20 +22,21 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::topoAction
Foam::polyPatchIntersection
Description
A virtual base class for topological actions
Intersection between two polyPatches
SourceFiles
polyPatchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef topoAction_H
#define topoAction_H
#ifndef polyPatchIntersection_H
#define polyPatchIntersection_H
#include "typeInfo.H"
#include "autoPtr.H"
#include "PatchIntersection.H"
#include "polyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -43,33 +44,32 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topoAction Declaration
Class polyPatchIntersection Declaration
\*---------------------------------------------------------------------------*/
class topoAction
class polyPatchIntersection
:
public PatchIntersection<polyPatch, polyPatch>
{
public:
// Static Data Members
// Runtime type information
//- Runtime type information
TypeName("topoAction");
virtual word type() const
{
return "primitive" + patchIntersection::typeName.capitalise();
}
// Constructors
//- Construct null
topoAction()
{}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const = 0;
//- Destructor
virtual ~topoAction()
{}
//- Construct from a source and a target patch
polyPatchIntersection
(
const polyPatch& srcPatch,
const polyPatch& tgtPatch,
const scalar snapTol
);
};

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "primitivePatchIntersection.H"
#include "FacePatchIntersection.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::primitivePatchIntersection::primitivePatchIntersection
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const scalar snapTol
)
:
PatchIntersection<primitivePatch, primitivePatch>
(
FacePatchIntersection<primitivePatch, primitivePatch>
(
srcPatch,
tgtPatch,
snapTol
)
)
{}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,18 +22,20 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyRemovePoint
Foam::primitivePatchIntersection
Description
Class containing data for point removal.
Intersection between two primitive patches
SourceFiles
primitivePatchIntersection.C
\*---------------------------------------------------------------------------*/
#ifndef polyRemovePoint_H
#define polyRemovePoint_H
#ifndef primitivePatchIntersection_H
#define primitivePatchIntersection_H
#include "label.H"
#include "topoAction.H"
#include "PatchIntersection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -41,68 +43,32 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyRemovePoint Declaration
Class primitivePatchIntersection Declaration
\*---------------------------------------------------------------------------*/
class polyRemovePoint
class primitivePatchIntersection
:
public topoAction
public PatchIntersection<primitivePatch, primitivePatch>
{
// Private Data
//- Point ID
label pointIndex_;
//- Merge point ID or -1
label mergePointIndex_;
public:
// Static Data Members
// Runtime type information
//- Runtime type information
TypeName("removePoint");
virtual word type() const
{
return "primitive" + patchIntersection::typeName.capitalise();
}
// Constructors
//- Construct null. Used for constructing lists
polyRemovePoint()
:
pointIndex_(-1),
mergePointIndex_(-1)
{}
//- Construct from components
polyRemovePoint(const label pointID, const label mergePointID = -1)
:
pointIndex_(pointID),
mergePointIndex_(mergePointID)
{}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyRemovePoint(*this));
}
// Default Destructor
// Member Functions
//- Return point ID
label pointIndex() const
{
return pointIndex_;
}
label mergePointIndex() const
{
return mergePointIndex_;
}
//- Construct from a source and a target patch
primitivePatchIntersection
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch,
const scalar snapTol
);
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,31 +21,25 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::primitiveTriPatch
\*---------------------------------------------------------------------------*/
#ifndef ZoneDynamicIDs_H
#define ZoneDynamicIDs_H
#ifndef primitiveTriPatch_H
#define primitiveTriPatch_H
#include "DynamicID.H"
#include "meshCellZonesFwd.H"
#include "meshFaceZonesFwd.H"
#include "meshPointZonesFwd.H"
#include "PrimitivePatch.H"
#include "triFace.H"
#include "List.H"
#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Foam::cellZoneDynamicID
typedef DynamicID<meshCellZones> cellZoneDynamicID;
//- Foam::faceZoneDynamicID
typedef DynamicID<meshFaceZones> faceZoneDynamicID;
//- Foam::pointZoneDynamicID
typedef DynamicID<meshPointZones> pointZoneDynamicID;
typedef PrimitivePatch<SubList<triFace>, const pointField&>
primitiveTriPatch;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,114 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021-2023 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 "star.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::star::reset()
{
forAll(starFaceFaces_, starFacei)
{
const label facei = starFaceFaces_[starFacei];
if (facei != -1)
{
faceStarFaces_[starFaceFaces_[starFacei]] = -1;
}
}
starFaceFaces_.clear();
forAll(starEdgeEdges_, starEdgei)
{
const label edgei = starEdgeEdges_[starEdgei].edgei_;
if (edgei != -1)
{
edgeStarEdges_[edgei] = -1;
}
}
starEdgeEdges_.clear();
}
void Foam::star::swapStarEdges(const label starEdgeiA, const label starEdgeiB)
{
const label starEdgeiA0 = starEdgeEdges_[starEdgeiA].starEdgei0_;
const label starEdgeiA1 = starEdgeEdges_[starEdgeiA].starEdgei1_;
const label starEdgeiB0 = starEdgeEdges_[starEdgeiB].starEdgei0_;
const label starEdgeiB1 = starEdgeEdges_[starEdgeiB].starEdgei1_;
if (starEdgeiA0 != -1)
{
starEdgeEdges_[starEdgeiA0].starEdgei1_ = starEdgeiB;
}
if (starEdgeiA1 != -1)
{
starEdgeEdges_[starEdgeiA1].starEdgei0_ = starEdgeiB;
}
if (starEdgeiB0 != -1)
{
starEdgeEdges_[starEdgeiB0].starEdgei1_ = starEdgeiA;
}
if (starEdgeiB1 != -1)
{
starEdgeEdges_[starEdgeiB1].starEdgei0_ = starEdgeiA;
}
Swap(starEdgeEdges_[starEdgeiA], starEdgeEdges_[starEdgeiB]);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::star::star()
:
starFaceFaces_(),
faceStarFaces_(),
starEdgeEdges_(),
edgeStarEdges_(),
work_()
{}
Foam::star::context::context(star& s)
:
star_(s)
{}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::star::~star()
{}
Foam::star::context::~context()
{
star_.reset();
}
// ************************************************************************* //

View File

@ -0,0 +1,274 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021-2023 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::star
Description
Engine for constructing a star-shaped domain by walking
SourceFiles
star.C
\*---------------------------------------------------------------------------*/
#ifndef star_H
#define star_H
#include "DynamicList.H"
#include "labelPair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class star Declaration
\*---------------------------------------------------------------------------*/
class star
{
public:
// Public Classes
//- Star-edge structure. Poor man's linked list link.
struct edgeLink
{
label starEdgei0_;
label edgei_;
label starEdgei1_;
};
//- Context class. Returned by populate and resets the star when it
// goes out of scope.
class context
{
private:
// Private Data
//- Reference to the star
star& star_;
public:
// Constructors
//- Construct from components
context(star& s);
// Destructor
~context();
};
private:
// Private Data
//- For each star face, the external face index, or -1
DynamicList<label> starFaceFaces_;
//- For each external face, the star face index, or -1
DynamicList<label> faceStarFaces_;
//- For each star edge, the external edge index and the two adjacent
// star indices, or {-1, -1, -1}
DynamicList<edgeLink> starEdgeEdges_;
//- For each external edge, the star edge index, or -1
DynamicList<label> edgeStarEdges_;
//- Generic work array
DynamicList<label> work_;
// Private Member Functions
//- Reset the contents of the star, ready for another populate
void reset();
//- Swap the given star edges. The star connectivity remains the same.
// Used to change order of iteration.
void swapStarEdges(const label starEdgeiA, const label starEdgeiB);
//- Expand the star by crossing a given edge into a given face
template<class FaceEdges>
void cross
(
const label edgei,
const label facei,
const UList<FaceEdges>& faceEdges
);
public:
// Constructors
//- Construct null
star();
// Destructor
~star();
// Member Functions
// Access
//- Access the star-face-faces
const DynamicList<label>& starFaceFaces() const
{
return starFaceFaces_;
}
//- Access the face-star-faces
const DynamicList<label>& faceStarFaces() const
{
return faceStarFaces_;
}
//- Access the star-edge-edges
const DynamicList<edgeLink>& starEdgeEdges() const
{
return starEdgeEdges_;
}
//- Access the edge-star-edges
const DynamicList<label>& edgeStarEdges() const
{
return edgeStarEdges_;
}
//- Populate the star given a seed face or edge, a function
// `canCross(edgei, facei)` to determine whether an edge may be cross,
// and face-edge and edge-face addressing. After this is called, the
// forAll... macros below can be user to iterate through the star
// faces, or around the edges of the star perimeter.
template<class CanCross, class FaceEdges>
context populate
(
const label faceOrEdgei,
const bool isFace,
CanCross canCross,
const UList<FaceEdges>& faceEdges,
const UList<labelPair>& edgeFaces,
const label maxIterations = labelMax
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline bool operator==(const star::edgeLink& elA, const star::edgeLink& elB)
{
return
elA.starEdgei0_ == elB.starEdgei0_
&& elA.edgei_ == elB.edgei_
&& elA.starEdgei1_ == elB.starEdgei1_;
}
inline bool operator!=(const star::edgeLink& elA, const star::edgeLink& elB)
{
return !(elA == elB);
}
inline Ostream& operator<<(Ostream& os, const star::edgeLink& el)
{
return os
<< '(' << el.starEdgei0_ << ')'
<< ' ' << el.edgei_ << ' '
<< '(' << el.starEdgei1_ << ')';
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define forAllStarFaces(star, starFacei, facei) \
for \
( \
label starFacei = 0, \
facei = \
( \
star.starFaceFaces().size() > 0 \
? (star).starFaceFaces()[starFacei] \
: -1 \
); \
starFacei < star.starFaceFaces().size(); \
++ starFacei, \
facei = \
( \
starFacei < (star).starFaceFaces().size() \
? (star).starFaceFaces()[starFacei] \
: -1 \
) \
)
#define forAllStarEdges(star, i, starEdgei, edgei) \
for \
( \
label i = 0, \
starEdgei = 0, \
edgei = \
( \
(star).starEdgeEdges().size() > 0 \
? (star).starEdgeEdges()[0].edgei_ \
: -1 \
); \
i < (star).starEdgeEdges().size(); \
++ i, \
starEdgei = (star).starEdgeEdges()[starEdgei].starEdgei1_, \
edgei = \
( \
starEdgei != -1 \
? (star).starEdgeEdges()[starEdgei].edgei_ \
: -1 \
) \
)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "starTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,324 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021-2023 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 "OFstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class FaceEdges>
void Foam::star::cross
(
const label edgei,
const label facei,
const UList<FaceEdges>& faceEdges
)
{
// Add the new face to the star
starFaceFaces_.append(facei);
faceStarFaces_[facei] = starFaceFaces_.size() - 1;
// Store information about this edge, then remove it from the star
const label faceEdgei = findIndex(faceEdges[facei], edgei);
const label starEdgei = edgeStarEdges_[faceEdges[facei][faceEdgei]];
const label starEdgei0 = starEdgeEdges_[starEdgei].starEdgei0_;
const label starEdgei1 = starEdgeEdges_[starEdgei].starEdgei1_;
edgeStarEdges_[starEdgeEdges_[starEdgei].edgei_] = -1;
starEdgeEdges_[starEdgei] = {-1, -1, -1};
// Walk backwards and forwards around the face until we find edges not
// already in the star. Remove edges in the star as we go.
label faceEdgeiR = faceEdgei, starEdgeiR1 = starEdgei1;
while (true)
{
const label faceEdgeiRR = faceEdges[facei].rcIndex(faceEdgeiR);
const label edgeiRR = faceEdges[facei][faceEdgeiRR];
const label starEdgeiRR = edgeStarEdges_[edgeiRR];
if (starEdgeiRR == -1 || starEdgeiRR != starEdgeiR1) break;
faceEdgeiR = faceEdgeiRR;
starEdgeiR1 = starEdgeEdges_[starEdgeiRR].starEdgei1_;
edgeStarEdges_[starEdgeEdges_[starEdgeiRR].edgei_] = -1;
starEdgeEdges_[starEdgeiRR] = {-1, -1, -1};
}
label faceEdgeiF = faceEdgei, starEdgeiF0 = starEdgei0;
while (true)
{
const label faceEdgeiFF = faceEdges[facei].fcIndex(faceEdgeiF);
const label edgeiFF = faceEdges[facei][faceEdgeiFF];
const label starEdgeiFF = edgeStarEdges_[edgeiFF];
if (starEdgeiFF == -1 || starEdgeiFF != starEdgeiF0) break;
faceEdgeiF = faceEdgeiFF;
starEdgeiF0 = starEdgeEdges_[starEdgeiFF].starEdgei0_;
edgeStarEdges_[starEdgeEdges_[starEdgeiFF].edgei_] = -1;
starEdgeEdges_[starEdgeiFF] = {-1, -1, -1};
}
// Get the first edge after the forwards edge
label faceEdgej = faceEdges[facei].fcIndex(faceEdgeiF);
// If there are no face edges not yet in the star, then all edges are
// to be removed. Just connect up the adjacent edges.
if (faceEdgej == faceEdgeiR)
{
starEdgeEdges_[starEdgeiF0].starEdgei1_ = starEdgeiR1;
starEdgeEdges_[starEdgeiR1].starEdgei0_ = starEdgeiF0;
}
// If there are face edges not yet in the star, then loop over them and
// add them into the star
while (faceEdgej != faceEdgeiR)
{
const bool isFirst =
faceEdgej == faceEdges[facei].fcIndex(faceEdgeiF);
const bool isLast =
faceEdgej == faceEdges[facei].rcIndex(faceEdgeiR);
const label starEdgej0 =
isFirst ? starEdgeiF0 : starEdgeEdges_.size() - 1;
const label starEdgej = starEdgeEdges_.size();
const label starEdgej1 =
isLast ? starEdgeiR1 : starEdgeEdges_.size() + 1;
const label edgej = faceEdges[facei][faceEdgej];
if (isFirst && starEdgeiF0 != -1)
{
starEdgeEdges_[starEdgeiF0].starEdgei1_ = starEdgej;
}
starEdgeEdges_.append({starEdgej0, edgej, starEdgej1});
edgeStarEdges_[edgej] = starEdgej;
if (isLast && starEdgeiR1 != -1)
{
starEdgeEdges_[starEdgeiR1].starEdgei0_ = starEdgej;
}
faceEdgej = faceEdges[facei].fcIndex(faceEdgej);
}
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class CanCross, class FaceEdges>
Foam::star::context Foam::star::populate
(
const label faceOrEdgei,
const bool isFace,
CanCross canCross,
const UList<FaceEdges>& faceEdges,
const UList<labelPair>& edgeFaces,
const label maxIterations
)
{
// Expand to fit supplied data
faceStarFaces_.resize(faceEdges.size(), -1);
edgeStarEdges_.resize(edgeFaces.size(), -1);
// Initialise
if (isFace)
{
// Add one edge from the to the star and then cross into the face
const label edgei = faceEdges[faceOrEdgei][0];
starEdgeEdges_.append({-1, edgei, -1});
edgeStarEdges_[edgei] = 0;
cross(edgei, faceOrEdgei, faceEdges);
// Add the edge back again as it was removed by the crossing. Reuse the
// first element of the starEdgeEdges array. Join up to close the star.
const label starEdgei0 = starEdgeEdges_.size() - 1;
const label starEdgei1 = 1;
starEdgeEdges_[0] = {starEdgei0, edgei, starEdgei1};
edgeStarEdges_[edgei] = 0;
starEdgeEdges_[starEdgei0].starEdgei1_ = 0;
starEdgeEdges_[starEdgei1].starEdgei0_ = 0;
}
else
{
bool first = true;
// Loop the adjacent faces
forAll(edgeFaces[faceOrEdgei], edgeFacei)
{
const label facei = edgeFaces[faceOrEdgei][edgeFacei];
if (facei == -1) continue;
// If the first face, add the edge to the star
if (first)
{
starEdgeEdges_.append({-1, faceOrEdgei, -1});
edgeStarEdges_[faceOrEdgei] = 0;
first = false;
}
// If the second face, add the edge into the star again as it was
// removed by the previous crossing. Reuse the first element of the
// starEdgeEdges array. Join up to close the star.
else
{
const label starEdgei0 = starEdgeEdges_.size() - 1;
const label starEdgei1 = 1;
starEdgeEdges_[0] = {starEdgei0, faceOrEdgei, starEdgei1};
edgeStarEdges_[faceOrEdgei] = 0;
starEdgeEdges_[starEdgei0].starEdgei1_ = 0;
starEdgeEdges_[starEdgei1].starEdgei0_ = 0;
}
// Cross the edge into the face
cross(faceOrEdgei, facei, faceEdges);
}
}
// Walk the star to get intersected faces and their perimeter
label iterations = 0;
forAll(starEdgeEdges_, starEdgei)
{
if (iterations == maxIterations) break;
++ iterations;
const label edgei = starEdgeEdges_[starEdgei].edgei_;
if (edgei == -1) continue;
// Get the adjacent non-star face
label facei = -1;
forAll(edgeFaces[edgei], edgeFacei)
{
facei = edgeFaces[edgei][edgeFacei];
if (facei != -1 && faceStarFaces_[facei] == -1)
{
break;
}
facei = -1;
}
// If the crossing condition is satisfied then expand the star
if (facei != -1 && canCross(edgei, facei))
{
cross(edgei, facei, faceEdges);
}
}
// Remove spikes
forAll(starEdgeEdges_, starEdgei)
{
label starEdgei0 = starEdgei;
label starEdgei1 = starEdgeEdges_[starEdgei0].starEdgei1_;
if (starEdgei1 == -1) continue;
label edgei0 = starEdgeEdges_[starEdgei0].edgei_;
label edgei1 = starEdgeEdges_[starEdgei1].edgei_;
while (edgei0 != -1 && edgei1 != -1 && edgei0 == edgei1)
{
const label starEdgei00 = starEdgeEdges_[starEdgei0].starEdgei0_;
const label starEdgei11 = starEdgeEdges_[starEdgei1].starEdgei1_;
edgeStarEdges_[edgei0] = -1;
starEdgeEdges_[starEdgei0] = {-1, -1, -1};
edgeStarEdges_[edgei1] = -1;
starEdgeEdges_[starEdgei1] = {-1, -1, -1};
if (starEdgei00 != -1)
{
starEdgeEdges_[starEdgei00].starEdgei1_ = starEdgei11;
starEdgei0 = starEdgei00;
edgei0 = starEdgeEdges_[starEdgei00].edgei_;
}
else
{
starEdgei0 = edgei0 = -1;
}
if (starEdgei11 != -1)
{
starEdgeEdges_[starEdgei11].starEdgei0_ = starEdgei00;
starEdgei1 = starEdgei11;
edgei1 = starEdgeEdges_[starEdgei11].edgei_;
}
else
{
starEdgei1 = edgei0 = -1;
}
}
}
// Allocate map storage
DynamicList<label>& oldStarEdgeNewStarEdges = work_;
oldStarEdgeNewStarEdges.resize(starEdgeEdges_.size());
oldStarEdgeNewStarEdges = -1;
// Remove empty edges by shuffling up
label starEdgei = 0;
forAll(starEdgeEdges_, starEdgej)
{
if (starEdgeEdges_[starEdgej].edgei_ != -1)
{
oldStarEdgeNewStarEdges[starEdgej] = starEdgei;
starEdgeEdges_[starEdgei] = starEdgeEdges_[starEdgej];
++ starEdgei;
}
}
starEdgeEdges_.resize(starEdgei);
// Map
forAll(starEdgeEdges_, starEdgei)
{
label& starEdgei0 = starEdgeEdges_[starEdgei].starEdgei0_;
if (starEdgei0 != -1)
{
starEdgei0 = oldStarEdgeNewStarEdges[starEdgei0];
}
label& starEdgei1 = starEdgeEdges_[starEdgei].starEdgei1_;
if (starEdgei1 != -1)
{
starEdgei1 = oldStarEdgeNewStarEdges[starEdgei1];
}
}
/*
// Randomise swapping
static Random rndGen(12345);
const label starEdgeiA = rndGen.sampleAB<label>(0, starEdgeEdges_.size());
const label starEdgeiB = rndGen.sampleAB<label>(0, starEdgeEdges_.size());
swapStarEdges(starEdgeiA, starEdgeiB);
*/
// If the star is open, put the first edge first, so that forAllStarEdges
// always starts in the right place
forAll(starEdgeEdges_, starEdgei)
{
if (starEdgeEdges_[starEdgei].starEdgei0_ == -1)
{
swapStarEdges(starEdgei, 0);
}
}
return context(*this);
}
// ************************************************************************* //

View File

@ -1,5 +1,3 @@
polyTopoChange/topoAction/topoActions.C
polyTopoChange/polyTopoChange.C
polyTopoChange/addPatchCellLayer.C
polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
@ -17,26 +15,8 @@ $(hexRef8)/hexRef8.C
$(hexRef8)/hexRef8Data.C
$(hexRef8)/refinementHistory.C
polyTopoChanger/polyTopoChanger/polyTopoChanger.C
polyTopoChanger/polyMeshModifier/polyMeshModifier.C
polyTopoChanger/attachPolyTopoChanger/attachPolyTopoChanger.C
enrichedPatch = polyTopoChanger/slidingInterface/enrichedPatch
$(enrichedPatch)/enrichedPatch.C
$(enrichedPatch)/enrichedPatchPointMap.C
$(enrichedPatch)/enrichedPatchFaces.C
$(enrichedPatch)/enrichedPatchPointPoints.C
$(enrichedPatch)/enrichedPatchCutFaces.C
$(enrichedPatch)/enrichedPatchMasterPoints.C
polyTopoChanger/slidingInterface/slidingInterface.C
polyTopoChanger/slidingInterface/slidingInterfaceProjectPoints.C
polyTopoChanger/slidingInterface/coupleSlidingInterface.C
polyTopoChanger/slidingInterface/slidingInterfaceAttachedAddressing.C
polyTopoChanger/slidingInterface/slidingInterfaceClearCouple.C
polyTopoChanger/slidingInterface/decoupleSlidingInterface.C
polyTopoChanger/perfectInterface/perfectInterface.C
perfectInterface/perfectInterface.C
mergePatchPairs/mergePatchPairs.C
repatchMesh/repatchMesh.C
repatchMesh/repatchPatch.C

View File

@ -1,4 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/meshCheck/lnInclude \
@ -6,6 +7,7 @@ EXE_INC = \
-I$(LIB_SRC)/mesh/extrudeModel/lnInclude
LIB_LIBS = \
-lfileFormats \
-lfiniteVolume \
-lmeshTools \
-lmeshCheck \

View File

@ -0,0 +1,754 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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 "mergePatchPairs.H"
#include "fvMesh.H"
#include "polyPatchIntersection.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineDebugSwitch(mergePatchPairs, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::mergePatchPairs::findPatchIndex(const word& patchName) const
{
label patchIndex = mesh_.boundaryMesh().findIndex(patchName);
if (patchIndex == -1)
{
FatalErrorInFunction
<< "Cannot find patch " << patchName << exit(FatalError);
}
return patchIndex;
}
Foam::Pair<Foam::label> Foam::mergePatchPairs::findPatchIndices
(
const Pair<word>& patchNamePair
) const
{
return Pair<label>
{
findPatchIndex(patchNamePair.first()),
findPatchIndex(patchNamePair.second())
};
}
Foam::List<Foam::Pair<Foam::label>> Foam::mergePatchPairs::patchPairs
(
const List<Pair<word>>& patchNamePairs
) const
{
List<Pair<label>> patchPairs(patchNamePairs.size());
forAll(patchNamePairs, pairi)
{
patchPairs[pairi] = findPatchIndices(patchNamePairs[pairi]);
}
return patchPairs;
}
void Foam::mergePatchPairs::removePoints
(
polyTopoChange& meshMod,
const List<Pair<label>>& patchPairs
) const
{
boolList removedPoints(mesh_.nPoints(), false);
forAll(patchPairs, ppi)
{
UIndirectList<bool>
(
removedPoints,
mesh_.boundaryMesh()[patchPairs[ppi].first()].meshPoints()
) = true;
UIndirectList<bool>
(
removedPoints,
mesh_.boundaryMesh()[patchPairs[ppi].second()].meshPoints()
) = true;
}
forAll(removedPoints, pi)
{
if (removedPoints[pi])
{
meshMod.removePoint(pi, -1);
}
}
}
void Foam::mergePatchPairs::addPoints
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const
{
const pointField& intersectionPoints = intersection.points();
newPoints_.setSize(intersectionPoints.size(), -1);
forAll(intersectionPoints, pi)
{
const label pointi = meshMod.addPoint
(
intersectionPoints[pi],
-1,
-1,
false
);
if (debug)
{
Info<< "Adding point "
<< pointi << " "
<< intersectionPoints[pi]
<< endl;
}
newPoints_[pi] = pointi;
}
}
void Foam::mergePatchPairs::removeFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const
{
const polyPatch& srcPatch = intersection.srcPatch();
const polyPatch& tgtPatch = intersection.tgtPatch();
// Remove existing source face
forAll(srcPatch, fi)
{
meshMod.removeFace(srcPatch.start() + fi, -1);
}
// Remove existing target face
forAll(tgtPatch, fi)
{
meshMod.removeFace(tgtPatch.start() + fi, -1);
}
}
Foam::face Foam::mergePatchPairs::mapFace(const face& f) const
{
face newFace(f);
forAll(newFace, fpi)
{
newFace[fpi] = newPoints_[newFace[fpi]];
}
return newFace;
}
void Foam::mergePatchPairs::addFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const
{
const faceList& faces = intersection.faces();
const DynamicList<label>& faceSrcFaces = intersection.faceSrcFaces();
const DynamicList<label>& faceTgtFaces = intersection.faceTgtFaces();
const label srcPatchi = intersection.srcPatch().index();
const label srcPatchStart = intersection.srcPatch().start();
const label tgtPatchi = intersection.tgtPatch().index();
const label tgtPatchStart = intersection.tgtPatch().start();
forAll(faces, fi)
{
const face f = mapFace(faces[fi]);
if (faceSrcFaces[fi] != -1 && faceTgtFaces[fi] != -1)
{
const label srcFacei = srcPatchStart + faceSrcFaces[fi];
const label tgtFacei = tgtPatchStart + faceTgtFaces[fi];
const label srcOwn = mesh_.faceOwner()[srcFacei];
const label tgtOwn = mesh_.faceOwner()[tgtFacei];
if (srcOwn < tgtOwn)
{
meshMod.addFace
(
f, // Face to add
srcOwn, // Owner cell
tgtOwn, // Neighbour cell
-1, // Master point index
-1, // Master edge index
srcFacei, // Master face index
false, // Flip
-1, // Patch index
-1, // Zone index
false // Zone sign
);
if (debug)
{
Info<< "Adding internal face " << f
<< " owner celli:" << srcOwn
<< " neighbour celli:" << tgtOwn << endl;
}
}
else
{
meshMod.addFace
(
f.reverseFace(), // Face to add
tgtOwn, // Owner cell
srcOwn, // Neighbour cell
-1, // Master point index
-1, // Master edge index
tgtFacei, // Master face index
false, // Flip
-1, // Patch index
-1, // Zone index
false // Zone sign
);
if (debug)
{
Info<< "Adding internal face " << f
<< " owner celli:" << tgtOwn
<< " neighbour celli:" << srcOwn << endl;
}
}
}
else if (faceSrcFaces[fi] != -1 && faceTgtFaces[fi] == -1)
{
const label srcFacei = srcPatchStart + faceSrcFaces[fi];
const label srcOwn = mesh_.faceOwner()[srcFacei];
// Get source face zone info
const label zoneIndex = mesh_.faceZones().whichZone(srcFacei);
bool zoneFlip = false;
if (zoneIndex >= 0)
{
const faceZone& fZone = mesh_.faceZones()[zoneIndex];
zoneFlip = fZone.flipMap()[fZone.whichFace(srcFacei)];
}
meshMod.addFace
(
f, // Face to add
srcOwn, // Owner cell
-1, // Neighbour cell
-1, // Master point index
-1, // Master edge index
srcFacei, // Master face index
false, // Flip
srcPatchi, // Patch index
zoneIndex, // Zone index
zoneFlip // Zone sign
);
if (debug)
{
Info<< "Adding patch face " << f
<< " owner celli:" << srcOwn
<< " patchi:" << srcPatchi << endl;
}
}
else if (faceSrcFaces[fi] == -1 && faceTgtFaces[fi] != -1)
{
const label tgtFacei = tgtPatchStart + faceTgtFaces[fi];
const label tgtOwn = mesh_.faceOwner()[tgtFacei];
// Get target face zone info
const label zoneIndex = mesh_.faceZones().whichZone(tgtFacei);
bool zoneFlip = false;
if (zoneIndex >= 0)
{
const faceZone& fZone = mesh_.faceZones()[zoneIndex];
zoneFlip = fZone.flipMap()[fZone.whichFace(tgtFacei)];
}
meshMod.addFace
(
f, // Face to add
tgtOwn, // Owner cell
-1, // Neighbour cell
-1, // Master point index
-1, // Master edge index
tgtFacei, // Master face index
false, // Flip
tgtPatchi, // Patch index
zoneIndex, // Zone index
zoneFlip // Zone sign
);
if (debug)
{
Info<< "Adding patch face " << f
<< " owner celli:" << tgtOwn
<< " patchi:" << tgtPatchi << endl;
}
}
else
{
FatalErrorInFunction
<< "Both faceSrcFaces and faceTgtFaces are -1 for face " << fi
<< exit(FatalError);
}
}
}
void Foam::mergePatchPairs::addEdgeAddedPoints
(
HashTable<labelList, edge, Hash<edge>>& edgeAddedPoints,
const primitivePatch& patch,
const List<DynamicList<label>>& patchEdgeAddedPoints
) const
{
const edgeList& patchEdges = patch.edges();
const labelList& pmp = patch.meshPoints();
forAll(patchEdgeAddedPoints, ei)
{
if (patchEdgeAddedPoints[ei].size())
{
labelList newEdgePoints(patchEdgeAddedPoints[ei].size());
forAll(newEdgePoints, npi)
{
newEdgePoints[npi] = newPoints_[patchEdgeAddedPoints[ei][npi]];
}
edgeAddedPoints.insert
(
edge(pmp[patchEdges[ei].start()], pmp[patchEdges[ei].end()]),
newEdgePoints
);
}
}
}
void Foam::mergePatchPairs::updatePoints
(
labelList& pointMap,
const primitivePatch& patch,
const labelList& pointPoints
) const
{
const labelList& pmp = patch.meshPoints();
forAll(pointPoints, ppi)
{
pointMap[pmp[ppi]] = newPoints_[pointPoints[ppi]];
}
}
void Foam::mergePatchPairs::modifyFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const
{
HashTable<labelList, edge, Hash<edge>> edgeAddedPoints;
addEdgeAddedPoints
(
edgeAddedPoints,
intersection.srcPatch(),
intersection.srcEdgePoints()
);
addEdgeAddedPoints
(
edgeAddedPoints,
intersection.tgtPatch(),
intersection.tgtEdgePoints()
);
labelList pointMap(identityMap(mesh_.nPoints()));
updatePoints
(
pointMap,
intersection.srcPatch(),
intersection.srcPointPoints()
);
updatePoints
(
pointMap,
intersection.tgtPatch(),
intersection.tgtPointPoints()
);
DynamicList<label> modifiedFace;
const faceList& meshFaces = mesh_.faces();
forAll(meshFaces, fi)
{
if (meshMod.faceRemoved(fi))
{
continue;
}
const face& f = meshFaces[fi];
bool addNext = true;
bool modified = false;
forAll(f, pi)
{
const edge e(f[pi], f[f.fcIndex(pi)]);
const HashTable<labelList, edge, Hash<edge>>::const_iterator iter =
edgeAddedPoints.find(e);
if (iter != edgeAddedPoints.end())
{
const labelList& addedPoints = iter();
if (edge::compare(e, iter.key()) == 1)
{
modifiedFace.append(addedPoints);
}
else
{
forAllReverse(addedPoints, i)
{
modifiedFace.append(addedPoints[i]);
}
}
if (f[f.fcIndex(pi)] == f[0])
{
// modifiedFace.first() = modifiedFace.last();
// modifiedFace.setSize(modifiedFace.size() - 1);
modifiedFace.erase(modifiedFace.begin());
}
addNext = false;
modified = true;
}
else
{
if (addNext)
{
modifiedFace.append(pointMap[f[pi]]);
}
else
{
addNext = true;
}
}
}
// Update points of point-connected faces if necessary
if (!modified)
{
forAll(f, pi)
{
if (pointMap[f[pi]] != f[pi])
{
if (!modified)
{
modifiedFace = f;
modified = true;
}
modifiedFace[pi] = pointMap[f[pi]];
}
}
}
if (modified)
{
// Get current zone info
const label zoneIndex = mesh_.faceZones().whichZone(fi);
bool zoneFlip = false;
if (zoneIndex >= 0)
{
const faceZone& fZone = mesh_.faceZones()[zoneIndex];
zoneFlip = fZone.flipMap()[fZone.whichFace(fi)];
}
if (mesh_.isInternalFace(fi))
{
if (debug)
{
Info<< "Modifying internal face " << fi
<< " newface:" << modifiedFace
<< " oldFace:" << f
<< " centre:" << mesh_.faceCentres()[fi]
<< " owner:" << mesh_.faceOwner()[fi]
<< " neighbour:" << mesh_.faceNeighbour()[fi]
<< endl;
}
meshMod.modifyFace
(
face(modifiedFace), // Modified face
fi, // index of face being modified
mesh_.faceOwner()[fi], // Owner cell
mesh_.faceNeighbour()[fi], // Neighbour cell
false, // Face flip
-1, // Patch index
zoneIndex, // Zone index
zoneFlip // Zone flip
);
}
else
{
if (debug)
{
Info<< "Modifying boundary face " << fi
<< " newface:" << modifiedFace
<< " oldFace:" << f
<< " centre:" << mesh_.faceCentres()[fi]
<< " owner:" << mesh_.faceOwner()[fi]
<< " patch:" << mesh_.boundaryMesh().whichPatch(fi)
<< endl;
}
meshMod.modifyFace
(
face(modifiedFace), // Modified face
fi, // index of face being modified
mesh_.faceOwner()[fi], // Owner cell
-1, // Neighbour cell
false, // Face flip
mesh_.boundaryMesh().whichPatch(fi), // Patch index
zoneIndex, // Zone index
zoneFlip // Zone flip
);
}
}
modifiedFace.clear();
}
}
void Foam::mergePatchPairs::intersectPatchPair
(
polyTopoChange& meshMod,
const polyPatch& srcPatch,
const polyPatch& tgtPatch
) const
{
polyPatchIntersection intersection(srcPatch, tgtPatch, snapTol_);
removeFaces(meshMod, intersection);
addPoints(meshMod, intersection);
addFaces(meshMod, intersection);
modifyFaces(meshMod, intersection);
}
Foam::autoPtr<Foam::polyTopoChangeMap> Foam::mergePatchPairs::merge
(
const List<Pair<label>>& patchPairs
) const
{
// Construct the mesh modifier for merging patch-pairs
polyTopoChange meshMod(mesh_);
removePoints(meshMod, patchPairs);
// Accumulate the mesh modifications for merging the patch-pairs
forAll(patchPairs, ppi)
{
// Intersect patch-pair and update the mesh modifier
intersectPatchPair
(
meshMod,
mesh_.boundaryMesh()[patchPairs[ppi].first()],
mesh_.boundaryMesh()[patchPairs[ppi].second()]
);
}
// Change mesh and return map
return meshMod.changeMesh
(
mesh_,
false // Update mesh without moving points and calculating meshPhi
);
}
bool Foam::mergePatchPairs::connected
(
boolList& patchPoints,
const label patchi
) const
{
const labelList& pmp = mesh_.boundaryMesh()[patchi].meshPoints();
forAll(pmp, pi)
{
if (patchPoints[pmp[pi]])
{
return true;
}
patchPoints[pmp[pi]] = true;
}
return false;
}
template<class Mesh>
inline void Foam::mergePatchPairs::merge
(
Mesh& mesh,
const List<Pair<word>>& patchPairNames
) const
{
const List<Pair<label>> patchPairs(this->patchPairs(patchPairNames));
bool connected = false;
if (patchPairs.size() > 1)
{
boolList patchPoints(mesh_.nPoints(), false);
forAll(patchPairs, ppi)
{
if
(
this->connected(patchPoints, patchPairs[ppi].first())
|| this->connected(patchPoints, patchPairs[ppi].second())
)
{
connected = true;
break;
}
}
}
if (connected)
{
// Merge patch-pairs and update fields
forAll(patchPairNames, ppi)
{
Info<< "Merging patch pair " << patchPairNames[ppi] << endl;
mesh.topoChange
(
merge(List<Pair<label>>(patchPairs[ppi]))
);
}
}
else
{
Info<< "Merging patch pairs " << patchPairNames << endl;
mesh.topoChange(merge(patchPairs));
}
// Find the stitched patches that are now zero-sized
labelHashSet zeroSizedMergedPatches;
forAll(patchPairs, ppi)
{
if (mesh.boundaryMesh()[patchPairs[ppi].first()].size() == 0)
{
zeroSizedMergedPatches.insert(patchPairs[ppi].first());
}
if (mesh.boundaryMesh()[patchPairs[ppi].second()].size() == 0)
{
zeroSizedMergedPatches.insert(patchPairs[ppi].second());
}
}
// Create a list of the remaining patch old patch labels
labelList remainingPatches
(
mesh.boundaryMesh().size() - zeroSizedMergedPatches.size()
);
label rpi = 0;
forAll(mesh.boundaryMesh(), patchi)
{
if (!zeroSizedMergedPatches.found(patchi))
{
remainingPatches[rpi++] = patchi;
}
}
// Remove the zero-sized stitched patches
mesh_.reorderPatches(remainingPatches, true);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::mergePatchPairs::mergePatchPairs
(
polyMesh& mesh,
const List<Pair<word>>& patchPairNames,
const scalar snapTol
)
:
mesh_(mesh),
snapTol_(snapTol)
{
// Merge patch-pairs and update fields
merge(mesh, patchPairNames);
}
Foam::mergePatchPairs::mergePatchPairs
(
fvMesh& mesh,
const List<Pair<word>>& patchPairNames,
const scalar snapTol
)
:
mesh_(mesh),
snapTol_(snapTol)
{
// Merge patch-pairs and update fields
merge(mesh, patchPairNames);
}
// ************************************************************************* //

View File

@ -0,0 +1,206 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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::mergePatchPairs
Description
Class to stitch mesh by merging patch-pairs
SourceFiles
mergePatchPairs.C
\*---------------------------------------------------------------------------*/
#ifndef mergePatchPairs_H
#define mergePatchPairs_H
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class fvMesh;
class polyTopoChange;
class polyPatchIntersection;
/*---------------------------------------------------------------------------*\
Class mergePatchPairs Declaration
\*---------------------------------------------------------------------------*/
class mergePatchPairs
{
// Private data
polyMesh& mesh_;
const scalar snapTol_;
mutable labelList newPoints_;
// Private member functions
//- Return patch index of given patch name
// Throw a fatal error if the patch name is not found
label findPatchIndex(const word& patchName) const;
//- Return the patch index pairs
// corresponding to the give patch name pairs
Pair<label> findPatchIndices
(
const Pair<word>& patchNamePair
) const;
//- Set the patchPairs_ patch index pair list
// from the patchNamePairs_ patch name pair list
List<Pair<label>> patchPairs
(
const List<Pair<word>>& patchNamePairs
) const;
//- Remove all the points corresponding to the patch pairs to be merged
void removePoints
(
polyTopoChange& meshMod,
const List<Pair<label>>& patchPairs
) const;
//- Add the points generated by the patch intersection
void addPoints
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const;
//- Remove all the faces of the source and target patches to me merged
void removeFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const;
//- Map face from patch-local point indices to global mesh point indices
face mapFace(const face& f) const;
//- Add a new internal and patch face to the mesh
void addFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const;
void addEdgeAddedPoints
(
HashTable<labelList, edge, Hash<edge>>& edgeAddedPoints,
const primitivePatch& patch,
const List<DynamicList<label>>& patchEdgeAddedPoints
) const;
//- Update the point indices of the changed points in the pointMap
void updatePoints
(
labelList& pointMap,
const primitivePatch& patch,
const labelList& pointPoints
) const;
//- Modify the edges of existing edge-connected faces
// by changing the end-point and adding any additional points
// and point-connected faces by updating the point indices
void modifyFaces
(
polyTopoChange& meshMod,
const polyPatchIntersection& intersection
) const;
//- Intersect the given pair of patches and insert all face changes
// into meshMod
void intersectPatchPair
(
polyTopoChange& meshMod,
const polyPatch& srcPatch,
const polyPatch& tgtPatch
) const;
//- Merge all the given patch pairs
// returning the polyTopoChangeMap for subsequent field mapping
autoPtr<polyTopoChangeMap> merge
(
const List<Pair<label>>& patchPairIndices
) const;
//- Return true if the points of patchi are connected to any
// of the patches already tested
bool connected
(
boolList& patchPoints,
const label patchi
) const;
//- Merge all the given patch pairs
// and update all the fields cached in the given mesh
template<class Mesh>
inline void merge
(
Mesh& mesh,
const List<Pair<word>>& patchPairNames
) const;
public:
// Static Data Members
ClassName("mergePatchPairs");
// Constructors
mergePatchPairs
(
polyMesh& mesh,
const List<Pair<word>>& patchPairNames,
const scalar snapTol
);
mergePatchPairs
(
fvMesh& mesh,
const List<Pair<word>>& patchPairNames,
const scalar snapTol
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -74,10 +74,10 @@ Foam::perfectInterface::perfectInterface
const word& slavePatchName
)
:
polyMeshModifier(name, mesh),
faceZoneIndex_(mesh.faceZones().findIndex(faceZoneName)),
masterPatchIndex_(mesh.boundaryMesh().findIndex(masterPatchName)),
slavePatchIndex_(mesh.boundaryMesh().findIndex(slavePatchName))
mesh_(mesh),
faceZoneIndex_(mesh_.faceZones().findIndex(faceZoneName)),
masterPatchIndex_(mesh_.boundaryMesh().findIndex(masterPatchName)),
slavePatchIndex_(mesh_.boundaryMesh().findIndex(slavePatchName))
{}
@ -89,12 +89,6 @@ Foam::perfectInterface::~perfectInterface()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::perfectInterface::changeTopology() const
{
return true;
}
void Foam::perfectInterface::setRefinement
(
const indirectPrimitivePatch& pp0,
@ -102,9 +96,7 @@ void Foam::perfectInterface::setRefinement
polyTopoChange& ref
) const
{
const polyMesh& mesh = this->mesh();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Some aliases
const edgeList& edges0 = pp0.edges();
@ -135,7 +127,7 @@ void Foam::perfectInterface::setRefinement
// Determine pointMapping in mesh point labels. Uses geometric
// comparison to find correspondence between patch points.
labelList renumberPoints(mesh.points().size());
labelList renumberPoints(mesh_.points().size());
forAll(renumberPoints, i)
{
renumberPoints[i] = i;
@ -202,7 +194,7 @@ void Foam::perfectInterface::setRefinement
if (meshPointi != renumberPoints[meshPointi])
{
const labelList& pFaces = mesh.pointFaces()[meshPointi];
const labelList& pFaces = mesh_.pointFaces()[meshPointi];
forAll(pFaces, pFacei)
{
@ -226,7 +218,7 @@ void Foam::perfectInterface::setRefinement
{
WarningInFunction
<< "Found face " << facei << " vertices "
<< mesh.faces()[facei] << " whose points are"
<< mesh_.faces()[facei] << " whose points are"
<< " used both by master patch and slave patch" << endl;
}
}
@ -236,7 +228,7 @@ void Foam::perfectInterface::setRefinement
forAllConstIter(labelHashSet, affectedFaces, iter)
{
const label facei = iter.key();
const face& f = mesh.faces()[facei];
const face& f = mesh_.faces()[facei];
face newFace(f.size());
@ -249,22 +241,22 @@ void Foam::perfectInterface::setRefinement
label patchi = -1;
if (mesh.isInternalFace(facei))
if (mesh_.isInternalFace(facei))
{
nbr = mesh.faceNeighbour()[facei];
nbr = mesh_.faceNeighbour()[facei];
}
else
{
patchi = patches.whichPatch(facei);
}
const label zoneID = mesh.faceZones().whichZone(facei);
const label zoneID = mesh_.faceZones().whichZone(facei);
bool zoneFlip = false;
if (zoneID >= 0)
{
const faceZone& fZone = mesh.faceZones()[zoneID];
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
}
@ -272,7 +264,7 @@ void Foam::perfectInterface::setRefinement
(
newFace, // modified face
facei, // label of face being modified
mesh.faceOwner()[facei], // owner
mesh_.faceOwner()[facei], // owner
nbr, // neighbour
false, // face flip
patchi, // patch for face
@ -306,12 +298,12 @@ void Foam::perfectInterface::setRefinement
// comment above about patch1 and patch0 never sharing points) and
// becoming internal.
const boolList& mfFlip =
mesh.faceZones()[faceZoneIndex_].flipMap();
mesh_.faceZones()[faceZoneIndex_].flipMap();
forAll(pp0, i)
{
const label facei = pp0.addressing()[i];
const face& f = mesh.faces()[facei];
const face& f = mesh_.faces()[facei];
face newFace(f.size());
@ -320,9 +312,9 @@ void Foam::perfectInterface::setRefinement
newFace[fp] = renumberPoints[f[fp]];
}
const label own = mesh.faceOwner()[facei];
const label own = mesh_.faceOwner()[facei];
const label pp1Facei = pp1.addressing()[from0To1Faces[i]];
const label nbr = mesh.faceOwner()[pp1Facei];
const label nbr = mesh_.faceOwner()[pp1Facei];
if (own < nbr)
{
@ -361,30 +353,27 @@ void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
if (debug)
{
Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
<< "for object " << name() << " : "
<< "masterPatchIndex_:" << masterPatchIndex_
<< " slavePatchIndex_:" << slavePatchIndex_
<< " faceZoneIndex_:" << faceZoneIndex_ << endl;
}
const polyMesh& mesh = this->mesh();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const polyPatch& patch0 = patches[masterPatchIndex_];
const polyPatch& patch1 = patches[slavePatchIndex_];
labelList pp0Labels(identityMap(patch0.size())+patch0.start());
indirectPrimitivePatch pp0
(
IndirectList<face>(mesh.faces(), pp0Labels),
mesh.points()
IndirectList<face>(mesh_.faces(), pp0Labels),
mesh_.points()
);
labelList pp1Labels(identityMap(patch1.size())+patch1.start());
indirectPrimitivePatch pp1
(
IndirectList<face>(mesh.faces(), pp1Labels),
mesh.points()
IndirectList<face>(mesh_.faces(), pp1Labels),
mesh_.points()
);
setRefinement(pp0, pp1, ref);

View File

@ -36,7 +36,7 @@ SourceFiles
#ifndef perfectInterface_H
#define perfectInterface_H
#include "polyMeshModifier.H"
#include "polyMesh.H"
#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,16 +44,18 @@ SourceFiles
namespace Foam
{
class polyTopoChange;
/*---------------------------------------------------------------------------*\
Class perfectInterface Declaration
\*---------------------------------------------------------------------------*/
class perfectInterface
:
public polyMeshModifier
{
// Private Data
const polyMesh& mesh_;
//- Master face zone ID
label faceZoneIndex_;
@ -101,9 +103,6 @@ public:
// Member Functions
//- Check for topology change
virtual bool changeTopology() const;
//- Insert the layer addition/removal instructions
// into the topological change
virtual void setRefinement(polyTopoChange&) const;
@ -119,14 +118,6 @@ public:
polyTopoChange&
) const;
//- Modify motion points to comply with the topological change
virtual void modifyMotionPoints(pointField& motionPoints) const
{}
//- Force recalculation of locally stored data on topological change
virtual void topoChange(const polyTopoChangeMap&)
{}
// Member Operators

View File

@ -1,350 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::polyAddFace
Description
A face addition data class. A face can be inflated either from a
point or from another face and can either be in internal or a
boundary face.
\*---------------------------------------------------------------------------*/
#ifndef polyAddFace_H
#define polyAddFace_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "label.H"
#include "face.H"
#include "topoAction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyAddFace Declaration
\*---------------------------------------------------------------------------*/
class polyAddFace
:
public topoAction
{
// Private Data
//- Face identifier
face face_;
//- Face owner
label owner_;
//- Face neighbour
label neighbour_;
//- Master point ID for faces blown up from points
label masterPointIndex_;
//- Master edge ID for faces blown up from edges
label masterEdgeIndex_;
//- Master face ID for faces blown up from faces
label masterFaceIndex_;
//- Does the face flux need to be flipped
bool flipFaceFlux_;
//- Boundary patch ID
label patchIndex_;
//- Face zone ID
label zoneIndex_;
//- Face zone flip
bool zoneFlip_;
public:
// Static Data Members
//- Runtime type information
TypeName("addFace");
// Constructors
//- Construct null. Used for constructing lists
polyAddFace()
:
face_(0),
owner_(-1),
neighbour_(-1),
masterPointIndex_(-1),
masterEdgeIndex_(-1),
masterFaceIndex_(-1),
flipFaceFlux_(false),
patchIndex_(-1),
zoneIndex_(-1),
zoneFlip_(false)
{}
//- Construct from components
polyAddFace
(
const face& f,
const label owner,
const label neighbour,
const label masterPointID,
const label masterEdgeID,
const label masterFaceID,
const bool flipFaceFlux,
const label patchID,
const label zoneID,
const bool zoneFlip
)
:
face_(f),
owner_(owner),
neighbour_(neighbour),
masterPointIndex_(masterPointID),
masterEdgeIndex_(masterEdgeID),
masterFaceIndex_(masterFaceID),
flipFaceFlux_(flipFaceFlux),
patchIndex_(patchID),
zoneIndex_(zoneID),
zoneFlip_(zoneFlip)
{
if (face_.size() < 3)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (min(face_) < 0)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (min(owner_, neighbour_) >= 0 && owner_ == neighbour_)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (neighbour_ >= 0 && patchID >= 0)
{
FatalErrorInFunction
<< ". This is not allowed.\n"
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (owner_ < 0 && zoneID < 0)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (zoneIndex_ == -1 && zoneFlip)
{
FatalErrorInFunction
<< "belong to zone. This is not allowed.\n"
<< "Face: " << face_
<< " masterPointID:" << masterPointIndex_
<< " masterEdgeID:" << masterEdgeIndex_
<< " masterFaceID:" << masterFaceIndex_
<< " patchID:" << patchIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyAddFace(*this));
}
// Default Destructor
// Member Functions
//- Return face
const face& newFace() const
{
return face_;
}
//- Return owner cell
label owner() const
{
return owner_;
}
//- Return neighbour cell
label neighbour() const
{
return neighbour_;
}
//- Is the face mastered by a point
bool isPointMaster() const
{
return masterPointIndex_ >= 0;
}
//- Is the face mastered by an edge
bool isEdgeMaster() const
{
return masterEdgeIndex_ >= 0;
}
//- Is the face mastered by another face
bool isFaceMaster() const
{
return masterFaceIndex_ >= 0;
}
//- Is the face appended with no master
bool appended() const
{
return !isPointMaster() && !isEdgeMaster() && !isFaceMaster();
}
//- Return master point ID
label masterPointIndex() const
{
return masterPointIndex_;
}
//- Return master edge ID
label masterEdgeIndex() const
{
return masterEdgeIndex_;
}
//- Return master face ID
label masterFaceIndex() const
{
return masterFaceIndex_;
}
//- Does the face flux need to be flipped
bool flipFaceFlux() const
{
return flipFaceFlux_;
}
//- Does the face belong to a boundary patch?
bool isInPatch() const
{
return patchIndex_ >= 0;
}
//- Boundary patch ID
label patchIndex() const
{
return patchIndex_;
}
//- Does the face belong to a zone?
bool isInZone() const
{
return zoneIndex_ >= 0;
}
//- Is the face only a zone face (i.e. not belonging to a cell)
bool onlyInZone() const
{
return zoneIndex_ >= 0 && owner_ < 0 && neighbour_ < 0;
}
//- Face zone ID
label zoneIndex() const
{
return zoneIndex_;
}
//- Face zone flip
label zoneFlip() const
{
return zoneFlip_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,170 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::polyAddPoint
Description
Class containing data for point addition.
\*---------------------------------------------------------------------------*/
#ifndef polyAddPoint_H
#define polyAddPoint_H
#include "label.H"
#include "point.H"
#include "topoAction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyAddPoint Declaration
\*---------------------------------------------------------------------------*/
class polyAddPoint
:
public topoAction
{
// Private Data
//- Point to add
point p_;
//- Master point
label masterPointIndex_;
//- Point zone ID
label zoneIndex_;
//- Does the point support a cell
bool inCell_;
public:
// Static Data Members
//- Runtime type information
TypeName("addPoint");
// Constructors
//- Construct null. Used only for list construction
polyAddPoint()
:
p_(Zero),
masterPointIndex_(-1),
zoneIndex_(-1),
inCell_(false)
{}
//- Construct from components
polyAddPoint
(
const point& p,
const label masterPointID,
const label zoneID,
const bool inCell
)
:
p_(p),
masterPointIndex_(masterPointID),
zoneIndex_(zoneID),
inCell_(inCell)
{
if (zoneIndex_ < 0 && !inCell)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "point: " << p
<< " master: " << masterPointIndex_
<< " zone: " << zoneIndex_
<< abort(FatalError);
}
}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyAddPoint(*this));
}
// Default Destructor
// Member Functions
//- Point location
const point& newPoint() const
{
return p_;
}
//- Master point label
label masterPointIndex() const
{
return masterPointIndex_;
}
//- Is the point appended with no master
bool appended() const
{
return masterPointIndex_ < 0;
}
//- Does the point belong to a zone?
bool isInZone() const
{
return zoneIndex_ >= 0;
}
//- Point zone ID
label zoneIndex() const
{
return zoneIndex_;
}
//- Does the point support a cell
bool inCell() const
{
return inCell_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,279 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::polyModifyFace
Description
Class describing modification of a face.
\*---------------------------------------------------------------------------*/
#ifndef polyModifyFace_H
#define polyModifyFace_H
#include "label.H"
#include "face.H"
#include "topoAction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyModifyFace Declaration
\*---------------------------------------------------------------------------*/
class polyModifyFace
:
public topoAction
{
// Private Data
//- Face
face face_;
//- Master face ID
label faceIndex_;
//- Face owner
label owner_;
//- Face neighbour
label neighbour_;
//- Does the face flux need to be flipped
bool flipFaceFlux_;
//- Boundary patch ID
label patchIndex_;
//- Remove from current zone
bool removeFromZone_;
//- Face zone ID
label zoneIndex_;
//- Face zone flip
bool zoneFlip_;
public:
// Static Data Members
//- Runtime type information
TypeName("modifyFace");
// Constructors
//- Construct null. Used in constructing lists
polyModifyFace()
:
face_(0),
faceIndex_(-1),
owner_(-1),
neighbour_(-1),
flipFaceFlux_(false),
patchIndex_(-1),
removeFromZone_(false),
zoneIndex_(-1),
zoneFlip_(false)
{}
//- Construct from components
polyModifyFace
(
const face& f,
const label faceID,
const label owner,
const label neighbour,
const bool flipFaceFlux,
const label patchID,
const bool removeFromZone,
const label zoneID,
const bool zoneFlip
)
:
face_(f),
faceIndex_(faceID),
owner_(owner),
neighbour_(neighbour),
flipFaceFlux_(flipFaceFlux),
patchIndex_(patchID),
removeFromZone_(removeFromZone),
zoneIndex_(zoneID),
zoneFlip_(zoneFlip)
{
if (face_.size() < 3)
{
FatalErrorInFunction
<< "Invalid face: less than 3 points. This is not allowed\n"
<< "Face: " << face_
<< " faceID:" << faceIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (min(face_) < 0)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< " faceID:" << faceIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (min(owner_, neighbour_) >= 0 && owner_ == neighbour_)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< " faceID:" << faceIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
if (neighbour_ >= 0 && patchIndex_ >= 0)
{
FatalErrorInFunction
<< "This is not allowed.\n"
<< "Face: " << face_
<< " faceID:" << faceIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< " patchID:" << patchIndex_
<< abort(FatalError);
}
if (zoneIndex_ < 0 && zoneFlip )
{
FatalErrorInFunction
<< "belong to zone. This is not allowed.\n"
<< "Face: " << face_
<< " faceID:" << faceIndex_
<< " owner:" << owner_
<< " neighbour:" << neighbour_
<< abort(FatalError);
}
}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyModifyFace(*this));
}
// Default Destructor
// Member Functions
//- Return face
const face& newFace() const
{
return face_;
}
//- Return master face ID
label faceIndex() const
{
return faceIndex_;
}
//- Return owner cell ID
label owner() const
{
return owner_;
}
//- Return owner cell ID
label neighbour() const
{
return neighbour_;
}
//- Does the face flux need to be flipped
bool flipFaceFlux() const
{
return flipFaceFlux_;
}
//- Does the face belong to a boundary patch?
bool isInPatch() const
{
return patchIndex_ >= 0;
}
//- Boundary patch ID
label patchIndex() const
{
return patchIndex_;
}
//- Does the face belong to a zone?
bool isInZone() const
{
return zoneIndex_ >= 0;
}
//- Is the face only a zone face (i.e. not belonging to a cell)
bool onlyInZone() const
{
return zoneIndex_ >= 0 && owner_ < 0 && neighbour_ < 0;
}
bool removeFromZone() const
{
return removeFromZone_;
}
//- Face zone ID
label zoneIndex() const
{
return zoneIndex_;
}
//- Face zone flip
label zoneFlip() const
{
return zoneFlip_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,164 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::polyModifyPoint
Description
Class describing modification of a point.
\*---------------------------------------------------------------------------*/
#ifndef polyModifyPoint_H
#define polyModifyPoint_H
#include "label.H"
#include "point.H"
#include "topoAction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyModifyPoint Declaration
\*---------------------------------------------------------------------------*/
class polyModifyPoint
:
public topoAction
{
// Private Data
//- Point ID
label pointIndex_;
//- New point location
point location_;
//- Remove from current zone
bool removeFromZone_;
//- New zone ID
label zoneIndex_;
//- Does the point support a cell
bool inCell_;
public:
// Static Data Members
//- Runtime type information
TypeName("modifyPoint");
// Constructors
//- Construct null. Used only for list construction
polyModifyPoint()
:
pointIndex_(-1),
location_(Zero),
removeFromZone_(false),
zoneIndex_(-1),
inCell_(false)
{}
//- Construct from components
polyModifyPoint
(
const label pointID,
const point& newP,
const bool removeFromZone,
const label newZoneID,
const bool inCell
)
:
pointIndex_(pointID),
location_(newP),
removeFromZone_(removeFromZone),
zoneIndex_(newZoneID),
inCell_(inCell)
{}
//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
return autoPtr<topoAction>(new polyModifyPoint(*this));
}
// Default Destructor
// Member Functions
//- Point ID
label pointIndex() const
{
return pointIndex_;
}
//- New point location
const point& newPoint() const
{
return location_;
}
//- Does the point belong to a zone?
bool isInZone() const
{
return zoneIndex_ >= 0;
}
//- Should the point be removed from current zone
bool removeFromZone() const
{
return removeFromZone_;
}
//- Point zone ID
label zoneIndex() const
{
return zoneIndex_;
}
//- Does the point support a cell
bool inCell() const
{
return inCell_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -26,12 +26,6 @@ License
#include "polyTopoChange.H"
#include "SortableList.H"
#include "polyMesh.H"
#include "polyAddPoint.H"
#include "polyModifyPoint.H"
#include "polyRemovePoint.H"
#include "polyAddFace.H"
#include "polyModifyFace.H"
#include "polyRemoveFace.H"
#include "objectMap.H"
#include "processorPolyPatch.H"
#include "fvMesh.H"
@ -2490,98 +2484,6 @@ void Foam::polyTopoChange::setCapacity
}
Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
{
if (isType<polyAddPoint>(action))
{
const polyAddPoint& pap = refCast<const polyAddPoint>(action);
return addPoint
(
pap.newPoint(),
pap.masterPointIndex(),
pap.zoneIndex(),
pap.inCell()
);
}
else if (isType<polyModifyPoint>(action))
{
const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
modifyPoint
(
pmp.pointIndex(),
pmp.newPoint(),
pmp.zoneIndex(),
pmp.inCell()
);
return -1;
}
else if (isType<polyRemovePoint>(action))
{
const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
removePoint(prp.pointIndex(), prp.mergePointIndex());
return -1;
}
else if (isType<polyAddFace>(action))
{
const polyAddFace& paf = refCast<const polyAddFace>(action);
return addFace
(
paf.newFace(),
paf.owner(),
paf.neighbour(),
paf.masterPointIndex(),
paf.masterEdgeIndex(),
paf.masterFaceIndex(),
paf.flipFaceFlux(),
paf.patchIndex(),
paf.zoneIndex(),
paf.zoneFlip()
);
}
else if (isType<polyModifyFace>(action))
{
const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
modifyFace
(
pmf.newFace(),
pmf.faceIndex(),
pmf.owner(),
pmf.neighbour(),
pmf.flipFaceFlux(),
pmf.patchIndex(),
pmf.zoneIndex(),
pmf.zoneFlip()
);
return -1;
}
else if (isType<polyRemoveFace>(action))
{
const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
removeFace(prf.faceIndex(), prf.mergeFaceIndex());
return -1;
}
else
{
FatalErrorInFunction
<< "Unknown type of topoChange: " << action.type()
<< abort(FatalError);
// Dummy return to keep compiler happy
return -1;
}
}
Foam::label Foam::polyTopoChange::addPoint
(
const point& pt,

View File

@ -499,9 +499,6 @@ public:
//- Move all points. Incompatible with other topology changes.
void movePoints(const pointField& newPoints);
//- For compatibility with polyTopoChange: set topological action.
label setAction(const topoAction& action);
//- Add point. Return new point label.
// Notes:
// - masterPointID can be < 0 (appended points)

View File

@ -1,115 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "attachPolyTopoChanger.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::attachPolyTopoChanger::attachPolyTopoChanger
(
polyMesh& mesh
)
:
polyTopoChanger(mesh)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::attachPolyTopoChanger::attach(const bool removeEmptyPatches)
{
if (debug)
{
Pout<< "void attachPolyTopoChanger::attach(): "
<< "Attaching mesh" << endl;
}
// Save current file instance
const fileName oldInst = mesh_.facesInstance();
// Execute all polyMeshModifiers
changeMesh(false); // no inflation
const pointField p = mesh_.oldPoints();
mesh_.movePoints(p);
if (debug)
{
Pout<< "Clearing mesh." << endl;
}
if (removeEmptyPatches)
{
// Re-do the boundary patches, removing the ones with zero size
const polyBoundaryMesh& oldPatches = mesh_.boundaryMesh();
List<polyPatch*> newPatches(oldPatches.size());
label nNewPatches = 0;
forAll(oldPatches, patchi)
{
if (oldPatches[patchi].size())
{
newPatches[nNewPatches] = oldPatches[patchi].clone
(
mesh_.boundaryMesh(),
nNewPatches,
oldPatches[patchi].size(),
oldPatches[patchi].start()
).ptr();
nNewPatches++;
}
else
{
if (debug)
{
Pout<< "Removing zero-sized patch " << patchi
<< " named " << oldPatches[patchi].name() << endl;
}
}
}
newPatches.setSize(nNewPatches);
mesh_.removeBoundary();
mesh_.addPatches(newPatches);
}
// Reset the file instance to overwrite the original mesh
mesh_.setInstance(oldInst);
if (debug)
{
Pout<< "void attachPolyTopoChanger::attach(): "
<< "Finished attaching mesh" << endl;
}
}
// ************************************************************************* //

View File

@ -1,102 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::attachPolyTopoChanger
Description
This class is derived from polyMesh and serves as a tool for
statically connecting pieces of a mesh by executing the mesh
modifiers and cleaning the mesh.
The idea is that a mesh can be built from pieces and put together
using various mesh modifiers (mainly sliding interfaces) which are
not needed during the run. Therefore, once the mesh is assembled
and mesh modification triggered, the newly created point, face and
cell zones can be cleared together with the mesh modifiers thus
creating a singly connected static mesh.
Note:
All point, face and cell zoning will be lost! Do it after
attaching the parts of the mesh, as the point, face and cell
numbering changes.
\*---------------------------------------------------------------------------*/
#ifndef attachPolyTopoChanger_H
#define attachPolyTopoChanger_H
#include "polyTopoChanger.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class attachPolyTopoChanger Declaration
\*---------------------------------------------------------------------------*/
class attachPolyTopoChanger
:
public polyTopoChanger
{
public:
// Constructors
//- Read constructor for given polyMesh
explicit attachPolyTopoChanger(polyMesh&);
//- Disallow default bitwise copy construction
attachPolyTopoChanger(const attachPolyTopoChanger&) = delete;
//- Destructor
virtual ~attachPolyTopoChanger()
{}
// Member Functions
//- Attach mesh. By default filter out empty patches.
void attach(const bool removeEmptyPatches = true);
// Member Operators
//- Disallow default bitwise assignment
void operator=(const attachPolyTopoChanger&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,131 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::polyMeshModifier
Description
Virtual base class for mesh modifiers.
SourceFiles
polyMeshModifier.C
\*---------------------------------------------------------------------------*/
#ifndef polyMeshModifier_H
#define polyMeshModifier_H
#include "typeInfo.H"
#include "pointFieldFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class polyMesh;
class polyTopoChange;
class polyTopoChangeMap;
/*---------------------------------------------------------------------------*\
Class polyMeshModifier Declaration
\*---------------------------------------------------------------------------*/
class polyMeshModifier
{
// Private Data
//- Name of modifier
word name_;
//- Reference to mesh
const polyMesh& mesh_;
public:
// Static Data Members
//- Runtime type information
TypeName("meshModifier");
// Constructors
//- Construct from components
polyMeshModifier
(
const word& name,
const polyMesh& mme
);
//- Disallow default bitwise copy construction
polyMeshModifier(const polyMeshModifier&) = delete;
//- Destructor
virtual ~polyMeshModifier();
// Member Functions
//- Return name of this modifier
const word& name() const
{
return name_;
}
//- Return reference to mesh
const polyMesh& mesh() const;
//- Check for topology change
virtual bool changeTopology() const = 0;
//- Insert the topological change instructions
virtual void setRefinement(polyTopoChange&) const = 0;
//- Modify motion points to comply with the topological change
virtual void modifyMotionPoints(pointField& motionPoints) const = 0;
//- Force recalculation of locally stored data on topological change
virtual void topoChange(const polyTopoChangeMap&) = 0;
// Member Operators
//- Disallow default bitwise assignment
void operator=(const polyMeshModifier&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,143 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "polyTopoChanger.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polyTopoChanger, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyTopoChanger::polyTopoChanger(polyMesh& mesh)
:
PtrList<polyMeshModifier>(),
mesh_(mesh)
{}
bool Foam::polyTopoChanger::changeTopology() const
{
// Go through all mesh modifiers and accumulate the morphing information
const PtrList<polyMeshModifier>& topoChanges = *this;
bool triggerChange = false;
forAll(topoChanges, morphI)
{
bool curTriggerChange = topoChanges[morphI].changeTopology();
if (debug)
{
Info<< "Modifier " << morphI << " named "
<< topoChanges[morphI].name();
if (curTriggerChange)
{
Info<< " morphing" << endl;
}
else
{
Info<< " unchanged" << endl;
}
}
triggerChange = triggerChange || curTriggerChange;
}
return triggerChange;
}
Foam::autoPtr<Foam::polyTopoChange>
Foam::polyTopoChanger::topoChangeRequest() const
{
// Collect changes from all modifiers
const PtrList<polyMeshModifier>& topoChanges = *this;
polyTopoChange* refPtr(new polyTopoChange(mesh()));
polyTopoChange& ref = *refPtr;
forAll(topoChanges, morphI)
{
topoChanges[morphI].setRefinement(ref);
}
return autoPtr<polyTopoChange>(refPtr);
}
void Foam::polyTopoChanger::update(const polyTopoChangeMap& m)
{
// Go through all mesh modifiers and accumulate the morphing information
PtrList<polyMeshModifier>& topoChanges = *this;
forAll(topoChanges, morphI)
{
topoChanges[morphI].topoChange(m);
}
}
Foam::autoPtr<Foam::polyTopoChangeMap> Foam::polyTopoChanger::changeMesh
(
const bool inflate,
const bool syncParallel,
const bool orderCells,
const bool orderPoints
)
{
if (changeTopology())
{
autoPtr<polyTopoChange> ref = topoChangeRequest();
autoPtr<polyTopoChangeMap> topoChangeMap = ref().changeMesh
(
mesh_,
inflate,
syncParallel,
orderCells,
orderPoints
);
update(topoChangeMap());
mesh_.topoChange(topoChangeMap());
return topoChangeMap;
}
else
{
return autoPtr<polyTopoChangeMap>(nullptr);
}
}
// ************************************************************************* //

View File

@ -1,423 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "slidingInterface.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "polyTopoChange.H"
#include "polyModifyFace.H"
#include "polyModifyPoint.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::slidingInterface::decoupleInterface
(
polyTopoChange& ref
) const
{
if (debug)
{
Pout<< "void slidingInterface::decoupleInterface("
<< "polyTopoChange& ref) const : "
<< "Decoupling sliding interface " << name() << endl;
}
if (!attached_)
{
if (debug)
{
Pout<< "void slidingInterface::decoupleInterface("
<< "polyTopoChange& ref) const : "
<< "Interface already decoupled." << endl;
}
return;
}
// Clear previous couple
clearCouple(ref);
const polyMesh& mesh = this->mesh();
const faceList& faces = mesh.faces();
const cellList& cells = mesh.cells();
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
// Master side
const primitiveFacePatch& masterPatch =
mesh.faceZones()[masterFaceZoneIndex_.index()]();
const labelList& masterPatchAddr =
mesh.faceZones()[masterFaceZoneIndex_.index()];
const boolList& masterPatchFlip =
mesh.faceZones()[masterFaceZoneIndex_.index()].flipMap();
const labelList& masterFc = masterFaceCells();
// Recover faces in master patch
forAll(masterPatchAddr, facei)
{
// Make a copy of the face and turn it if necessary
face newFace = faces[masterPatchAddr[facei]];
if (masterPatchFlip[facei])
{
newFace.flip();
}
ref.setAction
(
polyModifyFace
(
newFace, // new face
masterPatchAddr[facei], // master face index
masterFc[facei], // owner
-1, // neighbour
false, // flux flip
masterPatchIndex_.index(), // patch ID
false, // remove from zone
masterFaceZoneIndex_.index(), // zone ID
false // zone flip. Face corrected
)
);
// Pout<< "Modifying master patch face no "
// << masterPatchAddr[facei]
// << " face: " << faces[masterPatchAddr[facei]]
// << " old owner: " << own[masterPatchAddr[facei]]
// << " new owner: " << masterFc[facei]
// << endl;
}
// Slave side
const primitiveFacePatch& slavePatch =
mesh.faceZones()[slaveFaceZoneIndex_.index()]();
const labelList& slavePatchAddr =
mesh.faceZones()[slaveFaceZoneIndex_.index()];
const boolList& slavePatchFlip =
mesh.faceZones()[slaveFaceZoneIndex_.index()].flipMap();
const labelList& slaveFc = slaveFaceCells();
// Grab retired point mapping
const Map<label>& rpm = retiredPointMap();
// Recover faces in slave patch
forAll(slavePatchAddr, facei)
{
// Make a copy of face and turn it if necessary
face newFace = faces[slavePatchAddr[facei]];
if (slavePatchFlip[facei])
{
newFace.flip();
}
// Recover retired points on the slave side
forAll(newFace, pointi)
{
Map<label>::const_iterator rpmIter = rpm.find(newFace[pointi]);
if (rpmIter != rpm.end())
{
// Master of retired point; grab its original
// Pout<< "Reinstating retired point: " << newFace[pointi]
// << " with old: " << rpm.find(newFace[pointi])()
// << endl;
newFace[pointi] = rpmIter();
}
}
ref.setAction
(
polyModifyFace
(
newFace, // new face
slavePatchAddr[facei], // master face index
slaveFc[facei], // owner
-1, // neighbour
false, // flux flip
slavePatchIndex_.index(), // patch ID
false, // remove from zone
slaveFaceZoneIndex_.index(), // zone ID
false // zone flip. Face corrected
)
);
}
// Re-create the master stick-out faces
// Grab the list of faces in the layer
const labelList& masterStickOuts = masterStickOutFaces();
forAll(masterStickOuts, facei)
{
// Renumber the face and remove additional points
const label curFaceID = masterStickOuts[facei];
const face& oldFace = faces[curFaceID];
DynamicList<label> newFaceLabels(oldFace.size());
bool changed = false;
forAll(oldFace, pointi)
{
// Check if the point is removed
if (ref.pointRemoved(oldFace[pointi]))
{
// Point removed; skip it
changed = true;
}
else
{
newFaceLabels.append(oldFace[pointi]);
}
}
if (changed)
{
if (newFaceLabels.size() < 3)
{
FatalErrorInFunction
<< "Face " << curFaceID << " reduced to less than "
<< "3 points. Topological/cutting error." << nl
<< "Old face: " << oldFace << " new face: " << newFaceLabels
<< abort(FatalError);
}
// Get face zone and its flip
label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
bool modifiedFaceZoneFlip = false;
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
mesh.faceZones()[modifiedFaceZone].flipMap()
[
mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
];
}
face newFace;
newFace.transfer(newFaceLabels);
// Pout<< "Modifying master stick-out face " << curFaceID
// << " old face: " << oldFace
// << " new face: " << newFace
// << endl;
// Modify the face
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
nei[curFaceID], // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
}
}
// Re-create the slave stick-out faces
labelHashSet slaveLayerCellFaceMap
(
primitiveMesh::facesPerCell_*(masterPatch.size() + slavePatch.size())
);
forAll(slaveFc, facei)
{
const labelList& curFaces = cells[slaveFc[facei]];
forAll(curFaces, facei)
{
// Check if the face belongs to the slave face zone; and
// if it has been removed; if not add it
if
(
mesh.faceZones().whichZone(curFaces[facei])
!= slaveFaceZoneIndex_.index()
&& !ref.faceRemoved(curFaces[facei])
)
{
slaveLayerCellFaceMap.insert(curFaces[facei]);
}
}
}
// Grab the list of faces in the layer
const labelList& slaveStickOuts = slaveStickOutFaces();
// Grab master point mapping
const Map<label>& masterPm = masterPatch.meshPointMap();
forAll(slaveStickOuts, facei)
{
// Renumber the face and remove additional points
const label curFaceID = slaveStickOuts[facei];
const face& oldFace = faces[curFaceID];
DynamicList<label> newFaceLabels(oldFace.size());
bool changed = false;
forAll(oldFace, pointi)
{
// Check if the point is removed or retired
if (rpm.found(oldFace[pointi]))
{
// Master of retired point; grab its original
changed = true;
// Pout<< "Reinstating retired point: " << oldFace[pointi]
// << " with old: " << rpm.find(oldFace[pointi])()
// << endl;
newFaceLabels.append(rpm.find(oldFace[pointi])());
}
else if (ref.pointRemoved(oldFace[pointi]))
{
// Point removed; skip it
changed = true;
}
else if (masterPm.found(oldFace[pointi]))
{
// Point from master patch only; skip it
changed = true;
}
else
{
newFaceLabels.append(oldFace[pointi]);
}
}
if (changed)
{
if (newFaceLabels.size() < 3)
{
FatalErrorInFunction
<< "Face " << curFaceID << " reduced to less than "
<< "3 points. Topological/cutting error." << nl
<< "Old face: " << oldFace << " new face: " << newFaceLabels
<< abort(FatalError);
}
// Get face zone and its flip
label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
bool modifiedFaceZoneFlip = false;
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
mesh.faceZones()[modifiedFaceZone].flipMap()
[
mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
];
}
face newFace;
newFace.transfer(newFaceLabels);
// Pout<< "Modifying slave stick-out face " << curFaceID
// << " old face: " << oldFace
// << " new face: " << newFace
// << endl;
// Modify the face
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
nei[curFaceID], // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
}
}
// Bring all slave patch points back to life
const pointField& points = mesh.points();
const labelList& slaveMeshPoints =
mesh.faceZones()[slaveFaceZoneIndex_.index()]().meshPoints();
forAll(slaveMeshPoints, pointi)
{
ref.setAction
(
polyModifyPoint
(
slaveMeshPoints[pointi], // point ID
points[slaveMeshPoints[pointi]], // point
false, // remove from zone
mesh.pointZones().whichZone(slaveMeshPoints[pointi]), // zone
true // in a cell
)
);
}
// Clear the retired point numbering
retiredPointMapPtr_->clear();
// Finished decoupling
attached_ = false;
if (debug)
{
Pout<< "void slidingInterface::coupleInterface("
<< "polyTopoChange& ref) const : "
<< "Finished decoupling sliding interface " << name() << endl;
}
}
// ************************************************************************* //

View File

@ -1,287 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "enrichedPatch.H"
#include "demandDrivenData.H"
#include "OFstream.H"
#include "meshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(enrichedPatch, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::enrichedPatch::calcMeshPoints() const
{
if (meshPointsPtr_)
{
FatalErrorInFunction
<< "Mesh points already calculated."
<< abort(FatalError);
}
meshPointsPtr_ = new labelList(pointMap().toc());
labelList& mp = *meshPointsPtr_;
sort(mp);
}
void Foam::enrichedPatch::calcLocalFaces() const
{
if (localFacesPtr_)
{
FatalErrorInFunction
<< "Local faces already calculated."
<< abort(FatalError);
}
// Invert mesh points and renumber faces using it
const labelList& mp = meshPoints();
Map<label> mpLookup(2*mp.size());
forAll(mp, mpI)
{
mpLookup.insert(mp[mpI], mpI);
}
const faceList& faces = enrichedFaces();
localFacesPtr_ = new faceList(faces.size());
faceList& lf = *localFacesPtr_;
forAll(faces, facei)
{
const face& f = faces[facei];
face& curlf = lf[facei];
curlf.setSize(f.size());
forAll(f, pointi)
{
curlf[pointi] = mpLookup.find(f[pointi])();
}
}
}
void Foam::enrichedPatch::calcLocalPoints() const
{
if (localPointsPtr_)
{
FatalErrorInFunction
<< "Local points already calculated."
<< abort(FatalError);
}
const labelList& mp = meshPoints();
localPointsPtr_ = new pointField(mp.size());
pointField& lp = *localPointsPtr_;
forAll(lp, i)
{
lp[i] = pointMap().find(mp[i])();
}
}
void Foam::enrichedPatch::clearOut()
{
deleteDemandDrivenData(enrichedFacesPtr_);
deleteDemandDrivenData(meshPointsPtr_);
deleteDemandDrivenData(localFacesPtr_);
deleteDemandDrivenData(localPointsPtr_);
deleteDemandDrivenData(pointPointsPtr_);
deleteDemandDrivenData(masterPointFacesPtr_);
clearCutFaces();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::enrichedPatch::enrichedPatch
(
const primitiveFacePatch& masterPatch,
const primitiveFacePatch& slavePatch,
const labelList& slavePointPointHits,
const labelList& slavePointEdgeHits,
const List<objectHit>& slavePointFaceHits
)
:
masterPatch_(masterPatch),
slavePatch_(slavePatch),
pointMap_
(
masterPatch_.meshPoints().size()
+ slavePatch_.meshPoints().size()
),
pointMapComplete_(false),
pointMergeMap_(2*slavePatch_.meshPoints().size()),
slavePointPointHits_(slavePointPointHits),
slavePointEdgeHits_(slavePointEdgeHits),
slavePointFaceHits_(slavePointFaceHits),
enrichedFacesPtr_(nullptr),
meshPointsPtr_(nullptr),
localFacesPtr_(nullptr),
localPointsPtr_(nullptr),
pointPointsPtr_(nullptr),
masterPointFacesPtr_(nullptr),
cutFacesPtr_(nullptr),
cutFaceMasterPtr_(nullptr),
cutFaceSlavePtr_(nullptr)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::enrichedPatch::~enrichedPatch()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::labelList& Foam::enrichedPatch::meshPoints() const
{
if (!meshPointsPtr_)
{
calcMeshPoints();
}
return *meshPointsPtr_;
}
const Foam::faceList& Foam::enrichedPatch::localFaces() const
{
if (!localFacesPtr_)
{
calcLocalFaces();
}
return *localFacesPtr_;
}
const Foam::pointField& Foam::enrichedPatch::localPoints() const
{
if (!localPointsPtr_)
{
calcLocalPoints();
}
return *localPointsPtr_;
}
const Foam::labelListList& Foam::enrichedPatch::pointPoints() const
{
if (!pointPointsPtr_)
{
calcPointPoints();
}
return *pointPointsPtr_;
}
bool Foam::enrichedPatch::checkSupport() const
{
const faceList& faces = enrichedFaces();
bool error = false;
forAll(faces, facei)
{
const face& curFace = faces[facei];
forAll(curFace, pointi)
{
if (!pointMap().found(curFace[pointi]))
{
WarningInFunction
<< "Point " << pointi << " of face " << facei
<< " global point index: " << curFace[pointi]
<< " not supported in point map. This is not allowed."
<< endl;
error = true;
}
}
}
return error;
}
void Foam::enrichedPatch::writeOBJ(const fileName& fName) const
{
OFstream str(fName);
const pointField& lp = localPoints();
forAll(lp, pointi)
{
meshTools::writeOBJ(str, lp[pointi]);
}
const faceList& faces = localFaces();
forAll(faces, facei)
{
const face& f = faces[facei];
str << 'f';
forAll(f, fp)
{
str << ' ' << f[fp]+1;
}
str << nl;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -1,296 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::enrichedPatch
Description
The enriched patch contains a double set of faces from the two
sides of the sliding interface before the cutting.
The patch basically consists of two over-lapping sets of faces sitting
on a common point support, where every edge may be shared by more than
2 faces. The patch points are collected in a map. Additional
information needed for cutting is the point insertion into every edge
of master and slave.
Note:
If new points are created during master-slave edge cutting, they
should be registered with the pointMap.
SourceFiles
enrichedPatch.C
enrichedPatchCutFaces.C
enrichedPatchFaces.C
enrichedPatchPointMap.C
enrichedPatchPointMergeMap.C
enrichedPatchPointPoints.C
\*---------------------------------------------------------------------------*/
#ifndef enrichedPatch_H
#define enrichedPatch_H
#include "primitiveFacePatch.H"
#include "Map.H"
#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class enrichedPatch Declaration
\*---------------------------------------------------------------------------*/
class enrichedPatch
{
// Private Data
//- Reference to master patch
const primitiveFacePatch& masterPatch_;
//- Reference to slave patch
const primitiveFacePatch& slavePatch_;
//- Map of points supporting patch faces
mutable Map<point> pointMap_;
//- Has the point map been completed?
mutable bool pointMapComplete_;
//- Map of point merges
mutable Map<label> pointMergeMap_;
//- Slave point point hits
const labelList& slavePointPointHits_;
//- Slave point edge hits
const labelList& slavePointEdgeHits_;
//- Slave point face hits
const List<objectHit>& slavePointFaceHits_;
// Demand-driven private data
//- Enriched patch
mutable faceList* enrichedFacesPtr_;
//- Mesh points
mutable labelList* meshPointsPtr_;
//- Local faces
mutable faceList* localFacesPtr_;
//- Local points
mutable pointField* localPointsPtr_;
//- Point-point addressing
mutable labelListList* pointPointsPtr_;
// Master point addressing
mutable Map<labelList>* masterPointFacesPtr_;
// Cut faces and addressing
//- Cut faces
mutable faceList* cutFacesPtr_;
//- Cut face master
// - the face on the master patch for internal faces
// - the creator face for boundary face
mutable labelList* cutFaceMasterPtr_;
//- Cut face slave
// - the face on the slave patch for internal faces
// - -1 for boundary face
mutable labelList* cutFaceSlavePtr_;
// Private Member Functions
// Creation of demand-driven private data
//- Calculate point merge map
void calcPointMergeMap() const;
//- Complete point map
void completePointMap() const;
//- Calculate mesh points
void calcMeshPoints() const;
//- Calculate local points
void calcLocalPoints() const;
//- Calculate local faces
void calcLocalFaces() const;
//- Calculate point-point addressing
void calcPointPoints() const;
//- Calculate master point addressing
void calcMasterPointFaces() const;
//- Calculate cut faces
void calcCutFaces() const;
//- Clear cut faces
void clearCutFaces();
//- Clear out demand-driven data
void clearOut();
// Static Data Members
//- Estimated ratio of original-to-enriched face size
static const label enrichedFaceRatio_;
//- Estimated number of master face hits by slave points
static const label nFaceHits_;
//- Size of face on which the check is forced
static const label maxFaceSizeDebug_;
public:
// Static Data Members
ClassName("enrichedPatch");
// Constructors
//- Construct from components
enrichedPatch
(
const primitiveFacePatch& masterPatch,
const primitiveFacePatch& slavePatch,
const labelList& slavePointPointHits,
// -1 or common point snapped to
const labelList& slavePointEdgeHits,
// -1 or common edge snapped to
const List<objectHit>& slavePointFaceHits
// master face snapped to
);
//- Disallow default bitwise copy construction
enrichedPatch(const enrichedPatch&) = delete;
//- Destructor
~enrichedPatch();
// Member Functions
// Access
//- Return non-const access to point map to add points
Map<point>& pointMap();
//- Return map of points
const Map<point>& pointMap() const;
//- Return map of point merges
Map<label>& pointMergeMap()
{
return pointMergeMap_;
}
//- Return map of point merges
const Map<label>& pointMergeMap() const
{
return pointMergeMap_;
}
// Topological data
//- Calculate enriched faces
void calcEnrichedFaces
(
const labelListList& pointsIntoMasterEdges,
const labelListList& pointsIntoSlaveEdges,
const pointField& projectedSlavePoints
);
//- Return enriched faces
const faceList& enrichedFaces() const;
//- Return mesh points
const labelList& meshPoints() const;
//- Return local faces
const faceList& localFaces() const;
//- Return local points
const pointField& localPoints() const;
//- Return point-point addressing
const labelListList& pointPoints() const;
//- Master point face addressing
const Map<labelList>& masterPointFaces() const;
// Cut faces
//- Return list of cut faces
const faceList& cutFaces() const;
//- Return cut face master list
const labelList& cutFaceMaster() const;
//- Return cut face slave list
const labelList& cutFaceSlave() const;
//- Check if the patch is fully supported
bool checkSupport() const;
//- Debugging: dump graphical representation to obj format file
void writeOBJ(const fileName&) const;
// Member Operators
//- Disallow default bitwise assignment
void operator=(const enrichedPatch&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,701 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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/>.
Description
Calculating cut faces of the enriched patch, together with the addressing
into master and slave patch.
\*---------------------------------------------------------------------------*/
#include "enrichedPatch.H"
#include "boolList.H"
#include "DynamicList.H"
#include "labelPair.H"
#include "primitiveMesh.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// If the cut face gets more than this number of points, it will be checked
const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Index of debug signs:
// x - skip a point
// l - left turn
// r - right turn
void Foam::enrichedPatch::calcCutFaces() const
{
if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
{
FatalErrorInFunction
<< "Cut faces addressing already calculated."
<< abort(FatalError);
}
const faceList& enFaces = enrichedFaces();
const labelList& mp = meshPoints();
const faceList& lf = localFaces();
const pointField& lp = localPoints();
const labelListList& pp = pointPoints();
// Pout<< "enFaces: " << enFaces << endl;
// Pout<< "lf: " << lf << endl;
// Pout<< "lp: " << lp << endl;
// Pout<< "pp: " << pp << endl;
const Map<labelList>& masterPointFaceAddr = masterPointFaces();
// Prepare the storage
DynamicList<face> cf(2*lf.size());
DynamicList<label> cfMaster(2*lf.size());
DynamicList<label> cfSlave(2*lf.size());
// Algorithm
// Go through all the faces
// 1) For each face, start at point zero and grab the edge zero.
// Both points are added into the cut face. Mark the first edge
// as used and position the current point as the end of the first
// edge and previous point as point zero.
// 2) Grab all possible points from the current point. Excluding
// the previous point find the point giving the biggest right
// turn. Add the point to the face and mark the edges as used.
// Continue doing this until the face is complete, i.e. the next point
// to add is the first point of the face.
// 3) Once the facet is completed, register its owner from the face
// it has been created from (remember that the first lot of faces
// inserted into the enriched faces list are the slaves, to allow
// matching of the other side).
// 4) If the facet is created from an original slave face, go
// through the master patch and try to identify the master face
// this facet belongs to. This is recognised by the fact that the
// faces consists exclusively of the points of the master face and
// the points projected onto the face.
// Create a set of edge usage parameters
HashSet<edge, Hash<edge>> edgesUsedOnce(pp.size());
HashSet<edge, Hash<edge>> edgesUsedTwice
(pp.size()*primitiveMesh::edgesPerPoint_);
forAll(lf, facei)
{
const face& curLocalFace = lf[facei];
const face& curGlobalFace = enFaces[facei];
// Pout<< "Doing face " << facei
// << " local: " << curLocalFace
// << " or " << curGlobalFace
// << endl;
// if (facei < slavePatch_.size())
// {
// Pout<< "original slave: " << slavePatch_[facei]
// << " local: " << slavePatch_.localFaces()[facei] << endl;
// }
// else
// {
// Pout<< "original master: "
// << masterPatch_[facei - slavePatch_.size()] << " "
// << masterPatch_.localFaces()[facei - slavePatch_.size()]
// << endl;
// }
// {
// pointField facePoints = curLocalFace.points(lp);
// forAll(curLocalFace, pointi)
// {
// Pout<< "v " << facePoints[pointi].x() << " "
// << facePoints[pointi].y() << " "
// << facePoints[pointi].z() << endl;
// }
// }
// Track the usage of face edges. When all edges are used, the
// face decomposition is complete.
// The seed edges include all the edges of the original face + all edges
// of other faces that have been used in the construction of the
// facet. Edges from other faces can be considered as
// internal to the current face if used only once.
// Track the edge usage to avoid duplicate faces and reset it to unused
boolList usedFaceEdges(curLocalFace.size(), false);
SLList<edge> edgeSeeds;
// Insert the edges of current face into the seed list.
edgeList cfe = curLocalFace.edges();
forAll(curLocalFace, edgeI)
{
edgeSeeds.append(cfe[edgeI]);
}
// Grab face normal
const vector normal = curLocalFace.normal(lp);
while (edgeSeeds.size())
{
// Pout<< "edgeSeeds.size(): "
// << edgeSeeds.size()
// << endl;
const edge curEdge = edgeSeeds.removeHead();
// Locate the edge in current face
const label curEdgeWhich = curLocalFace.which(curEdge.start());
// Check if the edge is in current face and if it has been
// used already. If so, skip it
if
(
curEdgeWhich > -1
&& curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
)
{
// Edge is in the starting face.
// If unused, mark it as used; if used, skip it
if (usedFaceEdges[curEdgeWhich]) continue;
usedFaceEdges[curEdgeWhich] = true;
}
// If the edge has already been used twice, skip it
if (edgesUsedTwice.found(curEdge)) continue;
// Pout<< "Trying new edge (" << mp[curEdge.start()]
// << ", " << mp[curEdge.end()]
// << ") seed: " << curEdge
// << " used: " << edgesUsedTwice.found(curEdge)
// << endl;
// Estimate the size of cut face as twice the size of original face
DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
// Found unused edge.
label prevPointLabel = curEdge.start();
cutFaceGlobalPoints.append(mp[prevPointLabel]);
cutFaceLocalPoints.append(prevPointLabel);
// Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;
// Grab current point and append it to the list
label curPointLabel = curEdge.end();
point curPoint = lp[curPointLabel];
cutFaceGlobalPoints.append(mp[curPointLabel]);
cutFaceLocalPoints.append(curPointLabel);
bool completedCutFace = false;
label faceSizeDebug = maxFaceSizeDebug_;
do
{
// Grab the next point options
// Pout<< "curPointLabel: " << mp[curPointLabel] << endl;
const labelList& nextPoints = pp[curPointLabel];
// Pout<< "nextPoints: "
// << UIndirectList<label>(mp, nextPoints)
// << endl;
// Get the vector along the edge and the right vector
vector ahead = curPoint - lp[prevPointLabel];
ahead -= normal*(normal & ahead);
ahead /= mag(ahead);
vector right = normal ^ ahead;
right /= mag(right);
// Pout<< "normal: " << normal
// << " ahead: " << ahead
// << " right: " << right
// << endl;
scalar atanTurn = -great;
label bestAtanPoint = -1;
forAll(nextPoints, nextI)
{
// Exclude the point we are coming from; there will always
// be more than one edge, so this is safe
if (nextPoints[nextI] != prevPointLabel)
{
// Pout<< "cur point: " << curPoint
// << " trying for point: "
// << mp[nextPoints[nextI]]
// << " " << lp[nextPoints[nextI]];
vector newDir = lp[nextPoints[nextI]] - curPoint;
// Pout<< " newDir: " << newDir
// << " mag: " << mag(newDir) << flush;
newDir -= normal*(normal & newDir);
scalar magNewDir = mag(newDir);
// Pout<< " corrected: " << newDir
// << " mag: " << mag(newDir) << flush;
if (magNewDir < small)
{
FatalErrorInFunction
<< "projection error: slave patch probably "
<< "does not project onto master. "
<< "Please switch on "
<< "enriched patch debug for more info"
<< abort(FatalError);
}
newDir /= magNewDir;
scalar curAtanTurn =
atan2(newDir & right, newDir & ahead);
// Pout<< " atan: " << curAtanTurn << endl;
if (curAtanTurn > atanTurn)
{
bestAtanPoint = nextPoints[nextI];
atanTurn = curAtanTurn;
}
} // end of prev point skip
} // end of next point selection
// Pout<< " bestAtanPoint: " << bestAtanPoint << " or "
// << mp[bestAtanPoint]
// << endl;
// Selected next best point.
// Pout<< "cutFaceGlobalPoints: "
// << cutFaceGlobalPoints
// << endl;
// Check if the edge about to be added has been used
// in the current face or twice in other faces. If
// so, the face is bad.
const label whichNextPoint = curLocalFace.which(curPointLabel);
if
(
whichNextPoint > -1
&& curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
&& usedFaceEdges[whichNextPoint]
)
{
// This edge is already used in current face
// face cannot be good; start on a new one
// Pout<< "Double usage in current face, cannot be good"
// << endl;
completedCutFace = true;
}
if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
{
// This edge is already used -
// face cannot be good; start on a new one
completedCutFace = true;
// Pout<< "Double usage elsewhere, cannot be good" << endl;
}
if (completedCutFace) continue;
// Insert the next best point into the list
if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
{
// Face is completed. Save it and move on to the next
// available edge
completedCutFace = true;
if (debug)
{
Pout<< " local: " << cutFaceLocalPoints
<< " one side: " << facei;
}
// Append the face
face cutFaceGlobal;
cutFaceGlobal.transfer(cutFaceGlobalPoints);
cf.append(cutFaceGlobal);
face cutFaceLocal;
cutFaceLocal.transfer(cutFaceLocalPoints);
// Pout<< "\ncutFaceLocal: "
// << cutFaceLocal.points(lp)
// << endl;
// Go through all edges of the cut faces.
// If the edge corresponds to a starting face edge,
// mark the starting face edge as true
forAll(cutFaceLocal, cutI)
{
const edge curCutFaceEdge
(
cutFaceLocal[cutI],
cutFaceLocal.nextLabel(cutI)
);
// Increment the usage count using two hash sets
HashSet<edge, Hash<edge>>::iterator euoIter =
edgesUsedOnce.find(curCutFaceEdge);
if (euoIter == edgesUsedOnce.end())
{
// Pout<< "Found edge not used before: "
// << curCutFaceEdge
// << endl;
edgesUsedOnce.insert(curCutFaceEdge);
}
else
{
// Pout<< "Found edge used once: "
// << curCutFaceEdge
// << endl;
edgesUsedOnce.erase(euoIter);
edgesUsedTwice.insert(curCutFaceEdge);
}
const label curCutFaceEdgeWhich = curLocalFace.which
(
curCutFaceEdge.start()
);
if
(
curCutFaceEdgeWhich > -1
&& curLocalFace.nextLabel(curCutFaceEdgeWhich)
== curCutFaceEdge.end()
)
{
// Found edge in original face
// Pout<< "Found edge in orig face: "
// << curCutFaceEdge << ": "
// << curCutFaceEdgeWhich
// << endl;
usedFaceEdges[curCutFaceEdgeWhich] = true;
}
else
{
// Edge not in original face. Add it to seeds
// Pout<< "Found new edge seed: "
// << curCutFaceEdge
// << endl;
edgeSeeds.append(curCutFaceEdge.reverseEdge());
}
}
// Find out what the other side is
// Algorithm
// If the face is in the slave half of the
// enrichedFaces lists, it may be matched against
// the master face. It can be recognised by the
// fact that all its points belong to one master
// face. For this purpose:
// 1) Grab the master faces around the point zero
// of the cut face and collect all master faces
// around it.
// 2) Loop through the rest of cut face points and
// if the point refers to one of the faces marked
// by point zero, increment its count.
// 3) When completed, the face whose count is
// equal to the number of points in the cut face
// is the other side. If this is not the case, there is no
// face on the other side.
if (facei < slavePatch_.size())
{
Map<labelList>::const_iterator mpfAddrIter =
masterPointFaceAddr.find(cutFaceGlobal[0]);
bool otherSideFound = false;
if (mpfAddrIter != masterPointFaceAddr.end())
{
bool miss = false;
// Create the label count using point zero
const labelList& masterFacesOfPZero = mpfAddrIter();
labelList hits(masterFacesOfPZero.size(), 1);
for
(
label pointi = 1;
pointi < cutFaceGlobal.size();
pointi++
)
{
Map<labelList>::const_iterator
mpfAddrPointIter =
masterPointFaceAddr.find
(
cutFaceGlobal[pointi]
);
if
(
mpfAddrPointIter
== masterPointFaceAddr.end()
)
{
// Point is off the master patch. Skip
miss = true;
break;
}
const labelList& curMasterFaces =
mpfAddrPointIter();
// For every current face, try to find it in the
// zero-list
forAll(curMasterFaces, i)
{
forAll(masterFacesOfPZero, j)
{
if
(
curMasterFaces[i]
== masterFacesOfPZero[j]
)
{
hits[j]++;
break;
}
}
}
}
// If all point are found attempt matching
if (!miss)
{
forAll(hits, pointi)
{
if (hits[pointi] == cutFaceGlobal.size())
{
// Found other side.
otherSideFound = true;
cfMaster.append
(
masterFacesOfPZero[pointi]
);
cfSlave.append(facei);
// Reverse the face such that it
// points out of the master patch
cf.last().flip();
if (debug)
{
Pout<< " other side: "
<< masterFacesOfPZero[pointi]
<< endl;
}
} // end of hits
} // end of for all hits
} // end of not miss
if (!otherSideFound || miss)
{
if (debug)
{
Pout<< " solo slave A" << endl;
}
cfMaster.append(-1);
cfSlave.append(facei);
}
}
else
{
// First point not in master patch
if (debug)
{
Pout<< " solo slave B" << endl;
}
cfMaster.append(-1);
cfSlave.append(facei);
}
}
else
{
if (debug)
{
Pout<< " master side" << endl;
}
cfMaster.append(facei - slavePatch_.size());
cfSlave.append(-1);
}
}
else
{
// Face not completed. Prepare for the next point search
prevPointLabel = curPointLabel;
curPointLabel = bestAtanPoint;
curPoint = lp[curPointLabel];
cutFaceGlobalPoints.append(mp[curPointLabel]);
cutFaceLocalPoints.append(curPointLabel);
if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
{
faceSizeDebug *= 2;
// Check for duplicate points in the face
forAll(cutFaceGlobalPoints, checkI)
{
for
(
label checkJ = checkI + 1;
checkJ < cutFaceGlobalPoints.size();
checkJ++
)
{
if
(
cutFaceGlobalPoints[checkI]
== cutFaceGlobalPoints[checkJ]
)
{
// Shrink local points for debugging
cutFaceLocalPoints.shrink();
face origFace;
face origFaceLocal;
if (facei < slavePatch_.size())
{
origFace = slavePatch_[facei];
origFaceLocal =
slavePatch_.localFaces()[facei];
}
else
{
origFace =
masterPatch_
[facei - slavePatch_.size()];
origFaceLocal =
masterPatch_.localFaces()
[facei - slavePatch_.size()];
}
FatalErrorInFunction
<< "Duplicate point found in cut face. "
<< "Error in the face cutting "
<< "algorithm for global face "
<< origFace << " local face "
<< origFaceLocal << nl
<< "Slave size: " << slavePatch_.size()
<< " Master size: "
<< masterPatch_.size()
<< " index: " << facei << ".\n"
<< "Face: " << curGlobalFace << nl
<< "Cut face: " << cutFaceGlobalPoints
<< " local: " << cutFaceLocalPoints
<< nl << "Points: "
<< face(cutFaceLocalPoints).points(lp)
<< abort(FatalError);
}
}
}
} // end of debug
}
} while (!completedCutFace);
} // end of used edges
if (debug)
{
Pout<< " Finished face " << facei << endl;
}
} // end of local faces
// Re-pack the list into compact storage
cutFacesPtr_ = new faceList();
cutFacesPtr_->transfer(cf);
cutFaceMasterPtr_ = new labelList();
cutFaceMasterPtr_->transfer(cfMaster);
cutFaceSlavePtr_ = new labelList();
cutFaceSlavePtr_->transfer(cfSlave);
}
void Foam::enrichedPatch::clearCutFaces()
{
deleteDemandDrivenData(cutFacesPtr_);
deleteDemandDrivenData(cutFaceMasterPtr_);
deleteDemandDrivenData(cutFaceSlavePtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::faceList& Foam::enrichedPatch::cutFaces() const
{
if (!cutFacesPtr_)
{
calcCutFaces();
}
return *cutFacesPtr_;
}
const Foam::labelList& Foam::enrichedPatch::cutFaceMaster() const
{
if (!cutFaceMasterPtr_)
{
calcCutFaces();
}
return *cutFaceMasterPtr_;
}
const Foam::labelList& Foam::enrichedPatch::cutFaceSlave() const
{
if (!cutFaceSlavePtr_)
{
calcCutFaces();
}
return *cutFaceSlavePtr_;
}
// ************************************************************************* //

View File

@ -1,398 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "enrichedPatch.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::label Foam::enrichedPatch::enrichedFaceRatio_ = 3;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::enrichedPatch::calcEnrichedFaces
(
const labelListList& pointsIntoMasterEdges,
const labelListList& pointsIntoSlaveEdges,
const pointField& projectedSlavePoints
)
{
if (enrichedFacesPtr_)
{
FatalErrorInFunction
<< "Enriched faces already calculated."
<< abort(FatalError);
}
// Create a list of enriched faces
// Algorithm:
// 1) Grab the original face and start from point zero.
// 2) If the point has been merged away, grab the merge label;
// otherwise, keep the original label.
// 3) Go to the next edge. Collect all the points to be added along
// the edge; order them in the necessary direction and insert onto the
// face.
// 4) Grab the next point and return on step 2.
enrichedFacesPtr_ = new faceList(masterPatch_.size() + slavePatch_.size());
faceList& enrichedFaces = *enrichedFacesPtr_;
label nEnrichedFaces = 0;
const pointField& masterLocalPoints = masterPatch_.localPoints();
const faceList& masterLocalFaces = masterPatch_.localFaces();
const labelListList& masterFaceEdges = masterPatch_.faceEdges();
const faceList& slaveLocalFaces = slavePatch_.localFaces();
const labelListList& slaveFaceEdges = slavePatch_.faceEdges();
// For correct functioning of the enrichedPatch class, the slave
// faces need to be inserted first. See comments in
// enrichedPatch.H
// Get reference to the point merge map
const Map<label>& pmm = pointMergeMap();
// Add slave faces into the enriched faces list
forAll(slavePatch_, facei)
{
const face oldFace = slavePatch_[facei];
const face oldLocalFace = slaveLocalFaces[facei];
// Info<< "old slave face " << facei << ": " << oldFace << endl;
const labelList& curEdges = slaveFaceEdges[facei];
DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
// Note: The number of points and edges in a face is always identical
// so both can be done is the same loop
forAll(oldFace, i)
{
// Add the point
Map<label>::const_iterator mpIter =
pmm.find(oldFace[i]);
if (mpIter == pmm.end())
{
// Point not mapped
newFace.append(oldFace[i]);
// Add the projected point into the patch support
pointMap().insert
(
oldFace[i], // Global label of point
projectedSlavePoints[oldLocalFace[i]] // Projected position
);
}
else
{
// Point mapped
newFace.append(mpIter());
// Add the projected point into the patch support
pointMap().insert
(
mpIter(), // Merged global label of point
projectedSlavePoints[oldLocalFace[i]] // Projected position
);
}
// Grab the edge points
const labelList& slavePointsOnEdge =
pointsIntoSlaveEdges[curEdges[i]];
// Info<< "slavePointsOnEdge for "
// << curEdges[i] << ": " << slavePointsOnEdge
// << endl;
// If there are no points on the edge, skip everything
// If there is only one point, no need for sorting
if (slavePointsOnEdge.size())
{
// Sort edge points in order
scalarField edgePointWeights(slavePointsOnEdge.size());
const point& startPoint = projectedSlavePoints[oldLocalFace[i]];
vector e =
projectedSlavePoints[oldLocalFace.nextLabel(i)]
- startPoint;
scalar magSqrE = magSqr(e);
if (magSqrE > small)
{
e /= magSqrE;
}
else
{
FatalErrorInFunction
<< "Zero length edge in slave patch for face " << i
<< ". This is not allowed."
<< abort(FatalError);
}
pointField slavePosOnEdge(slavePointsOnEdge.size());
forAll(slavePointsOnEdge, edgePointi)
{
slavePosOnEdge[edgePointi] =
pointMap().find(slavePointsOnEdge[edgePointi])();
edgePointWeights[edgePointi] =
(e & (slavePosOnEdge[edgePointi] - startPoint));
}
if (debug)
{
// Check weights: all new points should be on the edge
if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
{
FatalErrorInFunction
<< " not on the edge for edge " << curEdges[i]
<< " of face " << facei << " in slave patch." << nl
<< "Min weight: " << min(edgePointWeights)
<< " Max weight: " << max(edgePointWeights)
<< abort(FatalError);
}
}
// Go through the points and collect them based on
// weights from lower to higher. This gives the
// correct order of points along the edge.
forAll(edgePointWeights, passI)
{
// Max weight can only be one, so the sorting is
// done by elimination.
label nextPoint = -1;
scalar dist = 2;
forAll(edgePointWeights, wI)
{
if (edgePointWeights[wI] < dist)
{
dist = edgePointWeights[wI];
nextPoint = wI;
}
}
// Insert the next point and reset its weight to exclude it
// from future picks
newFace.append(slavePointsOnEdge[nextPoint]);
edgePointWeights[nextPoint] = great;
// Add the point into patch support
pointMap().insert
(
slavePointsOnEdge[nextPoint],
slavePosOnEdge[nextPoint]
);
}
}
}
// Info<< "New slave face " << facei << ": " << newFace << endl;
// Add the new face to the list
enrichedFaces[nEnrichedFaces].transfer(newFace);
nEnrichedFaces++;
}
// Add master faces into the enriched faces list
forAll(masterPatch_, facei)
{
const face& oldFace = masterPatch_[facei];
const face& oldLocalFace = masterLocalFaces[facei];
// Info<< "old master face: " << oldFace << endl;
const labelList& curEdges = masterFaceEdges[facei];
DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
// Note: The number of points and edges in a face is always identical
// so both can be done is the same loop
forAll(oldFace, i)
{
// Add the point
Map<label>::const_iterator mpIter =
pmm.find(oldFace[i]);
if (mpIter == pmm.end())
{
// Point not mapped
newFace.append(oldFace[i]);
// Add the point into patch support
pointMap().insert
(
oldFace[i],
masterLocalPoints[oldLocalFace[i]]
);
}
else
{
// Point mapped
newFace.append(mpIter());
// Add the point into support
pointMap().insert(mpIter(), masterLocalPoints[oldLocalFace[i]]);
}
// Grab the edge points
const labelList& masterPointsOnEdge =
pointsIntoMasterEdges[curEdges[i]];
// If there are no points on the edge, skip everything
// If there is only one point, no need for sorting
if (masterPointsOnEdge.size())
{
// Sort edge points in order
scalarField edgePointWeights(masterPointsOnEdge.size());
const point& startPoint = masterLocalPoints[oldLocalFace[i]];
vector e =
masterLocalPoints[oldLocalFace.nextLabel(i)]
- startPoint;
scalar magSqrE = magSqr(e);
if (magSqrE > small)
{
e /= magSqrE;
}
else
{
FatalErrorInFunction
<< "Zero length edge in master patch for face " << i
<< ". This is not allowed."
<< abort(FatalError);
}
pointField masterPosOnEdge(masterPointsOnEdge.size());
forAll(masterPointsOnEdge, edgePointi)
{
masterPosOnEdge[edgePointi] =
pointMap().find(masterPointsOnEdge[edgePointi])();
edgePointWeights[edgePointi] =
(e & (masterPosOnEdge[edgePointi] - startPoint));
}
if (debug)
{
// Check weights: all new points should be on the edge
if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
{
FatalErrorInFunction
<< " not on the edge for edge " << curEdges[i]
<< " of face " << facei << " in master patch." << nl
<< "Min weight: " << min(edgePointWeights)
<< " Max weight: " << max(edgePointWeights)
<< abort(FatalError);
}
}
// Go through the points and collect them based on
// weights from lower to higher. This gives the
// correct order of points along the edge.
forAll(edgePointWeights, passI)
{
// Max weight can only be one, so the sorting is
// done by elimination.
label nextPoint = -1;
scalar dist = 2;
forAll(edgePointWeights, wI)
{
if (edgePointWeights[wI] < dist)
{
dist = edgePointWeights[wI];
nextPoint = wI;
}
}
// Insert the next point and reset its weight to exclude it
// from future picks
newFace.append(masterPointsOnEdge[nextPoint]);
edgePointWeights[nextPoint] = great;
// Add the point into patch support
pointMap().insert
(
masterPointsOnEdge[nextPoint],
masterPosOnEdge[nextPoint]
);
}
}
}
// Info<< "New master face: " << newFace << endl;
// Add the new face to the list
enrichedFaces[nEnrichedFaces].transfer(newFace);
nEnrichedFaces++;
}
// Check the support for the enriched patch
if (debug)
{
if (!checkSupport())
{
Info<< "Enriched patch support OK. Slave faces: "
<< slavePatch_.size() << " Master faces: "
<< masterPatch_.size() << endl;
}
else
{
FatalErrorInFunction
<< "Error in enriched patch support"
<< abort(FatalError);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::faceList& Foam::enrichedPatch::enrichedFaces() const
{
if (!enrichedFacesPtr_)
{
FatalErrorInFunction
<< "void enrichedPatch::calcEnrichedFaces\n"
<< "(\n"
<< " const labelListList& pointsIntoMasterEdges,\n"
<< " const labelListList& pointsIntoSlaveEdges,\n"
<< " const pointField& projectedSlavePoints\n"
<< ")"
<< " before trying to access faces."
<< abort(FatalError);
}
return *enrichedFacesPtr_;
}
// ************************************************************************* //

View File

@ -1,159 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "enrichedPatch.H"
#include "primitiveMesh.H"
#include "demandDrivenData.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::label Foam::enrichedPatch::nFaceHits_ = 4;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::enrichedPatch::calcMasterPointFaces() const
{
if (masterPointFacesPtr_)
{
FatalErrorInFunction
<< "Master point face addressing already calculated."
<< abort(FatalError);
}
// Note:
// Master point face addressing lists the master faces for all points
// in the enriched patch support (if there are no master faces, which is
// normal, the list will be empty). The index represents the index of
// the master face rather than the index from the enriched patch
// Master face points lists the points of the enriched master face plus
// points projected into the master face
Map<DynamicList<label>> mpf(meshPoints().size());
const faceList& ef = enrichedFaces();
// Add the original face points
forAll(masterPatch_, facei)
{
const face& curFace = ef[facei + slavePatch_.size()];
// Pout<< "Cur face in pfAddr: " << curFace << endl;
forAll(curFace, pointi)
{
Map<DynamicList<label>>::iterator mpfIter =
mpf.find(curFace[pointi]);
if (mpfIter == mpf.end())
{
// Not found, add new dynamic list
mpf.insert
(
curFace[pointi],
DynamicList<label>(primitiveMesh::facesPerPoint_)
);
// Iterator is invalidated - have to find again
mpf.find(curFace[pointi])().append(facei);
}
else
{
mpfIter().append(facei);
}
}
}
// Add the projected points which hit the face
const labelList& slaveMeshPoints = slavePatch_.meshPoints();
forAll(slavePointFaceHits_, pointi)
{
if
(
slavePointPointHits_[pointi] < 0
&& slavePointEdgeHits_[pointi] < 0
&& slavePointFaceHits_[pointi].hit()
)
{
// Get the index of projected point corresponding to this slave
// point
const label mergedSmp =
pointMergeMap().find(slaveMeshPoints[pointi])();
Map<DynamicList<label>>::iterator mpfIter =
mpf.find(mergedSmp);
if (mpfIter == mpf.end())
{
// Not found, add new dynamic list
mpf.insert
(
mergedSmp,
DynamicList<label>(primitiveMesh::facesPerPoint_)
);
// Iterator is invalidated - have to find again
mpf.find(mergedSmp)().append
(
slavePointFaceHits_[pointi].hitObject()
);
}
else
{
mpfIter().append(slavePointFaceHits_[pointi].hitObject());
}
}
}
// Re-pack dynamic lists into normal lists
const labelList mpfToc = mpf.toc();
masterPointFacesPtr_ = new Map<labelList>(2*mpfToc.size());
Map<labelList>& masterPointFaceAddr = *masterPointFacesPtr_;
forAll(mpfToc, mpfTocI)
{
labelList l;
l.transfer(mpf.find(mpfToc[mpfTocI])());
masterPointFaceAddr.insert(mpfToc[mpfTocI], l);
}
// Pout<< "masterPointFaceAddr: " << masterPointFaceAddr << endl;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::Map<Foam::labelList>& Foam::enrichedPatch::masterPointFaces() const
{
if (!masterPointFacesPtr_)
{
calcMasterPointFaces();
}
return *masterPointFacesPtr_;
}
// ************************************************************************* //

View File

@ -1,105 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "enrichedPatch.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::enrichedPatch::completePointMap() const
{
if (pointMapComplete_)
{
FatalErrorInFunction
<< "Point map already completed"
<< abort(FatalError);
}
pointMapComplete_ = true;
const Map<label>& pmm = pointMergeMap();
// Get the mesh points for both patches. If the point has not been
// merged away, add it to the map
// Do master patch
const labelList& masterMeshPoints = masterPatch_.meshPoints();
const pointField& masterLocalPoints = masterPatch_.localPoints();
forAll(masterMeshPoints, pointi)
{
if (!pmm.found(masterMeshPoints[pointi]))
{
pointMap_.insert
(
masterMeshPoints[pointi],
masterLocalPoints[pointi]
);
}
}
// Do slave patch
const labelList& slaveMeshPoints = slavePatch_.meshPoints();
const pointField& slaveLocalPoints = slavePatch_.localPoints();
forAll(slaveMeshPoints, pointi)
{
if (!pmm.found(slaveMeshPoints[pointi]))
{
pointMap_.insert
(
slaveMeshPoints[pointi],
slaveLocalPoints[pointi]
);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap()
{
if (!pointMapComplete_)
{
completePointMap();
}
return pointMap_;
}
const Foam::Map<Foam::point>& Foam::enrichedPatch::pointMap() const
{
if (!pointMapComplete_)
{
completePointMap();
}
return pointMap_;
}
// ************************************************************************* //

View File

@ -1,115 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "enrichedPatch.H"
#include "primitiveMesh.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::enrichedPatch::calcPointPoints() const
{
// Calculate point-point addressing
if (pointPointsPtr_)
{
FatalErrorInFunction
<< "Point-point addressing already calculated."
<< abort(FatalError);
}
// Algorithm:
// Go through all faces and add the previous and next point as the
// neighbour for each point. While inserting points, reject the
// duplicates (as every internal edge will be visited twice).
List<DynamicList<label, primitiveMesh::edgesPerPoint_>>
pp(meshPoints().size());
const faceList& lf = localFaces();
bool found = false;
forAll(lf, facei)
{
const face& curFace = lf[facei];
forAll(curFace, pointi)
{
DynamicList<label, primitiveMesh::edgesPerPoint_>&
curPp = pp[curFace[pointi]];
// Do next label
label next = curFace.nextLabel(pointi);
found = false;
forAll(curPp, i)
{
if (curPp[i] == next)
{
found = true;
break;
}
}
if (!found)
{
curPp.append(next);
}
// Do previous label
label prev = curFace.prevLabel(pointi);
found = false;
forAll(curPp, i)
{
if (curPp[i] == prev)
{
found = true;
break;
}
}
if (!found)
{
curPp.append(prev);
}
}
}
// Re-pack the list
pointPointsPtr_ = new labelListList(pp.size());
labelListList& ppAddr = *pointPointsPtr_;
forAll(pp, pointi)
{
ppAddr[pointi].transfer(pp[pointi]);
}
}
// ************************************************************************* //

View File

@ -1,639 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "slidingInterface.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
#include "addToRunTimeSelectionTable.H"
#include "plane.H"
// Index of debug signs:
// p - adjusting a projection point
// * - adjusting edge intersection
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(slidingInterface, 0);
template<>
const char* Foam::NamedEnum
<
Foam::slidingInterface::typeOfMatch,
2
>::names[] =
{
"integral",
"partial"
};
}
const Foam::NamedEnum<Foam::slidingInterface::typeOfMatch, 2>
Foam::slidingInterface::typeOfMatchNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::slidingInterface::checkDefinition()
{
const polyMesh& mesh = this->mesh();
if
(
!masterFaceZoneIndex_.active()
|| !slaveFaceZoneIndex_.active()
|| !cutPointZoneIndex_.active()
|| !cutFaceZoneIndex_.active()
|| !masterPatchIndex_.active()
|| !slavePatchIndex_.active()
)
{
FatalErrorInFunction
<< "Not all zones and patches needed in the definition "
<< "have been found. Please check your mesh definition."
<< abort(FatalError);
}
// Check the sizes and set up state
if
(
mesh.faceZones()[masterFaceZoneIndex_.index()].empty()
|| mesh.faceZones()[slaveFaceZoneIndex_.index()].empty()
)
{
FatalErrorInFunction
<< "Please check your mesh definition."
<< abort(FatalError);
}
if (debug)
{
Pout<< "Sliding interface object " << name() << " :" << nl
<< " master face zone: " << masterFaceZoneIndex_.index() << nl
<< " slave face zone: " << slaveFaceZoneIndex_.index() << endl;
}
}
void Foam::slidingInterface::clearOut() const
{
clearPointProjection();
clearAttachedAddressing();
clearAddressing();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::slidingInterface::slidingInterface
(
const word& name,
const polyMesh& mesh,
const word& masterFaceZoneName,
const word& slaveFaceZoneName,
const word& cutPointZoneName,
const word& cutFaceZoneName,
const word& masterPatchName,
const word& slavePatchName,
const typeOfMatch tom,
const bool coupleDecouple,
const intersection::algorithm algo
)
:
polyMeshModifier(name, mesh),
masterFaceZoneIndex_
(
masterFaceZoneName,
mesh.faceZones()
),
slaveFaceZoneIndex_
(
slaveFaceZoneName,
mesh.faceZones()
),
cutPointZoneIndex_
(
cutPointZoneName,
mesh.pointZones()
),
cutFaceZoneIndex_
(
cutFaceZoneName,
mesh.faceZones()
),
masterPatchIndex_
(
masterPatchName,
mesh.boundaryMesh()
),
slavePatchIndex_
(
slavePatchName,
mesh.boundaryMesh()
),
matchType_(tom),
coupleDecouple_(coupleDecouple),
attached_(false),
projectionAlgo_(algo),
trigger_(false),
pointMergeTol_(pointMergeTolDefault_),
edgeMergeTol_(edgeMergeTolDefault_),
nFacesPerSlaveEdge_(nFacesPerSlaveEdgeDefault_),
edgeFaceEscapeLimit_(edgeFaceEscapeLimitDefault_),
integralAdjTol_(integralAdjTolDefault_),
edgeMasterCatchFraction_(edgeMasterCatchFractionDefault_),
edgeCoPlanarTol_(edgeCoPlanarTolDefault_),
edgeEndCutoffTol_(edgeEndCutoffTolDefault_),
cutFaceMasterPtr_(nullptr),
cutFaceSlavePtr_(nullptr),
masterFaceCellsPtr_(nullptr),
slaveFaceCellsPtr_(nullptr),
masterStickOutFacesPtr_(nullptr),
slaveStickOutFacesPtr_(nullptr),
retiredPointMapPtr_(nullptr),
cutPointEdgePairMapPtr_(nullptr),
slavePointPointHitsPtr_(nullptr),
slavePointEdgeHitsPtr_(nullptr),
slavePointFaceHitsPtr_(nullptr),
masterPointEdgeHitsPtr_(nullptr),
projectedSlavePointsPtr_(nullptr)
{
checkDefinition();
if (attached_)
{
FatalErrorInFunction
<< "Creation of a sliding interface from components "
<< "in attached state not supported."
<< abort(FatalError);
}
else
{
calcAttachedAddressing();
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::slidingInterface::~slidingInterface()
{
clearOut();
}
void Foam::slidingInterface::clearAddressing() const
{
deleteDemandDrivenData(cutFaceMasterPtr_);
deleteDemandDrivenData(cutFaceSlavePtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::faceZoneDynamicID&
Foam::slidingInterface::masterFaceZoneIndex() const
{
return masterFaceZoneIndex_;
}
const Foam::faceZoneDynamicID&
Foam::slidingInterface::slaveFaceZoneIndex() const
{
return slaveFaceZoneIndex_;
}
bool Foam::slidingInterface::changeTopology() const
{
if (coupleDecouple_)
{
// Always changes. If not attached, project points
if (debug)
{
Pout<< "bool slidingInterface::changeTopology() const "
<< "for object " << name() << " : "
<< "Couple-decouple mode." << endl;
}
if (!attached_)
{
projectPoints();
}
else
{
}
return true;
}
if
(
attached_
&& !mesh().changing()
)
{
// If the mesh is not moving or morphing and the interface is
// already attached, the topology will not change
return false;
}
else
{
// Check if the motion changes point projection
return projectPoints();
}
}
void Foam::slidingInterface::setRefinement(polyTopoChange& ref) const
{
if (coupleDecouple_)
{
if (attached_)
{
// Attached, coupling
decoupleInterface(ref);
}
else
{
// Detached, coupling
coupleInterface(ref);
}
return;
}
if (trigger_)
{
if (attached_)
{
// Clear old coupling data
clearCouple(ref);
}
coupleInterface(ref);
trigger_ = false;
}
}
void Foam::slidingInterface::modifyMotionPoints(pointField& motionPoints) const
{
if (debug)
{
Pout<< "void slidingInterface::modifyMotionPoints("
<< "pointField& motionPoints) const for object " << name() << " : "
<< "Adjusting motion points." << endl;
}
const polyMesh& mesh = this->mesh();
// Get point from the cut zone
const labelList& cutPoints = mesh.pointZones()[cutPointZoneIndex_.index()];
if (cutPoints.size() && !projectedSlavePointsPtr_)
{
return;
}
else
{
const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
const Map<label>& rpm = retiredPointMap();
const Map<Pair<edge>>& cpepm = cutPointEdgePairMap();
const Map<label>& slaveZonePointMap =
mesh.faceZones()[slaveFaceZoneIndex_.index()]().meshPointMap();
const primitiveFacePatch& masterPatch =
mesh.faceZones()[masterFaceZoneIndex_.index()]();
const edgeList& masterEdges = masterPatch.edges();
const pointField& masterLocalPoints = masterPatch.localPoints();
const primitiveFacePatch& slavePatch =
mesh.faceZones()[slaveFaceZoneIndex_.index()]();
const edgeList& slaveEdges = slavePatch.edges();
const pointField& slaveLocalPoints = slavePatch.localPoints();
const vectorField& slavePointNormals = slavePatch.pointNormals();
forAll(cutPoints, pointi)
{
// Try to find the cut point in retired points
Map<label>::const_iterator rpmIter = rpm.find(cutPoints[pointi]);
if (rpmIter != rpm.end())
{
if (debug)
{
Pout<< "p";
}
// Cut point is a retired point
motionPoints[cutPoints[pointi]] =
projectedSlavePoints[slaveZonePointMap.find(rpmIter())()];
}
else
{
// A cut point is not a projected slave point. Therefore, it
// must be an edge-to-edge intersection.
Map<Pair<edge>>::const_iterator cpepmIter =
cpepm.find(cutPoints[pointi]);
if (cpepmIter != cpepm.end())
{
// Pout<< "Need to re-create hit for point "
// << cutPoints[pointi]
// << " lookup: " << cpepmIter()
// << endl;
// Note.
// The edge cutting code is repeated in
// slidingInterface::coupleInterface. This is done for
// efficiency reasons and avoids multiple creation of
// cutting planes. Please update both simultaneously.
//
const edge& globalMasterEdge = cpepmIter().first();
const label curMasterEdgeIndex =
masterPatch.whichEdge
(
edge
(
masterPatch.whichPoint
(
globalMasterEdge.start()
),
masterPatch.whichPoint
(
globalMasterEdge.end()
)
)
);
const edge& cme = masterEdges[curMasterEdgeIndex];
// Pout<< "curMasterEdgeIndex: " << curMasterEdgeIndex
// << " cme: " << cme
// << endl;
const edge& globalSlaveEdge = cpepmIter().second();
const label curSlaveEdgeIndex =
slavePatch.whichEdge
(
edge
(
slavePatch.whichPoint
(
globalSlaveEdge.start()
),
slavePatch.whichPoint
(
globalSlaveEdge.end()
)
)
);
const edge& curSlaveEdge = slaveEdges[curSlaveEdgeIndex];
// Pout<< "curSlaveEdgeIndex: " << curSlaveEdgeIndex
// << " curSlaveEdge: " << curSlaveEdge
// << endl;
const point& a = projectedSlavePoints[curSlaveEdge.start()];
const point& b = projectedSlavePoints[curSlaveEdge.end()];
point c =
0.5*
(
slaveLocalPoints[curSlaveEdge.start()]
+ slavePointNormals[curSlaveEdge.start()]
+ slaveLocalPoints[curSlaveEdge.end()]
+ slavePointNormals[curSlaveEdge.end()]
);
// Create the plane
plane cutPlane(a, b, c);
linePointRef curSlaveLine =
curSlaveEdge.line(slaveLocalPoints);
const scalar curSlaveLineMag = curSlaveLine.mag();
scalar cutOnMaster =
cutPlane.lineIntersect
(
cme.line(masterLocalPoints)
);
if
(
cutOnMaster > edgeEndCutoffTol_
&& cutOnMaster < 1.0 - edgeEndCutoffTol_
)
{
// Master is cut, check the slave
point masterCutPoint =
masterLocalPoints[cme.start()]
+ cutOnMaster*cme.vec(masterLocalPoints);
pointHit slaveCut =
curSlaveLine.nearestDist(masterCutPoint);
if (slaveCut.hit())
{
// Strict checking of slave cut to avoid capturing
// end points.
scalar cutOnSlave =
(
(
slaveCut.hitPoint()
- curSlaveLine.start()
) & curSlaveLine.vec()
)/sqr(curSlaveLineMag);
// Calculate merge tolerance from the
// target edge length
scalar mergeTol =
edgeCoPlanarTol_*mag(b - a);
if
(
cutOnSlave > edgeEndCutoffTol_
&& cutOnSlave < 1.0 - edgeEndCutoffTol_
&& slaveCut.distance() < mergeTol
)
{
// Cut both master and slave.
motionPoints[cutPoints[pointi]] =
masterCutPoint;
}
}
else
{
Pout<< "Missed slave edge!!! This is an error. "
<< "Master edge: "
<< cme.line(masterLocalPoints)
<< " slave edge: " << curSlaveLine
<< " point: " << masterCutPoint
<< " weight: " <<
(
(
slaveCut.missPoint()
- curSlaveLine.start()
) & curSlaveLine.vec()
)/sqr(curSlaveLineMag)
<< endl;
}
}
else
{
Pout<< "Missed master edge!!! This is an error"
<< endl;
}
}
else
{
FatalErrorInFunction
<< "Cut point " << cutPoints[pointi]
<< " not recognised as either the projected "
<< "or as intersection point. Error in point "
<< "projection or data mapping"
<< abort(FatalError);
}
}
}
if (debug)
{
Pout<< endl;
}
}
}
void Foam::slidingInterface::topoChange(const polyTopoChangeMap& m)
{
if (debug)
{
Pout<< "void slidingInterface::topoChange(const polyTopoChangeMap& m)"
<< " const for object " << name() << " : "
<< "Updating topology." << endl;
}
// Mesh has changed topologically. Update local topological data
const polyMesh& mesh = this->mesh();
masterFaceZoneIndex_.update(mesh.faceZones());
slaveFaceZoneIndex_.update(mesh.faceZones());
cutPointZoneIndex_.update(mesh.pointZones());
cutFaceZoneIndex_.update(mesh.faceZones());
masterPatchIndex_.update(mesh.boundaryMesh());
slavePatchIndex_.update(mesh.boundaryMesh());
//MJ:Disabled updating
// if (!attached())
// {
// calcAttachedAddressing();
// }
// else
// {
// renumberAttachedAddressing(m);
// }
}
const Foam::pointField& Foam::slidingInterface::pointProjection() const
{
if (!projectedSlavePointsPtr_)
{
projectPoints();
}
return *projectedSlavePointsPtr_;
}
void Foam::slidingInterface::setTolerances(const dictionary&dict, bool report)
{
pointMergeTol_ = dict.lookupOrDefault<scalar>
(
"pointMergeTol",
pointMergeTol_
);
edgeMergeTol_ = dict.lookupOrDefault<scalar>
(
"edgeMergeTol",
edgeMergeTol_
);
nFacesPerSlaveEdge_ = dict.lookupOrDefault<label>
(
"nFacesPerSlaveEdge",
nFacesPerSlaveEdge_
);
edgeFaceEscapeLimit_ = dict.lookupOrDefault<label>
(
"edgeFaceEscapeLimit",
edgeFaceEscapeLimit_
);
integralAdjTol_ = dict.lookupOrDefault<scalar>
(
"integralAdjTol",
integralAdjTol_
);
edgeMasterCatchFraction_ = dict.lookupOrDefault<scalar>
(
"edgeMasterCatchFraction",
edgeMasterCatchFraction_
);
edgeCoPlanarTol_ = dict.lookupOrDefault<scalar>
(
"edgeCoPlanarTol",
edgeCoPlanarTol_
);
edgeEndCutoffTol_ = dict.lookupOrDefault<scalar>
(
"edgeEndCutoffTol",
edgeEndCutoffTol_
);
if (report)
{
Info<< "Sliding interface parameters:" << nl
<< "pointMergeTol : " << pointMergeTol_ << nl
<< "edgeMergeTol : " << edgeMergeTol_ << nl
<< "nFacesPerSlaveEdge : " << nFacesPerSlaveEdge_ << nl
<< "edgeFaceEscapeLimit : " << edgeFaceEscapeLimit_ << nl
<< "integralAdjTol : " << integralAdjTol_ << nl
<< "edgeMasterCatchFraction : " << edgeMasterCatchFraction_ << nl
<< "edgeCoPlanarTol : " << edgeCoPlanarTol_ << nl
<< "edgeEndCutoffTol : " << edgeEndCutoffTol_ << endl;
}
}
// ************************************************************************* //

View File

@ -1,382 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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::slidingInterface
Description
Sliding interface mesh modifier. Given two face zones, couple the
master and slave side using a cutting procedure.
The coupled faces are collected into the "coupled" zone and can become
either internal or placed into a master and slave coupled zone. The
remaining faces (uncovered master or slave) are placed into the master
and slave patch.
The definition of the sliding interface can be either integral or partial.
Integral interface implies that the slave side completely covers
the master (i.e. no faces are uncovered); partial interface
implies that the uncovered part of master/slave face zone should
become boundary faces.
SourceFiles
slidingInterface.C
coupleSlidingInterface.C
decoupleSlidingInterface.C
slidingInterfaceProjectPoints.C
slidingInterfaceAttachedAddressing.C
slidingInterfaceClearCouple.C
\*---------------------------------------------------------------------------*/
#ifndef slidingInterface_H
#define slidingInterface_H
#include "polyMeshModifier.H"
#include "primitiveFacePatch.H"
#include "polyPatchDynamicID.H"
#include "ZoneDynamicIDs.H"
#include "intersection.H"
#include "Pair.H"
#include "objectHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class slidingInterface Declaration
\*---------------------------------------------------------------------------*/
class slidingInterface
:
public polyMeshModifier
{
public:
// Public enumerations
//- Type of match
enum typeOfMatch
{
INTEGRAL,
PARTIAL
};
//- Direction names
static const NamedEnum<typeOfMatch, 2> typeOfMatchNames_;
private:
// Private Data
//- Master face zone ID
faceZoneDynamicID masterFaceZoneIndex_;
//- Slave face zone ID
faceZoneDynamicID slaveFaceZoneIndex_;
//- Cut point zone ID
pointZoneDynamicID cutPointZoneIndex_;
//- Cut face zone ID
faceZoneDynamicID cutFaceZoneIndex_;
//- Master patch ID
polyPatchDynamicID masterPatchIndex_;
//- Slave patch ID
polyPatchDynamicID slavePatchIndex_;
//- Type of match
const typeOfMatch matchType_;
//- Couple-decouple operation.
// If the interface is coupled, decouple it and vice versa.
// Used in conjunction with automatic mesh motion
mutable Switch coupleDecouple_;
//- State of the modifier
mutable Switch attached_;
//- Point projection algorithm
intersection::algorithm projectionAlgo_;
//- Trigger topological change
mutable bool trigger_;
// Tolerances. Initialised to static ones below.
//- Point merge tolerance
scalar pointMergeTol_;
//- Edge merge tolerance
scalar edgeMergeTol_;
//- Estimated number of faces an edge goes through
label nFacesPerSlaveEdge_;
//- Edge-face interaction escape limit
label edgeFaceEscapeLimit_;
//- Integral match point adjustment tolerance
scalar integralAdjTol_;
//- Edge intersection master catch fraction
scalar edgeMasterCatchFraction_;
//- Edge intersection co-planar tolerance
scalar edgeCoPlanarTol_;
//- Edge end cut-off tolerance
scalar edgeEndCutoffTol_;
// Private addressing data.
//- Cut face master face. Gives the index of face in master patch
// the cut face has been created from. For a slave-only face
// this will be -1
mutable labelList* cutFaceMasterPtr_;
//- Cut face slave face. Gives the index of face in slave patch
// the cut face has been created from. For a master-only face
// this will be -1
mutable labelList* cutFaceSlavePtr_;
//- Master zone faceCells
mutable labelList* masterFaceCellsPtr_;
//- Slave zone faceCells
mutable labelList* slaveFaceCellsPtr_;
//- Master stick-out faces
mutable labelList* masterStickOutFacesPtr_;
//- Slave stick-out faces
mutable labelList* slaveStickOutFacesPtr_;
//- Retired point mapping.
// For every retired slave side point, gives the label of the
// master point that replaces it
mutable Map<label>* retiredPointMapPtr_;
//- Cut edge pairs
// For cut points created by intersection two edges,
// store the master-slave edge pair used in creation
mutable Map<Pair<edge>>* cutPointEdgePairMapPtr_;
//- Slave point hit. The index of master point hit by the
// slave point in projection. For no point hit, set to -1
mutable labelList* slavePointPointHitsPtr_;
//- Slave edge hit. The index of master edge hit by the
// slave point in projection. For point or no edge hit, set to -1
mutable labelList* slavePointEdgeHitsPtr_;
//- Slave face hit. The index of master face hit by the
// slave point in projection.
mutable List<objectHit>* slavePointFaceHitsPtr_;
//- Master point edge hit. The index of slave edge hit by
// a master point. For no hit set to -1
mutable labelList* masterPointEdgeHitsPtr_;
//- Projected slave points
mutable pointField* projectedSlavePointsPtr_;
// Private Member Functions
//- Clear out
void clearOut() const;
//- Check validity of construction data
void checkDefinition();
//- Calculate attached addressing
void calcAttachedAddressing() const;
//- Calculate decoupled zone face-cell addressing
void renumberAttachedAddressing(const polyTopoChangeMap&) const;
//- Clear attached addressing
void clearAttachedAddressing() const;
// Topological changes
//- Master faceCells
const labelList& masterFaceCells() const;
//- Slave faceCells
const labelList& slaveFaceCells() const;
//- Master stick-out faces
const labelList& masterStickOutFaces() const;
//- Slave stick-out faces
const labelList& slaveStickOutFaces() const;
//- Retired point map
const Map<label>& retiredPointMap() const;
//- Cut point edge pair map
const Map<Pair<edge>>& cutPointEdgePairMap() const;
//- Clear addressing
void clearAddressing() const;
//- Project slave points and compare with the current projection.
// If the projection has changed, the sliding interface
// changes topologically
bool projectPoints() const;
//- Couple sliding interface
void coupleInterface(polyTopoChange& ref) const;
//- Clear projection
void clearPointProjection() const;
//- Clear old couple
void clearCouple(polyTopoChange& ref) const;
//- Decouple interface (returns it to decoupled state)
// Note: this should not be used in normal operation of the
// sliding mesh, but only to return the mesh to its
// original state
void decoupleInterface(polyTopoChange& ref) const;
// Static Data Members
//- Point merge tolerance
static const scalar pointMergeTolDefault_;
//- Edge merge tolerance
static const scalar edgeMergeTolDefault_;
//- Estimated number of faces an edge goes through
static const label nFacesPerSlaveEdgeDefault_;
//- Edge-face interaction escape limit
static const label edgeFaceEscapeLimitDefault_;
//- Integral match point adjustment tolerance
static const scalar integralAdjTolDefault_;
//- Edge intersection master catch fraction
static const scalar edgeMasterCatchFractionDefault_;
//- Edge intersection co-planar tolerance
static const scalar edgeCoPlanarTolDefault_;
//- Edge end cut-off tolerance
static const scalar edgeEndCutoffTolDefault_;
public:
//- Runtime type information
TypeName("slidingInterface");
// Constructors
//- Construct from components
slidingInterface
(
const word& name,
const polyMesh& mesh,
const word& masterFaceZoneName,
const word& slaveFaceZoneName,
const word& cutPointZoneName,
const word& cutFaceZoneName,
const word& masterPatchName,
const word& slavePatchName,
const typeOfMatch tom,
const bool coupleDecouple = false,
const intersection::algorithm algo =
intersection::algorithm::visible
);
//- Disallow default bitwise copy construction
slidingInterface(const slidingInterface&) = delete;
//- Destructor
virtual ~slidingInterface();
// Member Functions
//- Return master face zone ID
const faceZoneDynamicID& masterFaceZoneIndex() const;
//- Return slave face zone ID
const faceZoneDynamicID& slaveFaceZoneIndex() const;
//- Return true if attached
bool attached() const
{
return attached_;
}
//- Check for topology change
virtual bool changeTopology() const;
//- Insert the layer addition/removal instructions
// into the topological change
virtual void setRefinement(polyTopoChange&) const;
//- Modify motion points to comply with the topological change
virtual void modifyMotionPoints(pointField& motionPoints) const;
//- Force recalculation of locally stored data on topological change
virtual void topoChange(const polyTopoChangeMap&);
//- Return projected points for a slave patch
const pointField& pointProjection() const;
//- Set the tolerances from the values in a dictionary
void setTolerances(const dictionary&, bool report=false);
// Member Operators
//- Disallow default bitwise assignment
void operator=(const slidingInterface&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,555 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "slidingInterface.H"
#include "polyMesh.H"
#include "polyTopoChangeMap.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::slidingInterface::calcAttachedAddressing() const
{
if (debug)
{
Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
<< " for object " << name() << " : "
<< "Calculating zone face-cell addressing."
<< endl;
}
if (!attached_)
{
// Clear existing addressing
clearAttachedAddressing();
const polyMesh& mesh = this->mesh();
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const meshFaceZones& faceZones = mesh.faceZones();
// Master side
const primitiveFacePatch& masterPatch =
faceZones[masterFaceZoneIndex_.index()]();
const labelList& masterPatchFaces =
faceZones[masterFaceZoneIndex_.index()];
const boolList& masterFlip =
faceZones[masterFaceZoneIndex_.index()].flipMap();
masterFaceCellsPtr_ = new labelList(masterPatchFaces.size());
labelList& mfc = *masterFaceCellsPtr_;
forAll(masterPatchFaces, facei)
{
if (masterFlip[facei])
{
mfc[facei] = nei[masterPatchFaces[facei]];
}
else
{
mfc[facei] = own[masterPatchFaces[facei]];
}
}
// Slave side
const primitiveFacePatch& slavePatch =
faceZones[slaveFaceZoneIndex_.index()]();
const labelList& slavePatchFaces =
faceZones[slaveFaceZoneIndex_.index()];
const boolList& slaveFlip =
faceZones[slaveFaceZoneIndex_.index()].flipMap();
slaveFaceCellsPtr_ = new labelList(slavePatchFaces.size());
labelList& sfc = *slaveFaceCellsPtr_;
forAll(slavePatchFaces, facei)
{
if (slaveFlip[facei])
{
sfc[facei] = nei[slavePatchFaces[facei]];
}
else
{
sfc[facei] = own[slavePatchFaces[facei]];
}
}
// Check that the addressing is valid
if (min(mfc) < 0 || min(sfc) < 0)
{
if (debug)
{
forAll(mfc, facei)
{
if (mfc[facei] < 0)
{
Pout<< "No cell next to master patch face " << facei
<< ". Global face no: " << mfc[facei]
<< " own: " << own[masterPatchFaces[facei]]
<< " nei: " << nei[masterPatchFaces[facei]]
<< " flip: " << masterFlip[facei] << endl;
}
}
forAll(sfc, facei)
{
if (sfc[facei] < 0)
{
Pout<< "No cell next to slave patch face " << facei
<< ". Global face no: " << sfc[facei]
<< " own: " << own[slavePatchFaces[facei]]
<< " nei: " << nei[slavePatchFaces[facei]]
<< " flip: " << slaveFlip[facei] << endl;
}
}
}
FatalErrorInFunction
<< "decoupled mesh or sliding interface definition."
<< abort(FatalError);
}
// Calculate stick-out faces
const labelListList& pointFaces = mesh.pointFaces();
// Master side
labelHashSet masterStickOutFaceMap
(
primitiveMesh::facesPerCell_*(masterPatch.size())
);
const labelList& masterMeshPoints = masterPatch.meshPoints();
forAll(masterMeshPoints, pointi)
{
const labelList& curFaces = pointFaces[masterMeshPoints[pointi]];
forAll(curFaces, facei)
{
// Check if the face belongs to the master face zone;
// if not add it
if
(
faceZones.whichZone(curFaces[facei])
!= masterFaceZoneIndex_.index()
)
{
masterStickOutFaceMap.insert(curFaces[facei]);
}
}
}
masterStickOutFacesPtr_ = new labelList(masterStickOutFaceMap.toc());
// Slave side
labelHashSet slaveStickOutFaceMap
(
primitiveMesh::facesPerCell_*(slavePatch.size())
);
const labelList& slaveMeshPoints = slavePatch.meshPoints();
forAll(slaveMeshPoints, pointi)
{
const labelList& curFaces = pointFaces[slaveMeshPoints[pointi]];
forAll(curFaces, facei)
{
// Check if the face belongs to the slave face zone;
// if not add it
if
(
faceZones.whichZone(curFaces[facei])
!= slaveFaceZoneIndex_.index()
)
{
slaveStickOutFaceMap.insert(curFaces[facei]);
}
}
}
slaveStickOutFacesPtr_ = new labelList(slaveStickOutFaceMap.toc());
// Retired point addressing does not exist at this stage.
// It will be filled when the interface is coupled.
retiredPointMapPtr_ =
new Map<label>
(
2*faceZones[slaveFaceZoneIndex_.index()]().nPoints()
);
// Ditto for cut point edge map. This is a rough guess of its size
//
cutPointEdgePairMapPtr_ =
new Map<Pair<edge>>
(
faceZones[slaveFaceZoneIndex_.index()]().nEdges()
);
}
else
{
FatalErrorInFunction
<< "cannot be assembled for object " << name()
<< abort(FatalError);
}
if (debug)
{
Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
<< " for object " << name() << " : "
<< "Finished calculating zone face-cell addressing."
<< endl;
}
}
void Foam::slidingInterface::clearAttachedAddressing() const
{
deleteDemandDrivenData(masterFaceCellsPtr_);
deleteDemandDrivenData(slaveFaceCellsPtr_);
deleteDemandDrivenData(masterStickOutFacesPtr_);
deleteDemandDrivenData(slaveStickOutFacesPtr_);
deleteDemandDrivenData(retiredPointMapPtr_);
deleteDemandDrivenData(cutPointEdgePairMapPtr_);
}
void Foam::slidingInterface::renumberAttachedAddressing
(
const polyTopoChangeMap& m
) const
{
// Get reference to reverse cell renumbering
// The renumbering map is needed the other way around, i.e. giving
// the new cell number for every old cell next to the interface.
const labelList& reverseCellMap = m.reverseCellMap();
const labelList& mfc = masterFaceCells();
const labelList& sfc = slaveFaceCells();
// Master side
labelList* newMfcPtr = new labelList(mfc.size(), -1);
labelList& newMfc = *newMfcPtr;
const labelList& mfzRenumber =
m.faceZoneFaceMap()[masterFaceZoneIndex_.index()];
forAll(mfc, facei)
{
label newCelli = reverseCellMap[mfc[mfzRenumber[facei]]];
if (newCelli >= 0)
{
newMfc[facei] = newCelli;
}
}
// Slave side
labelList* newSfcPtr = new labelList(sfc.size(), -1);
labelList& newSfc = *newSfcPtr;
const labelList& sfzRenumber =
m.faceZoneFaceMap()[slaveFaceZoneIndex_.index()];
forAll(sfc, facei)
{
label newCelli = reverseCellMap[sfc[sfzRenumber[facei]]];
if (newCelli >= 0)
{
newSfc[facei] = newCelli;
}
}
if (debug)
{
// Check if all the mapped cells are live
if (min(newMfc) < 0 || min(newSfc) < 0)
{
FatalErrorInFunction
<< "Error in cell renumbering for object " << name()
<< ". Some of master cells next "
<< "to the interface have been removed."
<< abort(FatalError);
}
}
// Renumber stick-out faces
const labelList& reverseFaceMap = m.reverseFaceMap();
// Master side
const labelList& msof = masterStickOutFaces();
labelList* newMsofPtr = new labelList(msof.size(), -1);
labelList& newMsof = *newMsofPtr;
forAll(msof, facei)
{
label newFacei = reverseFaceMap[msof[facei]];
if (newFacei >= 0)
{
newMsof[facei] = newFacei;
}
}
// Pout<< "newMsof: " << newMsof << endl;
// Slave side
const labelList& ssof = slaveStickOutFaces();
labelList* newSsofPtr = new labelList(ssof.size(), -1);
labelList& newSsof = *newSsofPtr;
forAll(ssof, facei)
{
label newFacei = reverseFaceMap[ssof[facei]];
if (newFacei >= 0)
{
newSsof[facei] = newFacei;
}
}
// Pout<< "newSsof: " << newSsof << endl;
if (debug)
{
// Check if all the mapped cells are live
if (min(newMsof) < 0 || min(newSsof) < 0)
{
FatalErrorInFunction
<< "Error in face renumbering for object " << name()
<< ". Some of stick-out next "
<< "to the interface have been removed."
<< abort(FatalError);
}
}
// Renumber the retired point map. Need to take a copy!
const Map<label> rpm = retiredPointMap();
Map<label>* newRpmPtr = new Map<label>(rpm.size());
Map<label>& newRpm = *newRpmPtr;
const labelList rpmToc = rpm.toc();
// Get reference to point renumbering
const labelList& reversePointMap = m.reversePointMap();
label key, value;
forAll(rpmToc, rpmTocI)
{
key = reversePointMap[rpmToc[rpmTocI]];
value = reversePointMap[rpm.find(rpmToc[rpmTocI])()];
if (debug)
{
// Check if all the mapped cells are live
if (key < 0 || value < 0)
{
FatalErrorInFunction
<< "Error in retired point numbering for object "
<< name() << ". Some of master "
<< "points have been removed."
<< abort(FatalError);
}
}
newRpm.insert(key, value);
}
// Renumber the cut point edge pair map. Need to take a copy!
const Map<Pair<edge>> cpepm = cutPointEdgePairMap();
Map<Pair<edge>>* newCpepmPtr = new Map<Pair<edge>>(cpepm.size());
Map<Pair<edge>>& newCpepm = *newCpepmPtr;
const labelList cpepmToc = cpepm.toc();
forAll(cpepmToc, cpepmTocI)
{
key = reversePointMap[cpepmToc[cpepmTocI]];
const Pair<edge>& oldPe = cpepm.find(cpepmToc[cpepmTocI])();
// Re-do the edges in global addressing
const label ms = reversePointMap[oldPe.first().start()];
const label me = reversePointMap[oldPe.first().end()];
const label ss = reversePointMap[oldPe.second().start()];
const label se = reversePointMap[oldPe.second().end()];
if (debug)
{
// Check if all the mapped cells are live
if (key < 0 || ms < 0 || me < 0 || ss < 0 || se < 0)
{
FatalErrorInFunction
<< "Error in cut point edge pair map numbering for object "
<< name() << ". Some of master points have been removed."
<< abort(FatalError);
}
}
newCpepm.insert(key, Pair<edge>(edge(ms, me), edge(ss, se)));
}
if (!projectedSlavePointsPtr_)
{
FatalErrorInFunction
<< "Error in projected point numbering for object " << name()
<< abort(FatalError);
}
// Renumber the projected slave zone points
const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
pointField* newProjectedSlavePointsPtr
(
new pointField(projectedSlavePoints.size())
);
pointField& newProjectedSlavePoints = *newProjectedSlavePointsPtr;
const labelList& sfzPointRenumber =
m.faceZonePointMap()[slaveFaceZoneIndex_.index()];
forAll(newProjectedSlavePoints, pointi)
{
if (sfzPointRenumber[pointi] > -1)
{
newProjectedSlavePoints[pointi] =
projectedSlavePoints[sfzPointRenumber[pointi]];
}
}
// Re-set the lists
clearAttachedAddressing();
deleteDemandDrivenData(projectedSlavePointsPtr_);
masterFaceCellsPtr_ = newMfcPtr;
slaveFaceCellsPtr_ = newSfcPtr;
masterStickOutFacesPtr_ = newMsofPtr;
slaveStickOutFacesPtr_ = newSsofPtr;
retiredPointMapPtr_ = newRpmPtr;
cutPointEdgePairMapPtr_ = newCpepmPtr;
projectedSlavePointsPtr_ = newProjectedSlavePointsPtr;
}
const Foam::labelList& Foam::slidingInterface::masterFaceCells() const
{
if (!masterFaceCellsPtr_)
{
FatalErrorInFunction
<< "Master zone face-cell addressing not available for object "
<< name()
<< abort(FatalError);
}
return *masterFaceCellsPtr_;
}
const Foam::labelList& Foam::slidingInterface::slaveFaceCells() const
{
if (!slaveFaceCellsPtr_)
{
FatalErrorInFunction
<< "Slave zone face-cell addressing not available for object "
<< name()
<< abort(FatalError);
}
return *slaveFaceCellsPtr_;
}
const Foam::labelList& Foam::slidingInterface::masterStickOutFaces() const
{
if (!masterStickOutFacesPtr_)
{
FatalErrorInFunction
<< "Master zone stick-out face addressing not available for object "
<< name()
<< abort(FatalError);
}
return *masterStickOutFacesPtr_;
}
const Foam::labelList& Foam::slidingInterface::slaveStickOutFaces() const
{
if (!slaveStickOutFacesPtr_)
{
FatalErrorInFunction
<< "Slave zone stick-out face addressing not available for object "
<< name()
<< abort(FatalError);
}
return *slaveStickOutFacesPtr_;
}
const Foam::Map<Foam::label>& Foam::slidingInterface::retiredPointMap() const
{
if (!retiredPointMapPtr_)
{
FatalErrorInFunction
<< "Retired point map not available for object " << name()
<< abort(FatalError);
}
return *retiredPointMapPtr_;
}
const Foam::Map<Foam::Pair<Foam::edge>>&
Foam::slidingInterface::cutPointEdgePairMap() const
{
if (!cutPointEdgePairMapPtr_)
{
FatalErrorInFunction
<< "Retired point map not available for object " << name()
<< abort(FatalError);
}
return *cutPointEdgePairMapPtr_;
}
// ************************************************************************* //

View File

@ -1,76 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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 "slidingInterface.H"
#include "polyTopoChange.H"
#include "polyMesh.H"
#include "polyRemovePoint.H"
#include "polyRemoveFace.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::slidingInterface::clearCouple
(
polyTopoChange& ref
) const
{
if (debug)
{
Pout<< "void slidingInterface::clearCouple("
<< "polyTopoChange& ref) const for object " << name() << " : "
<< "Clearing old couple points and faces." << endl;
}
// Remove all points from the point zone
const polyMesh& mesh = this->mesh();
const labelList& cutPointZoneLabels =
mesh.pointZones()[cutPointZoneIndex_.index()];
forAll(cutPointZoneLabels, pointi)
{
ref.setAction(polyRemovePoint(cutPointZoneLabels[pointi]));
}
// Remove all faces from the face zone
const labelList& cutFaceZoneLabels =
mesh.faceZones()[cutFaceZoneIndex_.index()];
forAll(cutFaceZoneLabels, facei)
{
ref.setAction(polyRemoveFace(cutFaceZoneLabels[facei]));
}
if (debug)
{
Pout<< "void slidingInterface::clearCouple("
<< "polyTopoChange& ref) const for object " << name() << " : "
<< "Finished clearing old couple points and faces." << endl;
}
}
// ************************************************************************* //