Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2013-09-17 15:58:41 +01:00
29 changed files with 2591 additions and 132 deletions

View File

@ -1312,8 +1312,16 @@ int main(int argc, char *argv[])
// Use the resolveFeatureAngle from the refinement parameters
scalar curvature = refineParams.curvature();
scalar planarAngle = refineParams.planarAngle();
snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
snapDriver.doSnap
(
snapDict,
motionDict,
curvature,
planarAngle,
snapParams
);
writeMesh
(

View File

@ -194,10 +194,13 @@ castellatedMeshControls
// Planar angle:
// - used to determine if surface normals
// are roughly the same or opposite. Used e.g. in proximity refinement
// and to decide when to merge free-standing baffles
// If e.g. running in surfaceSimplify mode set this to 180 to
// merge all baffles
// are roughly the same or opposite. Used
// - in proximity refinement
// - to decide when to merge free-standing baffles
// (if e.g. running in surfaceSimplify mode set this to 180 to
// merge all baffles)
// - in snapping to avoid snapping to nearest on 'wrong' side
// of thin gap
//
// If not specified same as resolveFeatureAngle
planarAngle 30;
@ -284,6 +287,10 @@ snapControls
// Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
multiRegionFeatureSnap false;
// wip: disable snapping to opposite near surfaces (revert to 22x behaviour)
// detectNearSurfacesSnap false;
}
// Settings for the layer addition.

View File

@ -0,0 +1,4 @@
patchFaceOrientation.C
orientFaceZone.C
EXE = $(FOAM_APPBIN)/orientFaceZone

View File

@ -0,0 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools \
-ltriSurface

View File

@ -0,0 +1,22 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Get application name
application=`getApplication`
runApplication blockMesh
# Make random cell numbering to get faceZone orientation random
runApplication renumberMesh -dict system/renumberMeshDict
# Construct faceZone
runApplication topoSet
runApplication decomposePar -force -cellDist
runParallel orientFaceZone 2 f0 '(10 0 0)'
runApplication reconstructParMesh -latestTime -mergeTol 1e-6
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1 @@
2013-09-16 To test orientFaceZone. Make random orientation. Orient afterwards.

View File

@ -0,0 +1,75 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.1;
vertices
(
(0 0 0)
(1 0 0)
(1 1 0)
(0 1 0)
(0 0 0.1)
(1 0 0.1)
(1 1 0.1)
(0 1 0.1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
movingWall
{
type wall;
faces
(
(3 7 6 2)
);
}
fixedWalls
{
type wall;
faces
(
(0 4 7 3)
(2 6 5 1)
(1 5 4 0)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(4 5 6 7)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

View File

@ -0,0 +1,41 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3
(
movingWall
{
type wall;
nFaces 20;
startFace 760;
}
fixedWalls
{
type wall;
nFaces 60;
startFace 780;
}
frontAndBack
{
type empty;
inGroups 1(empty);
nFaces 800;
startFace 840;
}
)
// ************************************************************************* //

View File

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu nu [ 0 2 -1 0 0 0 0 ] 0.01;
// ************************************************************************* //

View File

@ -0,0 +1,49 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 0.5;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //

View File

@ -0,0 +1,135 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
note "mesh decomposition control dictionary";
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 2;
//- Keep owner and neighbour on same processor for faces in zones:
// preserveFaceZones (heater solid1 solid3);
//- Keep owner and neighbour on same processor for faces in patches:
// (makes sense only for cyclic patches)
//preservePatches (cyclic_half0 cyclic_half1);
//- Keep all of faceSet on a single processor. This puts all cells
// connected with a point, edge or face on the same processor.
// (just having face connected cells might not guarantee a balanced
// decomposition)
// The processor can be -1 (the decompositionMethod chooses the processor
// for a good load balance) or explicitly provided (upsets balance).
//singleProcessorFaceSets ((f0 -1));
//- Use the volScalarField named here as a weight for each cell in the
// decomposition. For example, use a particle population field to decompose
// for a balanced number of particles in a lagrangian simulation.
// weightField dsmcRhoNMean;
//method scotch;
method hierarchical;
// method simple;
// method metis;
// method manual;
// method multiLevel;
// method structured; // does 2D decomposition of structured mesh
multiLevelCoeffs
{
// Decomposition methods to apply in turn. This is like hierarchical but
// fully general - every method can be used at every level.
level0
{
numberOfSubdomains 64;
//method simple;
//simpleCoeffs
//{
// n (2 1 1);
// delta 0.001;
//}
method scotch;
}
level1
{
numberOfSubdomains 4;
method scotch;
}
}
// Desired output
simpleCoeffs
{
n (2 1 1);
delta 0.001;
}
hierarchicalCoeffs
{
n (2 1 1);
delta 0.001;
order xyz;
}
metisCoeffs
{
/*
processorWeights
(
1
1
1
1
);
*/
}
scotchCoeffs
{
//processorWeights
//(
// 1
// 1
// 1
// 1
//);
//writeGraph true;
//strategy "b";
}
manualCoeffs
{
dataFile "decompositionData";
}
structuredCoeffs
{
// Patches to do 2D decomposition on. Structured mesh only; cells have
// to be in 'columns' on top of patches.
patches (bottomPatch);
}
//// Is the case distributed? Note: command-line argument -roots takes
//// precedence
//distributed yes;
//// Per slave (so nProcs-1 entries) the directory above the case.
//roots
//(
// "/tmp"
// "/tmp"
//);
// ************************************************************************* //

View File

@ -0,0 +1,60 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear orthogonal;
laplacian((1|A(U)),p) Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradSchemes
{
default orthogonal;
}
fluxRequired
{
default no;
p ;
}
// ************************************************************************* //

View File

@ -0,0 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
}
U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,110 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
note "mesh renumbering dictionary";
object renumberMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Write maps from renumbered back to original mesh
writeMaps true;
// Optional entry: sort cells on coupled boundaries to last for use with
// e.g. nonBlockingGaussSeidel.
sortCoupledFaceCells false;
// Optional entry: renumber on a block-by-block basis. It uses a
// blockCoeffs dictionary to construct a decompositionMethod to do
// a block subdivision) and then applies the renumberMethod to each
// block in turn. This can be used in large cases to keep the blocks
// fitting in cache with all the the cache misses bunched at the end.
// This number is the approximate size of the blocks - this gets converted
// to a number of blocks that is the input to the decomposition method.
//blockSize 1000;
// Optional entry: sort points into internal and boundary points
//orderPoints false;
//method CuthillMcKee;
//method Sloan;
//method manual;
method random;
//method structured;
//method spring;
//method zoltan; // only if compiled with zoltan support
//CuthillMcKeeCoeffs
//{
// // Reverse CuthillMcKee (RCM) or plain
// reverse true;
//}
manualCoeffs
{
// In system directory: new-to-original (i.e. order) labelIOList
dataFile "cellMap";
}
// For extruded (i.e. structured in one direction) meshes
structuredCoeffs
{
// Patches that mesh was extruded from. These determine the starting
// layer of cells
patches (movingWall);
// Method to renumber the starting layer of cells
method random;
// Renumber in columns (depthFirst) or in layers
depthFirst true;
// Optional: reverse ordering
//reverse false;
}
springCoeffs
{
// Maximum jump of cell indices. Is fraction of number of cells
maxCo 0.01;
// Limit the amount of movement; the fraction maxCo gets decreased
// with every iteration
freezeFraction 0.999;
// Maximum number of iterations
maxIter 1000;
}
blockCoeffs
{
method scotch;
//method hierarchical;
//hierarchicalCoeffs
//{
// n (1 2 1);
// delta 0.001;
// order xyz;
//}
}
zoltanCoeffs
{
ORDER_METHOD LOCAL_HSFC;
}
// ************************************************************************* //

View File

@ -0,0 +1,64 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object topoSetDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
actions
(
// Get 'T' of faces
{
name f0;
type faceSet;
action new;
source boxToFace;
sourceInfo
{
boxes
(
(0.0499 -100 -100)(0.0501 100 100)
(0.03 0.0499 -100)(0.05 0.0501 100)
);
}
}
// // Pick up some cells
// {
// name c0;
// type cellSet;
// action new;
// source boxToCell;
// sourceInfo
// {
// boxes
// (
// (-100 -100 -100)(0.05 0.05 100)
// (0.05 0.05 -100)(100 100 100)
// );
// }
// }
// Convert faceSet to faceZone
{
name f0;
type faceZoneSet;
action new;
source setToFaceZone;
sourceInfo
{
faceSet f0;
}
}
);
// ************************************************************************* //

View File

@ -0,0 +1,390 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
Application
orientFaceZone
Description
Corrects orientation of faceZone.
- correct in parallel - excludes coupled faceZones from walk
- correct for non-manifold faceZones - restarts walk
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "syncTools.H"
#include "patchFaceOrientation.H"
#include "PatchEdgeFaceWave.H"
#include "orientedSurface.H"
#include "globalIndex.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::validArgs.append("faceZone");
argList::validArgs.append("outsidePoint");
# include "setRootCase.H"
# include "createTime.H"
# include "createPolyMesh.H"
const word zoneName = args[1];
const point outsidePoint = args.argRead<point>(2);
Info<< "Orienting faceZone " << zoneName
<< " such that " << outsidePoint << " is outside"
<< nl << endl;
const faceZone& fZone = mesh.faceZones()[zoneName];
if (fZone.checkParallelSync())
{
FatalErrorIn(args.executable())
<< "Face zone " << fZone.name()
<< " is not parallel synchronised."
<< " Any coupled face also needs its coupled version to be included"
<< " and with opposite flipMap."
<< exit(FatalError);
}
const labelList& faceLabels = fZone;
const indirectPrimitivePatch patch
(
IndirectList<face>(mesh.faces(), faceLabels),
mesh.points()
);
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
// Data on all edges and faces
List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
List<patchFaceOrientation> allFaceInfo(patch.size());
// Make sure we don't walk through
// - slaves of coupled faces
// - non-manifold edges
{
const polyBoundaryMesh& bm = mesh.boundaryMesh();
label nProtected = 0;
forAll(faceLabels, faceI)
{
const label meshFaceI = faceLabels[faceI];
const label patchI = bm.whichPatch(meshFaceI);
if
(
patchI != -1
&& bm[patchI].coupled()
&& !isMasterFace[meshFaceI]
)
{
// Slave side. Mark so doesn't get visited.
allFaceInfo[faceI] = orientedSurface::NOFLIP;
nProtected++;
}
}
Info<< "Protected from visiting "
<< returnReduce(nProtected, sumOp<label>())
<< " slaves of coupled faces" << nl << endl;
}
{
// Number of (master)faces per edge
labelList nMasterFaces(patch.nEdges(), 0);
forAll(faceLabels, faceI)
{
const label meshFaceI = faceLabels[faceI];
if (isMasterFace[meshFaceI])
{
const labelList& fEdges = patch.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
nMasterFaces[fEdges[fEdgeI]]++;
}
}
}
syncTools::syncEdgeList
(
mesh,
patch.meshEdges(mesh.edges(), mesh.pointEdges()),
nMasterFaces,
plusEqOp<label>(),
0
);
label nProtected = 0;
forAll(nMasterFaces, edgeI)
{
if (nMasterFaces[edgeI] > 2)
{
allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
nProtected++;
}
}
Info<< "Protected from visiting "
<< returnReduce(nProtected, sumOp<label>())
<< " non-manifold edges" << nl << endl;
}
DynamicList<label> changedEdges;
DynamicList<patchFaceOrientation> changedInfo;
const scalar tol = PatchEdgeFaceWave
<
indirectPrimitivePatch,
patchFaceOrientation
>::propagationTol();
int dummyTrackData;
globalIndex globalFaces(patch.size());
while (true)
{
// Pick an unset face
label unsetFaceI = labelMax;
forAll(allFaceInfo, faceI)
{
if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
{
unsetFaceI = globalFaces.toGlobal(faceI);
break;
}
}
reduce(unsetFaceI, minOp<label>());
if (unsetFaceI == labelMax)
{
break;
}
label procI = globalFaces.whichProcID(unsetFaceI);
label seedFaceI = globalFaces.toLocal(procI, unsetFaceI);
Info<< "Seeding from processor " << procI << " face " << seedFaceI
<< endl;
if (procI == Pstream::myProcNo())
{
// Determine orientation of seedFace
vector d = outsidePoint-patch.faceCentres()[seedFaceI];
const vector& fn = patch.faceNormals()[seedFaceI];
// Set information to correct orientation
patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];
faceInfo = orientedSurface::NOFLIP;
if ((fn&d) < 0)
{
faceInfo.flip();
Pout<< "Face " << seedFaceI << " at "
<< patch.faceCentres()[seedFaceI]
<< " with normal " << fn
<< " needs to be flipped." << endl;
}
else
{
Pout<< "Face " << seedFaceI << " at "
<< patch.faceCentres()[seedFaceI]
<< " with normal " << fn
<< " points in positive direction (cos = " << (fn&d)/mag(d)
<< ")" << endl;
}
const labelList& fEdges = patch.faceEdges()[seedFaceI];
forAll(fEdges, fEdgeI)
{
label edgeI = fEdges[fEdgeI];
patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
if
(
edgeInfo.updateEdge<int>
(
mesh,
patch,
edgeI,
seedFaceI,
faceInfo,
tol,
dummyTrackData
)
)
{
changedEdges.append(edgeI);
changedInfo.append(edgeInfo);
}
}
}
if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
{
break;
}
// Walk
PatchEdgeFaceWave
<
indirectPrimitivePatch,
patchFaceOrientation
> calc
(
mesh,
patch,
changedEdges,
changedInfo,
allEdgeInfo,
allFaceInfo,
returnReduce(patch.nEdges(), sumOp<label>())
);
}
// Push master zone info over to slave (since slave faces never visited)
{
const polyBoundaryMesh& bm = mesh.boundaryMesh();
labelList neiStatus
(
mesh.nFaces()-mesh.nInternalFaces(),
orientedSurface::UNVISITED
);
forAll(faceLabels, i)
{
const label meshFaceI = faceLabels[i];
if (!mesh.isInternalFace(meshFaceI))
{
neiStatus[meshFaceI-mesh.nInternalFaces()] =
allFaceInfo[i].flipStatus();
}
}
syncTools::swapBoundaryFaceList(mesh, neiStatus);
forAll(faceLabels, i)
{
const label meshFaceI = faceLabels[i];
const label patchI = bm.whichPatch(meshFaceI);
if
(
patchI != -1
&& bm[patchI].coupled()
&& !isMasterFace[meshFaceI]
)
{
// Slave side. Take flipped from neighbour
label bFaceI = meshFaceI-mesh.nInternalFaces();
if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
{
allFaceInfo[i] = orientedSurface::FLIP;
}
else if (neiStatus[bFaceI] == orientedSurface::FLIP)
{
allFaceInfo[i] = orientedSurface::NOFLIP;
}
else
{
FatalErrorIn(args.executable())
<< "Incorrect status for face " << meshFaceI
<< abort(FatalError);
}
}
}
}
// Convert to flipmap and adapt faceZones
boolList newFlipMap(allFaceInfo.size(), false);
label nChanged = 0;
forAll(allFaceInfo, faceI)
{
if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
{
newFlipMap[faceI] = false;
}
else if (allFaceInfo[faceI] == orientedSurface::FLIP)
{
newFlipMap[faceI] = true;
}
else
{
FatalErrorIn(args.executable())
<< "Problem : unvisited face " << faceI
<< " centre:" << mesh.faceCentres()[faceLabels[faceI]]
<< abort(FatalError);
}
if (fZone.flipMap()[faceI] != newFlipMap[faceI])
{
nChanged++;
}
}
reduce(nChanged, sumOp<label>());
if (nChanged > 0)
{
Info<< "Flipping " << nChanged << " out of "
<< globalFaces.size() << " faces." << nl << endl;
mesh.faceZones()[zoneName].resetAddressing(faceLabels, newFlipMap);
if (!mesh.faceZones().write())
{
FatalErrorIn(args.executable())
<< "Failed writing faceZones" << exit(FatalError);
}
}
Info<< "End." << endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "patchFaceOrientation.H"
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Foam::Ostream& os,
const Foam::patchFaceOrientation& wDist
)
{
return os << wDist.flipStatus_;
}
Foam::Istream& Foam::operator>>
(
Foam::Istream& is,
Foam::patchFaceOrientation& wDist
)
{
return is >> wDist.flipStatus_;
}
// ************************************************************************* //

View File

@ -0,0 +1,176 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::patchFaceOrientation
Description
Transport of orientation for use in PatchEdgeFaceWave.
SourceFiles
patchFaceOrientationI.H
patchFaceOrientation.C
\*---------------------------------------------------------------------------*/
#ifndef patchFaceOrientation_H
#define patchFaceOrientation_H
#include "tensor.H"
#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class polyMesh;
/*---------------------------------------------------------------------------*\
Class patchFaceOrientation Declaration
\*---------------------------------------------------------------------------*/
class patchFaceOrientation
{
// Private data
//- flip status
label flipStatus_;
public:
// Constructors
//- Construct null
inline patchFaceOrientation();
//- Construct from components
inline patchFaceOrientation(const label);
// Member Functions
//- Orientation
inline label flipStatus() const;
//- Reverse orientation
inline void flip();
// Needed by meshWave
//- Check whether origin has been changed at all or
// still contains original (invalid) value.
template<class TrackingData>
inline bool valid(TrackingData& td) const;
//- Apply rotation matrix
template<class TrackingData>
inline void transform
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const tensor& rotTensor,
const scalar tol,
TrackingData& td
);
//- Influence of face on edge
template<class TrackingData>
inline bool updateEdge
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const label edgeI,
const label faceI,
const patchFaceOrientation& faceInfo,
const scalar tol,
TrackingData& td
);
//- New information for edge (from e.g. coupled edge)
template<class TrackingData>
inline bool updateEdge
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const patchFaceOrientation& edgeInfo,
const bool sameOrientation,
const scalar tol,
TrackingData& td
);
//- Influence of edge on face.
template<class TrackingData>
inline bool updateFace
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const label faceI,
const label edgeI,
const patchFaceOrientation& edgeInfo,
const scalar tol,
TrackingData& td
);
//- Same (like operator==)
template<class TrackingData>
inline bool equal(const patchFaceOrientation&, TrackingData&) const;
// Member Operators
// Needed for List IO
inline bool operator==(const patchFaceOrientation&) const;
inline bool operator!=(const patchFaceOrientation&) const;
// IOstream Operators
friend Ostream& operator<<(Ostream&, const patchFaceOrientation&);
friend Istream& operator>>(Istream&, patchFaceOrientation&);
};
//- Data associated with patchFaceOrientation type are contiguous
template<>
inline bool contiguous<patchFaceOrientation>()
{
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "patchFaceOrientationI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,227 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "polyMesh.H"
#include "transform.H"
#include "orientedSurface.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline Foam::patchFaceOrientation::patchFaceOrientation()
:
flipStatus_(orientedSurface::UNVISITED)
{}
inline Foam::patchFaceOrientation::patchFaceOrientation
(
const label flipStatus
)
:
flipStatus_(flipStatus)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::patchFaceOrientation::flipStatus() const
{
return flipStatus_;
}
inline void Foam::patchFaceOrientation::flip()
{
if (flipStatus_ == orientedSurface::NOFLIP)
{
flipStatus_ = orientedSurface::FLIP;
}
else if (flipStatus_ == orientedSurface::FLIP)
{
flipStatus_ = orientedSurface::NOFLIP;
}
}
template<class TrackingData>
inline bool Foam::patchFaceOrientation::valid(TrackingData& td) const
{
return flipStatus_ != orientedSurface::UNVISITED;
}
template<class TrackingData>
inline void Foam::patchFaceOrientation::transform
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const tensor& rotTensor,
const scalar tol,
TrackingData& td
)
{}
template<class TrackingData>
inline bool Foam::patchFaceOrientation::updateEdge
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const label edgeI,
const label faceI,
const patchFaceOrientation& faceInfo,
const scalar tol,
TrackingData& td
)
{
if (valid(td))
{
return false;
}
const face& f = patch.localFaces()[faceI];
const edge& e = patch.edges()[edgeI];
//Pout<< "Updating edge:" << edgeI << " verts:" << e << nl
// << " start:" << patch.localPoints()[e[0]] << nl
// << " end:" << patch.localPoints()[e[1]] << endl;
patchFaceOrientation consistentInfo(faceInfo);
// Check how edge relates to face
if (f.edgeDirection(e) < 0)
{
// Create flipped version of faceInfo
consistentInfo.flip();
}
operator=(consistentInfo);
return true;
}
template<class TrackingData>
inline bool Foam::patchFaceOrientation::updateEdge
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const patchFaceOrientation& edgeInfo,
const bool sameOrientation,
const scalar tol,
TrackingData& td
)
{
if (valid(td))
{
return false;
}
// Create (flipped/unflipped) version of edgeInfo
patchFaceOrientation consistentInfo(edgeInfo);
if (!sameOrientation)
{
consistentInfo.flip();
}
operator=(consistentInfo);
return true;
}
template<class TrackingData>
inline bool Foam::patchFaceOrientation::updateFace
(
const polyMesh& mesh,
const indirectPrimitivePatch& patch,
const label faceI,
const label edgeI,
const patchFaceOrientation& edgeInfo,
const scalar tol,
TrackingData& td
)
{
if (valid(td))
{
return false;
}
// Transfer flip to face
const face& f = patch.localFaces()[faceI];
const edge& e = patch.edges()[edgeI];
//Pout<< "Updating face:" << faceI << nl
// << " verts:" << f << nl
// << " with edge:" << edgeI << nl
// << " start:" << patch.localPoints()[e[0]] << nl
// << " end:" << patch.localPoints()[e[1]] << endl;
// Create (flipped/unflipped) version of edgeInfo
patchFaceOrientation consistentInfo(edgeInfo);
if (f.edgeDirection(e) > 0)
{
consistentInfo.flip();
}
operator=(consistentInfo);
return true;
}
template<class TrackingData>
inline bool Foam::patchFaceOrientation::equal
(
const patchFaceOrientation& rhs,
TrackingData& td
) const
{
return operator==(rhs);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline bool Foam::patchFaceOrientation::operator==
(
const Foam::patchFaceOrientation& rhs
) const
{
return flipStatus() == rhs.flipStatus();
}
inline bool Foam::patchFaceOrientation::operator!=
(
const Foam::patchFaceOrientation& rhs
) const
{
return !(*this == rhs);
}
// ************************************************************************* //

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "TimePaths.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -41,12 +42,21 @@ Foam::TimePaths::TimePaths
system_(systemName),
constant_(constantName)
{
// Find out from case name whether a processor directory
std::string::size_type pos = caseName.find("processor");
if (pos != string::npos)
{
processorCase_ = true;
if (pos == 0)
{
globalCaseName_ = ".";
}
else
{
globalCaseName_ = caseName(pos-1);
}
}
else
{
globalCaseName_ = caseName;
@ -70,7 +80,29 @@ Foam::TimePaths::TimePaths
case_(caseName),
system_(systemName),
constant_(constantName)
{}
{
if (!processorCase)
{
// For convenience: find out from case name whether it is a
// processor directory and set processorCase flag so file searching
// goes up one level.
std::string::size_type pos = caseName.find("processor");
if (pos != string::npos)
{
processorCase_ = true;
if (pos == 0)
{
globalCaseName_ = ".";
}
else
{
globalCaseName_ = caseName(pos-1);
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

View File

@ -792,11 +792,636 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
}
Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avgCellCentres
(
const fvMesh& mesh,
const indirectPrimitivePatch& pp
)
{
const labelListList& pointFaces = pp.pointFaces();
tmp<pointField> tavgBoundary
(
new pointField(pointFaces.size(), vector::zero)
);
pointField& avgBoundary = tavgBoundary();
labelList nBoundary(pointFaces.size(), 0);
forAll(pointFaces, pointI)
{
const labelList& pFaces = pointFaces[pointI];
forAll(pFaces, pfI)
{
label faceI = pFaces[pfI];
label meshFaceI = pp.addressing()[faceI];
label own = mesh.faceOwner()[meshFaceI];
avgBoundary[pointI] += mesh.cellCentres()[own];
nBoundary[pointI]++;
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
avgBoundary,
plusEqOp<point>(), // combine op
vector::zero // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
nBoundary,
plusEqOp<label>(), // combine op
label(0) // null value
);
forAll(avgBoundary, i)
{
avgBoundary[i] /= nBoundary[i];
}
return tavgBoundary;
}
//Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::calcEdgeLen
//(
// const indirectPrimitivePatch& pp
//) const
//{
// // Get local edge length based on refinement level
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// // (Ripped from autoLayerDriver)
//
// tmp<scalarField> tedgeLen(new scalarField(pp.nPoints()));
// scalarField& edgeLen = tedgeLen();
// {
// const fvMesh& mesh = meshRefiner_.mesh();
// const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
// const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
//
// labelList maxPointLevel(pp.nPoints(), labelMin);
//
// forAll(pp, i)
// {
// label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
// const face& f = pp.localFaces()[i];
// forAll(f, fp)
// {
// maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
// }
// }
//
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// maxPointLevel,
// maxEqOp<label>(),
// labelMin // null value
// );
//
//
// forAll(maxPointLevel, pointI)
// {
// // Find undistorted edge size for this level.
// edgeLen[pointI] = edge0Len/(1<<maxPointLevel[pointI]);
// }
// }
// return tedgeLen;
//}
void Foam::autoSnapDriver::detectNearSurfaces
(
const scalar planarCos,
const indirectPrimitivePatch& pp,
const pointField& nearestPoint,
const vectorField& nearestNormal,
vectorField& disp
) const
{
Info<< "Detecting near surfaces ..." << endl;
const pointField& localPoints = pp.localPoints();
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
const fvMesh& mesh = meshRefiner_.mesh();
//// Get local edge length based on refinement level
//const scalarField edgeLen(calcEdgeLen(pp));
//
//// Generate rays for every surface point
//// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//{
// const scalar cos45 = Foam::cos(degToRad(45));
// vector n(cos45, cos45, cos45);
// n /= mag(n);
//
// pointField start(14*pp.nPoints());
// pointField end(start.size());
//
// label rayI = 0;
// forAll(localPoints, pointI)
// {
// const point& pt = localPoints[pointI];
//
// // Along coordinate axes
//
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() -= edgeLen[pointI];
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() += edgeLen[pointI];
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.y() -= edgeLen[pointI];
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.y() += edgeLen[pointI];
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.z() -= edgeLen[pointI];
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.z() += edgeLen[pointI];
// }
//
// // At 45 degrees
//
// const vector vec(edgeLen[pointI]*n);
//
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() += vec.x();
// endPt.y() += vec.y();
// endPt.z() += vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() -= vec.x();
// endPt.y() += vec.y();
// endPt.z() += vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() += vec.x();
// endPt.y() -= vec.y();
// endPt.z() += vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() -= vec.x();
// endPt.y() -= vec.y();
// endPt.z() += vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() += vec.x();
// endPt.y() += vec.y();
// endPt.z() -= vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() -= vec.x();
// endPt.y() += vec.y();
// endPt.z() -= vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() += vec.x();
// endPt.y() -= vec.y();
// endPt.z() -= vec.z();
// }
// {
// start[rayI] = pt;
// point& endPt = end[rayI++];
// endPt = pt;
// endPt.x() -= vec.x();
// endPt.y() -= vec.y();
// endPt.z() -= vec.z();
// }
// }
//
// labelList surface1;
// List<pointIndexHit> hit1;
// labelList region1;
// vectorField normal1;
//
// labelList surface2;
// List<pointIndexHit> hit2;
// labelList region2;
// vectorField normal2;
// surfaces.findNearestIntersection
// (
// unzonedSurfaces, // surfacesToTest,
// start,
// end,
//
// surface1,
// hit1,
// region1,
// normal1,
//
// surface2,
// hit2,
// region2,
// normal2
// );
//
// // All intersections
// {
// OBJstream str
// (
// mesh.time().path()
// / "surfaceHits_" + meshRefiner_.timeName() + ".obj"
// );
//
// Info<< "Dumping intersections with rays to " << str.name()
// << endl;
//
// forAll(hit1, i)
// {
// if (hit1[i].hit())
// {
// str.write(linePointRef(start[i], hit1[i].hitPoint()));
// }
// if (hit2[i].hit())
// {
// str.write(linePointRef(start[i], hit2[i].hitPoint()));
// }
// }
// }
//
// // Co-planar intersections
// {
// OBJstream str
// (
// mesh.time().path()
// / "coplanarHits_" + meshRefiner_.timeName() + ".obj"
// );
//
// Info<< "Dumping intersections with co-planar surfaces to "
// << str.name() << endl;
//
// forAll(localPoints, pointI)
// {
// bool hasNormal = false;
// point surfPointA;
// vector surfNormalA;
// point surfPointB;
// vector surfNormalB;
//
// bool isCoplanar = false;
//
// label rayI = 14*pointI;
// for (label i = 0; i < 14; i++)
// {
// if (hit1[rayI].hit())
// {
// const point& pt = hit1[rayI].hitPoint();
// const vector& n = normal1[rayI];
//
// if (!hasNormal)
// {
// hasNormal = true;
// surfPointA = pt;
// surfNormalA = n;
// }
// else
// {
// if
// (
// meshRefiner_.isGap
// (
// planarCos,
// surfPointA,
// surfNormalA,
// pt,
// n
// )
// )
// {
// isCoplanar = true;
// surfPointB = pt;
// surfNormalB = n;
// break;
// }
// }
// }
// if (hit2[rayI].hit())
// {
// const point& pt = hit2[rayI].hitPoint();
// const vector& n = normal2[rayI];
//
// if (!hasNormal)
// {
// hasNormal = true;
// surfPointA = pt;
// surfNormalA = n;
// }
// else
// {
// if
// (
// meshRefiner_.isGap
// (
// planarCos,
// surfPointA,
// surfNormalA,
// pt,
// n
// )
// )
// {
// isCoplanar = true;
// surfPointB = pt;
// surfNormalB = n;
// break;
// }
// }
// }
//
// rayI++;
// }
//
// if (isCoplanar)
// {
// str.write(linePointRef(surfPointA, surfPointB));
// }
// }
// }
//}
const pointField avgCc(avgCellCentres(mesh, pp));
// Construct rays from localPoints to beyond cell centre
pointField end(pp.nPoints());
forAll(localPoints, pointI)
{
const point& pt = localPoints[pointI];
end[pointI] = pt + 2*(avgCc[pointI]-pt);
}
label nOverride = 0;
// 1. All points to non-interface surfaces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
const labelList unzonedSurfaces =
surfaceZonesInfo::getUnnamedSurfaces
(
meshRefiner_.surfaces().surfZones()
);
// Do intersection test
labelList surface1;
List<pointIndexHit> hit1;
labelList region1;
vectorField normal1;
labelList surface2;
List<pointIndexHit> hit2;
labelList region2;
vectorField normal2;
surfaces.findNearestIntersection
(
unzonedSurfaces,
localPoints,
end,
surface1,
hit1,
region1,
normal1,
surface2,
hit2,
region2,
normal2
);
forAll(localPoints, pointI)
{
// Current location
const point& pt = localPoints[pointI];
bool override = false;
if (hit1[pointI].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
nearestPoint[pointI],
nearestNormal[pointI],
hit1[pointI].hitPoint(),
normal1[pointI]
)
)
{
disp[pointI] = hit1[pointI].hitPoint()-pt;
override = true;
}
}
if (hit2[pointI].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
nearestPoint[pointI],
nearestNormal[pointI],
hit2[pointI].hitPoint(),
normal2[pointI]
)
)
{
disp[pointI] = hit2[pointI].hitPoint()-pt;
override = true;
}
}
if (override)
{
nOverride++;
}
}
}
// 2. All points on zones to their respective surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
// Surfaces with zone information
const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
const labelList zonedSurfaces = surfaceZonesInfo::getNamedSurfaces
(
surfZones
);
forAll(zonedSurfaces, i)
{
label zoneSurfI = zonedSurfaces[i];
const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
const labelList surfacesToTest(1, zoneSurfI);
// Get indices of points both on faceZone and on pp.
labelList zonePointIndices
(
getZoneSurfacePoints
(
mesh,
pp,
faceZoneName
)
);
// Do intersection test
labelList surface1;
List<pointIndexHit> hit1;
labelList region1;
vectorField normal1;
labelList surface2;
List<pointIndexHit> hit2;
labelList region2;
vectorField normal2;
surfaces.findNearestIntersection
(
surfacesToTest,
pointField(localPoints, zonePointIndices),
pointField(end, zonePointIndices),
surface1,
hit1,
region1,
normal1,
surface2,
hit2,
region2,
normal2
);
label nOverride = 0;
forAll(hit1, i)
{
label pointI = zonePointIndices[i];
// Current location
const point& pt = localPoints[pointI];
bool override = false;
if (hit1[i].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
nearestPoint[pointI],
nearestNormal[pointI],
hit1[i].hitPoint(),
normal1[i]
)
)
{
disp[pointI] = hit1[i].hitPoint()-pt;
override = true;
}
}
if (hit2[i].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
nearestPoint[pointI],
nearestNormal[pointI],
hit2[i].hitPoint(),
normal2[i]
)
)
{
disp[pointI] = hit2[i].hitPoint()-pt;
override = true;
}
}
if (override)
{
nOverride++;
}
}
}
}
Info<< "Overriding nearest with intersection of close gaps at "
<< returnReduce(nOverride, sumOp<label>())
<< " out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " points." << endl;
}
Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
(
const meshRefinement& meshRefiner,
const scalarField& snapDist,
const indirectPrimitivePatch& pp
const indirectPrimitivePatch& pp,
pointField& nearestPoint,
vectorField& nearestNormal
)
{
Info<< "Calculating patchDisplacement as distance to nearest surface"
@ -815,12 +1440,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
labelList snapSurf(localPoints.size(), -1);
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces =
const labelList zonedSurfaces =
surfaceZonesInfo::getNamedSurfaces
(
meshRefiner.surfaces().surfZones()
);
labelList unzonedSurfaces =
const labelList unzonedSurfaces =
surfaceZonesInfo::getUnnamedSurfaces
(
meshRefiner.surfaces().surfZones()
@ -833,6 +1458,33 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
{
List<pointIndexHit> hitInfo;
labelList hitSurface;
if (nearestNormal.size() == localPoints.size())
{
labelList hitRegion;
vectorField hitNormal;
surfaces.findNearestRegion
(
unzonedSurfaces,
localPoints,
sqr(snapDist),
hitSurface,
hitInfo,
hitRegion,
hitNormal
);
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
nearestPoint[pointI] = hitInfo[pointI].hitPoint();
nearestNormal[pointI] = hitNormal[pointI];
}
}
}
else
{
surfaces.findNearest
(
unzonedSurfaces,
@ -841,6 +1493,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
hitSurface,
hitInfo
);
}
forAll(hitInfo, pointI)
{
@ -888,14 +1541,43 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
// Find nearest for points both on faceZone and pp.
List<pointIndexHit> hitInfo;
labelList hitSurface;
if (nearestNormal.size() == localPoints.size())
{
labelList hitRegion;
vectorField hitNormal;
surfaces.findNearestRegion
(
surfacesToTest,
pointField(localPoints, zonePointIndices),
sqr(scalarField(minSnapDist, zonePointIndices)),
hitSurface,
hitInfo,
hitRegion,
hitNormal
);
forAll(hitInfo, i)
{
if (hitInfo[i].hit())
{
label pointI = zonePointIndices[i];
nearestPoint[pointI] = hitInfo[i].hitPoint();
nearestNormal[pointI] = hitNormal[i];
}
}
}
else
{
surfaces.findNearest
(
labelList(1, zoneSurfI),
surfacesToTest,
pointField(localPoints, zonePointIndices),
sqr(scalarField(minSnapDist, zonePointIndices)),
hitSurface,
hitInfo
);
}
forAll(hitInfo, i)
{
@ -1529,6 +2211,7 @@ void Foam::autoSnapDriver::doSnap
const dictionary& snapDict,
const dictionary& motionDict,
const scalar featureCos,
const scalar planarAngle,
const snapParameters& snapParams
)
{
@ -1791,13 +2474,40 @@ void Foam::autoSnapDriver::doSnap
// Calculate displacement at every patch point. Insert into
// meshMover.
// Calculate displacement at every patch point
pointField nearestPoint;
vectorField nearestNormal;
if (snapParams.detectNearSurfacesSnap())
{
nearestPoint.setSize(pp.nPoints(), vector::max);
nearestNormal.setSize(pp.nPoints(), vector::zero);
}
vectorField disp = calcNearestSurface
(
meshRefiner_,
snapDist,
meshMover.patch()
pp,
nearestPoint,
nearestNormal
);
// Override displacement at thin gaps
if (snapParams.detectNearSurfacesSnap())
{
detectNearSurfaces
(
Foam::cos(degToRad(planarAngle)),// planar cos for gaps
pp,
nearestPoint, // surfacepoint from nearest test
nearestNormal, // surfacenormal from nearest test
disp
);
}
// Override displacement with feature edge attempt
if (doFeatures)
{

View File

@ -148,6 +148,7 @@ class autoSnapDriver
(
const label iter,
const indirectPrimitivePatch& pp,
const scalarField& faceSnapDist,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
labelList& faceSurfaceRegion,
@ -432,6 +433,23 @@ public:
const word& zoneName
);
//- Helper: calculate average cell centre per point
static tmp<pointField> avgCellCentres
(
const fvMesh& mesh,
const indirectPrimitivePatch&
);
//- Per patch point override displacement if in gap situation
void detectNearSurfaces
(
const scalar planarCos,
const indirectPrimitivePatch&,
const pointField& nearestPoint,
const vectorField& nearestNormal,
vectorField& disp
) const;
//- Per patch point calculate point on nearest surface. Set as
// boundary conditions of motionSmoother displacement field. Return
// displacement of patch points.
@ -439,7 +457,9 @@ public:
(
const meshRefinement& meshRefiner,
const scalarField& snapDist,
const indirectPrimitivePatch&
const indirectPrimitivePatch&,
pointField& nearestPoint,
vectorField& nearestNormal
);
////- Per patch point calculate point on nearest surface. Set as
@ -483,6 +503,7 @@ public:
const dictionary& snapDict,
const dictionary& motionDict,
const scalar featureCos,
const scalar planarAngle,
const snapParameters& snapParams
);

View File

@ -319,6 +319,7 @@ void Foam::autoSnapDriver::calcNearestFace
(
const label iter,
const indirectPrimitivePatch& pp,
const scalarField& faceSnapDist,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
labelList& faceSurfaceGlobalRegion,
@ -409,7 +410,7 @@ void Foam::autoSnapDriver::calcNearestFace
(
labelList(1, zoneSurfI),
fc,
sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
sqr(faceSnapDist),// sqr of attract dist
hitSurface,
hitInfo,
hitRegion,
@ -429,6 +430,14 @@ void Foam::autoSnapDriver::calcNearestFace
hitRegion[hitI]
);
}
else
{
WarningIn
(
"autoSnapDriver::calcNearestFace(..)"
) << "Did not find surface near face centre " << fc[hitI]
<< endl;
}
}
}
@ -467,7 +476,7 @@ void Foam::autoSnapDriver::calcNearestFace
(
unzonedSurfaces,
fc,
sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
sqr(faceSnapDist),// sqr of attract dist
hitSurface,
hitInfo,
hitRegion,
@ -487,6 +496,14 @@ void Foam::autoSnapDriver::calcNearestFace
hitRegion[hitI]
);
}
else
{
WarningIn
(
"autoSnapDriver::calcNearestFace(..)"
) << "Did not find surface near face centre " << fc[hitI]
<< endl;
}
}
@ -560,24 +577,44 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
forAll(pp.pointFaces(), pointI)
{
const labelList& pFaces = pp.pointFaces()[pointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(pFaces.size());
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(pFaces.size());
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(pFaces.size());
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(pFaces.size());
// Count valid face normals
label nFaces = 0;
forAll(pFaces, i)
{
label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI];
pDisp[i] = faceDisp[faceI];
pFc[i] = pp.faceCentres()[faceI];
if (faceSurfaceGlobalRegion[faceI] != -1)
{
nFaces++;
}
}
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(nFaces);
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(nFaces);
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(nFaces);
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(nFaces);
nFaces = 0;
forAll(pFaces, i)
{
label faceI = pFaces[i];
label globalRegionI = faceSurfaceGlobalRegion[faceI];
if (globalRegionI != -1)
{
pNormals[nFaces] = faceSurfaceNormal[faceI];
pDisp[nFaces] = faceDisp[faceI];
pFc[nFaces] = pp.faceCentres()[faceI];
//label meshFaceI = pp.addressing()[faceI];
//pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
pFid[i] = globalToMasterPatch_[faceSurfaceGlobalRegion[faceI]];
//pFid[nFaces] = mesh.boundaryMesh().whichPatch(meshFaceI);
pFid[nFaces] = globalToMasterPatch_[globalRegionI];
nFaces++;
}
}
}
@ -2123,9 +2160,85 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
}
// Collect candidate points for attraction
DynamicList<label> attractPoints(pp.nPoints());
{
const fvMesh& mesh = meshRefiner_.mesh();
boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
label nFeats = 0;
forAll(rawPatchConstraints, pointI)
{
if (rawPatchConstraints[pointI].first() >= 2)
{
isFeatureEdgeOrPoint[pointI] = true;
nFeats++;
}
}
Info<< "Intially selected " << returnReduce(nFeats, sumOp<label>())
<< " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " for reverse attraction." << endl;
// Make sure is synchronised (note: check if constraint is already
// synced in which case this is not needed here)
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
isFeatureEdgeOrPoint,
orEqOp<bool>(), // combine op
false
);
for (label nGrow = 0; nGrow < 1; nGrow++)
{
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
forAll(f, fp)
{
if (isFeatureEdgeOrPoint[f[fp]])
{
// Mark all points on face
forAll(f, fp)
{
isFeatureEdgeOrPoint[f[fp]] = true;
}
break;
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
isFeatureEdgeOrPoint,
orEqOp<bool>(), // combine op
false
);
}
// Collect attractPoints
forAll(isFeatureEdgeOrPoint, pointI)
{
if (isFeatureEdgeOrPoint[pointI])
{
attractPoints.append(pointI);
}
}
Info<< "Selected " << returnReduce(attractPoints.size(), sumOp<label>())
<< " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " for reverse attraction." << endl;
}
indexedOctree<treeDataPoint> ppTree
(
treeDataPoint(pp.localPoints()),
treeDataPoint(pp.localPoints(), attractPoints),
bb, // overall search domain
8, // maxLevel
10, // leafsize
@ -2159,7 +2272,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
if (nearInfo.hit())
{
label pointI = nearInfo.index();
label pointI =
ppTree.shapes().pointLabels()[nearInfo.index()];
const point attraction = featPt-pp.localPoints()[pointI];
// Check if this point is already being attracted. If so
@ -2220,7 +2334,9 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
if (nearInfo.hit())
{
label pointI = nearInfo.index();
label pointI =
ppTree.shapes().pointLabels()[nearInfo.index()];
const point& pt = pp.localPoints()[pointI];
const point attraction = featPt-pt;
@ -2681,6 +2797,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const fvMesh& mesh = meshRefiner_.mesh();
// Calculate attraction distance per face (from the attraction distance
// per point)
scalarField faceSnapDist(pp.size(), -GREAT);
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
forAll(f, fp)
{
faceSnapDist[faceI] = max(faceSnapDist[faceI], snapDist[f[fp]]);
}
}
// Displacement and orientation per pp face
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2695,6 +2824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
(
iter,
pp,
faceSnapDist,
faceDisp,
faceSurfaceNormal,
faceSurfaceGlobalRegion,

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,11 +40,12 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
multiRegionFeatureSnap_
(
dict.lookupOrDefault("multiRegionFeatureSnap", false)
),
detectNearSurfacesSnap_
(
dict.lookupOrDefault("detectNearSurfacesSnap", true)
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -70,6 +70,8 @@ class snapParameters
const Switch multiRegionFeatureSnap_;
const Switch detectNearSurfacesSnap_;
// Private Member Functions
@ -141,6 +143,11 @@ public:
return multiRegionFeatureSnap_;
}
Switch detectNearSurfacesSnap() const
{
return detectNearSurfacesSnap_;
}
};

View File

@ -289,16 +289,6 @@ private:
label& nRefine
) const;
//- Is local topology a small gap?
bool isGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Mark cell if local a gap topology or
bool checkProximity
(
@ -664,6 +654,26 @@ public:
// Refinement
//- Is local topology a small gap?
bool isGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Is local topology a small gap normal to the test vector
bool isNormalGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Calculate list of cells to refine.
labelList refineCandidates
(

View File

@ -1061,13 +1061,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
meshMover
);
pointField nearestPoint;
vectorField nearestNormal;
const vectorField disp
(
autoSnapDriver::calcNearestSurface
(
*this,
snapDist, // attraction
pp
pp,
nearestPoint,
nearestNormal
)
);

View File

@ -1214,7 +1214,6 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
}
// Mark small gaps
bool Foam::meshRefinement::isGap
(
const scalar planarCos,
@ -1225,12 +1224,53 @@ bool Foam::meshRefinement::isGap
const vector& normal1
) const
{
////- hits differ and angles are oppositeish
//return
// (mag(point0-point1) > mergeDistance())
// && ((normal0 & normal1) < (-1+planarCos));
//- hits differ and angles are oppositeish and
// hits have a normal distance
vector d = point1-point0;
scalar magD = mag(d);
if (magD > mergeDistance())
{
scalar cosAngle = (normal0 & normal1);
vector avg = vector::zero;
if (cosAngle < (-1+planarCos))
{
// Opposite normals
avg = 0.5*(normal0-normal1);
}
else if (cosAngle > (1-planarCos))
{
avg = 0.5*(normal0+normal1);
}
if (avg != vector::zero)
{
return true;
}
else
{
return false;
}
}
else
{
return false;
}
}
// Mark small gaps
bool Foam::meshRefinement::isNormalGap
(
const scalar planarCos,
const vector& point0,
const vector& normal0,
const vector& point1,
const vector& normal1
) const
{
//- hits differ and angles are oppositeish and
// hits have a normal distance
vector d = point1-point0;
@ -1316,7 +1356,7 @@ bool Foam::meshRefinement::checkProximity
// - different location
// - opposite surface
bool closeSurfaces = isGap
bool closeSurfaces = isNormalGap
(
planarCos,
cellMaxLocation,
@ -1565,7 +1605,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
// Have valid data on both sides. Check planarCos.
if
(
isGap
isNormalGap
(
planarCos,
cellMaxLocation[own],
@ -1635,7 +1675,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
// Have valid data on both sides. Check planarCos.
if
(
isGap
isNormalGap
(
planarCos,
cellMaxLocation[own],

View File

@ -12,12 +12,20 @@ FoamFile
class dictionary;
object snappyHexMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Which of the steps to run
castellatedMesh false;
snap false;
addLayers true;
castellatedMesh true;
snap true;
addLayers false;
foamyHexMesh
{
#include "foamyHexMeshDict"
}
// Geometry. Definition of all surfaces. All surfaces are of class
@ -28,20 +36,8 @@ addLayers true;
// - to 'snap' the mesh boundary to the surface
geometry
{
// motorBike.obj
// {
// type triSurfaceMesh;
// name motorBike;
// }
//
// refinementBox
// {
// type searchableBox;
// min (-1.0 -0.7 0.0);
// max ( 8.0 0.7 2.5);
// }
};
${:foamyHexMesh.geometry};
}
// Settings for the castellatedMesh generation.
@ -67,19 +63,12 @@ castellatedMeshControls
// few cells. This setting will cause refinement to stop if <= minimumRefine
// are selected for refinement. Note: it will at least do one iteration
// (unless the number of cells to refine is 0)
minRefinementCells 10;
// Allow a certain level of imbalance during refining
// (since balancing is quite expensive)
// Expressed as fraction of perfect balance (= overall number of cells /
// nProcs). 0=balance always.
maxLoadUnbalance 0.10;
minRefinementCells 0;
// Number of buffer layers between different levels.
// 1 means normal 2:1 refinement restriction, larger means slower
// refinement.
nCellsBetweenLevels 3;
nCellsBetweenLevels 1;
@ -90,10 +79,10 @@ castellatedMeshControls
// This is a featureEdgeMesh, read from constant/triSurface for now.
features
(
//{
// file "someLine.eMesh";
// level 2;
//}
{
file "flange.eMesh";
level 0;
}
);
@ -109,19 +98,11 @@ castellatedMeshControls
refinementSurfaces
{
// motorBike
// {
// // Surface-wise min and max refinement level
// level (5 6);
//
// // Optional specification of patch type (default is wall). No
// // constraint types (cyclic, symmetry) etc. are allowed.
// patchInfo
// {
// type wall;
// inGroups (motorBike);
// }
// }
flange
{
// Surface-wise min and max refinement level
level (2 2);
}
}
// Resolve sharp angles
@ -143,11 +124,11 @@ castellatedMeshControls
refinementRegions
{
// refinementBox
// {
// mode inside;
// levels ((1E15 4));
// }
refineHole
{
mode inside;
levels ((1E15 3));
}
}
@ -159,7 +140,7 @@ castellatedMeshControls
// section reachable from the locationInMesh is kept.
// NOTE: This point should never be on a face, always inside a cell, even
// after refinement.
locationInMesh (-1 0 0);
locationInMesh ${:foamyHexMesh.surfaceConformation.locationInMesh};
// Whether any faceZones (as specified in the refinementSurfaces)
@ -180,19 +161,31 @@ snapControls
//- Relative distance for points to be attracted by surface feature point
// or edge. True distance is this factor times local
// maximum edge length.
tolerance 4.0;
tolerance 1.0;
//- Number of mesh displacement relaxation iterations.
nSolveIter 0;
nSolveIter 300;
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
nRelaxIter 5;
//- Highly experimental and wip: number of feature edge snapping
// iterations. Leave out altogether to disable.
// Do not use here since mesh resolution too low and baffles present
//nFeatureSnapIter 10;
// Feature snapping
// Number of feature edge snapping iterations.
// Leave out altogether to disable.
nFeatureSnapIter 10;
// Detect (geometric only) features by sampling the surface
// (default=false).
implicitFeatureSnap true;
// Use castellatedMeshControls::features (default = true)
explicitFeatureSnap false;
// Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
multiRegionFeatureSnap true;
}
@ -207,14 +200,15 @@ addLayersControls
// Per final patch (so not geometry!) the layer information
layers
{
"flange.obj.*"
"flange_.*"
{
nSurfaceLayers 1;
}
}
// Expansion factor for layer mesh
expansionRatio 1.5;
expansionRatio 1.0;
// Wanted thickness of final added cell layer. If multiple layers
// is the thickness of the layer furthest away from the wall.
@ -231,9 +225,9 @@ addLayersControls
// If points get not extruded do nGrow layers of connected faces that are
// also not grown. This helps convergence of the layer addition process
// close to features.
// Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
nGrow 0;
// Advanced settings
// When not to extrude surface. 0 is flat surface, 90 is when two faces
@ -261,10 +255,8 @@ addLayersControls
maxThicknessToMedialRatio 0.3;
// Angle used to pick up medial axis points
// Note: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130 in 17x.
minMedianAxisAngle 90;
// Create buffer region for new layer terminations
nBufferCellsNoExtrude 0;
@ -273,16 +265,35 @@ addLayersControls
// if it reaches this number of iterations; possibly with an illegal
// mesh.
nLayerIter 50;
// Max number of iterations after which relaxed meshQuality controls
// get used. Up to nRelaxIter it uses the settings in meshQualityControls,
// after nRelaxIter it uses the values in meshQualityControls::relaxed.
nRelaxedIter 20;
}
// Generic mesh quality settings. At any undoable phase these determine
// where to undo.
meshQualityControls
{
// Generic mesh quality settings. At any undoable phase these determine
// where to undo.
#include "meshQualityDict"
maxNonOrtho 65;
// Optional : some meshing phases allow usage of relaxed rules.
// See e.g. addLayersControls::nRelaxedIter.
relaxed
{
//- Maximum non-orthogonality allowed. Set to 180 to disable.
maxNonOrtho 75;
}
// Advanced
//- Number of error distribution iterations
nSmoothScale 4;
//- amount to scale back displacement at error points
@ -302,7 +313,7 @@ debug 0;
// Merge tolerance. Is fraction of overall bounding box of initial mesh.
// Note: the write tolerance needs to be higher than this.
mergeTolerance 1e-6;
mergeTolerance 1E-6;
// ************************************************************************* //