zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Henry
2012-10-30 11:04:18 +00:00
parent 83c61b0510 389170bc20
commit 2682817789
33 changed files with 971 additions and 547 deletions
Download Patch File Download Diff File Expand all files Collapse all files

78
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
View File

@ -25,15 +25,10 @@ License
#include "jumpCyclicFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF
@ -44,7 +39,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
template<class Type>
jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
(
const jumpCyclicFvPatchField<Type>& ptf,
const fvPatch& p,
@ -57,7 +52,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
template<class Type>
jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
@ -72,7 +67,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
template<class Type>
jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
(
const jumpCyclicFvPatchField<Type>& ptf
)
@ -82,7 +77,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
template<class Type>
jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
(
const jumpCyclicFvPatchField<Type>& ptf,
const DimensionedField<Type, volMesh>& iF
@ -95,7 +90,8 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
tmp<Field<Type> > jumpCyclicFvPatchField<Type>::patchNeighbourField() const
Foam::tmp<Foam::Field<Type> >
Foam::jumpCyclicFvPatchField<Type>::patchNeighbourField() const
{
const Field<Type>& iField = this->internalField();
const labelUList& nbrFaceCells =
@ -133,7 +129,7 @@ tmp<Field<Type> > jumpCyclicFvPatchField<Type>::patchNeighbourField() const
template<class Type>
void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
@ -142,48 +138,22 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
const Pstream::commsTypes
) const
{
scalarField pnf(this->size());
const labelUList& nbrFaceCells =
this->cyclicPatch().neighbFvPatch().faceCells();
// for AMG solve - only apply jump to finest level
if (psiInternal.size() == this->internalField().size())
{
Field<scalar> jf(this->jump()().component(cmpt));
if (!this->cyclicPatch().owner())
{
jf *= -1.0;
}
forAll(*this, facei)
{
pnf[facei] = psiInternal[nbrFaceCells[facei]] - jf[facei];
}
}
else
{
forAll(*this, facei)
{
pnf[facei] = psiInternal[nbrFaceCells[facei]];
}
}
// Transform according to the transformation tensors
this->transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
const labelUList& faceCells = this->cyclicPatch().faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
notImplemented
(
"void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix"
"("
"scalarField&, "
"const scalarField&, "
"const scalarField&, "
"const direction, "
"const Pstream::commsTypes"
") const"
);
}
template<class Type>
void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>& psiInternal,
@ -196,8 +166,8 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
const labelUList& nbrFaceCells =
this->cyclicPatch().neighbFvPatch().faceCells();
// for AMG solve - only apply jump to finest level
if (psiInternal.size() == this->internalField().size())
// only apply jump to original field
if (&psiInternal == &this->internalField())
{
Field<Type> jf(this->jump());
@ -231,8 +201,4 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

11
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H
View File

@ -144,6 +144,17 @@ public:
) const;
};
//- Update result field based on interface functionality
template<>
void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

58
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,17 +27,63 @@ License
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
makePatchFieldsTypeName(jumpCyclic);
}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
makePatchFieldsTypeName(jumpCyclic);
template<>
void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
scalarField pnf(this->size());
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
const labelUList& nbrFaceCells =
this->cyclicPatch().neighbFvPatch().faceCells();
// only apply jump to original field
if (&psiInternal == &this->internalField())
{
Field<scalar> jf(this->jump());
if (!this->cyclicPatch().owner())
{
jf *= -1.0;
}
forAll(*this, facei)
{
pnf[facei] = psiInternal[nbrFaceCells[facei]] - jf[facei];
}
}
else
{
forAll(*this, facei)
{
pnf[facei] = psiInternal[nbrFaceCells[facei]];
}
}
// Transform according to the transformation tensors
this->transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
const labelUList& faceCells = this->cyclicPatch().faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
} // End namespace Foam
// ************************************************************************* //

49
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C
View File

@ -129,33 +129,17 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
const Pstream::commsTypes
) const
{
const labelUList& nbrFaceCells =
this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);
pnf = this->cyclicAMIPatch().interpolate(pnf);
// for AMG solve - only apply jump to finest level
if (psiInternal.size() == this->internalField().size())
{
tmp<Field<scalar> > tjf = jump()().component(cmpt);
if (!this->cyclicAMIPatch().owner())
{
tjf = -tjf;
}
pnf -= tjf;
}
// Transform according to the transformation tensors
this->transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
const labelUList& faceCells = this->cyclicAMIPatch().faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
notImplemented
(
"void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix"
"("
"scalarField&, "
"const scalarField&, "
"const scalarField& coeffs,"
"const direction, "
"const Pstream::commsTypes"
") const"
);
}
@ -175,15 +159,16 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
pnf = this->cyclicAMIPatch().interpolate(pnf);
// for AMG solve - only apply jump to finest level
if (psiInternal.size() == this->internalField().size())
// only apply jump to original field
if (&psiInternal == &this->internalField())
{
tmp<Field<Type> > tjf = jump();
Field<Type> jf(this->jump());
if (!this->cyclicAMIPatch().owner())
{
tjf = -tjf;
jf *= -1.0;
}
pnf -= tjf;
pnf -= jf;
}
// Transform according to the transformation tensors

12
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H
View File

@ -148,6 +148,18 @@ public:
};
//- Update result field based on interface functionality
template<>
void jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

44
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
View File

@ -36,6 +36,50 @@ namespace Foam
makePatchFieldsTypeName(jumpCyclicAMI);
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<>
void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
const labelUList& nbrFaceCells =
this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);
pnf = this->cyclicAMIPatch().interpolate(pnf);
// only apply jump to original field
if (&psiInternal == &this->internalField())
{
Field<scalar> jf(this->jump());
if (!this->cyclicAMIPatch().owner())
{
jf *= -1.0;
}
pnf -= jf;
}
// Transform according to the transformation tensors
this->transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
const labelUList& faceCells = this->cyclicAMIPatch().faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

29
src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
View File

@ -25,6 +25,18 @@ License
#include "fanFvPatchField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
void Foam::fanFvPatchField<Type>::calcFanJump()
{
if (this->cyclicPatch().owner())
{
this->jump_ = this->jumpTable_->value(this->db().time().value());
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
@ -84,4 +96,21 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::fanFvPatchField<Type>::updateCoeffs()
{
if (this->updated())
{
return;
}
calcFanJump();
// call fixedJump variant - uniformJump will overwrite the jump value
fixedJumpFvPatchField<scalar>::updateCoeffs();
}
// ************************************************************************* //

19
src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
View File

@ -39,7 +39,6 @@ Description
\table
Property | Description | Required | Default value
patchType | underlying patch type should be \c cyclic| yes |
jump | current jump value | yes |
jumpTable | jump data, e.g. \c csvFile | yes |
\endtable
@ -49,7 +48,6 @@ Description
{
type fan;
patchType cyclic;
jump uniform 0;
jumpTable csvFile;
csvFileCoeffs
{
@ -101,6 +99,12 @@ class fanFvPatchField
public uniformJumpFvPatchField<Type>
{
// Private Member Functions
//- Calculate the fan pressure jump
void calcFanJump();
public:
//- Runtime type information
@ -170,14 +174,15 @@ public:
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
};
//- Specialisation of the jump-condition for the pressure
template<>
void fanFvPatchField<scalar>::calcFanJump();
template<>
fanFvPatchField<scalar>::fanFvPatchField
(
@ -185,8 +190,6 @@ fanFvPatchField<scalar>::fanFvPatchField
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
template<>
void fanFvPatchField<scalar>::updateCoeffs();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

79
src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
View File

@ -30,18 +30,41 @@ License
#include "Tuple2.H"
#include "polynomial.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
makeTemplatePatchTypeField
(
fvPatchScalarField,
fanFvPatchScalarField
);
}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<>
void Foam::fanFvPatchField<Foam::scalar>::calcFanJump()
{
if (this->cyclicPatch().owner())
{
const surfaceScalarField& phi =
db().lookupObject<surfaceScalarField>("phi");
const fvsPatchField<scalar>& phip =
patch().patchField<surfaceScalarField, scalar>(phi);
scalarField Un(max(phip/patch().magSf(), scalar(0)));
if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
{
Un /= patch().lookupPatchField<volScalarField, scalar>("rho");
}
this->jump_ = this->jumpTable_->value(Un);
}
}
makeTemplatePatchTypeField
(
fvPatchScalarField,
fanFvPatchScalarField
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -93,11 +116,6 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
this->jumpTable_ = DataEntry<scalar>::New("jumpTable", dict);
}
}
else
{
// Dummy jump table
this->jumpTable_.reset(new DataEntry<scalar>("jumpTable"));
}
if (dict.found("value"))
{
@ -113,41 +131,4 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Specialisation of the jump-condition for the pressure
template<>
void Foam::fanFvPatchField<Foam::scalar>::updateCoeffs()
{
if (this->updated())
{
return;
}
if (this->cyclicPatch().owner())
{
const surfaceScalarField& phi =
db().lookupObject<surfaceScalarField>("phi");
const fvsPatchField<scalar>& phip =
patch().patchField<surfaceScalarField, scalar>(phi);
scalarField Un(max(phip/patch().magSf(), scalar(0)));
if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
{
Un /= patch().lookupPatchField<volScalarField, scalar>("rho");
}
this->jump_ = this->jumpTable_->value(Un);
}
uniformJumpFvPatchField<scalar>::updateCoeffs();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

5
src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -138,7 +138,8 @@ template<class Type>
void Foam::fixedJumpFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
os.writeKeyword("patchType") << "cyclic" << token::END_STATEMENT << nl;
os.writeKeyword("patchType") << this->interfaceFieldType()
<< token::END_STATEMENT << nl;
jump_.writeEntry("jump", os);
this->writeEntry("value", os);
}

3
src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
View File

@ -141,7 +141,8 @@ template<class Type>
void Foam::fixedJumpAMIFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
os.writeKeyword("patchType") << "cyclicAMI" << token::END_STATEMENT << nl;
os.writeKeyword("patchType") << this->interfaceFieldType()
<< token::END_STATEMENT << nl;
jump_.writeEntry("jump", os);
this->writeEntry("value", os);
}

27
src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C
View File

@ -35,7 +35,7 @@ Foam::uniformJumpFvPatchField<Type>::uniformJumpFvPatchField
)
:
fixedJumpFvPatchField<Type>(p, iF),
jumpTable_(0)
jumpTable_(new DataEntry<Type>("jumpTable"))
{}
@ -76,6 +76,10 @@ Foam::uniformJumpFvPatchField<Type>::uniformJumpFvPatchField
Field<Type>("value", dict, p.size())
);
}
else
{
this->evaluate(Pstream::blocking);
}
}
@ -105,24 +109,19 @@ Foam::uniformJumpFvPatchField<Type>::uniformJumpFvPatchField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::uniformJumpFvPatchField<Type>::jump() const
void Foam::uniformJumpFvPatchField<Type>::updateCoeffs()
{
if (this->updated())
{
return;
}
if (this->cyclicPatch().owner())
{
const Type value = jumpTable_->value(this->db().time().value());
return tmp<Field<Type> >(new Field<Type>(this->patch().size(), value));
this->jump_ = jumpTable_->value(this->db().time().value());
}
else
{
const uniformJumpFvPatchField& nbrPatch =
refCast<const uniformJumpFvPatchField<Type> >
(
this->neighbourPatchField()
);
return nbrPatch.jump();
}
fixedJumpFvPatchField<Type>::updateCoeffs();
}

7
src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H
View File

@ -165,11 +165,8 @@ public:
// Member functions
// Access
//- Return the "jump" across the patch.
virtual tmp<Field<Type> > jump() const;
//- Update the coefficients
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;

22
src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C
View File

@ -102,25 +102,19 @@ Foam::uniformJumpAMIFvPatchField<Type>::uniformJumpAMIFvPatchField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::uniformJumpAMIFvPatchField<Type>::jump() const
void Foam::uniformJumpAMIFvPatchField<Type>::updateCoeffs()
{
if (this->updated())
{
return;
}
if (this->cyclicAMIPatch().owner())
{
Type j = jumpTable_->value(this->db().time().value());
return tmp<Field<Type> >(new Field<Type>(this->size(), j));
this->jump_ = jumpTable_->value(this->db().time().value());
}
else
{
const uniformJumpAMIFvPatchField& nbrPatch =
refCast<const uniformJumpAMIFvPatchField<Type> >
(
this->neighbourPatchField()
);
return this->cyclicAMIPatch().interpolate(nbrPatch.jump());
}
fixedJumpAMIFvPatchField<Type>::updateCoeffs();
}

7
src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.H
View File

@ -165,11 +165,8 @@ public:
// Member functions
// Access
//- Return the "jump" across the patch.
virtual tmp<Field<Type> > jump() const;
//- Update the coefficients
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;

3
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
View File

@ -97,6 +97,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
refineParams.curvature(),
true, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
@ -207,6 +208,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
refineParams.curvature(),
false, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
true, // surfaceRefinement
true, // curvatureRefinement
@ -368,6 +370,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
refineParams.curvature(),
false, // featureRefinement
true, // featureDistanceRefinement
true, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement

11
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -1406,9 +1406,10 @@ void Foam::autoSnapDriver::doSnap
adaptPatchIDs
)
);
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
const scalarField snapDist(calcSnapDistance(snapParams, pp));
// Construct iterative mesh mover.
@ -1420,7 +1421,7 @@ void Foam::autoSnapDriver::doSnap
motionSmoother meshMover
(
mesh,
ppPtr(),
pp,
adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
motionDict
@ -1475,7 +1476,7 @@ void Foam::autoSnapDriver::doSnap
}
// Check for displacement being outwards.
outwardsDisplacement(ppPtr(), disp);
outwardsDisplacement(pp, disp);
// Set initial distribution of displacement field (on patches)
// from patchDisp and make displacement consistent with b.c.
@ -1489,8 +1490,8 @@ void Foam::autoSnapDriver::doSnap
(
mesh.time().path()
/ "patchDisplacement_" + name(iter) + ".obj",
ppPtr().localPoints(),
ppPtr().localPoints() + disp
pp.localPoints(),
pp.localPoints() + disp
);
}

26
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -141,14 +141,6 @@ class autoSnapDriver
vectorField& pointSurfaceNormal,
vectorField& pointRotation
) const;
//void calcNearestFace
//(
// const label iter,
// const indirectPrimitivePatch& pp,
// vectorField& faceDisp,
// vectorField& faceSurfaceNormal,
// vectorField& faceRotation
//) const;
void calcNearestFace
(
const label iter,
@ -158,15 +150,19 @@ class autoSnapDriver
labelList& faceSurfaceRegion,
vectorField& faceRotation
) const;
void interpolateFaceToPoint
void calcNearestFacePointProperties
(
const label iter,
const indirectPrimitivePatch& pp,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceRotation,
const List<List<point> >& pointFaceCentres,
vectorField& patchDisp
//vectorField& patchRotationDisp
const vectorField& faceDisp,
const vectorField& faceSurfaceNormal,
const labelList& faceSurfaceRegion,
List<List<point> >& pointFaceSurfNormals,
List<List<point> >& pointFaceDisp,
List<List<point> >& pointFaceCentres,
List<labelList>& pointFacePatchID
) const;
void correctAttraction
(
@ -276,6 +272,7 @@ class autoSnapDriver
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch&,
const scalarField&,
@ -339,6 +336,7 @@ class autoSnapDriver
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,

510
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
View File

@ -321,7 +321,7 @@ void Foam::autoSnapDriver::calcNearestFace
const indirectPrimitivePatch& pp,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
labelList& faceSurfaceRegion,
labelList& faceSurfaceGlobalRegion,
vectorField& faceRotation
) const
{
@ -333,8 +333,8 @@ void Foam::autoSnapDriver::calcNearestFace
faceDisp = vector::zero;
faceSurfaceNormal.setSize(pp.size());
faceSurfaceNormal = vector::zero;
faceSurfaceRegion.setSize(pp.size());
faceSurfaceRegion = -1;
faceSurfaceGlobalRegion.setSize(pp.size());
faceSurfaceGlobalRegion = -1;
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces = surfaces.getNamedSurfaces();
@ -419,7 +419,7 @@ void Foam::autoSnapDriver::calcNearestFace
label faceI = ppFaces[hitI];
faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
faceSurfaceNormal[faceI] = hitNormal[hitI];
faceSurfaceRegion[faceI] = surfaces.globalRegion
faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
(
hitSurface[hitI],
hitRegion[hitI]
@ -477,7 +477,11 @@ void Foam::autoSnapDriver::calcNearestFace
label faceI = ppFaces[hitI];
faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
faceSurfaceNormal[faceI] = hitNormal[hitI];
faceSurfaceRegion[faceI] = hitRegion[hitI];
faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
(
hitSurface[hitI],
hitRegion[hitI]
);
}
}
@ -517,6 +521,169 @@ void Foam::autoSnapDriver::calcNearestFace
}
// Collect (possibly remote) per point data of all surrounding faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal
// - faceDisp
// - faceCentres&faceNormal
void Foam::autoSnapDriver::calcNearestFacePointProperties
(
const label iter,
const indirectPrimitivePatch& pp,
const vectorField& faceDisp,
const vectorField& faceSurfaceNormal,
const labelList& faceSurfaceGlobalRegion,
List<List<point> >& pointFaceSurfNormals,
List<List<point> >& pointFaceDisp,
List<List<point> >& pointFaceCentres,
List<labelList>& pointFacePatchID
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
// For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only
// (see e.g PatchToolsNormals.C)
pointFaceSurfNormals.setSize(pp.nPoints());
pointFaceDisp.setSize(pp.nPoints());
pointFaceCentres.setSize(pp.nPoints());
pointFacePatchID.setSize(pp.nPoints());
// Fill local data
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());
forAll(pFaces, i)
{
label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI];
pDisp[i] = faceDisp[faceI];
pFc[i] = pp.faceCentres()[faceI];
//label meshFaceI = pp.addressing()[faceI];
//pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
pFid[i] = globalToPatch_[faceSurfaceGlobalRegion[faceI]];
}
}
// Collect additionally 'normal' boundary faces for boundaryPoints of pp
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// points on the boundary of pp should pick up non-pp normals
// as well for the feature-reconstruction to behave correctly.
// (the movement is already constrained outside correctly so it
// is only that the unconstrained attraction vector is calculated
// correctly)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList patchID(pbm.patchID());
// Unmark all non-coupled boundary faces
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label meshFaceI = pp.start()+i;
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
}
}
// Remove any meshed faces
forAll(pp.addressing(), i)
{
label meshFaceI = pp.addressing()[i];
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
// See if pp point uses any non-meshed boundary faces
const labelList& boundaryPoints = pp.boundaryPoints();
forAll(boundaryPoints, i)
{
label pointI = boundaryPoints[i];
label meshPointI = pp.meshPoints()[pointI];
const point& pt = mesh.points()[meshPointI];
const labelList& pFaces = mesh.pointFaces()[meshPointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
forAll(pFaces, i)
{
label meshFaceI = pFaces[i];
if (!mesh.isInternalFace(meshFaceI))
{
label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
if (patchI != -1)
{
vector fn = mesh.faceAreas()[meshFaceI];
pNormals.append(fn/mag(fn));
pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
pFc.append(mesh.faceCentres()[meshFaceI]);
pFid.append(patchI);
}
}
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceSurfNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceDisp,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFacePatchID,
listPlusEqOp<label>(),
List<label>()
);
}
// Gets passed in offset to nearest point on feature edge. Calculates
// if the point has a different number of faces on either side of the feature
// and if so attracts the point to that non-dominant plane.
@ -580,56 +747,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findMultiPatchPoint
}
return pointIndexHit(false, vector::zero, labelMax);
}
////XXXXXXXX
//void Foam::autoSnapDriver::attractMultiPatchPoint
//(
// const label iter,
// const scalar featureCos,
//
// const indirectPrimitivePatch& pp,
// const scalarField& snapDist,
// const label pointI,
//
// const List<List<point> >& pointFaceSurfNormals,
// const labelListList& pointFaceSurfaceRegion,
// const List<List<point> >& pointFaceDisp,
// const List<List<point> >& pointFaceCentres,
// const labelListList& pointFacePatchID,
//
// vector& patchAttraction,
// pointConstraint& patchConstraint
//) const
//{
// // Collect
//
// );
//
// if
// (
// (constraint.first() > patchConstraints[pointI].first())
// || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
// )
// {
// patchAttraction[pointI] = attraction;
// patchConstraints[pointI] = constraint;
//
// // Check the number of directions
// if (patchConstraints[pointI].first() == 1)
// {
// // Flat surface. Check for different patchIDs
// pointIndexHit multiPatchPt
// (
// findMultiPatchPoint
// (
// pt,
// pointFacePatchID[pointI],
// pointFaceCentres[pointI]
// )
// );
// if (multiPatchPt.hit())
// {
// // Behave like when having two surface normals so
////XXXXXXXX
void Foam::autoSnapDriver::binFeatureFace
(
@ -1361,6 +1479,7 @@ void Foam::autoSnapDriver::determineFeatures
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
@ -1464,69 +1583,72 @@ void Foam::autoSnapDriver::determineFeatures
if (patchConstraints[pointI].first() == 1)
{
// Flat surface. Check for different patchIDs
pointIndexHit multiPatchPt
(
findMultiPatchPoint
(
pt,
pointFacePatchID[pointI],
pointFaceCentres[pointI]
)
);
if (multiPatchPt.hit())
if (multiRegionFeatureSnap)
{
// Behave like when having two surface normals so
// attract to nearest feature edge (with a guess for
// the multipatch point as starting point)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
pointIndexHit multiPatchPt
(
pp,
snapDist,
pointI,
multiPatchPt.hitPoint(), //estimatedPt
featI,
edgeAttractors,
edgeConstraints,
patchAttraction,
patchConstraints
findMultiPatchPoint
(
pt,
pointFacePatchID[pointI],
pointFaceCentres[pointI]
)
);
if (multiPatchPt.hit())
{
// Behave like when having two surface normals so
// attract to nearest feature edge (with a guess for
// the multipatch point as starting point)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
(
pp,
snapDist,
pointI,
multiPatchPt.hitPoint(), //estimatedPt
if (nearInfo.hit())
{
// Dump
if (featureEdgeStr.valid())
featI,
edgeAttractors,
edgeConstraints,
patchAttraction,
patchConstraints
);
if (nearInfo.hit())
{
meshTools::writeOBJ(featureEdgeStr(), pt);
featureEdgeVertI++;
meshTools::writeOBJ
(
featureEdgeStr(),
nearInfo.hitPoint()
);
featureEdgeVertI++;
featureEdgeStr()
<< "l " << featureEdgeVertI-1 << ' '
<< featureEdgeVertI << nl;
// Dump
if (featureEdgeStr.valid())
{
meshTools::writeOBJ(featureEdgeStr(), pt);
featureEdgeVertI++;
meshTools::writeOBJ
(
featureEdgeStr(),
nearInfo.hitPoint()
);
featureEdgeVertI++;
featureEdgeStr()
<< "l " << featureEdgeVertI-1 << ' '
<< featureEdgeVertI << nl;
}
}
}
else
{
if (missedEdgeStr.valid())
else
{
meshTools::writeOBJ(missedEdgeStr(), pt);
missedVertI++;
meshTools::writeOBJ
(
missedEdgeStr(),
nearInfo.missPoint()
);
missedVertI++;
missedEdgeStr()
<< "l " << missedVertI-1 << ' '
<< missedVertI << nl;
if (missedEdgeStr.valid())
{
meshTools::writeOBJ(missedEdgeStr(), pt);
missedVertI++;
meshTools::writeOBJ
(
missedEdgeStr(),
nearInfo.missPoint()
);
missedVertI++;
missedEdgeStr()
<< "l " << missedVertI-1 << ' '
<< missedVertI << nl;
}
}
}
}
@ -1645,6 +1767,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
@ -1697,6 +1820,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
(
iter,
featureCos,
multiRegionFeatureSnap,
pp,
snapDist,
@ -2131,6 +2255,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
//MEJ: any faces that have multi-patch points only keep the multi-patch
// points. The other points on the face will be dragged along
// (hopefully)
if (multiRegionFeatureSnap)
{
autoPtr<OFstream> multiPatchStr;
if (debug&meshRefinement::OBJINTERSECTIONS)
@ -2536,10 +2661,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
{
const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
Info<< "Overriding displacement on features :" << nl
<< " implicit features : " << implicitFeatureAttraction << nl
<< " explicit features : " << explicitFeatureAttraction << nl
<< " implicit features : " << implicitFeatureAttraction << nl
<< " explicit features : " << explicitFeatureAttraction << nl
<< " multi-patch features : " << multiRegionFeatureSnap << nl
<< endl;
@ -2547,6 +2674,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const pointField& localPoints = pp.localPoints();
const fvMesh& mesh = meshRefiner_.mesh();
// Displacement and orientation per pp face
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2554,7 +2682,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
vectorField faceDisp(pp.size(), vector::zero);
// normal of surface at point on surface
vectorField faceSurfaceNormal(pp.size(), vector::zero);
labelList faceSurfaceRegion(pp.size(), -1);
labelList faceSurfaceGlobalRegion(pp.size(), -1);
vectorField faceRotation(pp.size(), vector::zero);
calcNearestFace
@ -2563,7 +2691,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
pp,
faceDisp,
faceSurfaceNormal,
faceSurfaceRegion,
faceSurfaceGlobalRegion,
faceRotation
);
@ -2587,158 +2715,25 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal
// - faceDisp
// - faceRotation
// - faceCentres&faceNormal
// For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only
// (see e.g PatchToolsNormals.C)
List<List<point> > pointFaceSurfNormals(pp.nPoints());
List<List<point> > pointFaceDisp(pp.nPoints());
//List<List<point> > pointFaceRotation(pp.nPoints());
List<List<point> > pointFaceCentres(pp.nPoints());
List<labelList> pointFacePatchID(pp.nPoints());
// Fill local data
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>& pRot = pointFaceRotation[pointI];
//pRot.setSize(pFaces.size());
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(pFaces.size());
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(pFaces.size());
forAll(pFaces, i)
{
label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI];
pDisp[i] = faceDisp[faceI];
//pRot[i] = faceRotation[faceI];
pFc[i] = pp.faceCentres()[faceI];
label meshFaceI = pp.addressing()[faceI];
pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
}
}
// Collect additionally 'normal' boundary faces for boundaryPoints of pp
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// points on the boundary of pp should pick up non-pp normals
// as well for the feature-reconstruction to behave correctly.
// (the movement is already constrained outside correctly so it
// is only that the unconstrained attraction vector is calculated
// correctly)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList patchID(pbm.patchID());
// Unmark all non-coupled boundary faces
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label meshFaceI = pp.start()+i;
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
}
}
// Remove any meshed faces
forAll(pp.addressing(), i)
{
label meshFaceI = pp.addressing()[i];
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
// See if pp point uses any non-meshed boundary faces
const labelList& boundaryPoints = pp.boundaryPoints();
forAll(boundaryPoints, i)
{
label pointI = boundaryPoints[i];
label meshPointI = pp.meshPoints()[pointI];
const point& pt = mesh.points()[meshPointI];
const labelList& pFaces = mesh.pointFaces()[meshPointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
forAll(pFaces, i)
{
label meshFaceI = pFaces[i];
if (!mesh.isInternalFace(meshFaceI))
{
label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
if (patchI != -1)
{
vector fn = mesh.faceAreas()[meshFaceI];
pNormals.append(fn/mag(fn));
pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
pFc.append(mesh.faceCentres()[meshFaceI]);
pFid.append(patchI);
}
}
}
}
}
syncTools::syncPointList
calcNearestFacePointProperties
(
mesh,
pp.meshPoints(),
iter,
pp,
faceDisp,
faceSurfaceNormal,
faceSurfaceGlobalRegion,
pointFaceSurfNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceDisp,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
//syncTools::syncPointList
//(
// mesh,
// pp.meshPoints(),
// pointFaceRotation,
// listPlusEqOp<point>(),
// List<point>(),
// listTransform()
//);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFacePatchID,
listPlusEqOp<label>(),
List<label>()
pointFacePatchID
);
@ -2793,6 +2788,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
(
iter,
featureCos,
multiRegionFeatureSnap,
pp,
snapDist,

6
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
View File

@ -36,7 +36,11 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
nSnap_(readLabel(dict.lookup("nRelaxIter"))),
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
multiRegionFeatureSnap_
(
dict.lookupOrDefault("multiRegionFeatureSnap", false)
)
{}

7
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
View File

@ -68,6 +68,8 @@ class snapParameters
const Switch implicitFeatureSnap_;
const Switch multiRegionFeatureSnap_;
// Private Member Functions
@ -134,6 +136,11 @@ public:
return implicitFeatureSnap_;
}
Switch multiRegionFeatureSnap() const
{
return multiRegionFeatureSnap_;
}
};

9
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -264,6 +264,14 @@ private:
label& nRefine
) const;
//- Mark cells for distance-to-feature based refinement.
label markInternalDistanceToFeatureRefinement
(
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const;
//- Mark cells for refinement-shells based refinement.
label markInternalRefinement
(
@ -656,6 +664,7 @@ public:
const scalar curvature,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,

119
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -290,7 +290,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
forAll(features_, featI)
{
const featureEdgeMesh& featureMesh = features_[featI];
const label featureLevel = features_.levels()[featI];
const label featureLevel = features_.levels()[featI][0];
const labelListList& pointEdges = featureMesh.pointEdges();
// Find regions on edgeMesh
@ -511,6 +511,90 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
}
// Mark cells for distance-to-feature based refinement.
Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
(
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
const pointField& cellCentres = mesh_.cellCentres();
// Detect if there are any distance shells
if (features_.maxDistance() <= 0.0)
{
return 0;
}
label oldNRefine = nRefine;
// Collect cells to test
pointField testCc(cellLevel.size()-nRefine);
labelList testLevels(cellLevel.size()-nRefine);
label testI = 0;
forAll(cellLevel, cellI)
{
if (refineCell[cellI] == -1)
{
testCc[testI] = cellCentres[cellI];
testLevels[testI] = cellLevel[cellI];
testI++;
}
}
// Do test to see whether cells is inside/outside shell with higher level
labelList maxLevel;
features_.findHigherLevel(testCc, testLevels, maxLevel);
// Mark for refinement. Note that we didn't store the original cellID so
// now just reloop in same order.
testI = 0;
forAll(cellLevel, cellI)
{
if (refineCell[cellI] == -1)
{
if (maxLevel[testI] > testLevels[testI])
{
bool reachedLimit = !markForRefine
(
maxLevel[testI], // mark with any positive value
nAllowRefine,
refineCell[cellI],
nRefine
);
if (reachedLimit)
{
if (debug)
{
Pout<< "Stopped refining internal cells"
<< " since reaching my cell limit of "
<< mesh_.nCells()+7*nRefine << endl;
}
break;
}
}
testI++;
}
}
if
(
returnReduce(nRefine, sumOp<label>())
> returnReduce(nAllowRefine, sumOp<label>())
)
{
Info<< "Reached refinement limit." << endl;
}
return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
// Mark cells for non-surface intersection based refinement.
Foam::label Foam::meshRefinement::markInternalRefinement
(
@ -1128,6 +1212,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
const scalar curvature,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
@ -1196,8 +1281,24 @@ Foam::labelList Foam::meshRefinement::refineCandidates
nRefine
);
Info<< "Marked for refinement due to explicit features : "
<< nFeatures << " cells." << endl;
Info<< "Marked for refinement due to explicit features "
<< ": " << nFeatures << " cells." << endl;
}
// Inside distance-to-feature shells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (featureDistanceRefinement)
{
label nShell = markInternalDistanceToFeatureRefinement
(
nAllowRefine,
refineCell,
nRefine
);
Info<< "Marked for refinement due to distance to explicit features "
": " << nShell << " cells." << endl;
}
// Inside refinement shells
@ -1212,8 +1313,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to refinement shells : "
<< nShell << " cells." << endl;
Info<< "Marked for refinement due to refinement shells "
<< ": " << nShell << " cells." << endl;
}
// Refinement based on intersection of surface
@ -1230,8 +1331,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to surface intersection : "
<< nSurf << " cells." << endl;
Info<< "Marked for refinement due to surface intersection "
<< ": " << nSurf << " cells." << endl;
}
// Refinement based on curvature of surface
@ -1254,8 +1355,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to curvature/regions : "
<< nCurv << " cells." << endl;
Info<< "Marked for refinement due to curvature/regions "
<< ": " << nCurv << " cells." << endl;
}
// Pack cells-to-refine

180
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
View File

@ -25,6 +25,7 @@ License
#include "refinementFeatures.H"
#include "Time.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -34,15 +35,15 @@ void Foam::refinementFeatures::read
const PtrList<dictionary>& featDicts
)
{
forAll(featDicts, i)
forAll(featDicts, featI)
{
const dictionary& dict = featDicts[i];
const dictionary& dict = featDicts[featI];
fileName featFileName(dict.lookup("file"));
set
(
i,
featI,
new featureEdgeMesh
(
IOobject
@ -58,15 +59,74 @@ void Foam::refinementFeatures::read
)
);
const featureEdgeMesh& eMesh = operator[](i);
const featureEdgeMesh& eMesh = operator[](featI);
//eMesh.mergePoints(meshRefiner_.mergeDistance());
levels_[i] = readLabel(dict.lookup("level"));
Info<< "Refinement level " << levels_[i]
<< " for all cells crossed by feature " << featFileName
<< " (" << eMesh.points().size() << " points, "
if (dict.found("levels"))
{
List<Tuple2<scalar, label> > distLevels(dict["levels"]);
if (dict.size() < 1)
{
FatalErrorIn
(
"refinementFeatures::read"
"(const objectRegistry&"
", const PtrList<dictionary>&)"
) << " : levels should be at least size 1" << endl
<< "levels : " << dict["levels"]
<< exit(FatalError);
}
distances_[featI].setSize(distLevels.size());
levels_[featI].setSize(distLevels.size());
forAll(distLevels, j)
{
distances_[featI][j] = distLevels[j].first();
levels_[featI][j] = distLevels[j].second();
// Check in incremental order
if (j > 0)
{
if
(
(distances_[featI][j] <= distances_[featI][j-1])
|| (levels_[featI][j] > levels_[featI][j-1])
)
{
FatalErrorIn
(
"refinementFeatures::read"
"(const objectRegistry&"
", const PtrList<dictionary>&)"
) << " : Refinement should be specified in order"
<< " of increasing distance"
<< " (and decreasing refinement level)." << endl
<< "Distance:" << distances_[featI][j]
<< " refinementLevel:" << levels_[featI][j]
<< exit(FatalError);
}
}
}
}
else
{
// Look up 'level' for single level
levels_[featI] = labelList(1, readLabel(dict.lookup("level")));
distances_[featI] = scalarField(1, 0.0);
}
Info<< "Refinement level according to distance to "
<< featFileName << " (" << eMesh.points().size() << " points, "
<< eMesh.edges().size() << " edges)." << endl;
forAll(levels_[featI], j)
{
Info<< " level " << levels_[featI][j]
<< " for all cells within " << distances_[featI][j]
<< " meter." << endl;
}
}
}
@ -127,6 +187,80 @@ void Foam::refinementFeatures::buildTrees
}
// Find maximum level of a shell.
void Foam::refinementFeatures::findHigherLevel
(
const pointField& pt,
const label featI,
labelList& maxLevel
) const
{
const labelList& levels = levels_[featI];
const scalarField& distances = distances_[featI];
// Collect all those points that have a current maxLevel less than
// (any of) the shell. Also collect the furthest distance allowable
// to any shell with a higher level.
pointField candidates(pt.size());
labelList candidateMap(pt.size());
scalarField candidateDistSqr(pt.size());
label candidateI = 0;
forAll(maxLevel, pointI)
{
forAllReverse(levels, levelI)
{
if (levels[levelI] > maxLevel[pointI])
{
candidates[candidateI] = pt[pointI];
candidateMap[candidateI] = pointI;
candidateDistSqr[candidateI] = sqr(distances[levelI]);
candidateI++;
break;
}
}
}
candidates.setSize(candidateI);
candidateMap.setSize(candidateI);
candidateDistSqr.setSize(candidateI);
// Do the expensive nearest test only for the candidate points.
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
List<pointIndexHit> nearInfo(candidates.size());
forAll(candidates, candidateI)
{
nearInfo[candidateI] = tree.findNearest
(
candidates[candidateI],
candidateDistSqr[candidateI]
);
}
// Update maxLevel
forAll(nearInfo, candidateI)
{
if (nearInfo[candidateI].hit())
{
// Check which level it actually is in.
label minDistI = findLower
(
distances,
mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
);
label pointI = candidateMap[candidateI];
// pt is inbetween shell[minDistI] and shell[minDistI+1]
maxLevel[pointI] = levels[minDistI+1];
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::refinementFeatures::refinementFeatures
@ -136,6 +270,7 @@ Foam::refinementFeatures::refinementFeatures
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
@ -175,6 +310,7 @@ Foam::refinementFeatures::refinementFeatures
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
@ -336,4 +472,32 @@ void Foam::refinementFeatures::findNearestPoint
}
void Foam::refinementFeatures::findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const
{
// Maximum level of any shell. Start off with level of point.
maxLevel = ptLevel;
forAll(*this, featI)
{
findHigherLevel(pt, featI, maxLevel);
}
}
Foam::scalar Foam::refinementFeatures::maxDistance() const
{
scalar overallMax = -GREAT;
forAll(distances_, featI)
{
overallMax = max(overallMax, max(distances_[featI]));
}
return overallMax;
}
// ************************************************************************* //

34
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
View File

@ -57,8 +57,11 @@ private:
// Private data
//- Refinement levels
labelList levels_;
//- Per shell the list of ranges
List<scalarField> distances_;
//- Per shell per distance the refinement level
labelListList levels_;
//- Edge
PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
@ -75,6 +78,13 @@ private:
//- Build edge tree and feature point tree
void buildTrees(const label, const labelList&);
//- Find shell level higher than ptLevel
void findHigherLevel
(
const pointField& pt,
const label featI,
labelList& maxLevel
) const;
public:
@ -101,11 +111,18 @@ public:
// Access
const labelList& levels() const
//- Per featureEdgeMesh the list of level
const labelListList& levels() const
{
return levels_;
}
//- Per featureEdgeMesh the list of ranges
const List<scalarField>& distances() const
{
return distances_;
}
const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
{
return edgeTrees_;
@ -119,6 +136,9 @@ public:
// Query
//- Highest distance of all features
scalar maxDistance() const;
//- Find nearest point on nearest feature edge
void findNearestEdge
(
@ -141,6 +161,14 @@ public:
labelList& nearIndex
) const;
//- Find shell level higher than ptLevel
void findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const;
};

95
src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
View File

@ -45,31 +45,56 @@ void Foam::fieldMinMax::calcMinMaxFields
const label procI = Pstream::myProcNo();
const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
const fvMesh& mesh = field.mesh();
const volVectorField::GeometricBoundaryField& CfBoundary =
mesh.C().boundaryField();
switch (mode)
{
case mdMag:
{
const scalarField magField(mag(field));
const volScalarField magField(mag(field));
const volScalarField::GeometricBoundaryField& magFieldBoundary =
magField.boundaryField();
labelList minIs(Pstream::nProcs());
scalarList minVs(Pstream::nProcs());
List<vector> minCs(Pstream::nProcs());
minIs[procI] = findMin(magField);
minVs[procI] = magField[minIs[procI]];
minCs[procI] = field.mesh().C()[minIs[procI]];
label minProcI = findMin(magField);
minVs[procI] = magField[minProcI];
minCs[procI] = field.mesh().C()[minProcI];
Pstream::gatherList(minIs);
Pstream::gatherList(minVs);
Pstream::gatherList(minCs);
labelList maxIs(Pstream::nProcs());
scalarList maxVs(Pstream::nProcs());
List<vector> maxCs(Pstream::nProcs());
maxIs[procI] = findMax(magField);
maxVs[procI] = magField[maxIs[procI]];
maxCs[procI] = field.mesh().C()[maxIs[procI]];
label maxProcI = findMax(magField);
maxVs[procI] = magField[maxProcI];
maxCs[procI] = field.mesh().C()[maxProcI];
forAll(magFieldBoundary, patchI)
{
const scalarField& mfp = magFieldBoundary[patchI];
const vectorField& Cfp = CfBoundary[patchI];
label minPI = findMin(mfp);
if (mfp[minPI] < minVs[procI])
{
minVs[procI] = mfp[minPI];
minCs[procI] = Cfp[minPI];
}
label maxPI = findMax(mfp);
if (mfp[maxPI] > maxVs[procI])
{
maxVs[procI] = mfp[maxPI];
maxCs[procI] = Cfp[maxPI];
}
}
Pstream::gatherList(minVs);
Pstream::gatherList(minCs);
Pstream::gatherList(maxIs);
Pstream::gatherList(maxVs);
Pstream::gatherList(maxCs);
@ -127,25 +152,47 @@ void Foam::fieldMinMax::calcMinMaxFields
}
case mdCmpt:
{
List<Type> minVs(Pstream::nProcs());
labelList minIs(Pstream::nProcs());
List<vector> minCs(Pstream::nProcs());
minIs[procI] = findMin(field);
minVs[procI] = field[minIs[procI]];
minCs[procI] = field.mesh().C()[minIs[procI]];
const typename fieldType::GeometricBoundaryField&
fieldBoundary = field.boundaryField();
List<Type> minVs(Pstream::nProcs());
List<vector> minCs(Pstream::nProcs());
label minProcI = findMin(field);
minVs[procI] = field[minProcI];
minCs[procI] = field.mesh().C()[minProcI];
Pstream::gatherList(minIs);
Pstream::gatherList(minVs);
Pstream::gatherList(minCs);
List<Type> maxVs(Pstream::nProcs());
labelList maxIs(Pstream::nProcs());
List<vector> maxCs(Pstream::nProcs());
maxIs[procI] = findMax(field);
maxVs[procI] = field[maxIs[procI]];
maxCs[procI] = field.mesh().C()[maxIs[procI]];
label maxProcI = findMax(field);
maxVs[procI] = field[maxProcI];
maxCs[procI] = field.mesh().C()[maxProcI];
forAll(fieldBoundary, patchI)
{
const Field<Type>& fp = fieldBoundary[patchI];
const vectorField& Cfp = CfBoundary[patchI];
label minPI = findMin(fp);
if (fp[minPI] < minVs[procI])
{
minVs[procI] = fp[minPI];
minCs[procI] = Cfp[minPI];
}
label maxPI = findMax(fp);
if (fp[maxPI] > maxVs[procI])
{
maxVs[procI] = fp[maxPI];
maxCs[procI] = Cfp[maxPI];
}
}
Pstream::gatherList(minVs);
Pstream::gatherList(minCs);
Pstream::gatherList(maxIs);
Pstream::gatherList(maxVs);
Pstream::gatherList(maxCs);

8
src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
View File

@ -113,9 +113,15 @@ void Foam::energyJumpFvPatchScalarField::updateCoeffs()
thermo.T().boundaryField()[patchID]
);
fixedJumpFvPatchScalarField& Tbp =
const_cast<fixedJumpFvPatchScalarField&>(TbPatch);
// force update of jump
Tbp.updateCoeffs();
const labelUList& faceCells = this->patch().faceCells();
jump_ = thermo.he(pp, TbPatch.jump(), faceCells);
jump_ = thermo.he(pp, Tbp.jump(), faceCells);
}
fixedJumpFvPatchField<scalar>::updateCoeffs();

8
src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
View File

@ -113,9 +113,15 @@ void Foam::energyJumpAMIFvPatchScalarField::updateCoeffs()
thermo.T().boundaryField()[patchID]
);
fixedJumpAMIFvPatchScalarField& Tbp =
const_cast<fixedJumpAMIFvPatchScalarField&>(TbPatch);
// force update of jump
Tbp.updateCoeffs();
const labelUList& faceCells = this->patch().faceCells();
jump_ = thermo.he(pp, TbPatch.jump(), faceCells);
jump_ = thermo.he(pp, Tbp.jump(), faceCells);
}
fixedJumpAMIFvPatchField<scalar>::updateCoeffs();