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

This commit is contained in:
Henry
2013-10-21 10:02:36 +01:00
44 changed files with 889 additions and 1024 deletions

View File

@ -34,6 +34,7 @@
fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime()) fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime())
*fvc::snGrad(alpha1)*mesh.magSf() *fvc::snGrad(alpha1)*mesh.magSf()
); );
phiP1.boundaryField() == 0;
// Phase-2 pressure flux (e.g. due to particle-particle pressure) // Phase-2 pressure flux (e.g. due to particle-particle pressure)
surfaceScalarField phiP2 surfaceScalarField phiP2
@ -42,6 +43,7 @@
fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime()) fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime())
*fvc::snGrad(alpha2)*mesh.magSf() *fvc::snGrad(alpha2)*mesh.magSf()
); );
phiP2.boundaryField() == 0;
surfaceScalarField phiHbyA1 surfaceScalarField phiHbyA1
( (
@ -126,9 +128,11 @@
); );
pEqnComp1 = pEqnComp1 =
fvc::ddt(rho1) (
+ fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1) fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
+ correction - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1/rho1)*correction
( (
psi1*fvm::ddt(p) psi1*fvm::ddt(p)
+ fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p) + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
@ -137,9 +141,11 @@
pEqnComp1().relax(); pEqnComp1().relax();
pEqnComp2 = pEqnComp2 =
fvc::ddt(rho2) (
+ fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2) fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
+ correction - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2/rho2)*correction
( (
psi2*fvm::ddt(p) psi2*fvm::ddt(p)
+ fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p) + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
@ -150,12 +156,18 @@
else else
{ {
pEqnComp1 = pEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p)) (
+ fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1); fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p));
pEqnComp2 = pEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p)) (
+ fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2); fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p));
} }
// Cache p prior to solve for density update // Cache p prior to solve for density update
@ -171,11 +183,7 @@
solve solve
( (
( pEqnComp1() + pEqnComp2() + pEqnIncomp,
(alpha1/rho1)*pEqnComp1()
+ (alpha2/rho2)*pEqnComp2()
)
+ pEqnIncomp,
mesh.solver(p.select(pimple.finalInnerIter())) mesh.solver(p.select(pimple.finalInnerIter()))
); );
@ -201,8 +209,8 @@
fluid.dgdt() = fluid.dgdt() =
( (
pos(alpha2)*(pEqnComp2 & p)/rho2 pos(alpha2)*(pEqnComp2 & p)/max(alpha2, scalar(1e-3))
- pos(alpha1)*(pEqnComp1 & p)/rho1 - pos(alpha1)*(pEqnComp1 & p)/max(alpha1, scalar(1e-3))
); );
p.relax(); p.relax();

View File

@ -1016,7 +1016,6 @@ int i;
// found somewhere in the 2nd addition // found somewhere in the 2nd addition
uint32_t stroustrup (const char * s, int len) { uint32_t stroustrup (const char * s, int len) {
uint32_t h; uint32_t h;
int i;
for (register int i=0; i < len; ++s, ++i) for (register int i=0; i < len; ++s, ++i)
{ {

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -137,7 +137,7 @@ int main(int argc, char *argv[])
Info<< " std::list constructed from Foam::List: "; Info<< " std::list constructed from Foam::List: ";
std::list<vector>::iterator it; std::list<vector>::iterator it;
for (it=stlList.begin(); it != stlList.end(); it++) for (it=stlList.begin(); it != stlList.end(); ++it)
{ {
Info<< *it << " "; Info<< *it << " ";
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -79,7 +79,7 @@ int main(int argc, char *argv[])
while (iter != iter2) while (iter != iter2)
{ {
Info<< *iter; Info<< *iter;
iter++; ++iter;
} }
Info<< "<\n"; Info<< "<\n";

View File

@ -477,7 +477,9 @@ meshQualityControls
// 0 : only write final meshes // 0 : only write final meshes
// 1 : write intermediate meshes // 1 : write intermediate meshes
// 2 : write volScalarField with cellLevel for postprocessing // 2 : write volScalarField with cellLevel for postprocessing
// 4 : write current intersections as .obj files // 4 : write current mesh intersections as .obj files
// 8 : write information about explicit feature edge refinement
// 16 : write information about layers
debug 0; debug 0;
// Merge tolerance. Is fraction of overall bounding box of initial mesh. // Merge tolerance. Is fraction of overall bounding box of initial mesh.

View File

@ -1,22 +0,0 @@
#!/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

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

View File

@ -1,75 +0,0 @@
/*--------------------------------*- 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

@ -1,41 +0,0 @@
/*--------------------------------*- 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

@ -1,21 +0,0 @@
/*--------------------------------*- 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

@ -1,49 +0,0 @@
/*--------------------------------*- 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

@ -1,135 +0,0 @@
/*--------------------------------*- 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

@ -1,60 +0,0 @@
/*--------------------------------*- 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

@ -1,46 +0,0 @@
/*--------------------------------*- 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

@ -1,110 +0,0 @@
/*--------------------------------*- 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

@ -1,64 +0,0 @@
/*--------------------------------*- 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

@ -368,6 +368,8 @@ FoamFile
// { // {
// faceSet f0; // name of faceSet // faceSet f0; // name of faceSet
// cellSet c0; // name of cellSet of slave side // cellSet c0; // name of cellSet of slave side
// flip false; // optional: flip the faceZone (so now the cellSet
// // is the master side)
// } // }
// //
// // Select based on surface. Orientation from normals on surface // // Select based on surface. Orientation from normals on surface

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -78,6 +78,26 @@ public:
hexes(nCells), hexes(nCells),
polys(nCells) polys(nCells)
{} {}
// Member Functions
void setSize(const label nCells)
{
nTets = 0;
nPyrs = 0;
nPrisms = 0;
nHexesWedges = 0;
nPolys = 0;
tets.setSize(nCells);
pyrs.setSize(nCells);
prisms.setSize(nCells);
wedges.setSize(nCells);
hexes.setSize(nCells);
polys.setSize(nCells);
}
}; };

View File

@ -47,7 +47,8 @@ License
void Foam::ensightMesh::correct() void Foam::ensightMesh::correct()
{ {
patchPartOffset_ = 2; patchPartOffset_ = 2;
meshCellSets_ = mesh_.nCells(); meshCellSets_.setSize(mesh_.nCells());
boundaryFaceSets_.setSize(mesh_.boundary().size()); boundaryFaceSets_.setSize(mesh_.boundary().size());
allPatchNames_.clear(); allPatchNames_.clear();
patchNames_.clear(); patchNames_.clear();
@ -193,8 +194,6 @@ void Foam::ensightMesh::correct()
if (!noPatches_) if (!noPatches_)
{ {
forAll(mesh_.boundary(), patchi) forAll(mesh_.boundary(), patchi)
{
if (mesh_.boundary()[patchi].size())
{ {
const polyPatch& p = mesh_.boundaryMesh()[patchi]; const polyPatch& p = mesh_.boundaryMesh()[patchi];
@ -233,7 +232,6 @@ void Foam::ensightMesh::correct()
polys.setSize(nPolys); polys.setSize(nPolys);
} }
} }
}
forAll(allPatchNames_, patchi) forAll(allPatchNames_, patchi)
{ {
@ -241,14 +239,11 @@ void Foam::ensightMesh::correct()
nFacePrimitives nfp; nFacePrimitives nfp;
if (patchNames_.empty() || patchNames_.found(patchName)) if (patchNames_.empty() || patchNames_.found(patchName))
{
if (mesh_.boundary()[patchi].size())
{ {
nfp.nTris = boundaryFaceSets_[patchi].tris.size(); nfp.nTris = boundaryFaceSets_[patchi].tris.size();
nfp.nQuads = boundaryFaceSets_[patchi].quads.size(); nfp.nQuads = boundaryFaceSets_[patchi].quads.size();
nfp.nPolys = boundaryFaceSets_[patchi].polys.size(); nfp.nPolys = boundaryFaceSets_[patchi].polys.size();
} }
}
reduce(nfp.nTris, sumOp<label>()); reduce(nfp.nTris, sumOp<label>());
reduce(nfp.nQuads, sumOp<label>()); reduce(nfp.nQuads, sumOp<label>());
@ -1148,6 +1143,7 @@ void Foam::ensightMesh::write
if (nfp.nTris || nfp.nQuads || nfp.nPolys) if (nfp.nTris || nfp.nQuads || nfp.nPolys)
{ {
const polyPatch& p = mesh_.boundaryMesh()[patchi]; const polyPatch& p = mesh_.boundaryMesh()[patchi];
const labelList& tris = boundaryFaceSets_[patchi].tris; const labelList& tris = boundaryFaceSets_[patchi].tris;
const labelList& quads = boundaryFaceSets_[patchi].quads; const labelList& quads = boundaryFaceSets_[patchi].quads;
const labelList& polys = boundaryFaceSets_[patchi].polys; const labelList& polys = boundaryFaceSets_[patchi].polys;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -441,8 +441,6 @@ void pqPV3FoamReaderPanel::IncludeZonesToggled()
void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled() void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
{ {
vtkSMProperty* prop;
vtkSMIntVectorProperty::SafeDownCast vtkSMIntVectorProperty::SafeDownCast
( (
this->proxy()->GetProperty("UiExtrapolatePatches") this->proxy()->GetProperty("UiExtrapolatePatches")
@ -454,8 +452,6 @@ void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
void pqPV3FoamReaderPanel::InterpolateVolFieldsToggled() void pqPV3FoamReaderPanel::InterpolateVolFieldsToggled()
{ {
vtkSMProperty* prop;
vtkSMIntVectorProperty::SafeDownCast vtkSMIntVectorProperty::SafeDownCast
( (
this->proxy()->GetProperty("UiInterpolateVolFields") this->proxy()->GetProperty("UiInterpolateVolFields")

View File

@ -358,6 +358,59 @@ Type sum(const UList<Type>& f)
TMP_UNARY_FUNCTION(Type, sum) TMP_UNARY_FUNCTION(Type, sum)
template<class Type>
Type maxMagSqr(const UList<Type>& f)
{
if (f.size())
{
Type Max(f[0]);
TFOR_ALL_S_OP_FUNC_F_S
(
Type,
Max,
=,
maxMagSqrOp<Type>(),
Type,
f,
Type,
Max
)
return Max;
}
else
{
return pTraits<Type>::min;
}
}
TMP_UNARY_FUNCTION(Type, maxMagSqr)
template<class Type>
Type minMagSqr(const UList<Type>& f)
{
if (f.size())
{
Type Min(f[0]);
TFOR_ALL_S_OP_FUNC_F_S
(
Type,
Min,
=,
minMagSqrOp<Type>(),
Type,
f,
Type,
Min
)
return Min;
}
else
{
return pTraits<Type>::max;
}
}
TMP_UNARY_FUNCTION(Type, minMagSqr)
template<class Type> template<class Type>
scalar sumProd(const UList<Type>& f1, const UList<Type>& f2) scalar sumProd(const UList<Type>& f1, const UList<Type>& f2)
@ -488,6 +541,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
G_UNARY_FUNCTION(Type, gMax, max, max) G_UNARY_FUNCTION(Type, gMax, max, max)
G_UNARY_FUNCTION(Type, gMin, min, min) G_UNARY_FUNCTION(Type, gMin, min, min)
G_UNARY_FUNCTION(Type, gSum, sum, sum) G_UNARY_FUNCTION(Type, gSum, sum, sum)
G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum) G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum) G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum) G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)

View File

@ -171,6 +171,16 @@ Type sum(const UList<Type>& f);
TMP_UNARY_FUNCTION(Type, sum) TMP_UNARY_FUNCTION(Type, sum)
template<class Type>
Type maxMagSqr(const UList<Type>& f);
TMP_UNARY_FUNCTION(Type, maxMagSqr)
template<class Type>
Type minMagSqr(const UList<Type>& f);
TMP_UNARY_FUNCTION(Type, minMagSqr)
template<class Type> template<class Type>
scalar sumProd(const UList<Type>& f1, const UList<Type>& f2); scalar sumProd(const UList<Type>& f1, const UList<Type>& f2);
@ -208,6 +218,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
G_UNARY_FUNCTION(Type, gMax, max, max) G_UNARY_FUNCTION(Type, gMax, max, max)
G_UNARY_FUNCTION(Type, gMin, min, min) G_UNARY_FUNCTION(Type, gMin, min, min)
G_UNARY_FUNCTION(Type, gSum, sum, sum) G_UNARY_FUNCTION(Type, gSum, sum, sum)
G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum) G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum) G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum) G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)

View File

@ -1963,7 +1963,7 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
pointField sharedPoints(mesh_.points(), sharedPointLabels()); pointField sharedPoints(mesh_.points(), sharedPointLabels());
// Append from all processors // Append from all processors
combineReduce(sharedPoints, plusEqOp<pointField>()); combineReduce(sharedPoints, ListPlusEqOp<pointField>());
// Merge tolerance // Merge tolerance
scalar tolDim = matchTol_ * mesh_.bounds().mag(); scalar tolDim = matchTol_ * mesh_.bounds().mag();

View File

@ -109,29 +109,6 @@ class globalMeshData
public processorTopology public processorTopology
{ {
// Private class
// To combineReduce a pointField. Just appends all lists.
template<class T>
class plusEqOp
{
public:
void operator()(T& x, const T& y) const
{
label n = x.size();
x.setSize(x.size() + y.size());
forAll(y, i)
{
x[n++] = y[i];
}
}
};
// Private data // Private data
//- Reference to mesh //- Reference to mesh
@ -337,6 +314,29 @@ class globalMeshData
public: public:
// Public class
// To combineReduce a List. Just appends all lists.
template<class T>
class ListPlusEqOp
{
public:
void operator()(T& x, const T& y) const
{
label n = x.size();
x.setSize(x.size() + y.size());
forAll(y, i)
{
x[n++] = y[i];
}
}
};
//- Runtime type information //- Runtime type information
ClassName("globalMeshData"); ClassName("globalMeshData");

View File

@ -1446,9 +1446,9 @@ void Foam::autoLayerDriver::syncPatchDisplacement
} }
} }
Info<< "Prevented extrusion on " //Info<< "Prevented extrusion on "
<< returnReduce(nChangedTotal, sumOp<label>()) // << returnReduce(nChangedTotal, sumOp<label>())
<< " coupled patch points during syncPatchDisplacement." << endl; // << " coupled patch points during syncPatchDisplacement." << endl;
} }

View File

@ -524,7 +524,7 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
label nOldPointCounter = nPointCounter; //label nOldPointCounter = nPointCounter;
// Detect situation where two featureedge-neighbouring faces are partly or // Detect situation where two featureedge-neighbouring faces are partly or
// not extruded and the edge itself is extruded. In this case unmark the // not extruded and the edge itself is extruded. In this case unmark the
@ -633,9 +633,9 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
} }
} }
Info<< "Added " //Info<< "Added "
<< returnReduce(nPointCounter-nOldPointCounter, sumOp<label>()) // << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
<< " point not to extrude." << endl; // << " point not to extrude." << endl;
} }

View File

@ -910,6 +910,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
Info<< "Detecting near surfaces ..." << endl; Info<< "Detecting near surfaces ..." << endl;
const pointField& localPoints = pp.localPoints(); const pointField& localPoints = pp.localPoints();
const labelList& meshPoints = pp.meshPoints();
const refinementSurfaces& surfaces = meshRefiner_.surfaces(); const refinementSurfaces& surfaces = meshRefiner_.surfaces();
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
@ -1220,6 +1221,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
} }
const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
label nOverride = 0; label nOverride = 0;
// 1. All points to non-interface surfaces // 1. All points to non-interface surfaces
@ -1332,7 +1334,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
} }
} }
if (override) if (override && isMasterPoint[meshPoints[pointI]])
{ {
nOverride++; nOverride++;
} }
@ -1399,8 +1401,6 @@ void Foam::autoSnapDriver::detectNearSurfaces
); );
label nOverride = 0;
forAll(hit1, i) forAll(hit1, i)
{ {
label pointI = zonePointIndices[i]; label pointI = zonePointIndices[i];
@ -1472,7 +1472,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
} }
} }
if (override) if (override && isMasterPoint[meshPoints[pointI]])
{ {
nOverride++; nOverride++;
} }
@ -1690,7 +1690,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
scalarField magDisp(mag(patchDisp)); scalarField magDisp(mag(patchDisp));
Info<< "Wanted displacement : average:" Info<< "Wanted displacement : average:"
<< gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>()) << meshRefinement::gAverage
(
mesh,
syncTools::getMasterPoints(mesh),
pp.meshPoints(),
magDisp
)
<< " min:" << gMin(magDisp) << " min:" << gMin(magDisp)
<< " max:" << gMax(magDisp) << endl; << " max:" << gMax(magDisp) << endl;
} }
@ -2548,25 +2554,30 @@ void Foam::autoSnapDriver::doSnap
adaptPatchIDs adaptPatchIDs
) )
); );
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface // Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp)); const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
// Construct iterative mesh mover. // Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl; Info<< "Constructing mesh displacer ..." << endl;
Info<< "Using mesh parameters " << motionDict << nl << endl; Info<< "Using mesh parameters " << motionDict << nl << endl;
const pointMesh& pMesh = pointMesh::New(mesh); autoPtr<motionSmoother> meshMoverPtr
(
motionSmoother meshMover new motionSmoother
( (
mesh, mesh,
pp, ppPtr(),
adaptPatchIDs, adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs), meshRefinement::makeDisplacementField
(
pointMesh::New(mesh),
adaptPatchIDs
),
motionDict motionDict
)
); );
@ -2595,16 +2606,95 @@ void Foam::autoSnapDriver::doSnap
snapParams, snapParams,
nInitErrors, nInitErrors,
baffles, baffles,
meshMover meshMoverPtr()
); );
//- Only if in feature attraction mode:
// Nearest feature
vectorField patchAttraction;
// Constraints at feature
List<pointConstraint> patchConstraints;
for (label iter = 0; iter < nFeatIter; iter++) for (label iter = 0; iter < nFeatIter; iter++)
{ {
Info<< nl Info<< nl
<< "Morph iteration " << iter << nl << "Morph iteration " << iter << nl
<< "-----------------" << endl; << "-----------------" << endl;
//if (iter > 0 && iter == nFeatIter/2)
//{
// Info<< "Splitting diagonal attractions" << endl;
// const labelList& bFaces = ppPtr().addressing();
//
// DynamicList<label> splitFaces(bFaces.size());
// DynamicList<labelPair> splits(bFaces.size());
//
// forAll(bFaces, faceI)
// {
// const labelPair split
// (
// findDiagonalAttraction
// (
// ppPtr(),
// patchAttraction,
// patchConstraints,
// faceI
// )
// );
//
// if (split != labelPair(-1, -1))
// {
// splitFaces.append(bFaces[faceI]);
// splits.append(split);
// }
// }
//
// autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
// (
// splitFaces,
// splits
// );
//
// const labelList& faceMap = mapPtr().faceMap();
// meshRefinement::updateList(faceMap, -1, duplicateFace);
// const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
// forAll(baffles, i)
// {
// labelPair& baffle = baffles[i];
// baffle.first() = reverseFaceMap[baffle.first()];
// baffle.second() = reverseFaceMap[baffle.second()];
// }
//
// meshMoverPtr.clear();
// ppPtr.clear();
//
// ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
// meshMoverPtr.reset
// (
// new motionSmoother
// (
// mesh,
// ppPtr(),
// adaptPatchIDs,
// meshRefinement::makeDisplacementField
// (
// pointMesh::New(mesh),
// adaptPatchIDs
// ),
// motionDict
// )
// );
//}
indirectPrimitivePatch& pp = ppPtr();
motionSmoother& meshMover = meshMoverPtr();
// Calculate displacement at every patch point. Insert into // Calculate displacement at every patch point. Insert into
// meshMover. // meshMover.
// Calculate displacement at every patch point // Calculate displacement at every patch point
@ -2652,7 +2742,9 @@ void Foam::autoSnapDriver::doSnap
scalar(iter+1)/nFeatIter, scalar(iter+1)/nFeatIter,
snapDist, snapDist,
disp, disp,
meshMover meshMover,
patchAttraction,
patchConstraints
); );
} }

View File

@ -125,7 +125,9 @@ class autoSnapDriver
void smoothAndConstrain void smoothAndConstrain
( (
const PackedBoolList& isMasterEdge,
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
const labelList& meshEdges,
const List<pointConstraint>& constraints, const List<pointConstraint>& constraints,
vectorField& disp vectorField& disp
) const; ) const;
@ -196,6 +198,16 @@ class autoSnapDriver
List<pointConstraint>& patchConstraints List<pointConstraint>& patchConstraints
) const; ) const;
//- Detect any diagonal attraction. Returns indices in face
// or (-1, -1) if none
labelPair findDiagonalAttraction
(
const indirectPrimitivePatch& pp,
const vectorField& patchAttraction,
const List<pointConstraint>& patchConstraints,
const label faceI
) const;
//- Return hit if on multiple points //- Return hit if on multiple points
pointIndexHit findMultiPatchPoint pointIndexHit findMultiPatchPoint
( (
@ -371,7 +383,9 @@ class autoSnapDriver
const scalar featureAttract, const scalar featureAttract,
const scalarField& snapDist, const scalarField& snapDist,
const vectorField& nearestDisp, const vectorField& nearestDisp,
motionSmoother& meshMover motionSmoother& meshMover,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const; ) const;

View File

@ -38,50 +38,18 @@ License
#include "treeDataPoint.H" #include "treeDataPoint.H"
#include "indexedOctree.H" #include "indexedOctree.H"
#include "snapParameters.H" #include "snapParameters.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam namespace Foam
{ {
class listTransform
{
public:
void operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<List<point> >& fld
) const
{
const tensor T
(
forward
? vt.R()
: vt.R().T()
);
forAll(fld, i)
{
List<point>& elems = fld[i];
forAll(elems, elemI)
{
elems[elemI] = transform(T, elems[elemI]);
}
}
}
};
template<class T> template<class T>
class listPlusEqOp class listPlusEqOp
{ {
public: public:
void operator() void operator()(List<T>& x, const List<T>& y) const
(
List<T>& x,
const List<T>& y
) const
{ {
label sz = x.size(); label sz = x.size();
x.setSize(sz+y.size()); x.setSize(sz+y.size());
@ -159,7 +127,9 @@ bool Foam::autoSnapDriver::isFeaturePoint
void Foam::autoSnapDriver::smoothAndConstrain void Foam::autoSnapDriver::smoothAndConstrain
( (
const PackedBoolList& isMasterEdge,
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
const labelList& meshEdges,
const List<pointConstraint>& constraints, const List<pointConstraint>& constraints,
vectorField& disp vectorField& disp
) const ) const
@ -196,6 +166,10 @@ void Foam::autoSnapDriver::smoothAndConstrain
if (nConstraints <= 1) if (nConstraints <= 1)
{ {
forAll(pEdges, i) forAll(pEdges, i)
{
label edgeI = pEdges[i];
if (isMasterEdge[meshEdges[edgeI]])
{ {
label nbrPointI = edges[pEdges[i]].otherVertex(pointI); label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
if (constraints[nbrPointI].first() >= nConstraints) if (constraints[nbrPointI].first() >= nConstraints)
@ -206,6 +180,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
} }
} }
} }
}
syncTools::syncPointList syncTools::syncPointList
( (
@ -564,6 +539,9 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
{ {
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
// For now just get all surrounding face data. Expensive - should just // For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only // store and sync data on coupled points only
// (see e.g PatchToolsNormals.C) // (see e.g PatchToolsNormals.C)
@ -583,7 +561,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
forAll(pFaces, i) forAll(pFaces, i)
{ {
label faceI = pFaces[i]; label faceI = pFaces[i];
if (faceSurfaceGlobalRegion[faceI] != -1) if (isMasterFace[faceI] && faceSurfaceGlobalRegion[faceI] != -1)
{ {
nFaces++; nFaces++;
} }
@ -605,7 +583,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
label faceI = pFaces[i]; label faceI = pFaces[i];
label globalRegionI = faceSurfaceGlobalRegion[faceI]; label globalRegionI = faceSurfaceGlobalRegion[faceI];
if (globalRegionI != -1) if (isMasterFace[faceI] && globalRegionI != -1)
{ {
pNormals[nFaces] = faceSurfaceNormal[faceI]; pNormals[nFaces] = faceSurfaceNormal[faceI];
pDisp[nFaces] = faceDisp[faceI]; pDisp[nFaces] = faceDisp[faceI];
@ -694,7 +672,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceSurfNormals, pointFaceSurfNormals,
listPlusEqOp<point>(), listPlusEqOp<point>(),
List<point>(), List<point>(),
listTransform() mapDistribute::transform()
); );
syncTools::syncPointList syncTools::syncPointList
( (
@ -703,7 +681,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceDisp, pointFaceDisp,
listPlusEqOp<point>(), listPlusEqOp<point>(),
List<point>(), List<point>(),
listTransform() mapDistribute::transform()
); );
syncTools::syncPointList syncTools::syncPointList
( (
@ -712,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceCentres, pointFaceCentres,
listPlusEqOp<point>(), listPlusEqOp<point>(),
List<point>(), List<point>(),
listTransform() mapDistribute::transformPosition()
); );
syncTools::syncPointList syncTools::syncPointList
( (
@ -722,6 +700,25 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
listPlusEqOp<label>(), listPlusEqOp<label>(),
List<label>() List<label>()
); );
// Sort the data according to the face centres. This is only so we get
// consistent behaviour serial and parallel.
labelList visitOrder;
forAll(pointFaceDisp, pointI)
{
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
sortedOrder(mag(pFc)(), visitOrder);
pNormals = List<point>(pNormals, visitOrder);
pDisp = List<point>(pDisp, visitOrder);
pFc = List<point>(pFc, visitOrder);
pFid = UIndirectList<label>(pFid, visitOrder);
}
} }
@ -1360,6 +1357,70 @@ void Foam::autoSnapDriver::stringFeatureEdges
} }
// If only two attractions and across face return the face indices
Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
(
const indirectPrimitivePatch& pp,
const vectorField& patchAttraction,
const List<pointConstraint>& patchConstraints,
const label faceI
) const
{
const face& f = pp.localFaces()[faceI];
// For now just detect any attraction. Improve this to look at
// actual attraction position and orientation
labelPair attractIndices(-1, -1);
if (f.size() >= 4)
{
forAll(f, fp)
{
label pointI = f[fp];
if (patchConstraints[pointI].first() >= 2)
{
// Attract to feature edge or point
if (attractIndices[0] == -1)
{
// First attraction. Store
attractIndices[0] = fp;
}
else if (attractIndices[1] == -1)
{
// Second attraction. Check if not consecutive to first
// attraction
label fp0 = attractIndices[0];
if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
{
// Consecutive. Skip.
attractIndices = labelPair(-1, -1);
break;
}
else
{
attractIndices[1] = fp;
}
}
else
{
// More than two attractions. Skip.
attractIndices = labelPair(-1, -1);
break;
}
}
}
if (attractIndices[1] == -1)
{
// Found only one attraction. Skip.
attractIndices = labelPair(-1, -1);
}
}
return attractIndices;
}
Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
( (
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
@ -1929,7 +1990,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
edgeFaceNormals, edgeFaceNormals,
listPlusEqOp<point>(), listPlusEqOp<point>(),
List<point>(), List<point>(),
listTransform() mapDistribute::transform()
); );
} }
@ -2778,7 +2839,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const scalar featureAttract, const scalar featureAttract,
const scalarField& snapDist, const scalarField& snapDist,
const vectorField& nearestDisp, const vectorField& nearestDisp,
motionSmoother& meshMover motionSmoother& meshMover,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const ) const
{ {
const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap(); const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
@ -2851,7 +2914,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal // - faceSurfaceNormal
// - faceDisp // - faceDisp
// - faceCentres&faceNormal // - faceCentres
List<List<point> > pointFaceSurfNormals(pp.nPoints()); List<List<point> > pointFaceSurfNormals(pp.nPoints());
List<List<point> > pointFaceDisp(pp.nPoints()); List<List<point> > pointFaceDisp(pp.nPoints());
List<List<point> > pointFaceCentres(pp.nPoints()); List<List<point> > pointFaceCentres(pp.nPoints());
@ -2884,10 +2947,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// here. // here.
// Nearest feature // Nearest feature
vectorField patchAttraction(localPoints.size(), vector::zero); patchAttraction.setSize(localPoints.size());
patchAttraction = vector::zero;
// Constraints at feature // Constraints at feature
List<pointConstraint> patchConstraints(localPoints.size()); patchConstraints.setSize(localPoints.size());
patchConstraints = pointConstraint();
if (implicitFeatureAttraction) if (implicitFeatureAttraction)
{ {
@ -2951,15 +3015,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
patchConstraints patchConstraints
); );
const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
{
vector avgPatchDisp = meshRefinement::gAverage
(
mesh,
isMasterPoint,
pp.meshPoints(),
patchDisp
);
vector avgPatchAttr = meshRefinement::gAverage
(
mesh,
isMasterPoint,
pp.meshPoints(),
patchAttraction
);
Info<< "Attraction:" << endl Info<< "Attraction:" << endl
<< " linear : max:" << gMax(patchDisp) << " linear : max:" << gMaxMagSqr(patchDisp)
<< " avg:" << gAverage(patchDisp) << " avg:" << avgPatchDisp << endl
<< endl << " feature : max:" << gMaxMagSqr(patchAttraction)
<< " feature : max:" << gMax(patchAttraction) << " avg:" << avgPatchAttr << endl;
<< " avg:" << gAverage(patchAttraction) }
<< endl;
// So now we have: // So now we have:
// - patchDisp : point movement to go to nearest point on surface // - patchDisp : point movement to go to nearest point on surface
@ -2986,12 +3064,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// Count // Count
{ {
const labelList& meshPoints = pp.meshPoints();
label nMasterPoints = 0;
label nPlanar = 0; label nPlanar = 0;
label nEdge = 0; label nEdge = 0;
label nPoint = 0; label nPoint = 0;
forAll(patchConstraints, pointI) forAll(patchConstraints, pointI)
{ {
if (isMasterPoint[meshPoints[pointI]])
{
nMasterPoints++;
if (patchConstraints[pointI].first() == 1) if (patchConstraints[pointI].first() == 1)
{ {
nPlanar++; nPlanar++;
@ -3005,18 +3090,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
nPoint++; nPoint++;
} }
} }
}
label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>()); reduce(nMasterPoints, sumOp<label>());
reduce(nPlanar, sumOp<label>()); reduce(nPlanar, sumOp<label>());
reduce(nEdge, sumOp<label>()); reduce(nEdge, sumOp<label>());
reduce(nPoint, sumOp<label>()); reduce(nPoint, sumOp<label>());
Info<< "Feature analysis : total points:" Info<< "Feature analysis : total master points:"
<< nTotPoints << nMasterPoints
<< " attraction to :" << nl << " attraction to :" << nl
<< " feature point : " << nPoint << nl << " feature point : " << nPoint << nl
<< " feature edge : " << nEdge << nl << " feature edge : " << nEdge << nl
<< " nearest surface : " << nPlanar << nl << " nearest surface : " << nPlanar << nl
<< " rest : " << nTotPoints-nPoint-nEdge-nPlanar << " rest : " << nMasterPoints-nPoint-nEdge-nPlanar
<< nl << nl
<< endl; << endl;
} }
@ -3035,11 +3121,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
if (featureAttract < 1-0.001) if (featureAttract < 1-0.001)
{ {
const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
const vectorField pointNormals
(
PatchTools::pointNormals
(
mesh,
pp
)
);
const labelList meshEdges
(
pp.meshEdges(mesh.edges(), mesh.pointEdges())
);
// 1. Smoothed all displacement // 1. Smoothed all displacement
vectorField smoothedPatchDisp = patchDisp; vectorField smoothedPatchDisp = patchDisp;
smoothAndConstrain smoothAndConstrain
( (
isMasterEdge,
pp, pp,
meshEdges,
patchConstraints, patchConstraints,
smoothedPatchDisp smoothedPatchDisp
); );
@ -3047,16 +3151,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// 2. Smoothed tangential component // 2. Smoothed tangential component
vectorField tangPatchDisp = patchDisp; vectorField tangPatchDisp = patchDisp;
tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals(); tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
smoothAndConstrain smoothAndConstrain
( (
isMasterEdge,
pp, pp,
meshEdges,
patchConstraints, patchConstraints,
tangPatchDisp tangPatchDisp
); );
// Re-add normal component // Re-add normal component
tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals(); tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
if (debug&meshRefinement::OBJINTERSECTIONS) if (debug&meshRefinement::OBJINTERSECTIONS)
{ {

View File

@ -259,23 +259,14 @@ private:
label& nRefine label& nRefine
) const; ) const;
//- Mark cell if local curvature > curvature or //- Helper: count number of normals1 that are in normals2
// markDifferingRegions = true and intersections with different label countMatches
// regions.
bool checkCurvature
( (
const scalar curvature, const List<point>& normals1,
const label nAllowRefine, const List<point>& normals2,
const label surfaceLevel, const scalar tol = 1e-6
const vector& surfaceNormal,
const label cellI,
label& cellMaxLevel,
vector& cellMaxNormal,
labelList& refineCell,
label& nRefine
) const; ) const;
//- Mark cells for surface curvature based refinement. Marks if //- Mark cells for surface curvature based refinement. Marks if
// local curvature > curvature or if on different regions // local curvature > curvature or if on different regions
// (markDifferingRegions) // (markDifferingRegions)

View File

@ -38,7 +38,74 @@ License
#include "featureEdgeMesh.H" #include "featureEdgeMesh.H"
#include "Cloud.H" #include "Cloud.H"
//#include "globalIndex.H" //#include "globalIndex.H"
//#include "OBJstream.H" #include "OBJstream.H"
#include "cellSet.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//- To compare normals
static bool less(const vector& x, const vector& y)
{
for (direction i = 0; i < vector::nComponents; i++)
{
if (x[i] < y[i])
{
return true;
}
else if (x[i] > y[i])
{
return false;
}
}
// All components the same
return false;
}
//- To compare normals
class normalLess
{
const vectorList& values_;
public:
normalLess(const vectorList& values)
:
values_(values)
{}
bool operator()(const label a, const label b) const
{
return less(values_[a], values_[b]);
}
};
//- template specialization for pTraits<labelList> so we can have fields
template<>
class pTraits<labelList>
{
public:
//- Component type
typedef labelList cmptType;
};
//- template specialization for pTraits<labelList> so we can have fields
template<>
class pTraits<vectorList>
{
public:
//- Component type
typedef vectorList cmptType;
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -954,64 +1021,33 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
} }
// Checks if multiple intersections of a cell (by a surface with a higher // Count number of matches of first argument in second argument
// max than the cell level) and if so if the normals at these intersections Foam::label Foam::meshRefinement::countMatches
// make a large angle.
// Returns false if the nRefine limit has been reached, true otherwise.
bool Foam::meshRefinement::checkCurvature
( (
const scalar curvature, const List<point>& normals1,
const label nAllowRefine, const List<point>& normals2,
const scalar tol
const label surfaceLevel, // current intersection max level
const vector& surfaceNormal,// current intersection normal
const label cellI,
label& cellMaxLevel, // cached max surface level for this cell
vector& cellMaxNormal, // cached surface normal for this cell
labelList& refineCell,
label& nRefine
) const ) const
{ {
const labelList& cellLevel = meshCutter_.cellLevel(); label nMatches = 0;
// Test if surface applicable forAll(normals1, i)
if (surfaceLevel > cellLevel[cellI])
{ {
if (cellMaxLevel == -1) const vector& n1 = normals1[i];
{
// First visit of cell. Store
cellMaxLevel = surfaceLevel;
cellMaxNormal = surfaceNormal;
}
else
{
// Second or more visit. Check curvature.
if ((cellMaxNormal & surfaceNormal) < curvature)
{
return markForRefine
(
surfaceLevel, // mark with any non-neg number.
nAllowRefine,
refineCell[cellI],
nRefine
);
}
// Set normal to that of highest surface. Not really necessary forAll(normals2, j)
// over here but we reuse cellMax info when doing coupled faces.
if (surfaceLevel > cellMaxLevel)
{ {
cellMaxLevel = surfaceLevel; const vector& n2 = normals2[j];
cellMaxNormal = surfaceNormal;
if (magSqr(n1-n2) < tol)
{
nMatches++;
break;
} }
} }
} }
// Did not reach refinement limit. return nMatches;
return true;
} }
@ -1039,6 +1075,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// on a different surface gets refined (if its current level etc.) // on a different surface gets refined (if its current level etc.)
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
// Collect candidate faces (i.e. intersecting any surface and // Collect candidate faces (i.e. intersecting any surface and
// owner/neighbour not yet refined. // owner/neighbour not yet refined.
labelList testFaces(getRefineCandidateFaces(refineCell)); labelList testFaces(getRefineCandidateFaces(refineCell));
@ -1068,6 +1107,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
start[i] = cellCentres[own]; start[i] = cellCentres[own];
end[i] = neiCc[bFaceI]; end[i] = neiCc[bFaceI];
if (!isMasterFace[faceI])
{
Swap(start[i], end[i]);
}
minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]); minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
} }
} }
@ -1081,10 +1126,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// Test for all intersections (with surfaces of higher max level than // Test for all intersections (with surfaces of higher max level than
// minLevel) and cache per cell the max surface level and the local normal // minLevel) and cache per cell the interesting inter
// on that surface. labelListList cellSurfLevels(mesh_.nCells());
labelList cellMaxLevel(mesh_.nCells(), -1); List<vectorList> cellSurfNormals(mesh_.nCells());
vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
{ {
// Per segment the normals of the surfaces // Per segment the normals of the surfaces
@ -1104,12 +1148,29 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
surfaceNormal, surfaceNormal,
surfaceLevel surfaceLevel
); );
// Sort the data according to surface location. This will guarantee
// that on coupled faces both sides visit the intersections in
// the same order so will decide the same
labelList visitOrder;
forAll(surfaceNormal, pointI)
{
vectorList& pNormals = surfaceNormal[pointI];
labelList& pLevel = surfaceLevel[pointI];
sortedOrder(pNormals, visitOrder, normalLess(pNormals));
pNormals = List<point>(pNormals, visitOrder);
pLevel = UIndirectList<label>(pLevel, visitOrder);
}
// Clear out unnecessary data // Clear out unnecessary data
start.clear(); start.clear();
end.clear(); end.clear();
minLevel.clear(); minLevel.clear();
// Extract per cell information on the surface with the highest max // Convert face-wise data to cell.
forAll(testFaces, i) forAll(testFaces, i)
{ {
label faceI = testFaces[i]; label faceI = testFaces[i];
@ -1118,91 +1179,87 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
const vectorList& fNormals = surfaceNormal[i]; const vectorList& fNormals = surfaceNormal[i];
const labelList& fLevels = surfaceLevel[i]; const labelList& fLevels = surfaceLevel[i];
forAll(fLevels, hitI) forAll(fNormals, hitI)
{ {
checkCurvature if (fLevels[hitI] > cellLevel[own])
( {
curvature, cellSurfLevels[own].append(fLevels[hitI]);
nAllowRefine, cellSurfNormals[own].append(fNormals[hitI]);
fLevels[hitI],
fNormals[hitI],
own,
cellMaxLevel[own],
cellMaxNormal[own],
refineCell,
nRefine
);
} }
if (mesh_.isInternalFace(faceI)) if (mesh_.isInternalFace(faceI))
{ {
label nei = mesh_.faceNeighbour()[faceI]; label nei = mesh_.faceNeighbour()[faceI];
if (fLevels[hitI] > cellLevel[nei])
forAll(fLevels, hitI)
{ {
checkCurvature cellSurfLevels[nei].append(fLevels[hitI]);
cellSurfNormals[nei].append(fNormals[hitI]);
}
}
}
}
}
// Bit of statistics
{
label nSet = 0;
label nNormals = 9;
forAll(cellSurfNormals, cellI)
{
const vectorList& normals = cellSurfNormals[cellI];
if (normals.size())
{
nSet++;
nNormals += normals.size();
}
}
reduce(nSet, sumOp<label>());
reduce(nNormals, sumOp<label>());
Info<< "markSurfaceCurvatureRefinement :"
<< " cells:" << mesh_.globalData().nTotalCells()
<< " of which with normals:" << nSet
<< " ; total normals stored:" << nNormals
<< endl;
}
bool reachedLimit = false;
// 1. Check normals per cell
// ~~~~~~~~~~~~~~~~~~~~~~~~~
for
( (
curvature, label cellI = 0;
nAllowRefine, !reachedLimit && cellI < cellSurfNormals.size();
cellI++
fLevels[hitI], )
fNormals[hitI],
nei,
cellMaxLevel[nei],
cellMaxNormal[nei],
refineCell,
nRefine
);
}
}
}
}
// 2. Find out a measure of surface curvature and region edges.
// Send over surface region and surface normal to neighbour cell.
labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
{ {
label bFaceI = faceI-mesh_.nInternalFaces(); const vectorList& normals = cellSurfNormals[cellI];
label own = mesh_.faceOwner()[faceI]; const labelList& levels = cellSurfLevels[cellI];
neiBndMaxLevel[bFaceI] = cellMaxLevel[own]; // n^2 comparison of all normals in a cell
neiBndMaxNormal[bFaceI] = cellMaxNormal[own]; for (label i = 0; !reachedLimit && i < normals.size(); i++)
}
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
// Loop over all faces. Could only be checkFaces.. except if they're coupled
// Internal faces
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{ {
label own = mesh_.faceOwner()[faceI]; for (label j = i+1; !reachedLimit && j < normals.size(); j++)
label nei = mesh_.faceNeighbour()[faceI]; {
if ((normals[i] & normals[j]) < curvature)
{
label maxLevel = max(levels[i], levels[j]);
if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1) if (cellLevel[cellI] < maxLevel)
{
// Have valid data on both sides. Check curvature.
if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
{
// See which side to refine
if (cellLevel[own] < cellMaxLevel[own])
{ {
if if
( (
!markForRefine !markForRefine
( (
cellMaxLevel[own], maxLevel,
nAllowRefine, nAllowRefine,
refineCell[own], refineCell[cellI],
nRefine nRefine
) )
) )
@ -1213,17 +1270,93 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
<< " limit of " << mesh_.nCells()+7*nRefine << " limit of " << mesh_.nCells()+7*nRefine
<< endl; << endl;
} }
reachedLimit = true;
break; break;
} }
} }
}
}
}
}
if (cellLevel[nei] < cellMaxLevel[nei])
// 2. Find out a measure of surface curvature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Look at normals between neighbouring surfaces
// Loop over all faces. Could only be checkFaces, except if they're coupled
// Internal faces
for
(
label faceI = 0;
!reachedLimit && faceI < mesh_.nInternalFaces();
faceI++
)
{
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const vectorList& neiNormals = cellSurfNormals[nei];
const labelList& neiLevels = cellSurfLevels[nei];
// Special case: owner normals set is a subset of the neighbour
// normals. Do not do curvature refinement since other cell's normals
// do not add any information. Avoids weird corner extensions of
// refinement regions.
bool ownIsSubset =
countMatches(ownNormals, neiNormals)
== ownNormals.size();
bool neiIsSubset =
countMatches(neiNormals, ownNormals)
== neiNormals.size();
if (!ownIsSubset && !neiIsSubset)
{
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{
// Have valid data on both sides. Check curvature.
if ((ownNormals[i] & neiNormals[j]) < curvature)
{
// See which side to refine.
if (cellLevel[own] < ownLevels[i])
{ {
if if
( (
!markForRefine !markForRefine
( (
cellMaxLevel[nei], ownLevels[i],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching"
<< " my cell limit of "
<< mesh_.nCells()+7*nRefine << endl;
}
reachedLimit = true;
break;
}
}
if (cellLevel[nei] < neiLevels[j])
{
if
(
!markForRefine
(
neiLevels[j],
nAllowRefine, nAllowRefine,
refineCell[nei], refineCell[nei],
nRefine nRefine
@ -1232,32 +1365,70 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
{ {
if (debug) if (debug)
{ {
Pout<< "Stopped refining since reaching my cell" Pout<< "Stopped refining since reaching"
<< " limit of " << mesh_.nCells()+7*nRefine << " my cell limit of "
<< endl; << mesh_.nCells()+7*nRefine << endl;
} }
reachedLimit = true;
break; break;
} }
} }
} }
} }
} }
}
}
// Send over surface normal to neighbour cell.
List<vectorList> neiSurfaceNormals;
syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
// Boundary faces // Boundary faces
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++) for
(
label faceI = mesh_.nInternalFaces();
!reachedLimit && faceI < mesh_.nFaces();
faceI++
)
{ {
label own = mesh_.faceOwner()[faceI]; label own = mesh_.faceOwner()[faceI];
label bFaceI = faceI - mesh_.nInternalFaces(); label bFaceI = faceI - mesh_.nInternalFaces();
if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1) const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
// Special case: owner normals set is a subset of the neighbour
// normals. Do not do curvature refinement since other cell's normals
// do not add any information. Avoids weird corner extensions of
// refinement regions.
bool ownIsSubset =
countMatches(ownNormals, neiNormals)
== ownNormals.size();
bool neiIsSubset =
countMatches(neiNormals, ownNormals)
== neiNormals.size();
if (!ownIsSubset && !neiIsSubset)
{
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{ {
// Have valid data on both sides. Check curvature. // Have valid data on both sides. Check curvature.
if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature) if ((ownNormals[i] & neiNormals[j]) < curvature)
{
if (cellLevel[own] < ownLevels[i])
{ {
if if
( (
!markForRefine !markForRefine
( (
cellMaxLevel[own], ownLevels[i],
nAllowRefine, nAllowRefine,
refineCell[own], refineCell[own],
nRefine nRefine
@ -1266,15 +1437,21 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
{ {
if (debug) if (debug)
{ {
Pout<< "Stopped refining since reaching my cell" Pout<< "Stopped refining since reaching"
<< " limit of " << mesh_.nCells()+7*nRefine << " my cell limit of "
<< mesh_.nCells()+7*nRefine
<< endl; << endl;
} }
reachedLimit = true;
break; break;
} }
} }
} }
} }
}
}
}
if if
( (

View File

@ -241,7 +241,7 @@ void meshRefinement::collectAndPrint
Pstream::blocking Pstream::blocking
); );
Field<T> allData; List<T> allData;
globalPoints.gather globalPoints.gather
( (
Pstream::worldComm, Pstream::worldComm,

View File

@ -720,7 +720,7 @@ void Foam::refinementSurfaces::findAllHigherIntersections
labelList pRegions; labelList pRegions;
vectorField pNormals; vectorField pNormals;
forAll(surfaces(), surfI) forAll(surfaces_, surfI)
{ {
const searchableSurface& surface = allGeometry_[surfaces_[surfI]]; const searchableSurface& surface = allGeometry_[surfaces_[surfI]];

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -56,12 +56,14 @@ Foam::setsToFaceZone::setsToFaceZone
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& faceSetName, const word& faceSetName,
const word& cellSetName const word& cellSetName,
const Switch& flip
) )
: :
topoSetSource(mesh), topoSetSource(mesh),
faceSetName_(faceSetName), faceSetName_(faceSetName),
cellSetName_(cellSetName) cellSetName_(cellSetName),
flip_(flip)
{} {}
@ -74,7 +76,8 @@ Foam::setsToFaceZone::setsToFaceZone
: :
topoSetSource(mesh), topoSetSource(mesh),
faceSetName_(dict.lookup("faceSet")), faceSetName_(dict.lookup("faceSet")),
cellSetName_(dict.lookup("cellSet")) cellSetName_(dict.lookup("cellSet")),
flip_(dict.lookupOrDefault("flip", false))
{} {}
@ -87,7 +90,8 @@ Foam::setsToFaceZone::setsToFaceZone
: :
topoSetSource(mesh), topoSetSource(mesh),
faceSetName_(checkIs(is)), faceSetName_(checkIs(is)),
cellSetName_(checkIs(is)) cellSetName_(checkIs(is)),
flip_(false)
{} {}
@ -136,7 +140,7 @@ void Foam::setsToFaceZone::applyToSet
if (!fzSet.found(faceI)) if (!fzSet.found(faceI))
{ {
bool flip = false; bool flipFace = false;
label own = mesh_.faceOwner()[faceI]; label own = mesh_.faceOwner()[faceI];
bool ownFound = cSet.found(own); bool ownFound = cSet.found(own);
@ -148,11 +152,11 @@ void Foam::setsToFaceZone::applyToSet
if (ownFound && !neiFound) if (ownFound && !neiFound)
{ {
flip = false; flipFace = false;
} }
else if (!ownFound && neiFound) else if (!ownFound && neiFound)
{ {
flip = true; flipFace = true;
} }
else else
{ {
@ -174,11 +178,17 @@ void Foam::setsToFaceZone::applyToSet
} }
else else
{ {
flip = !ownFound; flipFace = !ownFound;
}
if (flip_)
{
flipFace = !flipFace;
} }
newAddressing.append(faceI); newAddressing.append(faceI);
newFlipMap.append(flip); newFlipMap.append(flipFace);
} }
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -36,6 +36,7 @@ SourceFiles
#define setsToFaceZone_H #define setsToFaceZone_H
#include "topoSetSource.H" #include "topoSetSource.H"
#include "Switch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,10 +57,13 @@ class setsToFaceZone
static addToUsageTable usage_; static addToUsageTable usage_;
//- Name of set to use //- Name of set to use
word faceSetName_; const word faceSetName_;
//- Name of set to use //- Name of set to use
word cellSetName_; const word cellSetName_;
//- Whether cellSet is slave cells or master cells
const Switch flip_;
public: public:
@ -73,7 +77,8 @@ public:
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& faceSetName, const word& faceSetName,
const word& cellSetName const word& cellSetName,
const Switch& flip
); );
//- Construct from dictionary //- Construct from dictionary

View File

@ -66,7 +66,11 @@ void thixotropicViscosity::updateMu()
const volScalarField filmMass("filmMass", film.netMass() + mSMALL); const volScalarField filmMass("filmMass", film.netMass() + mSMALL);
// weighting field to blend new and existing mass contributions // weighting field to blend new and existing mass contributions
const volScalarField w("w", max(0.0, min(1.0, deltaMass/filmMass))); const volScalarField w
(
"w",
max(scalar(0.0), min(scalar(1.0), deltaMass/filmMass))
);
// evaluate thixotropic viscosity // evaluate thixotropic viscosity
volScalarField muThx("muThx", muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL)); volScalarField muThx("muThx", muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL));

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -17,7 +17,7 @@ FoamFile
solvers solvers
{ {
alpha alpha.air
{ {
nAlphaCorr 1; nAlphaCorr 1;
nAlphaSubCycles 2; nAlphaSubCycles 2;

View File

@ -36,7 +36,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

View File

@ -34,7 +34,7 @@ divSchemes
"div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,.*rho.*\)" Gauss linear; "div\(alphaPhi.*,.*rho.*\)" Gauss linear;
"div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;