diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index f0cc40a8e4867cc7d1de6408832d6360006e9985..4bef1736f01a061bb88177b6e3cbb0620e069171 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -96,20 +96,29 @@ castellatedMeshControls // refinement. nCellsBetweenLevels 1; + // Explicit feature edge refinement // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Specifies a level for any cell intersected by explicitly provided // edges. // This is a featureEdgeMesh, read from constant/triSurface for now. + // Specify 'levels' in the same way as the 'distance' mode in the + // refinementRegions (see below). The old specification + // level 2; + // is equivalent to + // levels ((0 2)); + features ( //{ // file "someLine.eMesh"; - // level 2; + // //level 2; + // levels ((0.0 2) (1.0 3)); //} ); + // Surface based refinement // ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -178,7 +187,7 @@ castellatedMeshControls // three modes // - distance. 'levels' specifies per distance to the surface the // wanted refinement level. The distances need to be specified in - // descending order. + // increasing order. // - inside. 'levels' is only one entry and only the level is used. All // cells inside the surface get refined up to the level. The surface // needs to be closed for this to be possible. diff --git a/bin/foamNew b/bin/foamNew index 4c5a996aad01147a3e75afb743ecc31bc4003aba..5d61de63331ed59260c4bec23178278901684d4a 100755 --- a/bin/foamNew +++ b/bin/foamNew @@ -52,11 +52,11 @@ USAGE case "$1" in (-s | -source | source) shift - $WM_PROJECT_DIR/etc/codeTemplates/source/foamNewSource $* + $WM_PROJECT_DIR/etc/codeTemplates/source/foamNewSource "$@" ;; (-t | -template | template) shift - $WM_PROJECT_DIR/etc/codeTemplates/template/foamNewTemplate $* + $WM_PROJECT_DIR/etc/codeTemplates/template/foamNewTemplate "$@" ;; (*) usage "unknown type '$1'" diff --git a/etc/config/settings.csh b/etc/config/settings.csh index 069624f513770e59da28f60cf8aadb975c5d82b8..059ea30e336bca5078baed3c2e343c97d678adfa 100644 --- a/etc/config/settings.csh +++ b/etc/config/settings.csh @@ -391,22 +391,12 @@ case SYSTEMOPENMPI: # Use the system installed openmpi, get library directory via mpicc setenv FOAM_MPI openmpi-system - # Set compilation flags here instead of in wmake/rules/../mplibSYSTEMOPENMPI - setenv PINC "`mpicc --showme:compile`" - setenv PLIBS "`mpicc --showme:link`" - set libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'` + set libDir=`mpicc --showme:link | sed -e 's/.*-L\([^ ]*\).*/\1/'` # Bit of a hack: strip off 'lib' and hope this is the path to openmpi # include files and libraries. setenv MPI_ARCH_PATH "${libDir:h}" - if ($?FOAM_VERBOSE && $?prompt) then - echo "Using system installed MPI:" - echo " compile flags : $PINC" - echo " link flags : $PLIBS" - echo " libmpi dir : $libDir" - endif - _foamAddLib $libDir unset libDir breaksw diff --git a/etc/config/settings.sh b/etc/config/settings.sh index 1fc2136530a5374a4d9cc3d9f8955c8858dbb91f..ed45429a650de2aa87ad202bf0634d8af40f32aa 100644 --- a/etc/config/settings.sh +++ b/etc/config/settings.sh @@ -418,23 +418,12 @@ SYSTEMOPENMPI) # Use the system installed openmpi, get library directory via mpicc export FOAM_MPI=openmpi-system - # Set compilation flags here instead of in wmake/rules/../mplibSYSTEMOPENMPI - export PINC="`mpicc --showme:compile`" - export PLIBS="`mpicc --showme:link`" - libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'` + libDir=`mpicc --showme:link | sed -e 's/.*-L\([^ ]*\).*/\1/'` # Bit of a hack: strip off 'lib' and hope this is the path to openmpi # include files and libraries. export MPI_ARCH_PATH="${libDir%/*}" - if [ "$FOAM_VERBOSE" -a "$PS1" ] - then - echo "Using system installed MPI:" 1>&2 - echo " compile flags : $PINC" 1>&2 - echo " link flags : $PLIBS" 1>&2 - echo " libmpi dir : $libDir" 1>&2 - fi - _foamAddLib $libDir unset libDir ;; diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C index db6849d26746be1dbdd6a5ffe56feed4e8dde12b..6313dab9c187afe3b82500ff29651415167cd0c0 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C @@ -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 - // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H index a5fd2343bfbf2fdfe13cc55cfbe5e12f4686abe6..42d0c1805f86a345e061b45faedcaffeed665c60 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H @@ -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; + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C index 02e25e31da795a6cc4aa4c2d06676ea708f48511..2001e1e2f4bdf3b04dca197402ab49d0141255c7 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C @@ -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 * * * * * * * * * * * * * // + +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]]; + } + } -makePatchFieldsTypeName(jumpCyclic); + // 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 // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C index 87d7c20efd5ba6f36d52eb22b8596efd211eb1bc..ff36ee9797cfb8b61b259f85048881ce6bf62d64 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C @@ -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 diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H index 6759b09a36c57ba5a49b75d48147e216f010126e..a2e77f0396e8e31008b7a9e172a16254b144df22 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H @@ -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 diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C index 323f460d9f7d99139c13f5d878c6b6e2094a9c05..b7349812dabde90d340cc9a55a29e40bbcc7066b 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C @@ -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 diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C index 3a150c6f70792102426a0c851279fe19c0838c66..775c0f1ced57a7af7e54e0df3454b75afe5fe913 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C @@ -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(); +} + + // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H index 463ebeb543602d0a1dd0e4291b626b6ef1101307..218b4b18626d114127da29094a9f5b0ae4458180 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H @@ -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(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C index 668d43483de22fd7e282f9d999bea437d36ee6a3..8535c571279124b0016ac8bf53901ae4c4ebea74 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C @@ -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 - // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C index e58270ec638df26b663319261278f6661d3f9b12..9ef4be21a407f055c7f5831627e1e14a4120e49f 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C @@ -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); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C index d9cdc383a3f56c64518e87484ee8a442e27eae72..4c3290cf2d7e1d5a8b4c5006cd9c69601bab82fe 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C @@ -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); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C index c4ee64eaa842271da8e38d9c61b2e13121019ee9..250442dad05d416d0ac011b9acb7dfba9ba132ac 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C @@ -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->cyclicPatch().owner()) + if (this->updated()) { - const Type value = jumpTable_->value(this->db().time().value()); - - return tmp<Field<Type> >(new Field<Type>(this->patch().size(), value)); + return; } - else - { - const uniformJumpFvPatchField& nbrPatch = - refCast<const uniformJumpFvPatchField<Type> > - ( - this->neighbourPatchField() - ); - return nbrPatch.jump(); + if (this->cyclicPatch().owner()) + { + this->jump_ = jumpTable_->value(this->db().time().value()); } + + fixedJumpFvPatchField<Type>::updateCoeffs(); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H index 633dc1752e932bc92a6f646b8b73be79cf268623..aba63cb9d0eb70f6a70f3fa7edd13b9468b0384f 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H @@ -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; diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C index a3e1fe659190cf604f8749567e9be9e22a21d501..06f4e3c13457a484a4d525cdb5e645818099d090 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C @@ -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->cyclicAMIPatch().owner()) + if (this->updated()) { - Type j = jumpTable_->value(this->db().time().value()); - - return tmp<Field<Type> >(new Field<Type>(this->size(), j)); + return; } - else - { - const uniformJumpAMIFvPatchField& nbrPatch = - refCast<const uniformJumpAMIFvPatchField<Type> > - ( - this->neighbourPatchField() - ); - return this->cyclicAMIPatch().interpolate(nbrPatch.jump()); + if (this->cyclicAMIPatch().owner()) + { + this->jump_ = jumpTable_->value(this->db().time().value()); } + + fixedJumpAMIFvPatchField<Type>::updateCoeffs(); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.H index e5eb2e6826883fb37c33b72c3af8d2118525af10..e00f9f9991dd43b6598a54e008240d031e16930b 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.H @@ -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; diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C index 985bcfd103b4289ac77ddff1410c72f66f24eb43..dfd769fd15002f3474084f4c40bcf11733593c88 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C @@ -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 diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C index e358ba4ffa4d04a5ea96f1571d9286feb790ab31..54a78f477c6a5e7170572da5c8867654c564c543 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C @@ -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 ); } diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H index bd5dd86188fa9e64b467d600257b569679e63a77..7cf1e79f6121a963dc5ccdf108ac05e0bd7f75d9 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H @@ -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, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C index 8bcd9c6270a22ae42ce171a63d8f50000919204e..f32d9ce4bd3a775e76bac478a361c1a96d6aab16 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C @@ -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 + 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 - featI, - edgeAttractors, - edgeConstraints, + featI, + edgeAttractors, + edgeConstraints, - patchAttraction, - patchConstraints - ); + patchAttraction, + patchConstraints + ); - if (nearInfo.hit()) - { - // Dump - if (featureEdgeStr.valid()) + 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()]; + calcNearestFacePointProperties + ( + iter, + pp, - 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); - } - } - } - } - } + faceDisp, + faceSurfaceNormal, + faceSurfaceGlobalRegion, - 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(), - // 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, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C index b9d6f3b0c6bcb2771130e4c3c3fcd962c1fde4a0..77bff601870e835f49db1f33c27017047c517484 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C @@ -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) + ) {} diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H index 8103cd4082dac4c468603a86eff1d4c5d9a320f6..294d0de55d5538df2e9557262be482eb1d9b1157 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H @@ -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_; + } + }; diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H index 918aae9000d4b6c3603f0c057824b102e3a91c8e..b98230e4122c925527c9158dbcc57b8bc8c8abc4 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H @@ -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, diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C index 4d9d10ab80106b2bd566b04dcc8774a83b782d01..2c15a847ead3d6571c75656d7352cc9f762d3bb0 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C @@ -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 diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C index 0cdd440d1a1d8a0cf40cfdd9c6579a58f687f41b..e72f13d72663e3683a878f7267c116d2428be75e 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C @@ -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; +} + + // ************************************************************************* // diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H index 5442e140ddfab7b4626bfad256239fe6e4749561..aa2807bac16f008d5ba5ca31c1a1967338d01e1f 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H +++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H @@ -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; + }; diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C index 64a03c43b1d61e5feca2ca66db5e3d28e8848670..33564024b1c13146cacebe7e0add6ccb3005b64e 100644 --- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C +++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C @@ -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: { + const typename fieldType::GeometricBoundaryField& + fieldBoundary = field.boundaryField(); + 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]]; + 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); diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C index 62ac6516efbe0c408c2fa518c6737d38019d4b81..71a097f55f6b9b51aec1498e4f96657f272e920d 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C @@ -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(); diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C index 9f596c3cdec5589f8de979438241b12b1749023b..81a992b9ab2a93a89d510d3c2b059e68fa71705c 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C @@ -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(); diff --git a/tutorials/mesh/cvMesh/simpleShapes/Allrun b/tutorials/mesh/cvMesh/simpleShapes/Allrun index f1259d32d27f3acd46796fd6c2afac3c00cf10ae..4c2ecc94734fd310575de005bb4df43b710db1bb 100755 --- a/tutorials/mesh/cvMesh/simpleShapes/Allrun +++ b/tutorials/mesh/cvMesh/simpleShapes/Allrun @@ -12,11 +12,11 @@ runApplication surfaceClean \ constant/triSurface/coneAndSphere_clean.obj mv log.surfaceClean log.surfaceClean.coneAndSphere # Orient so point to be meshed is inside surface -surfaceOrient \ +runApplication surfaceOrient \ constant/triSurface/coneAndSphere_clean.obj \ -inside '(0 -0.5 0)' \ - constant/triSurface/coneAndSphere_clean_orient.obj \ - > log.surfaceOrient.coneAndSphere 2>&1 + constant/triSurface/coneAndSphere_clean_orient.obj +mv log.surfaceOrient log.surfaceOrient.coneAndSphere # Same for outside runApplication surfaceClean \ @@ -24,11 +24,11 @@ runApplication surfaceClean \ 1e-4 1e-6 \ constant/triSurface/domain_clean.stl mv log.surfaceClean log.surfaceClean.domain -surfaceOrient \ +runApplication surfaceOrient \ constant/triSurface/domain_clean.stl \ -inside '(0 -0.5 0)' \ constant/triSurface/domain_clean_orient.stl - > log.surfaceOrient.domain 2>&1 +mv log.surfaceOrient log.surfaceOrient.domain runApplication surfaceFeatureExtract mv log.surfaceFeatureExtract log.surfaceFeatureExtract.coneAndSphere_clean