diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files index 16b1c4e44621b4735a104dd4ae38299f3d4f8635..f130c2166194b00298ea8e6d6fe63421e1212683 100644 --- a/src/fvOptions/Make/files +++ b/src/fvOptions/Make/files @@ -56,6 +56,7 @@ $(interRegion)/interRegionExplicitPorositySource/interRegionExplicitPorositySour /* Constraints */ generalConstraints=constraints/general $(generalConstraints)/fixedValueConstraint/fixedValueConstraints.C +$(generalConstraints)/mapFieldConstraint/mapFieldConstraints.C derivedConstraints=constraints/derived $(derivedConstraints)/fixedTemperatureConstraint/fixedTemperatureConstraint.C diff --git a/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.C b/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.C new file mode 100644 index 0000000000000000000000000000000000000000..831ac0b664d3882c8c6c9ab792cb7d578f1e8285 --- /dev/null +++ b/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.C @@ -0,0 +1,454 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2023 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "MapFieldConstraint.H" +#include "fvMatrices.H" +#include "meshToMesh.H" + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ +static inline tmp<volScalarField> createField +( + const fvMesh& mesh, + const scalar val +) +{ + return tmp<volScalarField>::New + ( + IOobject + ( + polyMesh::defaultRegion, + mesh.time().timeName(), + mesh, + IOobjectOption::NO_READ, + IOobject::NO_WRITE, + IOobject::NO_REGISTER + ), + mesh, + dimensionedScalar(dimless, val) + ); +} +} // End namespace fv +} // End namespace Foam + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template<class Type> +void Foam::fv::MapFieldConstraint<Type>::setSourceMesh +( + refPtr<fvMesh>& meshRef, + const autoPtr<Time>& runTimePtr +) +{ + const Time& runTime = runTimePtr(); + const word meshName(polyMesh::defaultRegion); + + // Fetch mesh from Time database + meshRef.cref + ( + runTime.cfindObject<fvMesh>(meshName) + ); + + if (!meshRef) + { + // Fallback: load mesh from disk and cache it + meshRef.reset + ( + new fvMesh + ( + IOobject + ( + meshName, + runTime.timeName(), + runTime, + IOobject::MUST_READ, + IOobject::NO_WRITE, + IOobject::REGISTER + ) + ) + ); + } +} + + +template<class Type> +void Foam::fv::MapFieldConstraint<Type>::createInterpolation +( + const fvMesh& srcMesh, + const fvMesh& tgtMesh +) +{ + if (consistent_) + { + interpPtr_.reset + ( + new meshToMesh + ( + srcMesh, + tgtMesh, + mapMethodName_, + patchMapMethodName_ + ) + ); + } + else + { + interpPtr_.reset + ( + new meshToMesh + ( + srcMesh, + tgtMesh, + mapMethodName_, + patchMapMethodName_, + patchMap_, + cuttingPatches_ + ) + ); + } +} + + +template<class Type> +template<class VolFieldType> +VolFieldType& Foam::fv::MapFieldConstraint<Type>::getOrReadField +( + const fvMesh& thisMesh, + const word& fieldName +) const +{ + auto* ptr = thisMesh.getObjectPtr<VolFieldType>(fieldName); + + if (!ptr) + { + ptr = new VolFieldType + ( + IOobject + ( + fieldName, + thisMesh.time().timeName(), + thisMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE, + IOobject::REGISTER + ), + thisMesh + ); + thisMesh.objectRegistry::store(ptr); + } + + return *ptr; +} + + +template<class Type> +Foam::labelList Foam::fv::MapFieldConstraint<Type>::tgtCellIDs() const +{ + const fvMesh& srcMesh = srcMeshPtr_(); + const fvMesh& tgtMesh = mesh_; + + // Create mask fields + const volScalarField srcFld(createField(srcMesh, 1)); + volScalarField tgtFld(createField(tgtMesh, 0)); + + // Map the mask field of 1s onto the mask field of 0s + interpPtr_->mapSrcToTgt(srcFld, plusEqOp<scalar>(), tgtFld); + + // Identify and collect cell indices whose values were changed from 0 to 1 + DynamicList<label> cells; + forAll(tgtFld, i) + { + if (tgtFld[i] != 0) + { + cells.append(i); + } + } + + return cells; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Type> +Foam::fv::MapFieldConstraint<Type>::transform::transform() +: + positionPtr_(), + directionPtr_(), + points_(), + origin_(), + normal_(), + active_(false) +{} + + +template<class Type> +Foam::fv::MapFieldConstraint<Type>::MapFieldConstraint +( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh +) +: + fv::cellSetOption(name, modelType, dict, mesh), + transform_(), + srcTimePtr_(), + srcMeshPtr_(), + interpPtr_(), + patchMap_(), + cuttingPatches_(), + mapMethodName_(), + patchMapMethodName_(), + consistent_(false) +{ + read(dict); + + setSourceMesh(srcMeshPtr_, srcTimePtr_); + + const fvMesh& srcMesh = srcMeshPtr_(); + const fvMesh& tgtMesh = mesh_; + createInterpolation(srcMesh, tgtMesh); + + // Overwrite the members of cellSetOption + cells_ = tgtCellIDs(); + V_ = interpPtr_->V(); + + transform_.initialize(srcMesh, dict); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +bool Foam::fv::MapFieldConstraint<Type>::transform::initialize +( + const fvMesh& srcMesh, + const dictionary& dict +) +{ + const dictionary* subDictPtr = dict.findDict("transform"); + + if (!subDictPtr) + { + return false; + } + + positionPtr_.reset + ( + Function1<point>::NewIfPresent + ( + "position", + *subDictPtr, + word::null, + &srcMesh + ) + ); + + directionPtr_.reset + ( + Function1<point>::NewIfPresent + ( + "direction", + *subDictPtr, + word::null, + &srcMesh + ) + ); + + if (positionPtr_) + { + subDictPtr->readIfPresent("origin", origin_); + } + + if (directionPtr_) + { + subDictPtr->readIfPresent("normal", normal_); + normal_.normalise(); + } + + points_ = srcMesh.points(); + + active_ = true; + + return true; +} + + +template<class Type> +void Foam::fv::MapFieldConstraint<Type>::transform::translate +( + refPtr<fvMesh>& srcMeshPtr, + const scalar t +) +{ + if (!positionPtr_) + { + return; + } + + const pointField translate + ( + points_ + (positionPtr_->value(t) - origin_) + ); + + fvMesh& srcMesh = srcMeshPtr.ref(); + srcMesh.movePoints(translate); +} + + +template<class Type> +void Foam::fv::MapFieldConstraint<Type>::transform::rotate +( + refPtr<fvMesh>& srcMeshPtr, + const scalar t +) +{ + if (!directionPtr_) + { + return; + } + + const vector dir(normalised(directionPtr_->value(t))); + + const tensor rot(rotationTensor(normal_, dir)); + + pointField rotate(points_); + + Foam::transform(rotate, rot, rotate); + + fvMesh& srcMesh = srcMeshPtr.ref(); + srcMesh.movePoints(rotate); +} + + +template<class Type> +bool Foam::fv::MapFieldConstraint<Type>::read(const dictionary& dict) +{ + if (!fv::cellSetOption::read(dict)) + { + return false; + } + + fieldNames_.resize(1, coeffs_.getWord("field")); + + fv::option::resetApplied(); + + // Load the time database for the source mesh once per simulation + if (!srcTimePtr_) + { + fileName srcMesh(coeffs_.get<fileName>("srcMesh").expand()); + srcMesh.clean(); + + srcTimePtr_.reset(Time::New(srcMesh)); + + // Set time-step of source database to an arbitrary yet safe value + srcTimePtr_().setDeltaT(1.0); + } + + coeffs_.readEntry("mapMethod", mapMethodName_); + if (!meshToMesh::interpolationMethodNames_.found(mapMethodName_)) + { + FatalIOErrorInFunction(coeffs_) + << type() << " " << name() << ": unknown map method " + << mapMethodName_ << nl + << "Available methods include: " + << meshToMesh::interpolationMethodNames_ + << exit(FatalIOError); + } + + coeffs_.readIfPresent("consistent", consistent_); + coeffs_.readIfPresent("patchMap", patchMap_); + coeffs_.readIfPresent("cuttingPatches", cuttingPatches_); + + if (!coeffs_.readIfPresent("patchMapMethod", patchMapMethodName_)) + { + meshToMesh::interpolationMethod mapMethod + ( + meshToMesh::interpolationMethodNames_[mapMethodName_] + ); + patchMapMethodName_ = meshToMesh::interpolationMethodAMI(mapMethod); + } + + return true; +} + + +template<class Type> +void Foam::fv::MapFieldConstraint<Type>::constrain +( + fvMatrix<Type>& eqn, + const label +) +{ + DebugInfo + << "MapFieldConstraint<" + << pTraits<Type>::typeName + << ">::constrain for source " << name_ << endl; + + // Translate and/or rotate source mesh if requested + if (transform_.isActive()) + { + // Use time from mesh_ since source mesh does not advance in time + const scalar t = mesh_.time().timeOutputValue(); + transform_.translate(srcMeshPtr_, t); + transform_.rotate(srcMeshPtr_, t); + } + + typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType; + + const word& fldName = fieldNames_[0]; + + const fvMesh& srcMesh = srcMeshPtr_(); + const fvMesh& tgtMesh = mesh_; + + // Fetch source and target fields + const VolFieldType& srcFld = getOrReadField<VolFieldType>(srcMesh, fldName); + VolFieldType& tgtFld = tgtMesh.lookupObjectRef<VolFieldType>(fldName); + + // When mesh/src changes, reinitilize mesh-to-mesh and cellSetOption members + if (tgtMesh.changing() || transform_.isActive()) + { + createInterpolation(srcMesh, tgtMesh); + cells_ = tgtCellIDs(); + V_ = interpPtr_->V(); + } + + // Map source-mesh field onto target-mesh field + interpPtr_->mapSrcToTgt(srcFld, plusEqOp<Type>(), tgtFld); + + // Constrain mapped field in target mesh to avoid overwrite by solver + eqn.setValues(cells_, UIndirectList<Type>(tgtFld, cells_)); +} + + +// ************************************************************************* // diff --git a/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.H b/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.H new file mode 100644 index 0000000000000000000000000000000000000000..df06a4ed2bb85a5bdda52263ece394533fac84a4 --- /dev/null +++ b/src/fvOptions/constraints/general/mapFieldConstraint/MapFieldConstraint.H @@ -0,0 +1,311 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2023 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::fv::MapFieldConstraint + +Group + grpFvOptionsConstraints + +Description + The \c MapFieldConstraint constrains values of given fields of \c Type + with a source field from an external mesh, where + \c \<Type\>={scalar,vector,sphericalTensor,symmTensor,tensor}. + Optionally, the source field can be translated and/or rotated as a function + of time. + +Usage + Minimal example by using \c constant/fvOptions: + \verbatim + \<Type\>MapFieldConstraint1 + { + // Mandatory entries + type \<Type\>MapFieldConstraint; + field <word>; + srcMesh <fileName>; + mapMethod <word>; + + // Optional entries + consistent <bool>; + patchMapMethod <word>; + transform + { + // Optional entries + position <Function1<vector>>; + origin <vector>; + + direction <Function1<vector>>; + normal <vector>; + } + + // Conditional entries + + // when consistent=false + patchMap <HashTable<word>>; // (<patchSrc> <patchTgt>); + cuttingPatches <wordList>; // (<patchTgt1> ... <patchTgtN>); + + // Inherited entries + ... + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + type | Type name: \<Type\>MapFieldConstraint | word | yes | - + field | Name of operand field | word | yes | - + srcMesh | Directory path to mesh to map from | fileName | yes | - + mapMethod | Mapping method | word | yes | - + consistent | Flag to determine if meshes have consistent boundaries <!-- + --> | bool | no | false + patchMapMethod | Name of patch-map method | word | no | - + patchMap | Coincident source/target patches in two cases <!-- + --> | wordHashTable | no | - + cuttingPatches | Target patches cutting the source domain <!-- + --> | wordList | no | - + transform | Transform settings for source mesh points <!-- + --> | dict | no | - + position | Position of source mesh as a function of time <!-- + --> | Function1\<vector\> | no | - + direction | Direction of source mesh as a function of time <!-- + --> | Function1\<vector\> | no | - + origin | Origin of source mesh | vector | no | - + normal | Normal of reference plane representing source mesh <!-- + --> | vector | no | - + \endtable + + The inherited entries are elaborated in: + - \link cellSetOption.H \endlink + - \link Function1.H \endlink + +SourceFiles + MapFieldConstraint.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_fv_MapFieldConstraint_H +#define Foam_fv_MapFieldConstraint_H + +#include "cellSetOption.H" +#include "Function1.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class meshToMesh; + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class MapFieldConstraint Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type> +class MapFieldConstraint +: + public fv::cellSetOption +{ + // Private Classes + + class transform + { + // Private Data + + //- Position of source mesh as a function of time + autoPtr<Function1<point>> positionPtr_; + + //- Direction of source mesh as a function of time + autoPtr<Function1<point>> directionPtr_; + + //- Cached points of source mesh + pointField points_; + + //- Origin of source mesh + point origin_; + + //- Normal of reference plane representing source mesh + vector normal_; + + //- Flag to deduce if transformation is active + bool active_; + + + public: + + // Constructors + + //- Default construct + transform(); + + //- No copy construct + transform(const transform&) = delete; + + //- No copy assignment + void operator=(const transform&) = delete; + + + // Member Functions + + // Access + + //- Return flag to deduce if transformation is active + bool isActive() const noexcept { return active_; } + + + // Evaluation + + //- Translate source mesh as a function of time + void translate(refPtr<fvMesh>& srcMeshPtr, const scalar time); + + //- Rotate source mesh as a function of time + void rotate(refPtr<fvMesh>& srcMeshPtr, const scalar time); + + + // I-O + + //- Initialize the class members + bool initialize(const fvMesh& srcMesh, const dictionary& dict); + }; + + + // Private Data + + //- Transformation settings for source mesh + transform transform_; + + //- Time database for source mesh to map from + autoPtr<Time> srcTimePtr_; + + //- Source mesh to map from + refPtr<fvMesh> srcMeshPtr_; + + //- Mesh-to-mesh interpolation from source mesh to target mesh + autoPtr<meshToMesh> interpPtr_; + + //- List of coincident source/target patches in two cases + HashTable<word> patchMap_; + + //- List of names of target patches cutting the source domain + wordList cuttingPatches_; + + //- Name of map method + word mapMethodName_; + + //- Name of patch-map method + word patchMapMethodName_; + + //- Flag to determine if meshes have consistent boundaries + bool consistent_; + + + // Private Member Functions + + //- Helper function to set source mesh + // Fetch fvMesh from a given Time database + // Otherwise, load it from disk and cache it to the database + void setSourceMesh + ( + refPtr<fvMesh>& meshRef, + const autoPtr<Time>& runTimePtr + ); + + //- Helper function to create the mesh-to-mesh interpolation + void createInterpolation + ( + const fvMesh& srcMesh, + const fvMesh& tgtMesh + ); + + //- Return requested field from object registry + //- otherwise read it from disk and register it to the object registry + template<class VolFieldType> + VolFieldType& getOrReadField + ( + const fvMesh& thisMesh, + const word& fieldName + ) const; + + //- Return the local cell indices of the target mesh + labelList tgtCellIDs() const; + + +public: + + //- Runtime type information + TypeName("MapFieldConstraint"); + + + // Constructors + + //- Construct from components + MapFieldConstraint + ( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + //- No copy construct + MapFieldConstraint(const MapFieldConstraint&) = delete; + + //- No copy assignment + void operator=(const MapFieldConstraint&) = delete; + + + //- Destructor + virtual ~MapFieldConstraint() = default; + + + // Member Functions + + //- Read source dictionary + virtual bool read(const dictionary& dict); + + //- Set value on field + virtual void constrain(fvMatrix<Type>& eqn, const label); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "MapFieldConstraint.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/fvOptions/constraints/general/mapFieldConstraint/mapFieldConstraints.C b/src/fvOptions/constraints/general/mapFieldConstraint/mapFieldConstraints.C new file mode 100644 index 0000000000000000000000000000000000000000..e9279a4d061462a63807e418f67b16b6bb9d8aca --- /dev/null +++ b/src/fvOptions/constraints/general/mapFieldConstraint/mapFieldConstraints.C @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2023 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "makeFvOption.H" +#include "MapFieldConstraint.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFvOption(MapFieldConstraint, scalar); +makeFvOption(MapFieldConstraint, vector); +makeFvOption(MapFieldConstraint, sphericalTensor); +makeFvOption(MapFieldConstraint, symmTensor); +makeFvOption(MapFieldConstraint, tensor); + +// ************************************************************************* //