diff --git a/ThirdParty/AllMake b/ThirdParty/AllMake index ab563f980..48f14a41c 100755 --- a/ThirdParty/AllMake +++ b/ThirdParty/AllMake @@ -77,10 +77,6 @@ echo # this stage should be called last in your compilation process ./AllMake.stage5 -# Running stage eigen -./AllMake.eigen - - echo ======================================== echo Done ThirdParty Allwmake echo ======================================== diff --git a/ThirdParty/AllMake.eigen b/ThirdParty/AllMake.eigen deleted file mode 100755 index 009c39b55..000000000 --- a/ThirdParty/AllMake.eigen +++ /dev/null @@ -1,77 +0,0 @@ -#!/bin/bash -#------------------------------------------------------------------------------ -# ========= | -# \\ / F ield | foam-extend: Open Source CFD -# \\ / O peration | -# \\ / A nd | For copyright notice see file Copyright -# \\/ M anipulation | -#------------------------------------------------------------------------------ -# License -# This file is part of foam-extend. -# -# foam-extend 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. -# -# foam-extend 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 foam-extend. If not, see . -# -# Script -# AllMake.stage3 -# -# Description -# Build script for ThirdParty packages: eigen -# -# The ThirdParty libraries -# -# Requirements: -# 1: Your foam-extend environment must be properly initialized -# 2: AllMake.stage1 if you are overriding your system compiler -# 3: AllMake.stage2 if you are overriding your system comm. libraries -# -# Author: -# Martin Beaudoin, Hydro-Quebec, (2015) -# -#------------------------------------------------------------------------------ -# run from third-party directory only -cd ${0%/*} || exit 1 - -wmakeCheckPwd "$WM_THIRD_PARTY_DIR" || { - echo "Error: Current directory is not \$WM_THIRD_PARTY_DIR" - echo " The environment variables are inconsistent with the installation." - echo " Check the foam-extend entries in your dot-files and source them." - exit 1 -} -. tools/makeThirdPartyFunctionsForRPM -#------------------------------------------------------------------------------ - -echo ======================================== -echo Starting ThirdParty AllMake: eigen -echo ======================================== -echo - -# Eigen -if [ -d $WM_THIRD_PARTY_DIR/packages/eigen3 ] -then - echo "eigen3 found." -else - # Download eigen using wget - echo "Downloading eigen3 to eigen3" - wget http://bitbucket.org/eigen/eigen/get/3.2.8.tar.gz - \tar xvzf 3.2.8.tar.gz - \rm -f 3.2.8.tar.gz - \mv eigen-eigen-07105f7124f9 ./packages/eigen3 -fi - -echo ======================================== -echo Done ThirdParty AllMake: Stage3 -echo ======================================== -echo - -# ----------------------------------------------------------------- end-of-file diff --git a/ThirdParty/AllMake.pre b/ThirdParty/AllMake.pre index f5e2a167e..cebcc2a6f 100755 --- a/ThirdParty/AllMake.pre +++ b/ThirdParty/AllMake.pre @@ -72,8 +72,6 @@ echo # Running stage 4 ./AllMake.stage4 -# Running stage eigen -./AllMake.eigen echo ======================================== echo Done ThirdParty Allwmake.pre diff --git a/applications/solvers/FSI/Allwclean b/applications/solvers/FSI/Allwclean deleted file mode 100755 index 502a97429..000000000 --- a/applications/solvers/FSI/Allwclean +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/sh -set -x - -wclean fluidSolidInteraction - -wclean solvers/fsiFoam -wclean solvers/ampFsiFoam -wclean solvers/weakFsiFoam -wclean solvers/fluidFoam -wclean solvers/solidFoam -wclean solvers/thermalSolidFoam - -wclean utilities/functionObjects/pointHistory -wclean utilities/functionObjects/patchAvgTractionHistory -wclean utilities/functionObjects/centrifugalBodyForce - -# Wipe out all lnInclude directories and re-link -wcleanLnIncludeAll diff --git a/applications/solvers/FSI/Allwmake b/applications/solvers/FSI/Allwmake deleted file mode 100755 index 91f924cf0..000000000 --- a/applications/solvers/FSI/Allwmake +++ /dev/null @@ -1,15 +0,0 @@ -#!/bin/sh -set -x - -wmake libso fluidSolidInteraction - -wmake solvers/fsiFoam -wmake solvers/ampFsiFoam -wmake solvers/weakFsiFoam -wmake solvers/fluidFoam -wmake solvers/solidFoam -wmake solvers/thermalSolidFoam - -wmake libso utilities/functionObjects/pointHistory -wmake libso utilities/functionObjects/patchAvgTractionHistory -wmake libso utilities/functionObjects/centrifugalBodyForce diff --git a/applications/solvers/FSI/fluidSolidInteraction/Make/files b/applications/solvers/FSI/fluidSolidInteraction/Make/files deleted file mode 100755 index cbf9c1995..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/Make/files +++ /dev/null @@ -1,141 +0,0 @@ - -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFCoarsening.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.C -fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.C - -contactModels = solidSolvers/solidModels/fvPatchFields/solidContact/contactModels -$(contactModels)/frictionContactModels/frictionContactModel/frictionContactModel.C -$(contactModels)/frictionContactModels/frictionContactModel/newFrictionContactModel.C -$(contactModels)/frictionContactModels/frictionless/frictionless.C -$(contactModels)/frictionContactModels/standardPenaltyFriction/standardPenaltyFriction.C -$(contactModels)/frictionContactModels/frictionLaws/frictionLaw/frictionLaw.C -$(contactModels)/frictionContactModels/frictionLaws/frictionLaw/newFrictionLaw.C -$(contactModels)/frictionContactModels/frictionLaws/coulombFriction/coulombFriction.C -$(contactModels)/frictionContactModels/frictionLaws/coulombOrowanFriction/coulombOrowanFriction.C -$(contactModels)/frictionContactModels/frictionLaws/timeVaryingCoulomb/timeVaryingCoulomb.C -$(contactModels)/normalContactModels/normalContactModel/normalContactModel.C -$(contactModels)/normalContactModels/normalContactModel/newNormalContactModel.C -$(contactModels)/normalContactModels/standardPenalty/standardPenalty.C - -fluidSolvers/fluidSolver/fluidSolver.C -fluidSolvers/fluidSolver/newFluidSolver.C -fluidSolvers/icoFluid/icoFluid.C -fluidSolvers/pisoFluid/pisoFluid.C -fluidSolvers/consistentIcoFluid/consistentIcoFluid.C - -fluidSolvers/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.C -fluidSolvers/fvPatchFields/extrapolated/extrapolatedFvPatchFields.C -fluidSolvers/fvPatchFields/extrapolatedPressure/extrapolatedPressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/transitionalParabolicVelocity/transitionalParabolicVelocityFvPatchVectorField.C -fluidSolvers/fvPatchFields/movingWallVelocity/movingWallVelocityFvPatchVectorField.C -fluidSolvers/fvPatchFields/fixedVelocityPressure/fixedVelocityPressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/normalInletVelocity/normalInletVelocityFvPatchVectorField.C -fluidSolvers/fvPatchFields/outflowPressure/outflowPressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/movingWallPressure/movingWallPressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/robin/robinFvPatchFields.C -fluidSolvers/fvPatchFields/elasticWallPressure/elasticWallPressureFvPatchScalarField.C -fluidSolvers/fvPatchFields/skewCorrectedLinearDC/skewCorrectedLinearDC.C -fluidSolvers/fvPatchFields/elasticSlipWallVelocity/elasticSlipWallVelocityFvPatchVectorField.C -fluidSolvers/fvPatchFields/elasticWallVelocity/elasticWallVelocityFvPatchVectorField.C - -solidSolvers/solidModels/fvPatchFields/tractionDisplacement/tractionDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/tractionDisplacementIncrement/tractionDisplacementIncrementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/symmetryDisplacement/symmetryDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/fixedDisplacement/fixedDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/fixedNormalDisplacement/fixedNormalDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/fixedNormalDisplacementIncrement/fixedNormalDisplacementIncrementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/componentMixed/componentMixedPointPatchFields.C -solidSolvers/solidModels/fvPatchFields/cohesiveZoneModeI/cohesiveZoneFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/cohesiveZoneIncrementalModeI/cohesiveZoneIncrementalFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/planeContactDisplacement/planeContactDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/directionMixedDisplacement/directionMixedDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/timeVaryingFixedNormalDisplacement/timeVaryingFixedNormalDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/timeVaryingFixedDisplacement/timeVaryingFixedDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/solidWedge/solidWedgeFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/timeVaryingPressureDisplacement/timeVaryingPressureDisplacementFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/fixedHydPressure/fixedHydPressureFvPatchScalarField.C -solidSolvers/solidModels/fvPatchFields/fixedHydPressureIncrement/fixedHydPressureIncrementFvPatchScalarField.C -solidSolvers/solidModels/fvPatchFields/convectiveFlux/convectiveFluxFvPatchScalarField.C -solidSolvers/solidModels/fvPatchFields/wedge/wedgeFvPatchFields.C -solidSolvers/solidModels/fvPatchFields/wedge/wedgeFvPatchScalarField.C -solidSolvers/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C -solidSolvers/solidModels/fvPatchFields/correctedSymmetry/correctedSymmetryFvPatchFields.C -solidSolvers/solidModels/fvPatchFields/velocityTractionDisplacement/velocityTractionDisplacementFvPatchVectorField.C - -solidSolvers/solidSolver/solidSolver.C -solidSolvers/solidSolver/newSolidSolver.C -solidSolvers/unsTotalLagrangianSolid/unsTotalLagrangianSolid.C -solidSolvers/unsTotalLagrangianSolid/unsTotalLagrangianSolidSolve.C -solidSolvers/unsIncrTotalLagrangianSolid/unsIncrTotalLagrangianSolid.C -solidSolvers/unsIncrTotalLagrangianSolid/unsIncrTotalLagrangianSolidSolve.C - -fluidSolvers/finiteVolume/ddtSchemes/EulerDdtScheme.C -fluidSolvers/finiteVolume/ddtSchemes/backwardDdtScheme.C - -solidSolvers/finiteVolume/fvc/fvcCellLimitedGrad.C - -solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/plasticityStressReturn/plasticityStressReturn.C -solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/plasticityStressReturn/newPlasticityStressReturn.C -solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/aravasMises/aravasMises.C -solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/newAravasMises/newAravasMises.C -solidSolvers/solidModels/constitutiveModel/constitutiveModel.C -solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/kirchhoffMises/kirchhoffMises.C - -rheologyLaws = solidSolvers/solidModels/constitutiveModel/rheologyLaws -$(rheologyLaws)/rheologyLaw/rheologyLaw.C -$(rheologyLaws)/rheologyLaw/newRheologyLaw.C -$(rheologyLaws)/linearElastic/linearElastic.C -$(rheologyLaws)/elasticPlastic/elasticPlastic.C -$(rheologyLaws)/linearElasticTabulatedPlastic/linearElasticTabulatedPlastic.C -$(rheologyLaws)/multiMaterial/multiMaterial.C -$(rheologyLaws)/orthotropicElasticPlastic/orthotropicElasticPlastic.C -$(rheologyLaws)/orthotropicLinearElastic/orthotropicLinearElastic.C -$(rheologyLaws)/PronyViscoelastic/PronyViscoelastic.C -$(rheologyLaws)/expElasticPlastic/expElasticPlastic.C - - -solidSolvers/solidModels/thermalModel/thermalModel.C - -thermalLaws = solidSolvers/solidModels/thermalModel/thermalLaws -$(thermalLaws)/thermalLaw/thermalLaw.C -$(thermalLaws)/thermalLaw/newThermalLaw.C -$(thermalLaws)/constantThermal/constantThermal.C -$(thermalLaws)/multiMaterialThermal/multiMaterialThermal.C - - -solidSolvers/solidModels/arbitraryCrack/crackerMovingFvMesh/crackerMovingFvMesh.C -solidSolvers/solidModels/arbitraryCrack/crackerFvMesh/crackerFvMesh.C -solidSolvers/solidModels/arbitraryCrack/faceCracker/faceCracker.C -solidSolvers/solidModels/arbitraryCrack/faceCracker/detachFaceCracker.C -solidSolvers/solidModels/arbitraryCrack/cohesive/cohesivePolyPatch.C -solidSolvers/solidModels/arbitraryCrack/cohesiveLaw/cohesiveLawFvPatchVectorField.C - -solidSolvers/solidModels/SolidInterfaces/solidInterfaceTL/solidInterfaceTL.C -solidSolvers/solidModels/SolidInterfaces/solidInterfaceITL/solidInterfaceITL.C -solidSolvers/solidModels/materialInterfaces/materialInterface/materialInterface.C -solidSolvers/solidModels/materialInterfaces/TLMaterialInterface/TLMaterialInterface.C -solidSolvers/solidModels/materialInterfaces/ITLMaterialInterface/ITLMaterialInterface.C -solidSolvers/solidModels/materialInterfaces/ULLSMaterialInterface/ULLSMaterialInterface.C - - -solidSolvers/solidModels/simpleCohesiveLaws/simpleCohesiveLaw/simpleCohesiveLaw.C -solidSolvers/solidModels/simpleCohesiveLaws/linear/linearSimpleCohesiveLaw.C -solidSolvers/solidModels/simpleCohesiveLaws/Dugdale/DugdaleSimpleCohesiveLaw.C -solidSolvers/solidModels/eig3/eig3.C -solidSolvers/solidModels/eig3/eig3Field.C - - -solidSolvers/solidModels/nonLinearGeometry/nonLinearGeometry.C -solidSolvers/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C - -fluidSolidInterface/fluidSolidInterface.C - -LIB = $(FOAM_LIBBIN)/libfluidSolidInteraction diff --git a/applications/solvers/FSI/fluidSolidInteraction/Make/options b/applications/solvers/FSI/fluidSolidInteraction/Make/options deleted file mode 100755 index 50e3dce8a..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/Make/options +++ /dev/null @@ -1,54 +0,0 @@ -c++WARN = -Wno-deprecated -Wall -Wextra -Wno-unused-parameter -Wnon-virtual-dtor -Wunused-local-typedefs -Werror -Wredundant-decls -Wcast-align -Wmissing-declarations -Wswitch-enum -Winvalid-pch -Wredundant-decls -Wformat=2 -Wmissing-format-attribute -Wformat-nonliteral - - -EXE_INC = -std=c++11 \ - -I./fluidSolvers/finiteVolume/RBFMeshMotionSolver \ - -I./fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions \ - -I./solidSolvers/solidModels/fvPatchFields/symmetry \ - -I./solidSolvers/solidModels/contactModels/frictionContactModels/frictionContactModel/frictionLaws/frictionLaw \ - -I./solidSolvers/solidModels/contactModels/frictionContactModels/frictionContactModel \ - -I./solidSolvers/solidModels/contactModels/normalContactModels/normalContactModel \ - -I./solidSolvers/solidModels/constitutiveModel/tractionBoundaryGradient \ - -I./solidSolvers/solidModels/nonLinearGeometry \ - -I./solidSolvers/solidModels/constitutiveModel/rheologyLaws/PronyViscoelastic \ - -I./solidSolvers/solidModels/arbitraryCrack/pointDistDiff \ - -I$(WM_THIRD_PARTY_DIR)/packages/eigen3 \ - -I./solidSolvers/solidModels/componentReference \ - -I./solidSolvers/solidModels/constitutiveModel/lnInclude \ - -I./solidSolvers/solidModels/constitutiveModel/plasticityStressReturnMethods/plasticityStressReturn \ - -I./solidSolvers/solidModels/constitutiveModel/rheologyLaws/rheologyLaw \ - -I./solidSolvers/solidModels/constitutiveModel \ - -I./fluidSolidInterface \ - -I./solidSolvers/finiteVolume/fvc \ - -I./solidSolvers/solidModels/fvPatchFields/tractionDisplacement \ - -I./solidSolvers/solidModels/fvPatchFields/tractionDisplacementIncrement \ - -I./solidSolvers/solidModels/fvPatchFields/fixedDisplacement \ - -I./solidSolvers/solidModels/fvPatchFields/symmetryDisplacement \ - -I./solidSolvers/solidSolver \ - -I./solidSolvers/solidModels/arbitraryCrack/crackerFvMesh \ - -I./fluidSolvers/finiteVolume/ddtSchemes \ - -I./fluidSolvers/fluidSolver \ - -I$(LIB_SRC)/foam/lnInclude \ - -I$(LIB_SRC)/finiteArea/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/meshMotion/fvMotionSolver/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/meshMotion/tetMotionSolver/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/topoChangerFvMesh/lnInclude \ - -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/tetFiniteElement/lnInclude \ - -EXE_LIBS = \ - -lincompressibleTurbulenceModel \ - -lincompressibleRASModels \ - -lincompressibleLESModels \ - -lincompressibleTransportModels \ - -lfiniteVolume \ - -ldynamicFvMesh \ - -ldynamicMesh \ - -lmeshTools diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.C deleted file mode 100644 index a5dbb7316..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.C +++ /dev/null @@ -1,2841 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "fluidSolidInterface.H" -#include "volFields.H" -#include "polyPatchID.H" -#include "ZoneIDs.H" -#include "RectangularMatrix.H" -#include "primitivePatchInterpolation.H" -#include "twoDPointCorrector.H" - -#include "tetPointFields.H" -#include "fixedValueTetPolyPatchFields.H" -#include "tetPolyPatchInterpolation.H" -#include "tetFemMatrices.H" - -#include "fixedValuePointPatchFields.H" -#include "ggiInterpolation.H" -#include "IOmanip.H" -#include "dynamicMotionSolverFvMesh.H" -#include "motionSolver.H" -#include "elasticWallPressureFvPatchScalarField.H" -#include "movingWallPressureFvPatchScalarField.H" - -// #include "faCFD.H" - - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(fluidSolidInterface, 0); -// defineRunTimeSelectionTable(fluidSolidInterface, dictionary); -} - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::fluidSolidInterface::calcCurrentSolidZonePoints() const -{ - // Find global face zones - if (currentSolidZonePointsPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcCurrentSolidZonePoints() const" - ) - << "Current solid zone points alarady exist" - << abort(FatalError); - } - - currentSolidZonePointsPtr_ = - new vectorField(stress().currentFaceZonePoints(solidZoneIndex())); -} - -void Foam::fluidSolidInterface::calcCurrentSolidZonePatch() const -{ - // Find global face zones - if (currentSolidZonePatchPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcCurrentSolidZonePatch() const" - ) - << "Current solid zone patch alarady exists" - << abort(FatalError); - } - - currentSolidZonePatchPtr_ = - new PrimitivePatch - ( - solidMesh().faceZones()[solidZoneIndex_]().localFaces(), - currentSolidZonePoints() - ); -} - -void Foam::fluidSolidInterface::calcFluidToSolidInterpolator() const -{ - // Find global face zones - if (fluidToSolidPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcFluidToSolidInterpolator() const" - ) - << "Fluid to solid interpolator already exists" - << abort(FatalError); - } - - std::shared_ptr rbfFunction; - rbfFunction = std::shared_ptr(new TPSFunction()); - - // std::shared_ptr rbf; - fluidToSolidPtr_ = - std::shared_ptr - ( - new RBFInterpolation(rbfFunction) - // new RBFInterpolation(rbfFunctionPtr_) - ); - - vectorField solidZoneFaceCentres = - currentSolidZonePatch().faceCentres(); - - vectorField fluidZoneFaceCentres = - fluidMesh().faceZones()[fluidZoneIndex_]().faceCentres(); - - matrix fluidX(fluidZoneFaceCentres.size(), 3); - matrix solidX(solidZoneFaceCentres.size(), 3); - - forAll(fluidZoneFaceCentres, faceI) - { - fluidX(faceI, 0) = fluidZoneFaceCentres[faceI].x(); - fluidX(faceI, 1) = fluidZoneFaceCentres[faceI].y(); - fluidX(faceI, 2) = fluidZoneFaceCentres[faceI].z(); - } - - forAll(solidZoneFaceCentres, faceI) - { - solidX(faceI, 0) = solidZoneFaceCentres[faceI].x(); - solidX(faceI, 1) = solidZoneFaceCentres[faceI].y(); - solidX(faceI, 2) = solidZoneFaceCentres[faceI].z(); - } - - fluidToSolidPtr_->compute(fluidX, solidX); - - Info << "Checking fluid-to-solid interpolator" << endl; - { - vectorField fluidPatchFaceCentres = - vectorField - ( - fluidMesh().boundaryMesh()[fluidPatchIndex_].faceCentres() - ); - - vectorField fluidZoneFaceCentres - ( - fluidMesh().faceZones()[fluidZoneIndex_].size(), - vector::zero - ); - - const label fluidPatchStart = - fluidMesh().boundaryMesh()[fluidPatchIndex_].start(); - - forAll (fluidPatchFaceCentres, i) - { - fluidZoneFaceCentres - [ - fluidMesh().faceZones()[fluidZoneIndex_].whichFace - ( - fluidPatchStart + i - ) - ] = - fluidPatchFaceCentres[i]; - } - - // Parallel data exchange: collect faceCentres field on all processors - reduce(fluidZoneFaceCentres, sumOp()); - - vectorField solidPatchFaceCentres = - vectorField - ( - solidMesh().boundaryMesh()[solidPatchIndex_].faceCentres() - ); - - matrix fluidX(fluidZoneFaceCentres.size(), 3); - // matrix solidX(solidPatchFaceCentres.size(), 3); - matrix fluidXsolid(solidPatchFaceCentres.size(), 3); - - forAll(fluidZoneFaceCentres, faceI) - { - fluidX(faceI, 0) = fluidZoneFaceCentres[faceI].x(); - fluidX(faceI, 1) = fluidZoneFaceCentres[faceI].y(); - fluidX(faceI, 2) = fluidZoneFaceCentres[faceI].z(); - } - - // forAll(solidPatchFaceCentres, faceI) - // { - // solidX(faceI, 0) = solidPatchFaceCentres[faceI].x(); - // solidX(faceI, 1) = solidPatchFaceCentres[faceI].y(); - // solidX(faceI, 2) = solidPatchFaceCentres[faceI].z(); - // } - - // fluidToSolidPtr_->compute(fluidX, solidX); - fluidToSolidPtr_->interpolate(fluidX, fluidXsolid); - - vectorField fluidPatchFaceCentresAtSolid - ( - solidPatchFaceCentres.size(), - vector::zero - ); - - forAll(fluidPatchFaceCentresAtSolid, faceI) - { - fluidPatchFaceCentresAtSolid[faceI].x() = fluidXsolid(faceI, 0); - fluidPatchFaceCentresAtSolid[faceI].y() = fluidXsolid(faceI, 1); - fluidPatchFaceCentresAtSolid[faceI].z() = fluidXsolid(faceI, 2); - } - - scalar maxDist = gMax - ( - mag - ( - fluidPatchFaceCentresAtSolid - - solidPatchFaceCentres - ) - ); - - Info << "Fluid-to-solid face interpolation error: " << maxDist - << endl; - } -} - -// void Foam::fluidSolidInterface::calcGgiFluidToSolidInterpolator() const -// { -// // Find global face zones -// if (ggiFluidToSolidPtr_) -// { -// FatalErrorIn -// ( -// "void fluidSolidInterface::" -// "calcGgiFluidToSolidInterpolator() const" -// ) -// << "Ggi fluid to solid interpolator already exists" -// << abort(FatalError); -// } - -// ggiFluidToSolidPtr_ = -// new ggiZoneInterpolation -// ( -// fluidMesh().faceZones()[fluidZoneIndex_](), -// solidMesh().faceZones()[solidZoneIndex_](), -// tensorField(0), -// tensorField(0), -// vectorField(0), // Slave-to-master separation. Bug fix -// 0, // Non-overlapping face tolerances -// 0, // HJ, 24/Oct/2008 -// true, // Rescale weighting factors. Bug fix, MB. -// ggiInterpolation::AABB //BB_OCTREE // Octree search, MB. -// ); - - -// Info << "Checking fluid-to-solid interpolator" << endl; -// { -// vectorField fluidPatchFaceCentres = -// vectorField -// ( -// fluidMesh().boundaryMesh()[fluidPatchIndex_].faceCentres() -// ); - -// vectorField fluidZoneFaceCentres -// ( -// fluidMesh().faceZones()[fluidZoneIndex_].size(), -// vector::zero -// ); - -// const label fluidPatchStart = -// fluidMesh().boundaryMesh()[fluidPatchIndex_].start(); - -// forAll (fluidPatchFaceCentres, i) -// { -// fluidZoneFaceCentres -// [ -// fluidMesh().faceZones()[fluidZoneIndex_].whichFace -// ( -// fluidPatchStart + i -// ) -// ] = -// fluidPatchFaceCentres[i]; -// } - -// // Parallel data exchange: collect faceCentres field on all processors -// reduce(fluidZoneFaceCentres, sumOp()); - -// vectorField solidZoneFaceCentres = -// ggiFluidToSolidPtr_->masterToSlave -// ( -// fluidZoneFaceCentres -// ); - -// vectorField solidPatchFaceCentres -// ( -// solidMesh().boundaryMesh()[solidPatchIndex_].size(), -// vector::zero -// ); - -// const label solidPatchStart = -// solidMesh().boundaryMesh()[solidPatchIndex_].start(); - -// forAll(solidPatchFaceCentres, i) -// { -// solidPatchFaceCentres[i] = -// solidZoneFaceCentres -// [ -// solidMesh().faceZones()[solidZoneIndex_] -// .whichFace(solidPatchStart + i) -// ]; -// } - -// scalar maxDist = gMax -// ( -// mag -// ( -// solidPatchFaceCentres -// - solidMesh().boundaryMesh()[solidPatchIndex_].faceCentres() -// ) -// ); - -// Info << "Fluid-to-solid face interpolation error: " << maxDist -// << endl; -// } - -// Info << "Number of uncovered master faces: " -// << ggiFluidToSolidPtr_->uncoveredMasterFaces().size() << endl; - -// Info << "Number of uncovered slave faces: " -// << ggiFluidToSolidPtr_->uncoveredSlaveFaces().size() << endl; -// } - - -void Foam::fluidSolidInterface::calcGgiInterpolator() const -{ - // Create extended ggi interpolation - if (ggiInterpolatorPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcGgiInterpolator() const" - ) - << "Ggi interpolator already exists" - << abort(FatalError); - } - - // Create copy of solid face zone primitive patch in current configuration - - deleteDemandDrivenData(currentSolidZonePatchPtr_); - deleteDemandDrivenData(currentSolidZonePointsPtr_); - -// currentSolidZonePatch().movePoints(currentSolidZonePoints()); - - Info << "Create extended GGI zone-to-zone interpolator" << endl; - - ggiInterpolatorPtr_ = - new ggiZoneInterpolation - ( - fluidMesh().faceZones()[fluidZoneIndex_](), - currentSolidZonePatch(), -// solidMesh().faceZones()[solidZoneIndex_](), - tensorField(0), - tensorField(0), - vectorField(0), // Slave-to-master separation. Bug fix - true, // Patch data is complete on all processors - SMALL, // Non-overlapping face tolerances - SMALL, // HJ, 24/Oct/2008 - true, // Rescale weighting factors. Bug fix, MB. - ggiInterpolation::BB_OCTREE - // BB_OCTREE AABB THREE_D_DISTANCE - // Octree search, MB. - ); - -// currentSolidZonePatch().writeVTK("solidZone"); -// fluidMesh().faceZones()[fluidZoneIndex_]().writeVTK("fluidZone"); -// fluidMesh().boundaryMesh()[fluidPatchIndex_].writeVTK -// ( -// "fluidPatch", -// fluidMesh().boundaryMesh()[fluidPatchIndex_], -// fluidMesh().points() -// ); - - Info << "Checking fluid-to-solid face interpolator" << endl; - { - vectorField fluidPatchFaceCentres = - vectorField - ( - fluidMesh().boundaryMesh()[fluidPatchIndex_].faceCentres() - ); - - vectorField fluidZoneFaceCentres - ( - fluidMesh().faceZones()[fluidZoneIndex_].size(), - vector::zero - ); - - const label fluidPatchStart = - fluidMesh().boundaryMesh()[fluidPatchIndex_].start(); - - forAll (fluidPatchFaceCentres, i) - { - fluidZoneFaceCentres - [ - fluidMesh().faceZones()[fluidZoneIndex_].whichFace - ( - fluidPatchStart + i - ) - ] = - fluidPatchFaceCentres[i]; - } - - // Parallel data exchange: collect faceCentres field on all processors - reduce(fluidZoneFaceCentres, sumOp()); - - vectorField solidZoneFaceCentres = - ggiInterpolatorPtr_->masterToSlave - ( - fluidZoneFaceCentres - ); - - vectorField solidPatchFaceCentres - ( - solidMesh().boundaryMesh()[solidPatchIndex_].size(), - vector::zero - ); - - const label solidPatchStart = - solidMesh().boundaryMesh()[solidPatchIndex_].start(); - - forAll(solidPatchFaceCentres, i) - { - solidPatchFaceCentres[i] = - solidZoneFaceCentres - [ - solidMesh().faceZones()[solidZoneIndex_] - .whichFace(solidPatchStart + i) - ]; - } - - scalar maxDist = gMax - ( - mag - ( - solidPatchFaceCentres - - solidMesh().boundaryMesh()[solidPatchIndex_].faceCentres() - ) - ); - - Info << "Fluid-to-solid face interpolation error: " << maxDist - << endl; - } - - Info << "Checking solid-to-fluid point interpolator (GGI)" << endl; - { - vectorField solidZonePoints_ = - currentSolidZonePoints(); -// solidMesh().faceZones()[solidZoneIndex_]().localPoints(); - - vectorField solidZonePoints = - ggiInterpolatorPtr_->slaveToMasterPointInterpolate - ( - solidZonePoints_ - ); - - vectorField fluidZonePoints = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - scalar maxDist = gMax - ( - mag - ( - fluidZonePoints - - solidZonePoints - ) - ); - - Info << "Solid-to-fluid point interpolation error (GGI): " << maxDist - << endl; - } - - Info << "Number of uncovered master faces: " - << ggiInterpolatorPtr_->uncoveredMasterFaces().size() << endl; - - Info << "Number of uncovered slave faces: " - << ggiInterpolatorPtr_ ->uncoveredSlaveFaces().size() << endl; - - ggiInterpolatorPtr_->slavePointDistanceToIntersection(); - ggiInterpolatorPtr_->masterPointDistanceToIntersection(); -} - - -void Foam::fluidSolidInterface::calcSolidToFluidInterpolator() const -{ - // Find global face zones - if (solidToFluidPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcSolidToFluidInterpolator() const" - ) - << "Solid to fluid interpolator already exists" - << abort(FatalError); - } - - std::shared_ptr rbfFunction; - rbfFunction = std::shared_ptr(new TPSFunction()); - - // std::shared_ptr rbf; - solidToFluidPtr_ = - std::shared_ptr - ( - new RBFInterpolation(rbfFunction) - // new RBFInterpolation(rbfFunctionPtr_) - ); - - vectorField solidZonePoints = - currentSolidZonePatch().localPoints(); - - vectorField fluidZonePoints = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - matrix fluidX(fluidZonePoints.size(), 3); - matrix solidX(solidZonePoints.size(), 3); - - forAll(fluidZonePoints, faceI) - { - fluidX(faceI, 0) = fluidZonePoints[faceI].x(); - fluidX(faceI, 1) = fluidZonePoints[faceI].y(); - fluidX(faceI, 2) = fluidZonePoints[faceI].z(); - } - - forAll(solidZonePoints, faceI) - { - solidX(faceI, 0) = solidZonePoints[faceI].x(); - solidX(faceI, 1) = solidZonePoints[faceI].y(); - solidX(faceI, 2) = solidZonePoints[faceI].z(); - } - - solidToFluidPtr_->compute(solidX, fluidX); - - Info << "Checking solid-to-fluid interpolator" << endl; - { - matrix fluidPoints(fluidZonePoints.size(), 3); - matrix solidPoints(solidZonePoints.size(), 3); - vectorField fluidZonePointsInterp(fluidZonePoints.size(), vector::zero); - - - forAll(solidZonePoints, faceI) - { - solidPoints(faceI, 0) = solidZonePoints[faceI].x(); - solidPoints(faceI, 1) = solidZonePoints[faceI].y(); - solidPoints(faceI, 2) = solidZonePoints[faceI].z(); - } - - solidToFluidPtr_->interpolate(solidPoints, fluidPoints); - - forAll(fluidZonePoints, faceI) - { - fluidZonePointsInterp[faceI].x() = fluidPoints(faceI, 0); - fluidZonePointsInterp[faceI].y() = fluidPoints(faceI, 1); - fluidZonePointsInterp[faceI].z() = fluidPoints(faceI, 2); - } - - scalar maxDist = gMax - ( - mag - ( - fluidZonePointsInterp - - fluidZonePoints - ) - ); - - Info << "Solid-to-fluid point interpolation error: " << maxDist - << endl; - } -} - - -void Foam::fluidSolidInterface:: -calcAccumulatedFluidInterfaceDisplacement() const -{ - // Read accumulated displacement - if (accumulatedFluidInterfaceDisplacementPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcAccumulatedFluidInterfaceDisplacement() const" - ) - << "Accumulated displacement field already exists" - << abort(FatalError); - } - - // Accumulated fluid interface displacement - IOobject accumulatedFluidInterfaceDisplacementHeader - ( - "accumulatedFluidInterfaceDisplacement", - flow().runTime().timeName(), - fluidMesh(), - IOobject::MUST_READ - ); - - if(accumulatedFluidInterfaceDisplacementHeader.headerOk()) - { - Pout << "Reading accumulated fluid interface displacement" << endl; - - accumulatedFluidInterfaceDisplacementPtr_ = - new vectorIOField - ( - IOobject - ( - "accumulatedFluidInterfaceDisplacement", - flow().runTime().timeName(), - fluidMesh(), - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ) - ); - } - else - { - Pout << "Creating accumulated fluid interface displacement" << endl; - - accumulatedFluidInterfaceDisplacementPtr_ = - new vectorIOField - ( - IOobject - ( - "accumulatedFluidInterfaceDisplacement", - flow().runTime().timeName(), - fluidMesh(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - vectorField - ( - fluidMesh().boundaryMesh()[fluidPatchIndex()].nPoints(), - vector::zero - ) - ); - } -} - - -void Foam::fluidSolidInterface::calcMinEdgeLength() const -{ - // Read accumulated displacement - if (minEdgeLengthPtr_) - { - FatalErrorIn - ( - "void fluidSolidInterface::" - "calcMinEdgeLength() const" - ) - << "Minimal edge lengths already exist" - << abort(FatalError); - } - - minEdgeLengthPtr_ = - new scalarField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - 0 - ); - scalarField& minEdgeLength = *minEdgeLengthPtr_; - - - const edgeList& edges = - fluidMesh().faceZones()[fluidZoneIndex_]().edges(); - - const vectorField& points = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const labelListList& pointEdges = - fluidMesh().faceZones()[fluidZoneIndex_]().pointEdges(); - - forAll(points, pointI) - { - const labelList& curPointEdges = pointEdges[pointI]; - - scalar minLength = GREAT; - - forAll(curPointEdges, edgeI) - { - const edge& curEdge = edges[curPointEdges[edgeI]]; - - scalar Le = curEdge.mag(points); - - if (Le < minLength) - { - minLength = Le; - } - } - - minEdgeLength[pointI] = minLength; - } - -// Pout << "Min edge length: " << min(minEdgeLength) << endl; -// Pout << "gMin edge length: " << gMin(minEdgeLength) << endl; -} - -Foam::vectorIOField& -Foam::fluidSolidInterface::accumulatedFluidInterfaceDisplacement() -{ - if (!accumulatedFluidInterfaceDisplacementPtr_) - { - calcAccumulatedFluidInterfaceDisplacement(); - } - - return *accumulatedFluidInterfaceDisplacementPtr_; -} - - -const Foam::scalarField& Foam::fluidSolidInterface::minEdgeLength() const -{ - if (!minEdgeLengthPtr_) - { - calcMinEdgeLength(); - } - - return *minEdgeLengthPtr_; -} - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::fluidSolidInterface::fluidSolidInterface -( - dynamicFvMesh& fMesh, - fvMesh& sMesh -) -: - IOdictionary - ( - IOobject - ( - "fsiProperties", - fMesh.time().constant(), - fMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ), - fluidMesh_(fMesh), - flow_(fluidSolver::New(fluidMesh_)), - solidMesh_(sMesh), - stress_(solidSolver::New(solidMesh_)), - solidPatchIndex_(-1), - solidZoneIndex_(-1), - fluidPatchIndex_(-1), - fluidZoneIndex_(-1), - currentSolidZonePointsPtr_(NULL), - currentSolidZonePatchPtr_(NULL), - fluidToSolidPtr_(NULL), -// ggiFluidToSolidPtr_(NULL), - ggiInterpolatorPtr_(NULL), - solidToFluidPtr_(NULL), - couplingScheme_(lookup("couplingScheme")), - relaxationFactor_(readScalar(lookup("relaxationFactor"))), - aitkenRelaxationFactor_(relaxationFactor_), - outerCorrTolerance_(readScalar(lookup("outerCorrTolerance"))), - nOuterCorr_(readInt(lookup("nOuterCorr"))), - coupled_(lookup("coupled")), - predictor_(false), //(lookup("predictor")), - rbfInterpolation_(false), //(lookup("predictor")), - couplingReuse_(readInt(lookup("couplingReuse"))), - interfaceDeformationLimit_ - ( - readScalar(lookup("interfaceDeformationLimit")) - ), - fluidZonePointsDispl_(), - fluidZonePointsDisplRef_(), - fluidZonePointsDisplPrev_(), - solidZonePointsDispl_(), - solidZonePointsDisplRef_(), - interfacePointsDispl_(), - interfacePointsDisplPrev_(), - solidZonePressure_(), - solidZoneTraction_ - ( - IOobject - ( - "solidZoneTraction", - fMesh.time().timeName(), - fMesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - vectorField() - ), - solidZoneTractionPrev_(), - predictedSolidZoneTraction_ - ( - IOobject - ( - "predictedSolidZoneTraction", - fMesh.time().timeName(), - fMesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - vectorField() - ), - residual_(), - residualPrev_(), - maxResidualNorm_(0), - maxIntDisplNorm_(0), - outerCorr_(0), -// closedFluidDomain_(lookup("closedFluidDomain")), -// refPressure_(0), -// refPressureIncrement_(0), -// compressibility_(readScalar(lookup("compressibility"))), - interpolatorUpdateFrequency_ - ( - readInt(lookup("interpolatorUpdateFrequency")) - ), - fluidPatchPointsV_(), - fluidPatchPointsW_(), - fluidPatchPointsT_(), - accumulatedFluidInterfaceDisplacementPtr_(NULL), - minEdgeLengthPtr_(NULL) -{ - // Solid patch index - - word solidPatchName(lookup("solidPatch")); - - polyPatchID solidPatch(solidPatchName, solidMesh().boundaryMesh()); - - if (!solidPatch.active()) - { - FatalErrorIn("fluidSolidInterface::fluidSolidInterface(...)") - << "Solid patch name " << solidPatchName << " not found." - << abort(FatalError); - } - - solidPatchIndex_ = solidPatch.index(); - - - // Solid face zone index - - word solidZoneName(lookup("solidZone")); - - faceZoneID solidZone - ( - solidZoneName, - solidMesh().faceZones() - ); - - if (!solidZone.active()) - { - FatalErrorIn("") - << "Solid face zone name " << solidZoneName - << " not found. Please check your face zone definition." - << abort(FatalError); - } - - solidZoneIndex_ = solidZone.index(); - - - // Fluid patch index - - word fluidPatchName(lookup("fluidPatch")); - - polyPatchID fluidPatch(fluidPatchName, fluidMesh().boundaryMesh()); - - if (!fluidPatch.active()) - { - FatalErrorIn("fluidSolidInterface::fluidSolidInterface(...)") - << "Fluid patch name " << fluidPatchName << " not found." - << abort(FatalError); - } - - fluidPatchIndex_ = fluidPatch.index(); - - - // Fluid face zone index - - word fluidZoneName(lookup("fluidZone")); - - faceZoneID fluidZone - ( - fluidZoneName, - fluidMesh().faceZones() - ); - - if (!fluidZone.active()) - { - FatalErrorIn("fluidSolidInterface::fluidSolidInterface(...)") - << "Fluid face zone name " << fluidZoneName - << " not found. Please check your face zone definition." - << abort(FatalError); - } - - fluidZoneIndex_ = fluidZone.index(); - - - // Check coupling scheme - if - ( - (couplingScheme_ == "IQN-ILS") - || (couplingScheme_ == "Aitken") - || (couplingScheme_ == "FixedRelaxation") - ) - { - Info<< "Selecting coupling scheme " << couplingScheme_ << endl; - } - else - { - FatalErrorIn - ( - "fluidSolidInterface::fluidSolidInterface(...)" - ) << "couplingScheme: " << couplingScheme_ - << " is not a valid choice. " - << "Options are: IQN-ILS, Aitken, FixedRelaxation" - << abort(FatalError); - } - - // Initialize solid zone pressure - solidZonePressure_ = - scalarField(solidMesh().faceZones()[solidZoneIndex()].size(), 0.0); - - if (solidZoneTraction_.size() == 0) - { - solidZoneTraction_ = - vectorField - ( - solidMesh().faceZones()[solidZoneIndex_]().size(), - vector::zero - ); - } - - solidZoneTractionPrev_ = - vectorField - ( - solidMesh().faceZones()[solidZoneIndex_]().size(), - vector::zero - ); - - if (predictedSolidZoneTraction_.size() == 0) - { - predictedSolidZoneTraction_ = - vectorField - ( - solidMesh().faceZones()[solidZoneIndex_]().size(), - vector::zero - ); - } - - residual_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - if (found("predictor")) - { - predictor_ = Switch(lookup("predictor")); - } - - if (found("rbfInterpolation")) - { - rbfInterpolation_ = Switch(lookup("rbfInterpolation")); - } -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::fluidSolidInterface::~fluidSolidInterface() -{ - deleteDemandDrivenData(currentSolidZonePointsPtr_); - deleteDemandDrivenData(currentSolidZonePatchPtr_); - // deleteDemandDrivenData(fluidToSolidPtr_); -// deleteDemandDrivenData(ggiFluidToSolidPtr_); - deleteDemandDrivenData(ggiInterpolatorPtr_); - // deleteDemandDrivenData(solidToFluidPtr_); - deleteDemandDrivenData(accumulatedFluidInterfaceDisplacementPtr_); - deleteDemandDrivenData(minEdgeLengthPtr_); -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -const Foam::vectorField& -Foam::fluidSolidInterface::currentSolidZonePoints() const -{ - if (!currentSolidZonePointsPtr_) - { - calcCurrentSolidZonePoints(); - } - - return *currentSolidZonePointsPtr_; -} - -const Foam::PrimitivePatch& -Foam::fluidSolidInterface::currentSolidZonePatch() const -{ - if (!currentSolidZonePatchPtr_) - { - calcCurrentSolidZonePatch(); - } - - return *currentSolidZonePatchPtr_; -} - -const std::shared_ptr& -Foam::fluidSolidInterface::fluidToSolid() const -{ - if (!fluidToSolidPtr_) - { - calcFluidToSolidInterpolator(); - } - - return fluidToSolidPtr_; -} - -// const Foam::ggiZoneInterpolation& -// Foam::fluidSolidInterface::ggiFluidToSolid() const -// { -// if (!ggiFluidToSolidPtr_) -// { -// calcGgiFluidToSolidInterpolator(); -// } - -// return *ggiFluidToSolidPtr_; -// } - -const Foam::ggiZoneInterpolation& -Foam::fluidSolidInterface::ggiInterpolator() const -{ - if (!ggiInterpolatorPtr_) - { - calcGgiInterpolator(); - } - - return *ggiInterpolatorPtr_; -} - -const std::shared_ptr& -Foam::fluidSolidInterface::solidToFluid() const -{ - if (!solidToFluidPtr_) - { - calcSolidToFluidInterpolator(); - } - - return solidToFluidPtr_; -} - -void Foam::fluidSolidInterface::initializeFields() -{ - fluidZonePointsDispl_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - fluidZonePointsDisplRef_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - fluidZonePointsDisplPrev_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - solidZonePointsDispl_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - solidZonePointsDisplRef_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - residualPrev_ = residual_; -// vectorField -// ( -// fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), -// vector::zero -// ); - - residual_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - maxResidualNorm_ = 0; - - outerCorr_ = 0; - - nOuterCorr_ = readInt(lookup("nOuterCorr")); - - outerCorrTolerance_ = readScalar(lookup("outerCorrTolerance")); - - coupled_ = Switch(lookup("coupled")); - - couplingReuse_ = readInt(lookup("couplingReuse")); - - relaxationFactor_ = readScalar(lookup("relaxationFactor")); - - - interfacePointsDispl_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - - interfacePointsDisplPrev_ = - vectorField - ( - fluidMesh().faceZones()[fluidZoneIndex_]().nPoints(), - vector::zero - ); - -// refPressure_ += refPressureIncrement_; - -// if (timeIndex_ < runTime().timeIndex()) -// { -// timeIndex_ = runTime().timeIndex(); - -// } -} - - -void Foam::fluidSolidInterface::updateInterpolator() -{ -// label interpolatorUpdateFrequency_ = 2; - -// Info << runTime().timeIndex() << endl; - - if (!ggiInterpolatorPtr_) - { - ggiInterpolator(); - } - else if (interpolatorUpdateFrequency_ != 0) - { - if (((runTime().timeIndex()-1)%interpolatorUpdateFrequency_) == 0) - { - deleteDemandDrivenData(ggiInterpolatorPtr_); - ggiInterpolator(); - } - } -// else -// { -// if ((runTime().timeIndex()-1) == 0) -// { -// deleteDemandDrivenData(ggiInterpolatorPtr_); -// ggiInterpolator(); -// } -// } -} - - -void Foam::fluidSolidInterface::updateDisplacement() -{ - Info << "\nTime = " << flow().runTime().timeName() - << ", iteration: " << outerCorr() << endl; - -// if (outerCorr_ == 1) -// { -// // Cancel residual from previous time step -// fluidZonePointsDisplPrev() = fluidZonePointsDispl(); -// fluidZonePointsDispl() += residualPrev(); -// } - - - if (couplingScheme() == "FixedRelaxation") - { - Info << "Current fsi under-relaxation factor: " - << relaxationFactor() << endl; - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - // if (outerCorr_ == 1) - // { - // // Cancel residual from previous time step - // fluidZonePointsDispl() += residualPrev(); - // } - - fluidZonePointsDispl() += relaxationFactor()*residual(); - } - else if (couplingScheme() == "Aitken") - { - if (outerCorr_ < 3) - { - Info << "Current fsi under-relaxation factor: " - << relaxationFactor() << endl; - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - // if (outerCorr_ == 1) - // { - // // Cancel residual from previous time step - // fluidZonePointsDispl() += residualPrev(); - // } - - if ((outerCorr_ == 1) && predictor()) - { - fluidZonePointsDispl() += residual(); - } - else - { - fluidZonePointsDispl() += relaxationFactor()*residual(); - } -// fluidZonePointsDispl() += relaxationFactor()*residual(); - } - else - { - aitkenRelaxationFactor() = - -aitkenRelaxationFactor() - *( - sum - ( - residualPrev() - & (residual() - residualPrev()) - ) - /( - sum - ( - (residual() - residualPrev()) - & (residual() - residualPrev()) - ) - ) - ); - - if (Pstream::parRun()) - { - if(!Pstream::master()) - { - aitkenRelaxationFactor() = 0.0; - } - - //- pass to all procs - reduce(aitkenRelaxationFactor(), sumOp()); - } - - aitkenRelaxationFactor() = mag(aitkenRelaxationFactor()); - - if (aitkenRelaxationFactor()>1) - { - aitkenRelaxationFactor() = relaxationFactor(); - } - - Info << "Current fsi under-relaxation factor (Aitken): " - << aitkenRelaxationFactor() << endl; - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - fluidZonePointsDispl() += - aitkenRelaxationFactor()*residual(); - } - } - else if (couplingScheme() == "IQN-ILS") - { -// A fluid structure interaction solver with IQN-ILS -// coupling algorithm (J. Degroote, K.-J. Bathe and J. Vierendeels. -// Performance of a new partitioned procedure versus a monolithic -// procedure in fluid-structure interaction. Computers & Structures - - // IQN-ILS - if (outerCorr_ == 1) - { - // Clean up data from old time steps - - Info << "Modes before clean-up : " << fluidPatchPointsT_.size(); - - while (true) - { - if (fluidPatchPointsT_.size()) - { - if - ( - flow().runTime().timeIndex()-couplingReuse() - > fluidPatchPointsT_[0] - ) - { - for (label i = 0; i < fluidPatchPointsT_.size()-1; i++) - { - fluidPatchPointsT_[i] = fluidPatchPointsT_[i+1]; - fluidPatchPointsV_[i] = fluidPatchPointsV_[i+1]; - fluidPatchPointsW_[i] = fluidPatchPointsW_[i+1]; - } - fluidPatchPointsT_.remove(); - fluidPatchPointsV_.remove(); - fluidPatchPointsW_.remove(); - } - else - { - break; - } - } - else - { - break; - } - } - - Info << ", modes after clean-up : " - << fluidPatchPointsT_.size() << endl; - } - else if (outerCorr_ == 2) - { - // Set reference in the first coupling iteration - solidZonePointsDisplRef() = solidZonePointsDispl(); - fluidZonePointsDisplRef() = fluidZonePointsDispl(); - } - else - { - // Reference has been set in the first coupling iteration - fluidPatchPointsV_.append - ( - ( - solidZonePointsDispl() - - fluidZonePointsDispl() - ) - - ( - solidZonePointsDisplRef() - - fluidZonePointsDisplRef() - ) - ); - - fluidPatchPointsW_.append - ( - solidZonePointsDispl() - - solidZonePointsDisplRef() - ); - - fluidPatchPointsT_.append - ( - flow().runTime().timeIndex() - ); - } - - if (fluidPatchPointsT_.size() > 1) - { - updateDisplacementUsingIQNILS(); - } - else - { - // Relax the interface displacement - Info << "Current fsi under-relaxation factor: " - << relaxationFactor() << endl; - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - // if (outerCorr_ == 1) - // { - // // Cancel residual from previous time step - // fluidZonePointsDispl() += residualPrev(); - // } - - if ((outerCorr_ == 1) && predictor()) - { - fluidZonePointsDispl() += residual(); - } - else - { - fluidZonePointsDispl() += relaxationFactor()*residual(); - } - } - } - - - // Set interface acceleration - if - ( - isA - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ) - { - Info << "Setting acceleration at fluid side of the interface" << endl; - - vectorField solidZoneAcceleration = - stress().faceZoneAcceleration - ( - solidZoneIndex(), - solidPatchIndex() - ); - - vectorField fluidZoneAcceleration = - ggiInterpolator().slaveToMaster - ( - solidZoneAcceleration - ); - - vectorField fluidPatchAcceleration - ( - fluidMesh().boundary()[fluidPatchIndex()].size(), - vector::zero - ); - - const label patchStart = - fluidMesh().boundaryMesh()[fluidPatchIndex()].start(); - - forAll(fluidPatchAcceleration, i) - { - fluidPatchAcceleration[i] = - fluidZoneAcceleration - [ - fluidMesh().faceZones()[fluidZoneIndex()] - .whichFace(patchStart + i) - ]; - } - - const_cast - ( - refCast - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ).prevAcceleration() = fluidPatchAcceleration; - } - - // Make sure that displacement on all processors is equal to one - // calculated on master processor - if (Pstream::parRun()) - { - if(!Pstream::master()) - { - fluidZonePointsDispl() *= 0.0; - } - - //- pass to all procs - reduce(fluidZonePointsDispl(), sumOp()); - - label globalFluidZoneIndex = - findIndex(flow().globalFaceZones(), fluidZoneIndex()); - - if (globalFluidZoneIndex == -1) - { - FatalErrorIn - ( - "fluidSolidInterface::updateDisplacement()" - ) << "global zone point map is not availabel" - << abort(FatalError); - } - - const labelList& map = - flow().globalToLocalFaceZonePointMap()[globalFluidZoneIndex]; - - if(!Pstream::master()) - { - vectorField fluidZonePointsDisplGlobal = - fluidZonePointsDispl(); - - forAll(fluidZonePointsDisplGlobal, globalPointI) - { - label localPoint = map[globalPointI]; - - fluidZonePointsDispl()[localPoint] = - fluidZonePointsDisplGlobal[globalPointI]; - } - } - } -} - - -Foam::scalar Foam::fluidSolidInterface::updateWeakDisplacement() -{ - vectorField solidZonePointsDisplAtSolid = - stress().faceZonePointDisplacementIncrement(solidZoneIndex()); - - solidZonePointsDispl() = - ggiInterpolator().slaveToMasterPointInterpolate - ( - solidZonePointsDisplAtSolid - ); - - vectorField solidZonePointsTotDisplAtSolid = - stress().faceZonePointDisplacement(solidZoneIndex()); - - vectorField solidZonePointsTotDispl = - ggiInterpolator().slaveToMasterPointInterpolate - ( - solidZonePointsTotDisplAtSolid - ); - - - residualPrev() = residual(); - - residual() = solidZonePointsDispl() - fluidZonePointsDispl(); - -// Info << "residual: " << average(mag(residual())) << endl; - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - fluidZonePointsDispl() += residual(); - - // Make sure that displacement on all processors is equal to one - // calculated on master processor - if (Pstream::parRun()) - { - if(!Pstream::master()) - { - fluidZonePointsDispl() *= 0.0; - } - - //- pass to all procs - reduce(fluidZonePointsDispl(), sumOp()); - - label globalFluidZoneIndex = - findIndex(flow().globalFaceZones(), fluidZoneIndex()); - - if (globalFluidZoneIndex == -1) - { - FatalErrorIn - ( - "fluidSolidInterface::updateDisplacement()" - ) << "global zone point map is not availabel" - << abort(FatalError); - } - - const labelList& map = - flow().globalToLocalFaceZonePointMap()[globalFluidZoneIndex]; - - if(!Pstream::master()) - { - vectorField fluidZonePointsDisplGlobal = - fluidZonePointsDispl(); - - forAll(fluidZonePointsDisplGlobal, globalPointI) - { - label localPoint = map[globalPointI]; - - fluidZonePointsDispl()[localPoint] = - fluidZonePointsDisplGlobal[globalPointI]; - } - } - } - - - // Set interface acceleration - if - ( - isA - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ) - { - vectorField solidZoneAcceleration = - stress().faceZoneAcceleration - ( - solidZoneIndex(), - solidPatchIndex() - ); - - vectorField fluidZoneAcceleration = - ggiInterpolator().slaveToMaster - ( - solidZoneAcceleration - ); - - vectorField fluidPatchAcceleration - ( - fluidMesh().boundary()[fluidPatchIndex()].size(), - vector::zero - ); - - const label patchStart = - fluidMesh().boundaryMesh()[fluidPatchIndex()].start(); - - forAll(fluidPatchAcceleration, i) - { - fluidPatchAcceleration[i] = - fluidZoneAcceleration - [ - fluidMesh().faceZones()[fluidZoneIndex()] - .whichFace(patchStart + i) - ]; - } - - const_cast - ( - refCast - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ).prevAcceleration() = fluidPatchAcceleration; - } - - - // Calculate residual norm - - scalar residualNorm = ::sqrt(gSum(magSqr(residual()))); - scalar residualNorm_2 = residualNorm; - - if (residualNorm > maxResidualNorm_) - { - maxResidualNorm_ = residualNorm; - } - - residualNorm /= maxResidualNorm_ + SMALL; - - Info << "Current fsi relative residual norm: " << residualNorm << endl; - - interfacePointsDisplPrev_ = interfacePointsDispl_; - - interfacePointsDispl_ = solidZonePointsDispl(); - - vectorField intTotDispl = - interfacePointsDispl_ + solidZonePointsTotDispl; - scalar intTotDisplNorm = ::sqrt(gSum(magSqr(intTotDispl))); - if (intTotDisplNorm > maxIntDisplNorm_) - { - maxIntDisplNorm_ = intTotDisplNorm; - } - - residualNorm_2 /= maxIntDisplNorm_ + SMALL; - - Info << "Alternative fsi residual: " << residualNorm_2 << endl; - - return min(residualNorm_2, residualNorm); -// return residualNorm; -} - - -void Foam::fluidSolidInterface::updateDisplacementUsingIQNILS() -{ - // Consider fluidPatchPointsV as a matrix V - // with as columns the items - // in the DynamicList and calculate the QR-decomposition of V - // with modified Gram-Schmidt - label cols = fluidPatchPointsV_.size(); - RectangularMatrix R(cols, cols, 0.0); - RectangularMatrix C(cols, 1); - RectangularMatrix Rcolsum(1, cols); - DynamicList Q; - - for (label i = 0; i < cols; i++) - { - Q.append(fluidPatchPointsV_[cols-1-i]); - } - - for (label i = 0; i < cols; i++) - { - // Normalize column i - R[i][i] = Foam::sqrt(sum(Q[i] & Q[i])); - Q[i] /= R[i][i]; - - // Orthogonalize columns to the right of column i - for (label j = i+1; j < cols; j++) - { - R[i][j] = sum(Q[i] & Q[j]); - Q[j] -= R[i][j]*Q[i]; - } - - // Project minus the residual vector on the Q - C[i][0] = sum - ( - Q[i] - & ( - fluidZonePointsDispl() - - solidZonePointsDispl() - ) - ); - } - - // Solve the upper triangular system - for (label j = 0; j < cols; j++) - { - Rcolsum[0][j] = 0.0; - for (label i = 0; i < j+1; i++) - { - Rcolsum[0][j] += cmptMag(R[i][j]); - } - } - scalar epsilon = 1.0E-10*max(Rcolsum); - for (label i = 0; i < cols; i++) - { - if (cmptMag(R[i][i]) > epsilon) - { - for (label j = i+1; j < cols; j++) - { - R[i][j] /= R[i][i]; - } - C[i][0] /= R[i][i]; - R[i][i] = 1.0; - } - } - for (label j = cols-1; j >= 0; j--) - { - if (cmptMag(R[j][j]) > epsilon) - { - for (label i = 0; i < j; i++) - { - C[i][0] -= C[j][0]*R[i][j]; - } - } - else - { - C[j][0] = 0.0; - } - } - - fluidZonePointsDisplPrev() = fluidZonePointsDispl(); - - fluidZonePointsDispl() = solidZonePointsDispl(); - - for (label i = 0; i < cols; i++) - { - fluidZonePointsDispl() += fluidPatchPointsW_[i]*C[cols-1-i][0]; - } -} - -void Foam::fluidSolidInterface::moveFluidMesh() -{ - // Get fluid patch displacement from fluid zone displacement - - vectorField fluidPatchPointsDispl - ( - fluidMesh().boundaryMesh()[fluidPatchIndex()].nPoints(), - vector::zero - ); - - vectorField fluidPatchPointsDisplPrev - ( - fluidMesh().boundaryMesh()[fluidPatchIndex()].nPoints(), - vector::zero - ); - - const labelList& fluidPatchMeshPoints = - fluidMesh().boundaryMesh()[fluidPatchIndex()].meshPoints(); - - forAll(fluidPatchPointsDispl, pointI) - { - label curMeshPointID = fluidPatchMeshPoints[pointI]; - - label curFluidZonePointID = - fluidMesh().faceZones()[fluidZoneIndex()]() - .whichPoint(curMeshPointID); - - fluidPatchPointsDispl[pointI] = - fluidZonePointsDispl()[curFluidZonePointID]; - - fluidPatchPointsDisplPrev[pointI] = - fluidZonePointsDisplPrev()[curFluidZonePointID]; - } - - // Move fluid mesh - const vectorField& n = - fluidMesh().boundaryMesh()[fluidPatchIndex()].pointNormals(); - - primitivePatchInterpolation patchInterpolator - ( - fluidMesh().boundaryMesh()[fluidPatchIndex()] - ); - - scalarField pointDeltaCoeffs = - patchInterpolator.faceToPointInterpolate - ( - fluidMesh().boundary()[fluidPatchIndex()].deltaCoeffs() - ); - - scalar delta = - gMax - ( - mag - ( - n - & ( - accumulatedFluidInterfaceDisplacement() - + fluidPatchPointsDispl - - fluidPatchPointsDisplPrev - ) - ) - *pointDeltaCoeffs - ); - - Info << "Maximal accumulated displacement of interface points: " - << delta << endl; - - if (false) -// if (delta < interfaceDeformationLimit()) - { - // Move only interface points - pointField newPoints = fluidMesh().allPoints(); - - const labelList& meshPoints = - fluidMesh().boundaryMesh()[fluidPatchIndex()].meshPoints(); - - forAll (fluidPatchPointsDispl, pointI) - { - newPoints[meshPoints[pointI]] += - fluidPatchPointsDispl[pointI] - - fluidPatchPointsDisplPrev[pointI]; - } - - twoDPointCorrector twoDCorrector(fluidMesh()); - - twoDCorrector.correctPoints(newPoints); - - fluidMesh_.movePoints(newPoints); - - // Accumulate interface points displacement - accumulatedFluidInterfaceDisplacement() += - fluidPatchPointsDispl - - fluidPatchPointsDisplPrev; - } - else - { -// // Move whole fluid mesh -// const pointField& oldAllPoints = fluidMesh().allPoints(); -// pointField newPoints = fluidMesh().allPoints(); - -// const labelList& meshPoints = -// fluidMesh().boundaryMesh()[fluidPatchIndex()].meshPoints(); - -// forAll (accumulatedFluidInterfaceDisplacement(), pointI) -// { -// newPoints[meshPoints[pointI]] -= -// accumulatedFluidInterfaceDisplacement()[pointI]; -// } - -// twoDPointCorrector twoDCorrector(fluidMesh()); - -// twoDCorrector.correctPoints(newPoints); - -// fluidMesh_.movePoints(newPoints); - -// accumulatedFluidInterfaceDisplacement() += -// fluidPatchPointsDispl -// - fluidPatchPointsDisplPrev; - - // Check mesh motion solver type - bool feMotionSolver = - fluidMesh().objectRegistry::foundObject - ( - "motionU" - ); - - bool fvMotionSolver = - fluidMesh().objectRegistry::foundObject - ( - "pointMotionU" - ); - - if (feMotionSolver) - { - tetPointVectorField& motionU = - const_cast - ( - fluidMesh().objectRegistry:: - lookupObject - ( - "motionU" - ) - ); - - fixedValueTetPolyPatchVectorField& motionUFluidPatch = - refCast - ( - motionU.boundaryField()[fluidPatchIndex()] - ); - - tetPolyPatchInterpolation tppi - ( - refCast(motionUFluidPatch.patch()) - ); - - motionUFluidPatch == - tppi.pointToPointInterpolate - ( - accumulatedFluidInterfaceDisplacement() - /flow().runTime().deltaT().value() - ); - } - else if (fvMotionSolver) - { - pointVectorField& motionU = - const_cast - ( - fluidMesh().objectRegistry:: - lookupObject - ( - "pointMotionU" - ) - ); - - fixedValuePointPatchVectorField& motionUFluidPatch = - refCast - ( - motionU.boundaryField()[fluidPatchIndex()] - ); - - motionUFluidPatch == - ( - fluidPatchPointsDispl - - fluidPatchPointsDisplPrev - )/flow().runTime().deltaT().value(); - } - else - { - FatalErrorIn("fluidSolidInterface::moveFluidMesh()") - << "Problem with fluid mesh motion solver selection" - << abort(FatalError); - } - - fluidMesh_.update(); - - accumulatedFluidInterfaceDisplacement() = - vectorField - ( - accumulatedFluidInterfaceDisplacement().size(), - vector::zero - ); - } - - - // Move unused fluid mesh points - { - vectorField newPoints = fluidMesh().allPoints(); - - const labelList& fluidZoneMeshPoints = - fluidMesh().faceZones()[fluidZoneIndex()]().meshPoints(); - - forAll(fluidZonePointsDispl(), pointI) - { - if (fluidZoneMeshPoints[pointI] >= fluidMesh().nPoints()) - { - newPoints[fluidZoneMeshPoints[pointI]] += - fluidZonePointsDispl()[pointI] - - fluidZonePointsDisplPrev()[pointI]; - } - } - - twoDPointCorrector twoDCorrector(fluidMesh()); - - twoDCorrector.correctPoints(newPoints); - - fluidMesh_.movePoints(newPoints); - } - -// // Evaluate interface velocity -// const_cast(flow().U()) -// .boundaryField()[fluidPatchIndex()].evaluate(); -} - - -void Foam::fluidSolidInterface::updateForce() -{ - Info << "Setting traction on solid patch" << endl; - - vectorField fluidZoneTraction = - flow().faceZoneViscousForce - ( - fluidZoneIndex(), - fluidPatchIndex() - ); - - // Info << "Viscous force, max: " << max(mag(fluidZoneTraction)) << endl; - - scalarField fluidZonePressure = - flow().faceZonePressureForce(fluidZoneIndex(), fluidPatchIndex()); - - - // Fluid zone face normals - - vectorField p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const labelList& mp = - fluidMesh().faceZones()[fluidZoneIndex_]().meshPoints(); - const vectorField& allPoints = fluidMesh().allPoints(); - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - vectorField n(f.size(), vector::zero); - forAll(n, faceI) - { - n[faceI] = f[faceI].normal(p); - n[faceI] /= mag(n[faceI]); - } - - // Fluid zone total traction - vectorField fluidZoneTotalTraction = - fluidZoneTraction - fluidZonePressure*n; - -// vectorField solidZoneTraction = -// ggiInterpolator().masterToSlave -// ( -// -fluidZoneTraction -// ); - - vectorField solidZoneTotalTraction - ( - solidMesh().faceZones()[solidZoneIndex()].size(), - vector::zero - ); - - if (rbfInterpolation_) - { - Info << "... using RBF interpolation" << endl; - - matrix fluidForce(fluidZoneTotalTraction.size(), 3); - matrix solidForce(solidZoneTotalTraction.size(), 3); - - forAll(fluidZoneTotalTraction, faceI) - { - fluidForce(faceI, 0) = fluidZoneTotalTraction[faceI].x(); - fluidForce(faceI, 1) = fluidZoneTotalTraction[faceI].y(); - fluidForce(faceI, 2) = fluidZoneTotalTraction[faceI].z(); - } - - fluidToSolid()->interpolate(fluidForce, solidForce); - - forAll(solidZoneTotalTraction, faceI) - { - solidZoneTotalTraction[faceI].x() = -solidForce(faceI, 0); - solidZoneTotalTraction[faceI].y() = -solidForce(faceI, 1); - solidZoneTotalTraction[faceI].z() = -solidForce(faceI, 2); - } - } - else - { - solidZoneTotalTraction = - ggiInterpolator().masterToSlave - ( - - fluidZoneTotalTraction - ); - } - - solidZonePressure_ = - ggiInterpolator().masterToSlave - ( - fluidZonePressure - ); - - if (false) - { - volVectorField fluidTraction - ( - IOobject - ( - "fluidTraction", - runTime().timeName(), - fluidMesh(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - fluidMesh(), - dimensionedVector("0", dimPressure, vector::zero) - ); - fluidTraction.boundaryField()[fluidPatchIndex_] = - fluidZoneTotalTraction; - fluidTraction.write(); - - volVectorField solidTraction - ( - IOobject - ( - "solidTraction", - runTime().timeName(), - solidMesh(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - solidMesh(), - dimensionedVector("0", dimPressure, vector::zero) - ); - solidTraction.boundaryField()[solidPatchIndex_] = - solidZoneTotalTraction; - solidTraction.write(); - } - - if (coupled()) - { -// stress().setPressure -// ( -// solidPatchIndex(), -// solidZoneIndex(), -// solidZonePressure_ -// ); - - stress().setTraction - ( - solidPatchIndex(), - solidZoneIndex(), - solidZoneTotalTraction - ); - - // Set interface pressure for elasticWallPressure - // boundary condition - if - ( - isA - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ) - { - const_cast - ( - refCast - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ).prevPressure() = flow().patchPressureForce(fluidPatchIndex_); - } - } - else - { - // Set interface pressure for elasticWallPressure - // boundary condition - if - ( - isA - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ) - { - const_cast - ( - refCast - ( - flow().p().boundaryField()[fluidPatchIndex_] - ) - ).prevPressure() = 0; - } - } - - // Total force at the fluid side of the interface - if (true) - { - vectorField p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const labelList& mp = - fluidMesh().faceZones()[fluidZoneIndex_]().meshPoints(); - const vectorField& allPoints = fluidMesh().allPoints(); - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); - - vectorField C(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); - C[faceI] = f[faceI].centre(p); - } - - vector totalTractionForce = - sum(fluidZoneTotalTraction*mag(S)); - -// vector totalPressureForce = sum(fluidZonePressure*S); - -// Info << setprecision(12); - Info << "Total force (fluid) = " - << totalTractionForce << endl; - - // Info << average(C) << ", " << sum(mag(S)) - // << ", " << sum(fluidZoneTotalTraction) << endl; - -// Info << "Total pressure force (fluid) = " -// << totalPressureForce << endl; - } - - // Totla force at the solid side of the interface - if (true) - { - vectorField p = - solidMesh().faceZones()[solidZoneIndex_]().localPoints(); - - const labelList& mp = - solidMesh().faceZones()[solidZoneIndex_]().meshPoints(); - const vectorField& allPoints = solidMesh().allPoints(); - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - solidMesh().faceZones()[solidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); - vectorField C(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); - C[faceI] = f[faceI].centre(p); - } - - vector totalTractionForce = - sum(solidZoneTotalTraction*mag(S)); - -// vector totalPressureForce = -// sum(solidZonePressure_*S); - -// Info << setprecision(12); - Info << "Total force (solid) = " - << totalTractionForce << endl; - - // Info << average(C) << ", " << sum(mag(S)) - // << ", " << sum(fluidZoneTotalTraction) << endl; - -// Info << sum(mag(S)) << ", " << sum(fluidZoneTotalTraction) << endl; - -// Info << "Total pressure force (solid) = " -// << totalPressureForce << endl; - } -} - - -void Foam::fluidSolidInterface::updateWeakForce() -{ - Info << "Setting weak traction on solid patch" << endl; - - vectorField fluidZoneTraction = - flow().faceZoneViscousForce - ( - fluidZoneIndex(), - fluidPatchIndex() - ); - - scalarField fluidZonePressure = - flow().faceZonePressureForce - ( - fluidZoneIndex(), - fluidPatchIndex() - ); - - - // Calculate current interface zone face normals - // (this is maybe not necessary but I want to be sure - // that normal is correct for whole zone) - vectorField n(fluidZonePressure.size(), vector::zero); - { - vectorField p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const labelList& mp = - fluidMesh().faceZones()[fluidZoneIndex_]().meshPoints(); - - const vectorField& allPoints = fluidMesh().allPoints(); - - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - forAll(n, faceI) - { - n[faceI] = f[faceI].normal(p); - n[faceI] /= mag(n[faceI]); - } - } - - - // Velocity at fluid and solid side of the interface - // needed for kinamtically coupled beta scheme - -// // After fluid model solution (Robin bc) -// vectorField fluidZoneVelocity = -// flow().faceZoneVelocity -// ( -// fluidZoneIndex(), -// fluidPatchIndex() -// ); - -// // Velocity calculated by solving structure model -// // in previous time step -// vectorField solidZoneVelocity = -// stress().faceZoneVelocity -// ( -// solidZoneIndex(), -// solidPatchIndex() -// ); - -// vectorField solidZoneVelocityAtFluid = -// ggiInterpolator().slaveToMaster -// ( -// solidZoneVelocity -// ); - - -// // Pressure force correction according to -// // kinematically coupled beta scheme -// scalar rhoSolid = stress().rheology().rho()()[0]; -// scalar mu = stress().rheology().mu()()[0]; -// scalar lambda = stress().rheology().lambda()()[0]; -// scalar ap = sqrt((lambda+2*mu)/rhoSolid); -// scalar hSolid = ap*flow().runTime().deltaT().value(); - -// Info << hSolid << endl; - -// // Fluid interface velocity from previous step -// vectorField fluidZoneVelocity = -// flow().U().boundaryField()[fluidPatchIndex_]; - -// vectorField fluidZoneVelocityOld = -// flow().U().oldTime().boundaryField()[fluidPatchIndex_]; - -// scalarField fluidZonePressureCorr = -// rhoSolid*hSolid* -// ((fluidZoneVelocity - fluidZoneVelocityOld)&n) -// /flow().runTime().deltaT().value(); -// scalarField fluidZonePressureCorr = -// rhoSolid*hSolid* -// ((fluidZoneVelocity - solidZoneVelocityAtFluid)&n) -// /flow().runTime().deltaT().value(); - - - - // Fluid zone total traction - -// Info << max(fluidZonePressureCorr) << ", " -// << average(fluidZonePressureCorr) << endl; - - vectorField fluidZoneTotalTraction = - fluidZoneTraction - fluidZonePressure*n; - - vectorField solidZoneTotalTraction = - ggiInterpolator().masterToSlave - ( - -fluidZoneTotalTraction - ); - -// vectorField solidZoneVelocity = -// ggiInterpolator().masterToSlave -// ( -// fluidZoneVelocity -// ); - - vectorField solidZoneNormal = - ggiInterpolator().masterToSlave(-n); - - solidZonePressure_ = - ggiInterpolator().masterToSlave - ( - fluidZonePressure - ); - - if (coupled()) - { -// stress().setVelocityAndTraction -// ( -// solidPatchIndex(), -// solidZoneIndex(), -// solidZoneTotalTraction, -// solidZoneVelocity, -// solidZoneNormal -// ); - - stress().setTraction - ( - solidPatchIndex(), - solidZoneIndex(), - solidZoneTotalTraction - ); - } - - // Total force at the fluid side of the interface - { - vectorField p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const labelList& mp = - fluidMesh().faceZones()[fluidZoneIndex_]().meshPoints(); - const vectorField& allPoints = fluidMesh().allPoints(); - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); -// vectorField C(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); -// C[faceI] = f[faceI].centre(p); - } - - vector totalTractionForce = sum(fluidZoneTotalTraction*mag(S)); - -// Info << setprecision(12); - Info << "Total force (fluid) = " - << totalTractionForce << endl; - } - - // Totla force at the solid side of the interface - { - vectorField p = - solidMesh().faceZones()[solidZoneIndex_]().localPoints(); - - const labelList& mp = - solidMesh().faceZones()[solidZoneIndex_]().meshPoints(); - const vectorField& allPoints = solidMesh().allPoints(); - forAll(mp, pI) - { - p[pI] = allPoints[mp[pI]]; - } - - const faceList& f = - solidMesh().faceZones()[solidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); -// vectorField C(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); -// C[faceI] = f[faceI].centre(p); - } - - vector totalTractionForce = - sum(solidZoneTotalTraction*mag(S)); - - // Info << setprecision(12); - Info << "Total force (solid) = " - << totalTractionForce << endl; - } -} - - -void Foam::fluidSolidInterface::updateWeakTraction() -{ - Info << "Update weak traction on solid patch" << endl; - - // Calc fluid traction - - const vectorField& p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - vectorField n(f.size(), vector::zero); - forAll(n, faceI) - { - n[faceI] = f[faceI].normal(p); - n[faceI] /= mag(n[faceI]); - } - - vectorField fluidZoneTraction = - flow().faceZoneViscousForce - ( - fluidZoneIndex(), - fluidPatchIndex() - ) - - flow().faceZonePressureForce(fluidZoneIndex(), fluidPatchIndex())*n; - - vectorField fluidZoneTractionAtSolid = - ggiInterpolator().masterToSlave - ( - -fluidZoneTraction - ); - -// scalar beta_ = 1.0; - scalar beta_ = relaxationFactor_; - - solidZoneTractionPrev_ = solidZoneTraction_; - - solidZoneTraction_ = - beta_*fluidZoneTractionAtSolid - + (1.0-beta_)*predictedSolidZoneTraction_; - - predictedSolidZoneTraction_ = - 2*solidZoneTraction_ - solidZoneTractionPrev_; - - if (coupled()) - { - Info << "Setting weak traction on solid patch" << endl; - - stress().setTraction - ( - solidPatchIndex(), - solidZoneIndex(), - predictedSolidZoneTraction_ - ); - } - - // Total force at the fluid side of the interface - { - const vectorField& p = - fluidMesh().faceZones()[fluidZoneIndex_]().localPoints(); - - const faceList& f = - fluidMesh().faceZones()[fluidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); - } - - vector totalTractionForce = sum(fluidZoneTraction*mag(S)); - - Info << "Total force (fluid) = " - << totalTractionForce << endl; - } - - // Totla force at the solid side of the interface - { - const vectorField& p = - solidMesh().faceZones()[solidZoneIndex_]().localPoints(); - - const faceList& f = - solidMesh().faceZones()[solidZoneIndex_]().localFaces(); - - vectorField S(f.size(), vector::zero); - - forAll(S, faceI) - { - S[faceI] = f[faceI].normal(p); - } - - vector totalTractionForce = - sum(fluidZoneTractionAtSolid*mag(S)); - - Info << "Total force (solid) = " - << totalTractionForce << endl; - } -} - - -void Foam::fluidSolidInterface::predictAndUpdateForce() -{ - if (coupled()) - { - Info << "Setting traction on solid patch using prediction" << endl; - - stress().setPressure - ( - solidPatchIndex(), - solidZoneIndex(), - stress().predictPressure(solidPatchIndex(), solidZoneIndex()) - ); - - stress().setTraction - ( - solidPatchIndex(), - solidZoneIndex(), - stress().predictTraction(solidPatchIndex(), solidZoneIndex()) - ); - } -} - - -void Foam::fluidSolidInterface::evolveStress() -{ -// if (closedFluidDomain()) -// { -// DynamicList p0; -// DynamicList dV; - -// scalar requiredVolumeIncrement = 0; -// forAll(flow().U().boundaryField(), patchI) -// { -// if (patchI != fluidPatchIndex()) -// { -// requiredVolumeIncrement += -// sum -// ( -// flow().U().boundaryField()[patchI] -// & fluidMesh().Sf().boundaryField()[patchI] -// ) -// *runTime().deltaT().value(); -// } -// } - -// scalar volumeIncrement = 0; - -// label nIter = 0; - -// do -// { -// // Calc ref. pressure increment - -// if (dV.size() == 1) -// { -// refPressureIncrement_ += compressibility_*dV[0]; -// } -// else if (dV.size() >= 2) -// { -// label i = p0.size() - 1; - -// refPressureIncrement_ = -// p0[i-1] -// - dV[i-1]*(p0[i] - p0[i-1]) -// /(dV[i] - dV[i-1] + SMALL); -// } - -// p0.append(refPressureIncrement_); - -// stress().setPressure -// ( -// solidPatchIndex(), -// solidZoneIndex(), -// solidZonePressure_ + refPressure() + refPressureIncrement() -// ); - - -// // Solve solid -// stress().evolve(); - -// // Calculate volume increment - -// const labelList& meshPoints = -// solidMesh().faceZones()[solidZoneIndex()]().meshPoints(); - -// const faceList& localFaces = -// solidMesh().faceZones()[solidZoneIndex()]().localFaces(); - -// const vectorField& localPoints = -// solidMesh().faceZones()[solidZoneIndex()]().localPoints(); - - -// vectorField oldLocalPoints = -// localPoints -// + vectorField -// ( -// stress().pointD().oldTime(), -// meshPoints -// ); - -// vectorField newLocalPoints = -// localPoints -// + vectorField -// ( -// stress().pointD(), -// meshPoints -// ); - -// volumeIncrement = 0; - -// forAll(localFaces, faceI) -// { -// volumeIncrement += -// localFaces[faceI].sweptVol -// ( -// oldLocalPoints, -// newLocalPoints -// ); -// } - -// volumeIncrement -= requiredVolumeIncrement; - -// dV.append(volumeIncrement); -// } -// while(mag(volumeIncrement) > SMALL*10 && ++nIter < 10); - - -// Info << "Solid volume increment: " << volumeIncrement << endl; -// Info << "Ref pressure: " << refPressure_ << endl; -// Info << "Ref pressure increment: " << refPressureIncrement_ << endl; -// Info << "Calculated compressibility = " -// << (refPressureIncrement() - p0[0])/(dV[0] + SMALL) << endl; -// } -// else -// { -// stress().evolve(); -// } - - stress().evolve(); -} - - -Foam::scalar Foam::fluidSolidInterface::updateResidual() -{ - vectorField solidZonePointsDisplAtSolid = - stress().faceZonePointDisplacementIncrement(solidZoneIndex()); - - vectorField solidZonePointsTotDisplAtSolid = - stress().faceZonePointDisplacement(solidZoneIndex()); - - vectorField solidZonePointsTotDispl - ( - solidZonePointsDispl().size(), - vector::zero - ); - - if (rbfInterpolation_) - { - Info << "Displacement interpolation using RBF interpolation" << endl; - - matrix fluidDispl(solidZonePointsDispl().size(), 3); - matrix solidDispl(solidZonePointsDisplAtSolid.size(), 3); - - forAll(solidZonePointsDisplAtSolid, pointI) - { - solidDispl(pointI, 0) = solidZonePointsDisplAtSolid[pointI].x(); - solidDispl(pointI, 1) = solidZonePointsDisplAtSolid[pointI].y(); - solidDispl(pointI, 2) = solidZonePointsDisplAtSolid[pointI].z(); - } - - solidToFluid()->interpolate(solidDispl, fluidDispl); - - forAll(solidZonePointsDispl(), pointI) - { - solidZonePointsDispl()[pointI].x() = fluidDispl(pointI, 0); - solidZonePointsDispl()[pointI].y() = fluidDispl(pointI, 1); - solidZonePointsDispl()[pointI].z() = fluidDispl(pointI, 2); - } - - // Total displacement - forAll(solidZonePointsTotDisplAtSolid, pointI) - { - solidDispl(pointI, 0) = solidZonePointsTotDisplAtSolid[pointI].x(); - solidDispl(pointI, 1) = solidZonePointsTotDisplAtSolid[pointI].y(); - solidDispl(pointI, 2) = solidZonePointsTotDisplAtSolid[pointI].z(); - } - - solidToFluid()->interpolate(solidDispl, fluidDispl); - - forAll(solidZonePointsTotDispl, pointI) - { - solidZonePointsTotDispl[pointI].x() = fluidDispl(pointI, 0); - solidZonePointsTotDispl[pointI].y() = fluidDispl(pointI, 1); - solidZonePointsTotDispl[pointI].z() = fluidDispl(pointI, 2); - } - } - else - { - solidZonePointsDispl() = - ggiInterpolator().slaveToMasterPointInterpolate - ( - solidZonePointsDisplAtSolid - ); - - solidZonePointsTotDispl = - ggiInterpolator().slaveToMasterPointInterpolate - ( - solidZonePointsTotDisplAtSolid - ); - } - - residualPrev() = residual(); - - residual() = solidZonePointsDispl() - fluidZonePointsDispl(); - -// const scalarField& minLe = minEdgeLength(); - -// scalar residualNorm = gMax(mag(residual())/(minLe + SMALL)); - - scalar residualNorm = ::sqrt(gSum(magSqr(residual()))); - scalar residualNorm_2 = residualNorm; - -// Info << "Current fsi residual norm: " << residualNorm << endl; - - if (residualNorm > maxResidualNorm_) - { - maxResidualNorm_ = residualNorm; - } - -// Info << "Current fsi max residual norm: " << maxResidualNorm_ << endl; - - residualNorm /= maxResidualNorm_ + SMALL; - - Info << "Current fsi relative residual norm: " << residualNorm << endl; - - interfacePointsDisplPrev_ = interfacePointsDispl_; - - interfacePointsDispl_ = solidZonePointsDispl(); -// 0.5*(solidZonePointsDispl() + fluidZonePointsDispl()); - - vectorField intTotDispl = - interfacePointsDispl_ + solidZonePointsTotDispl; - scalar intTotDisplNorm = ::sqrt(gSum(magSqr(intTotDispl))); - if (intTotDisplNorm > maxIntDisplNorm_) - { - maxIntDisplNorm_ = intTotDisplNorm; - } - - residualNorm_2 /= maxIntDisplNorm_ + SMALL; - -// scalar alterResidual = -// max(mag(interfacePointsDispl_ - interfacePointsDisplPrev_)) -// /(max(mag(interfacePointsDispl_ + solidZonePointsTotDispl)) + SMALL); - - Info << "Alternative fsi residual: " << residualNorm_2 << endl; - - return min(residualNorm_2, residualNorm); -// return residualNorm; -} - - -// bool Foam::fluidSolidInterface::read() -// { -// if (regIOobject::read()) -// { -// fluidProperties_ = subDict(type() + "Coeffs"); - -// return true; -// } -// else -// { -// return false; -// } -// } - - -// ************************************************************************* // diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.H deleted file mode 100644 index 359166777..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolidInterface/fluidSolidInterface.H +++ /dev/null @@ -1,589 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -Class - fluidSolidInterface - -Description - Fluid-structure interface class - -Author - Zeljko Tukovic, FSB Zagreb. All rights reserved. - -SourceFiles - fluidSolidInterface.C - -\*---------------------------------------------------------------------------*/ - -#ifndef fluidSolidInterface_H -#define fluidSolidInterface_H - -#include "fluidSolver.H" -#include "solidSolver.H" -#include "IOdictionary.H" -#include "patchToPatchInterpolation.H" -#include "dynamicFvMesh.H" -#include "ggiInterpolation.H" -#include "movingWallVelocityFvPatchVectorField.H" - -#include "RBFInterpolation.H" -#include "TPSFunction.H" -using namespace rbf; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class fluidSolidInterface Declaration -\*---------------------------------------------------------------------------*/ - -class fluidSolidInterface -: - public IOdictionary -{ - // Private data - - //- Fluid (flow) mesh - dynamicFvMesh& fluidMesh_; - - //- Fluid model - autoPtr flow_; - - //- Solid (stress) mesh - fvMesh& solidMesh_; - - //- Stress model - autoPtr stress_; - - //- Solid patch index - label solidPatchIndex_; - - //- Solid face zone index - label solidZoneIndex_; - - //- Fluid patch index - label fluidPatchIndex_; - - //- Fluid face zone index - label fluidZoneIndex_; - - //- Solid face zone current points - mutable vectorField* currentSolidZonePointsPtr_; - - //- Solid face zone primitive patch in current configuration - mutable PrimitivePatch* - currentSolidZonePatchPtr_; - - //- Fluid zone to solid zone interpolator - mutable std::shared_ptr fluidToSolidPtr_; - -// //- Ggi zone-to-zone interpolation -// mutable ggiZoneInterpolation* ggiFluidToSolidPtr_; - - //- Ggi zone-to-zone interpolation - mutable ggiZoneInterpolation* ggiInterpolatorPtr_; - - //- Solid zone to fluid zone interpolator - mutable std::shared_ptr solidToFluidPtr_; - - //- Coupling scheme - word couplingScheme_; - - //- Relaxation factor - scalar relaxationFactor_; - - //- Relaxation factor - scalar aitkenRelaxationFactor_; - - //- Outer correction loop stoping tolerance - scalar outerCorrTolerance_; - - //- Maximal number of outer correctors - label nOuterCorr_; - - //- Fsi coupling on/off - Switch coupled_; - - //- Predictor on/off - Switch predictor_; - - //- Use rbf interpolation - Switch rbfInterpolation_; - - //- Coupling reuse - label couplingReuse_; - - //- Interface deformation limit - scalar interfaceDeformationLimit_; - - //- Fluid zone point displacement - vectorField fluidZonePointsDispl_; - - //- Fluid zone ref. point displacement - vectorField fluidZonePointsDisplRef_; - - //- Fluid zone previous point displacement - vectorField fluidZonePointsDisplPrev_; - - //- Solid zone point displacement - vectorField solidZonePointsDispl_; - - //- Solid zone ref. point displacement - vectorField solidZonePointsDisplRef_; - - //- Solid zone point displacement - vectorField interfacePointsDispl_; - - //- Solid zone ref. point displacement - vectorField interfacePointsDisplPrev_; - - //- Solid zone pressure - scalarField solidZonePressure_; - - //- Solid zone traction (pressure + vicous) - vectorIOField solidZoneTraction_; - - //- Solid zone traction (pressure + vicous) - vectorField solidZoneTractionPrev_; - - //- Solid zone traction (pressure + vicous) - vectorIOField predictedSolidZoneTraction_; - - //- Current fsi residual - vectorField residual_; - - //- Previous fsi residual - vectorField residualPrev_; - - //- Maximal resudual norm - scalar maxResidualNorm_; - - //- Maximal interface displacement norm - scalar maxIntDisplNorm_; - - //- Outer corrector - label outerCorr_; - -// //- Is it fluid domain pure Dirichlet (witout outlets)? -// Switch closedFluidDomain_; - -// //- Reference pressure -// scalar refPressure_; - -// //- Reference pressure -// scalar refPressureIncrement_; - -// //- Current time index -// label timeIndex_; - -// //- Copressibility -// scalar compressibility_; - - //- Interpolator update frequency - label interpolatorUpdateFrequency_; - - - //- IQN-ILS coupling fields - DynamicList fluidPatchPointsV_; - DynamicList fluidPatchPointsW_; - DynamicList fluidPatchPointsT_; - - //- Accumulated fluid side interface displacement - mutable vectorIOField* accumulatedFluidInterfaceDisplacementPtr_; - - //- Min edge length for interface points at fluid side - mutable scalarField* minEdgeLengthPtr_; - - // Private Member Functions - - //- Calculate current solid zone points - void calcCurrentSolidZonePoints() const; - - //- Calculate current solid zone primitive patch - void calcCurrentSolidZonePatch() const; - - //- Calculate fluid to solid interpolator - void calcFluidToSolidInterpolator() const; - -// //- Calculate fluid to solid ggi interpolator -// void calcGgiFluidToSolidInterpolator() const; - - //- Calculate fluid to solid ggi interpolator - void calcGgiInterpolator() const; - - //- Calculate fluid to solid interpolator - void calcSolidToFluidInterpolator() const; - - //- Accumulated fluid interface displacement - void calcAccumulatedFluidInterfaceDisplacement() const; - - //- Calculate minimal edge lengths - void calcMinEdgeLength() const; - - //- Return accumulated interface displacement - vectorIOField& accumulatedFluidInterfaceDisplacement(); - - //- Return minimal edge length - const scalarField& minEdgeLength() const; - - //- Disallow default bitwise copy construct - fluidSolidInterface(const fluidSolidInterface&); - - //- Disallow default bitwise assignment - void operator=(const fluidSolidInterface&); - - -protected: - - // Protected member functions - -public: - - //- Runtime type information - TypeName("fluidSolidInterface"); - - - // Declare run-time constructor selection table - - - // Constructors - - //- Construct from components - fluidSolidInterface - ( - dynamicFvMesh& fluidMesh, - fvMesh& solidMesh - ); - - - // Selectors - - - // Destructor - - virtual ~fluidSolidInterface(); - - - // Member Functions - - // Access - - //- Return fluid mesh - const dynamicFvMesh& fluidMesh() const - { - return fluidMesh_; - } - - //- Return solid mesh - const fvMesh& solidMesh() const - { - return solidMesh_; - } - - //- Return time - const Time& runTime() const - { - return fluidMesh_.time(); - } - - //- Return fluid solver - const fluidSolver& flow() const - { - return flow_(); - } - - //- Return fluid solver - fluidSolver& flow() - { - return flow_(); - } - - //- Return solid solver - const solidSolver& stress() const - { - return stress_(); - } - - //- Return solid solver - solidSolver& stress() - { - return stress_(); - } - - //- Return solid patch index - label solidPatchIndex() const - { - return solidPatchIndex_; - } - - //- Return solid face zone index - label solidZoneIndex() const - { - return solidZoneIndex_; - } - - //- Return fluid patch index - label fluidPatchIndex() const - { - return fluidPatchIndex_; - } - - //- Return fluid face zone index - label fluidZoneIndex() const - { - return fluidZoneIndex_; - } - - //- Return current solid zone points - const vectorField& currentSolidZonePoints() const; - - //- Return current solid zone patch - const PrimitivePatch& - currentSolidZonePatch() const; - - //- Return fluid to solid interpolator - const std::shared_ptr& fluidToSolid() const; - -// //- Return fluid to solid interpolator -// const ggiZoneInterpolation& ggiFluidToSolid() const; - - //- Return fluid to solid interpolator - const ggiZoneInterpolation& ggiInterpolator() const; - - //- Return fluid to solid interpolator - const std::shared_ptr& solidToFluid() const; - - //- Return coupling scheme - const word& couplingScheme() const - { - return couplingScheme_; - } - - //- Return relaxation factor - scalar relaxationFactor() const - { - return relaxationFactor_; - } - - //- Return relaxation factor - scalar& aitkenRelaxationFactor() - { - return aitkenRelaxationFactor_; - } - - //- Return outer corrector loop tolerance - scalar outerCorrTolerance() const - { - return outerCorrTolerance_; - } - - //- Return max numter of outer correctors - label nOuterCorr() const - { - return nOuterCorr_; - } - - //- Is it fluid and structure coupled - const Switch& coupled() const - { - return coupled_; - } - - //- Is it fluid and structure coupled - const Switch& predictor() const - { - return predictor_; - } - - //- Is it fluid and structure coupled - Switch& coupled() - { - return coupled_; - } - - //- Is it fluid and structure coupled - label couplingReuse() const - { - return couplingReuse_; - } - - //- Return relaxation factor - scalar interfaceDeformationLimit() const - { - return interfaceDeformationLimit_; - } - - //- Return fluid zone point displacement - vectorField& fluidZonePointsDispl() - { - return fluidZonePointsDispl_; - } - - //- Return fluid zone ref. point displacement - vectorField& fluidZonePointsDisplRef() - { - return fluidZonePointsDisplRef_; - } - - //- Return fluid zone previous point displacement - vectorField& fluidZonePointsDisplPrev() - { - return fluidZonePointsDisplPrev_; - } - - //- Return solid zone point displacement - vectorField& solidZonePointsDispl() - { - return solidZonePointsDispl_; - } - - //- Return solid zone ref. point displacement - vectorField& solidZonePointsDisplRef() - { - return solidZonePointsDisplRef_; - } - - //- Return solid zone pressure - scalarField& solidZonePressure() - { - return solidZonePressure_; - } - - //- Return solid zone total traction - vectorField& solidZoneTraction() - { - return solidZoneTraction_; - } - - //- Return solid zone total traction - vectorField& solidZoneTractionPrev() - { - return solidZoneTractionPrev_; - } - - //- Return current fsi residual - vectorField& residual() - { - return residual_; - } - - //- Return previous fsi residual - vectorField& residualPrev() - { - return residualPrev_; - } - - //- Return current outer iteration - label& outerCorr() - { - return outerCorr_; - } - - //- Return current outer iteration - const label& outerCorr() const - { - return outerCorr_; - } - -// //- Is it fluid domain pure Dirichlet -// Switch closedFluidDomain() const -// { -// return closedFluidDomain_; -// } - -// //- Ref. pressure used in case of closed fluid domain -// scalar refPressure() -// { -// return refPressure_; -// } - -// //- Ref. pressure used in case of closed fluid domain -// scalar refPressureIncrement() -// { -// return refPressureIncrement_; -// } - -// //- Compressibility used in case of closed fluid domain -// scalar compressibility() -// { -// return compressibility_; -// } - - // Edit - - //- Initialize fields - void initializeFields(); - - //- Initialize fields - void updateInterpolator(); - - //- Calculate interface displacement - void updateDisplacement(); - - //- Calculate interface displacement - scalar updateWeakDisplacement(); - - //- Calculate interface displacement - void updateDisplacementUsingIQNILS(); - - //- Move fluid mesh - void moveFluidMesh(); - - //- Update interface force - void updateForce(); - - //- Update interface force - void updateWeakForce(); - - //- Update interface force - void updateWeakTraction(); - - //- Update interface force - void predictAndUpdateForce(); - - //- Solve solid - void evolveStress(); - - //- Update interface force - scalar updateResidual(); - -// //- Evolve the fluid solver -// virtual void evolve() = 0; - -// //- Read dictionary -// virtual bool read(); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/calcPhi.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/calcPhi.H deleted file mode 100644 index 6e5271002..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/calcPhi.H +++ /dev/null @@ -1,21 +0,0 @@ -{ - phi() = - ( - fvc::interpolate(HU, "interpolate(U)") - /fvc::interpolate(AU, "interpolate(U)") - ) & mesh.Sf(); - - forAll(phi().boundaryField(), patchI) - { - if (!phi().boundaryField()[patchI].coupled()) - { - phi().boundaryField()[patchI] = - ( - U().boundaryField()[patchI] - & mesh.Sf().boundaryField()[patchI] - ); - } - } - - phi() += fvc::ddtPhiCorr(1.0/AU, U(), phi()); -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/consistentIcoFluid.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/consistentIcoFluid.C deleted file mode 100644 index 8aeb5dbd1..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/consistentIcoFluid/consistentIcoFluid.C +++ /dev/null @@ -1,289 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "consistentIcoFluid.H" -#include "volFields.H" -#include "fvm.H" -#include "fvc.H" -#include "fvMatrices.H" -#include "addToRunTimeSelectionTable.H" -#include "adjustPhi.H" - -#include "EulerDdtScheme.H" -#include "CrankNicolsonDdtScheme.H" -#include "backwardDdtScheme.H" - -#include "elasticSlipWallVelocityFvPatchVectorField.H" -#include "elasticWallVelocityFvPatchVectorField.H" -#include "elasticWallPressureFvPatchScalarField.H" - -#include "findRefCell.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ -namespace fluidSolvers -{ - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -defineTypeNameAndDebug(consistentIcoFluid, 0); -addToRunTimeSelectionTable(fluidSolver, consistentIcoFluid, dictionary); -// addToRunTimeSelectionTable(icoFluid, consistentIcoFluid, dictionary); - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void consistentIcoFluid::makeSf() const -{ - // Find global face zones - if (SfPtr_) - { - FatalErrorIn - ( - "void fluidSolver::makeSf() const" - ) - << "Face surface vectors alrady created" - << abort(FatalError); - } - -// const fvMesh& mesh = fluidSolver::mesh(); - - IOobject SfHeader - ( - "Sf", - runTime().timeName(), - mesh(), - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ); - - SfPtr_ = - new surfaceVectorField - ( - SfHeader, - mesh(), - dimensionedVector("0", dimArea, vector::zero) - ); - surfaceVectorField& Sf = *SfPtr_; - - if(!SfHeader.headerOk()) - { - const vectorField& allFaceAreas = mesh().faceAreas(); - - Sf.internalField() = - vectorField::subField(allFaceAreas, mesh().nInternalFaces()); - - const fvPatchList& patches = mesh().boundary(); - - forAll (patches, patchI) - { - Sf.boundaryField()[patchI] = - patches[patchI].patchSlice(allFaceAreas); - } - } -} - - -void consistentIcoFluid::updateSf() -{ - Sf().oldTime(); - - const vectorField& allFaceAreas = mesh().faceAreas(); - - Sf().internalField() = - vectorField::subField(allFaceAreas, mesh().nInternalFaces()); - - const fvPatchList& patches = mesh().boundary(); - - forAll (patches, patchI) - { - Sf().boundaryField()[patchI] = - patches[patchI].patchSlice(allFaceAreas); - } -} - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -consistentIcoFluid::consistentIcoFluid(const fvMesh& mesh) -: - icoFluid(mesh), - SfPtr_(NULL) -{ - phi().oldTime(); - updateSf(); -} - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void consistentIcoFluid::evolve() -{ - Info << "Evolving fluid solver: " << this->type() << endl; - - const fvMesh& mesh = fluidSolver::mesh(); - - updateSf(); - - int nCorr(readInt(fluidProperties().lookup("nCorrectors"))); - - int nNonOrthCorr = - readInt(fluidProperties().lookup("nNonOrthogonalCorrectors")); - - // Prepare for the pressure solution - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p(), fluidProperties(), pRefCell, pRefValue); - - if(mesh.moving()) - { - // Make the fluxes relative - phi() -= fvc::meshPhi(U()); - } - - // CourantNo - { - scalar CoNum = 0.0; - scalar meanCoNum = 0.0; - scalar velMag = 0.0; - - if (mesh.nInternalFaces()) - { - surfaceScalarField SfUfbyDelta = - mesh.surfaceInterpolation::deltaCoeffs()*mag(phi()); - - CoNum = max(SfUfbyDelta/mesh.magSf()) - .value()*runTime().deltaT().value(); - - meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf())) - .value()*runTime().deltaT().value(); - - velMag = max(mag(phi())/mesh.magSf()).value(); - } - - Info<< "Courant Number mean: " << meanCoNum - << " max: " << CoNum - << " velocity magnitude: " << velMag << endl; - } - - fvVectorMatrix UEqn - ( - fvm::ddt(U()) - + fvm::div(phi(), U()) - - fvm::laplacian(nu(), U()) - ); - - solve(UEqn == -gradp()); - - // --- PISO loop - - volScalarField AU = UEqn.A(); - surfaceScalarField rAUf("rAUf", 1.0/fvc::interpolate(AU, "interpolate(U)")); - - for (int corr=0; corr ( new RBFInterpolation() ) ), - rbfCoarse( std::shared_ptr ( new RBFInterpolation( rbf->rbfFunction, rbf->polynomialTerm, rbf->cpu ) ) ), - enabled( false ), - livePointSelection( false ), - livePointSelectionSumValues( false ), - tol( 0 ), - tolLivePointSelection( 0 ), - coarseningMinPoints( 0 ), - coarseningMaxPoints( 0 ), - twoPointSelection( false ), - surfaceCorrection( false ), - ratioRadiusError( 10 ), - exportTxt( false ), - selectedPositions(), - nbStaticFaceCentersRemove( 0 ), - positions(), - positionsInterpolation(), - values(), - errorInterpolationCoarse(), - closestBoundaryIndexCorrection(), - valuesCorrection(), - nbMovingFaceCenters( 0 ), - fileExportIndex( 0 ) - { - assert( rbf ); - } - - RBFCoarsening::RBFCoarsening( std::shared_ptr rbf ) - : - rbf( rbf ), - rbfCoarse( std::shared_ptr ( new RBFInterpolation( rbf->rbfFunction, rbf->polynomialTerm, rbf->cpu ) ) ), - enabled( false ), - livePointSelection( false ), - livePointSelectionSumValues( false ), - tol( 0 ), - tolLivePointSelection( 0 ), - coarseningMinPoints( 0 ), - coarseningMaxPoints( 0 ), - twoPointSelection( false ), - surfaceCorrection( false ), - ratioRadiusError( 10 ), - exportTxt( false ), - selectedPositions(), - nbStaticFaceCentersRemove( 0 ), - positions(), - positionsInterpolation(), - values(), - errorInterpolationCoarse(), - closestBoundaryIndexCorrection(), - valuesCorrection(), - nbMovingFaceCenters( 0 ), - fileExportIndex( 0 ) - { - assert( rbf ); - } - - RBFCoarsening::RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool exportTxt - ) - : - rbf( rbf ), - rbfCoarse( std::shared_ptr ( new RBFInterpolation( rbf->rbfFunction, rbf->polynomialTerm, rbf->cpu ) ) ), - enabled( enabled ), - livePointSelection( livePointSelection ), - livePointSelectionSumValues( livePointSelectionSumValues ), - tol( tol ), - tolLivePointSelection( tolLivePointSelection ), - coarseningMinPoints( coarseningMinPoints ), - coarseningMaxPoints( coarseningMaxPoints ), - twoPointSelection( false ), - surfaceCorrection( false ), - ratioRadiusError( 10.0 ), - exportTxt( exportTxt ), - selectedPositions(), - nbStaticFaceCentersRemove( 0 ), - positions(), - positionsInterpolation(), - values(), - errorInterpolationCoarse(), - closestBoundaryIndexCorrection(), - valuesCorrection(), - nbMovingFaceCenters( 0 ), - fileExportIndex( 0 ) - { - assert( rbf ); - assert( coarseningMinPoints <= coarseningMaxPoints ); - assert( coarseningMinPoints > 0 ); - assert( coarseningMaxPoints > 0 ); - assert( tol > 0 ); - assert( tol < 1 ); - assert( tolLivePointSelection > 0 ); - assert( tolLivePointSelection < 1 ); - - // If unit displacement do not use polynomial for selection - if ( enabled && !livePointSelection && rbf->polynomialTerm ) - { - WarningIn( "RBFCoarsening::RBFCoarsening" ) - << "Unit displacement is combined with polynomial addition into RBF interpolation. Could cause 'strange' results." << endl; - } - } - - RBFCoarsening::RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool twoPointSelection, - bool exportTxt - ) - : - rbf( rbf ), - rbfCoarse( std::shared_ptr ( new RBFInterpolation( rbf->rbfFunction, rbf->polynomialTerm, rbf->cpu ) ) ), - enabled( enabled ), - livePointSelection( livePointSelection ), - livePointSelectionSumValues( livePointSelectionSumValues ), - tol( tol ), - tolLivePointSelection( tolLivePointSelection ), - coarseningMinPoints( coarseningMinPoints ), - coarseningMaxPoints( coarseningMaxPoints ), - twoPointSelection( twoPointSelection ), - surfaceCorrection( false ), - ratioRadiusError( 10.0 ), - exportTxt( exportTxt ), - selectedPositions(), - nbStaticFaceCentersRemove( 0 ), - positions(), - positionsInterpolation(), - values(), - errorInterpolationCoarse(), - closestBoundaryIndexCorrection(), - valuesCorrection(), - nbMovingFaceCenters( 0 ), - fileExportIndex( 0 ) - { - assert( rbf ); - assert( coarseningMinPoints <= coarseningMaxPoints ); - assert( coarseningMinPoints > 0 ); - assert( coarseningMaxPoints > 0 ); - assert( tol > 0 ); - assert( tol < 1 ); - assert( tolLivePointSelection > 0 ); - assert( tolLivePointSelection < 1 ); - - // If unit displacement do not use polynomial for selection - if ( enabled && !livePointSelection && rbf->polynomialTerm ) - { - WarningIn( "RBFCoarsening::RBFCoarsening" ) - << "Unit displacement is combined with polynomial addition into RBF interpolation. Could cause 'strange' results." << endl; - } - } - - RBFCoarsening::RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool twoPointSelection, - bool surfaceCorrection, - scalar ratioRadiusError, - bool exportTxt - ) - : - rbf( rbf ), - rbfCoarse( std::shared_ptr ( new RBFInterpolation( rbf->rbfFunction, rbf->polynomialTerm, rbf->cpu ) ) ), - enabled( enabled ), - livePointSelection( livePointSelection ), - livePointSelectionSumValues( livePointSelectionSumValues ), - tol( tol ), - tolLivePointSelection( tolLivePointSelection ), - coarseningMinPoints( coarseningMinPoints ), - coarseningMaxPoints( coarseningMaxPoints ), - twoPointSelection( twoPointSelection ), - surfaceCorrection( surfaceCorrection ), - ratioRadiusError( ratioRadiusError ), - exportTxt( exportTxt ), - selectedPositions(), - nbStaticFaceCentersRemove( 0 ), - positions(), - positionsInterpolation(), - values(), - errorInterpolationCoarse(), - closestBoundaryIndexCorrection(), - valuesCorrection(), - nbMovingFaceCenters( 0 ), - fileExportIndex( 0 ) - { - assert( rbf ); - assert( coarseningMinPoints <= coarseningMaxPoints ); - assert( coarseningMinPoints > 0 ); - assert( coarseningMaxPoints > 0 ); - assert( tol > 0 ); - assert( tol < 1 ); - assert( tolLivePointSelection > 0 ); - assert( tolLivePointSelection < 1 ); - - // If unit displacement do not use polynomial for selection - if ( enabled && !livePointSelection && rbf->polynomialTerm ) - { - WarningIn( "RBFCoarsening::RBFCoarsening" ) - << "Unit displacement is combined with polynomial addition into RBF interpolation. Could cause 'strange' results." << endl; - } - } - - /* Select a subset of control point with a greedy algorithm. - * The selection of the points is based on a displacement/motion of - * a unit displacement at every control point. Based on a user specified tolerance, - * the control points are selected. - */ - void RBFCoarsening::compute( - const matrix & positions, - const matrix & positionsInterpolation - ) - { - this->positions = positions; - this->positionsInterpolation = positionsInterpolation; - } - - void RBFCoarsening::greedySelection( const matrix & values ) - { - assert( positions.cols() == positionsInterpolation.cols() ); - assert( positions.cols() > 0 ); - assert( positions.rows() > 0 ); - assert( positionsInterpolation.rows() > 0 ); - assert( positions.rows() == values.rows() ); - - matrix usedPositions = positions; - - if ( enabled ) - { - // Greedy algorithm - - rbf::vector errorList( positions.rows() ); - selectedPositions.resize( 2 ); - - for ( int i = 0; i < selectedPositions.rows(); i++ ) - selectedPositions( i ) = i; - - assert( positions.rows() >= selectedPositions.rows() ); - - rbf::matrix positionsInterpolationCoarse = positions; - - int maxNbPoints = std::min( coarseningMaxPoints, static_cast( positions.rows() ) ); - int minPoints = std::min( coarseningMinPoints, static_cast( positions.rows() ) ); - scalar error = 0; - scalar errorMax = 0; - - // Create RBF interpolator - - // Run the greedy algorithm - scalar runTimeInterpolate = 0.0; - scalar runTimeError = 0.0; - scalar runTimeConvergence = 0.0; - bool addedSecondPoint = false; - int counter = selectedPositions.rows(); - - while ( true ) - { - std::clock_t t = std::clock(); - - // Build the matrices used for the RBF interpolation - rbf::matrix positionsCoarse( counter, positions.cols() ); - rbf::matrix valuesCoarse( positionsCoarse.rows(), positionsCoarse.cols() ); - rbf::matrix valuesInterpolationCoarse( positionsInterpolationCoarse.rows(), positionsInterpolationCoarse.cols() ); - - for ( int j = 0; j < selectedPositions.rows(); j++ ) - { - positionsCoarse.row( j ) = positions.row( selectedPositions( j ) ); - valuesCoarse.row( j ) = values.row( selectedPositions( j ) ); - } - - // Perform the RBF interpolation. - rbfCoarse->interpolate( positionsCoarse, positionsInterpolationCoarse, valuesCoarse, valuesInterpolationCoarse ); - - if ( debug > 0 ) - { - t = std::clock() - t; - runTimeInterpolate += static_cast(t) / CLOCKS_PER_SEC; - t = std::clock(); - } - - // Evaluate the error - for ( int j = 0; j < valuesInterpolationCoarse.rows(); j++ ) - errorList( j ) = ( valuesInterpolationCoarse.row( j ) - values.row( j ) ).norm(); - - // Select the point with the largest error which is not already selected. - int index = -1; - scalar largestError = errorList.maxCoeff( &index ); - - // Additional function to check whether the largestError = 0 ( 0 ) - { - t = std::clock() - t; - runTimeError += static_cast(t) / CLOCKS_PER_SEC; - t = std::clock(); - } - - scalar epsilon = std::sqrt( SMALL ); - error = (errorList).norm() / (values.norm() + epsilon); - errorMax = largestError / ( ( values.rowwise().norm() ).maxCoeff() + epsilon ); - - // bool convergence = (error < tol && counter >= minPoints) || counter >= maxNbPoints; - bool convergence = (error < tol && errorMax < tol && counter >= minPoints) || counter >= maxNbPoints; - - if ( convergence ) - { - if ( livePointSelection ) - errorInterpolationCoarse = valuesInterpolationCoarse - values; - - break; - } - - // Break if maximum point are reached - if ( counter >= maxNbPoints ) - { - if ( livePointSelection ) - errorInterpolationCoarse = valuesInterpolationCoarse - values; - - break; - } - - selectedPositions.conservativeResize( selectedPositions.rows() + 1 ); - selectedPositions( selectedPositions.rows() - 1 ) = index; - counter++; - - // Add second point if possible - if ( twoPointSelection && index2 >= 0 && index != index2 ) - { - addedSecondPoint = true; - selectedPositions.conservativeResize( selectedPositions.rows() + 1 ); - selectedPositions( selectedPositions.rows() - 1 ) = index2; - counter++; - } - - if ( debug > 0 ) - { - t = std::clock() - t; - runTimeConvergence += static_cast(t) / CLOCKS_PER_SEC; - t = std::clock(); - } - } - - if ( debug > 0 ) - { - Info << "RBFCoarsening::debug 1. interpolate to surface = " << runTimeInterpolate << " s" << endl; - Info << "RBFCoarsening::debug 2. find largest error = " << runTimeError << " s" << ". Added second point = " << addedSecondPoint << endl; - Info << "RBFCoarsening::debug 3. convergence check = " << runTimeConvergence << " s" << endl; - } - - Info << "RBF interpolation coarsening: selected " << selectedPositions.rows() << "/" << positions.rows() << " points, 2-norm(error) = " - << error << ", max(error) = " << errorMax << ", tol = " << tol << endl; - - rbf::matrix positionsCoarse( selectedPositions.rows(), positions.cols() ); - - for ( int i = 0; i < selectedPositions.rows(); i++ ) - positionsCoarse.row( i ) = positions.row( selectedPositions( i ) ); - - usedPositions = positionsCoarse; - - if ( exportTxt ) - { - std::string fileNameTXT = "rbf-coarsening-greedy-selection-" + std::to_string( fileExportIndex ) + ".txt"; - std::ofstream fileTXT( fileNameTXT ); - - if ( fileTXT.is_open() ) - fileTXT << usedPositions; - - std::string fileNameCSV = "rbf-coarsening-greedy-selection-" + std::to_string( fileExportIndex ) + ".csv"; - std::ofstream fileCSV( fileNameCSV ); - Eigen::IOFormat CSVFmt( Eigen::FullPrecision, Eigen::DontAlignCols, "," ); - - if ( fileCSV.is_open() ) - fileCSV << usedPositions.format( CSVFmt ); - - fileExportIndex++; - } - } - - rbf->compute( usedPositions, positionsInterpolation ); - } - - void RBFCoarsening::interpolate( - const matrix & values, - matrix & valuesInterpolation - ) - { - matrix usedValues = values; - - if ( enabled ) - { - if ( livePointSelection ) - { - // For RBF mesh interpolation, the values to be interpolated need to be - // the total displacements. As input, the incremental displacements - // are given. - if ( livePointSelectionSumValues ) - { - if ( this->values.cols() != values.cols() ) - this->values = values; - else - this->values.array() += values.array(); - } - else - this->values = values; - - // Check if re-selection is necessary - bool reselection = true; - - if ( rbfCoarse->computed ) - { - rbf::matrix valuesCoarse( selectedPositions.rows(), values.cols() ); - rbf::matrix valuesInterpolationCoarse( positions.rows(), valuesInterpolation.cols() ); - - for ( int j = 0; j < selectedPositions.rows(); j++ ) - valuesCoarse.row( j ) = this->values.row( selectedPositions( j ) ); - - rbfCoarse->interpolate2( valuesCoarse, valuesInterpolationCoarse ); - - scalar epsilon = std::sqrt( SMALL ); - - errorInterpolationCoarse = valuesInterpolationCoarse - this->values; - scalar error = (errorInterpolationCoarse).matrix().norm() / (this->values.norm() + epsilon); - scalar errorMax = ( errorInterpolationCoarse.rowwise().norm() ).maxCoeff() / ( ( this->values.rowwise().norm() ).maxCoeff() + epsilon ); - - // bool convergence = error < tolLivePointSelection; - bool convergence = (error < tolLivePointSelection && errorMax < tolLivePointSelection); - - if ( convergence ) - reselection = false; - - Info << "RBF interpolation coarsening: 2-norm(error) = " << error << ", max(error) = " << errorMax << ", tol = " << tolLivePointSelection << ", reselection = "; - - if ( reselection ) - Info << "true"; - else - Info << "false"; - - Info << endl; - } - - if ( reselection ) - { - greedySelection( this->values ); - - rbf->Hhat.conservativeResize( rbf->Hhat.rows(), rbf->Hhat.cols() - nbStaticFaceCentersRemove ); - } - } - else - if ( !rbf->computed ) - { - // Unit displacement of control points - matrix unitDisplacement( positions.rows(), positions.cols() ); - unitDisplacement.setZero(); - - assert( unitDisplacement.rows() >= nbMovingFaceCenters ); - - if ( nbMovingFaceCenters == 0 ) - unitDisplacement.fill( 1 ); - else - for ( int i = 0; i < nbMovingFaceCenters; i++ ) - for ( int j = 0; j < unitDisplacement.cols(); j++ ) - unitDisplacement( i, j ) = 1; - - greedySelection( unitDisplacement ); - - rbf->Hhat.conservativeResize( rbf->Hhat.rows(), rbf->Hhat.cols() - nbStaticFaceCentersRemove ); - } - - rbf::matrix selectedValues( selectedPositions.rows(), values.cols() ); - - for ( int i = 0; i < selectedValues.rows(); i++ ) - selectedValues.row( i ) = values.row( selectedPositions( i ) ); - - usedValues = selectedValues; - } - else - { - if ( !rbf->computed ) - { - rbf->compute( positions, positionsInterpolation ); - rbf->Hhat.conservativeResize( rbf->Hhat.rows(), rbf->Hhat.cols() - nbStaticFaceCentersRemove ); - } - } - - usedValues.conservativeResize( usedValues.rows() - nbStaticFaceCentersRemove, usedValues.cols() ); - rbf->interpolate( usedValues, valuesInterpolation ); - - // start doing correction of surface is requested - if ( livePointSelection && surfaceCorrection ) - { - correctSurface( valuesInterpolation ); - } - } - - void RBFCoarsening::correctSurface( matrix & valuesInterpolation ) - { - if ( valuesCorrection.rows() == 0 ) - { - valuesCorrection.conservativeResize( valuesInterpolation.rows(), valuesInterpolation.cols() ); - valuesCorrection.setZero(); - } - - scalar R = ratioRadiusError * ( errorInterpolationCoarse.rowwise().norm() ).maxCoeff(); - - if ( debug > 0 ) - { - Info << nl << "RBFCoarsening::correctSurface::debug 0: ratioRadiusError = " << ratioRadiusError << ", R = " << R << endl; - } - - // Find nearest boundary point for each internal point. Do this only the first time - vector closestBoundaryRadius( positionsInterpolation.rows() ); - - if ( closestBoundaryIndexCorrection.rows() == 0 ) - { - closestBoundaryIndexCorrection.conservativeResize( positionsInterpolation.rows() ); - std::clock_t t = std::clock(); - scalar runTimeNN = 0; - - for ( int i = 0; i < positionsInterpolation.rows(); i++ ) - { - scalar smallestRadius = GREAT; - int boundaryIndex = -1; - - for ( int j = 0; j < positions.rows(); j++ ) - { - scalar radius = ( positions.row( j ) - positionsInterpolation.row( i ) ).norm(); - - if ( radius < smallestRadius ) - { - boundaryIndex = j; - smallestRadius = radius; - } - } - - closestBoundaryIndexCorrection( i ) = boundaryIndex; - closestBoundaryRadius( i ) = smallestRadius; - } - - if ( debug > 0 ) - { - t = std::clock() - t; - runTimeNN += static_cast(t) / CLOCKS_PER_SEC; - Info << "RBFCoarsening::correctSurface::debug 1. nearest neighbour selection = " << runTimeNN << " s" << endl; - } - } - else - { - for ( int i = 0; i < positionsInterpolation.rows(); i++ ) - { - closestBoundaryRadius( i ) = ( positions.row( closestBoundaryIndexCorrection( i ) ) - positionsInterpolation.row( i ) ).norm(); - } - } - - // Start doing the correction - std::clock_t t = std::clock(); - scalar runTimeCorr = 0; - std::shared_ptr rbfFunction = std::shared_ptr ( new rbf::WendlandC2Function( R ) ); - - for ( int i = 0; i < positionsInterpolation.rows(); i++ ) - { - matrix fEval = -( rbfFunction->evaluate( closestBoundaryRadius( i ) ) ) * errorInterpolationCoarse.row( closestBoundaryIndexCorrection( i ) ); - valuesInterpolation.row( i ) += ( fEval - valuesCorrection.row( i ) ); - valuesCorrection.row( i ) = fEval; - } - - if ( debug > 0 ) - { - t = std::clock() - t; - runTimeCorr += static_cast(t) / CLOCKS_PER_SEC; - Info << "RBFCoarsening::correctSurface::debug 2. correction evaluation = " << runTimeCorr << " s" << endl; - } - } - - void RBFCoarsening::setNbMovingAndStaticFaceCenters( - int nbMovingFaceCenters, - int nbStaticFaceCenters - ) - { - nbStaticFaceCentersRemove = nbStaticFaceCenters; - this->nbMovingFaceCenters = nbMovingFaceCenters; - - if ( enabled ) - { - // Determine the number of selected static face centers - nbStaticFaceCentersRemove = 0; - - for ( int i = 0; i < selectedPositions.rows(); i++ ) - if ( selectedPositions( i ) >= nbMovingFaceCenters ) - nbStaticFaceCentersRemove++; - } - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFCoarsening.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFCoarsening.H deleted file mode 100644 index 49af45211..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFCoarsening.H +++ /dev/null @@ -1,110 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef RBFCoarsening_H -#define RBFCoarsening_H - -#include "RBFInterpolation.H" -#include "fvCFD.H" -#include - -namespace rbf -{ - class RBFCoarsening - { - public: - RBFCoarsening(); - - explicit RBFCoarsening( std::shared_ptr rbf ); - - RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool exportTxt - ); - - RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool twoPointSelection, - bool exportTxt - ); - - RBFCoarsening( - std::shared_ptr rbf, - bool enabled, - bool livePointSelection, - bool livePointSelectionSumValues, - scalar tol, - scalar tolLivePointSelection, - int coarseningMinPoints, - int coarseningMaxPoints, - bool twoPointSelection, - bool surfaceCorrection, - scalar ratioRadiusError, - bool exportTxt - ); - - void greedySelection( const matrix & values ); - - void compute( - const matrix & positions, - const matrix & positionsInterpolation - ); - - void interpolate( - const matrix & values, - matrix & valuesInterpolation - ); - - void setNbMovingAndStaticFaceCenters( - int nbMovingFaceCenters, - int nbStaticFaceCenters - ); - - void correctSurface( matrix & valuesInterpolation ); - - std::shared_ptr rbf; - std::shared_ptr rbfCoarse; - bool enabled; - bool livePointSelection; - bool livePointSelectionSumValues; - scalar tol; - scalar tolLivePointSelection; - int coarseningMinPoints; - int coarseningMaxPoints; - bool twoPointSelection; - bool surfaceCorrection; - scalar ratioRadiusError; - bool exportTxt; - Eigen::VectorXi selectedPositions; - int nbStaticFaceCentersRemove; - matrix positions; - matrix positionsInterpolation; - matrix values; - matrix errorInterpolationCoarse; - Eigen::VectorXi closestBoundaryIndexCorrection; - matrix valuesCorrection; - int nbMovingFaceCenters; - int fileExportIndex; - - static debug::debugSwitch debug; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctionInterface.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctionInterface.H deleted file mode 100644 index c5863a62c..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctionInterface.H +++ /dev/null @@ -1,27 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef RBFFunctionInterface_H -#define RBFFunctionInterface_H - -#include -#include -#include -#include -#include "fvCFD.H" - -namespace rbf -{ - class RBFFunctionInterface - { - public: - virtual ~RBFFunctionInterface(){} - - virtual scalar evaluate( scalar value ) = 0; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.C deleted file mode 100644 index e21557af1..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.C +++ /dev/null @@ -1,21 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "LinearFunction.H" - -namespace rbf -{ - LinearFunction::LinearFunction() - {} - - LinearFunction::~LinearFunction() - {} - - scalar LinearFunction::evaluate( scalar value ) - { - return value; - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.H deleted file mode 100644 index f62ffee0a..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/LinearFunction.H +++ /dev/null @@ -1,25 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef LinearFunction_H -#define LinearFunction_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class LinearFunction : public RBFFunctionInterface - { - public: - LinearFunction(); - - virtual ~LinearFunction(); - - virtual scalar evaluate( scalar value ); - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.C deleted file mode 100644 index 6e410d205..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.C +++ /dev/null @@ -1,24 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "TPSFunction.H" - -namespace rbf -{ - TPSFunction::TPSFunction() - {} - - TPSFunction::~TPSFunction() - {} - - scalar TPSFunction::evaluate( scalar value ) - { - if ( value > 0 ) - return std::log( value ) * value * value; - - return 0L; - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.H deleted file mode 100644 index 531ccb6f3..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/TPSFunction.H +++ /dev/null @@ -1,25 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef TPSFunction_H -#define TPSFunction_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class TPSFunction : public RBFFunctionInterface - { - public: - TPSFunction(); - - virtual ~TPSFunction(); - - virtual scalar evaluate( scalar value ); - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.C deleted file mode 100644 index 9fa89c3c8..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.C +++ /dev/null @@ -1,30 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "WendlandC0Function.H" - -namespace rbf -{ - WendlandC0Function::WendlandC0Function( scalar radius ) - : - radius( radius ) - { - assert( radius > 0 ); - } - - WendlandC0Function::~WendlandC0Function() - {} - - scalar WendlandC0Function::evaluate( scalar value ) - { - value /= radius; - - if ( 1 - value < 0 ) - return 0; - - return std::pow( 1 - value, 2 ); - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.H deleted file mode 100644 index a9ae48d35..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC0Function.H +++ /dev/null @@ -1,27 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef WendlandC0Function_H -#define WendlandC0Function_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class WendlandC0Function : public RBFFunctionInterface - { - public: - explicit WendlandC0Function( scalar radius ); - - virtual ~WendlandC0Function(); - - virtual scalar evaluate( scalar value ); - - scalar radius; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.C deleted file mode 100644 index 76c2c0afe..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.C +++ /dev/null @@ -1,30 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "WendlandC2Function.H" - -namespace rbf -{ - WendlandC2Function::WendlandC2Function( scalar radius ) - : - radius( radius ) - { - assert( radius > 0 ); - } - - WendlandC2Function::~WendlandC2Function() - {} - - scalar WendlandC2Function::evaluate( scalar value ) - { - value /= radius; - - if ( 1 - value < 0 ) - return 0; - - return std::pow( 1 - value, 4 ) * (4 * value + 1); - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.H deleted file mode 100644 index 33899ee7e..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC2Function.H +++ /dev/null @@ -1,27 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef WendlandC2Function_H -#define WendlandC2Function_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class WendlandC2Function : public RBFFunctionInterface - { - public: - explicit WendlandC2Function( scalar radius ); - - virtual ~WendlandC2Function(); - - virtual scalar evaluate( scalar value ); - - scalar radius; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.C deleted file mode 100644 index 9621fac52..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.C +++ /dev/null @@ -1,30 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "WendlandC4Function.H" - -namespace rbf -{ - WendlandC4Function::WendlandC4Function( scalar radius ) - : - radius( radius ) - { - assert( radius > 0 ); - } - - WendlandC4Function::~WendlandC4Function() - {} - - scalar WendlandC4Function::evaluate( scalar value ) - { - value /= radius; - - if ( 1 - value < 0 ) - return 0; - - return std::pow( 1 - value, 6 ) * (35 * std::pow( value, 2 ) + 18 * value + 3); - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.H deleted file mode 100644 index 34cfbe944..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC4Function.H +++ /dev/null @@ -1,27 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef WendlandC4Function_H -#define WendlandC4Function_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class WendlandC4Function : public RBFFunctionInterface - { - public: - explicit WendlandC4Function( scalar radius ); - - virtual ~WendlandC4Function(); - - virtual scalar evaluate( scalar value ); - - scalar radius; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.C deleted file mode 100644 index eaae82e8b..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.C +++ /dev/null @@ -1,30 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "WendlandC6Function.H" - -namespace rbf -{ - WendlandC6Function::WendlandC6Function( scalar radius ) - : - radius( radius ) - { - assert( radius > 0 ); - } - - WendlandC6Function::~WendlandC6Function() - {} - - scalar WendlandC6Function::evaluate( scalar value ) - { - value /= radius; - - if ( 1 - value < 0 ) - return 0; - - return std::pow( 1 - value, 8 ) * (32 * std::pow( value, 3 ) + 25 * std::pow( value, 2 ) + 8 * value + 1); - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.H deleted file mode 100644 index c20784939..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFFunctions/WendlandC6Function.H +++ /dev/null @@ -1,27 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef WendlandC6Function_H -#define WendlandC6Function_H - -#include "RBFFunctionInterface.H" - -namespace rbf -{ - class WendlandC6Function : public RBFFunctionInterface - { - public: - explicit WendlandC6Function( scalar radius ); - - virtual ~WendlandC6Function(); - - virtual scalar evaluate( scalar value ); - - scalar radius; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.C deleted file mode 100644 index a924950b3..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.C +++ /dev/null @@ -1,399 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "RBFInterpolation.H" -#include "TPSFunction.H" - -namespace rbf -{ - RBFInterpolation::RBFInterpolation() - : - rbfFunction( std::shared_ptr ( new TPSFunction() ) ), - polynomialTerm( true ), - cpu( false ), - computed( false ), - n_A( 0 ), - n_B( 0 ), - dimGrid( 0 ), - Hhat(), - Phi(), - lu(), - positions(), - positionsInterpolation() - {} - - RBFInterpolation::RBFInterpolation( std::shared_ptr rbfFunction ) - : - rbfFunction( rbfFunction ), - polynomialTerm( true ), - cpu( false ), - computed( false ), - n_A( 0 ), - n_B( 0 ), - dimGrid( 0 ), - Hhat(), - Phi(), - lu(), - positions(), - positionsInterpolation() - { - assert( rbfFunction ); - } - - RBFInterpolation::RBFInterpolation( - std::shared_ptr rbfFunction, - bool polynomialTerm, - bool cpu - ) - : - rbfFunction( rbfFunction ), - polynomialTerm( polynomialTerm ), - cpu( cpu ), - computed( false ), - n_A( 0 ), - n_B( 0 ), - dimGrid( 0 ), - Hhat(), - Phi(), - lu(), - positions(), - positionsInterpolation() - { - assert( rbfFunction ); - } - - void RBFInterpolation::evaluateH( - const matrix & positions, - matrix & H - ) - { - // RBF function evaluation - - scalar r = 0; - - for ( int i = 0; i < n_A; i++ ) - { - for ( int j = i; j < n_A; j++ ) - { - r = ( positions.row( i ) - positions.row( j ) ).norm(); - H( j, i ) = rbfFunction->evaluate( r ); - } - } - } - - void RBFInterpolation::evaluatePhi( - const matrix & positions, - const matrix & positionsInterpolation, - matrix & Phi - ) - { - // Evaluate Phi which contains the evaluation of the radial basis function - - scalar r = 0; - - for ( int i = 0; i < n_A; i++ ) - { - for ( int j = 0; j < n_B; j++ ) - { - r = ( positions.row( i ) - positionsInterpolation.row( j ) ).norm(); - Phi( j, i ) = rbfFunction->evaluate( r ); - } - } - } - - void RBFInterpolation::compute( - const matrix & positions, - const matrix & positionsInterpolation - ) - { - // Verify input - - assert( positions.cols() == positionsInterpolation.cols() ); - assert( positions.rows() > 0 ); - assert( positions.cols() > 0 ); - assert( positionsInterpolation.rows() > 0 ); - - n_A = positions.rows(); - n_B = positionsInterpolation.rows(); - dimGrid = positions.cols(); - - // Radial basis function interpolation - // Initialize matrices H and Phi - matrix H( n_A, n_A ), Phi( n_B, n_A ); - - if ( polynomialTerm ) - { - H.resize( n_A + dimGrid + 1, n_A + dimGrid + 1 ); - } - - // Evaluate radial basis functions for matrix H - - evaluateH( positions, H ); - - // Include polynomial contributions - if ( polynomialTerm ) - { - for ( int i = 0; i < n_A; i++ ) - H( n_A, i ) = 1; - - H.bottomLeftCorner( dimGrid, n_A ) = positions.block( 0, 0, n_A, dimGrid ).transpose(); - - for ( int i = 0; i < dimGrid + 1; i++ ) - for ( int j = 0; j < dimGrid + 1; j++ ) - H( H.rows() - dimGrid - 1 + i, H.rows() - dimGrid - 1 + j ) = 0; - } - - if ( cpu ) - { - this->positions = positions; - this->positionsInterpolation = positionsInterpolation; - - lu.compute( H.selfadjointView() ); - } - - if ( not cpu ) - { - // Evaluate Phi which contains the evaluation of the radial basis function - - if ( polynomialTerm ) - Phi.resize( n_B, n_A + dimGrid + 1 ); - - evaluatePhi( positions, positionsInterpolation, Phi ); - - // Include polynomial contributions in matrix Phi - if ( polynomialTerm ) - { - for ( int i = 0; i < Phi.rows(); i++ ) - Phi( i, n_A ) = 1; - - Phi.topRightCorner( n_B, dimGrid ) = positionsInterpolation.block( 0, 0, n_B, dimGrid ); - } - - // Compute the LU decomposition of the matrix H - Eigen::FullPivLU lu( H.selfadjointView() ); - - // Compute interpolation matrix - - Hhat.noalias() = Phi * lu.inverse(); - - Hhat.conservativeResize( n_B, n_A ); - } - - computed = true; - } - - void RBFInterpolation::interpolate( - const matrix & values, - matrix & valuesInterpolation - ) - { - if ( cpu && not computed ) - compute( positions, positionsInterpolation ); - - assert( computed ); - - if ( cpu ) - { - matrix B, valuesLU( n_A, values.cols() ), Phi( n_B, n_A ); - - if ( polynomialTerm ) - { - Phi.resize( n_B, n_A + dimGrid + 1 ); - valuesLU.resize( n_A + dimGrid + 1, values.cols() ); - } - - valuesLU.setZero(); - valuesLU.topLeftCorner( values.rows(), values.cols() ) = values; - - B = lu.solve( valuesLU ); - - evaluatePhi( positions, positionsInterpolation, Phi ); - - if ( polynomialTerm ) - { - // Include polynomial contributions in matrix Phi - - for ( int i = 0; i < Phi.rows(); i++ ) - Phi( i, n_A ) = 1; - - Phi.topRightCorner( n_B, dimGrid ) = positionsInterpolation.block( 0, 0, n_B, dimGrid ); - } - - valuesInterpolation.noalias() = Phi * B; - } - - if ( not cpu ) - { - valuesInterpolation.noalias() = Hhat * values; - } - - assert( valuesInterpolation.rows() == n_B ); - assert( values.cols() == valuesInterpolation.cols() ); - } - - /* - * Compute interpolation matrix and directly interpolate the values. - * The algorithms solves for the coefficients, and explicitly - * uses the coefficients to interpolate the data to the new positions. - */ - void RBFInterpolation::interpolate( - const matrix & positions, - const matrix & positionsInterpolation, - const matrix & values, - matrix & valuesInterpolation - ) - { - // Verify input - - assert( positions.cols() == positionsInterpolation.cols() ); - assert( positions.rows() > 0 ); - assert( positions.cols() > 0 ); - assert( positionsInterpolation.rows() > 0 ); - - n_A = positions.rows(); - n_B = positionsInterpolation.rows(); - dimGrid = positions.cols(); - - // Radial basis function interpolation - - // Initialize matrices - - matrix H( n_A, n_A ); - - if ( polynomialTerm ) - H.resize( n_A + dimGrid + 1, n_A + dimGrid + 1 ); - - // RBF function evaluation - - evaluateH( positions, H ); - - // Include polynomial contributions - - if ( polynomialTerm ) - { - // THIJS: initialize Phi if empty - if ( Phi.cols() == 0 ) - { - Phi.conservativeResize( n_B, dimGrid + 1 ); - } - - for ( int i = 0; i < n_A; i++ ) - H( n_A, i ) = 1; - - H.bottomLeftCorner( dimGrid, n_A ) = positions.block( 0, 0, n_A, dimGrid ).transpose(); - - for ( int i = 0; i < dimGrid + 1; i++ ) - for ( int j = 0; j < dimGrid + 1; j++ ) - H( H.rows() - dimGrid - 1 + i, H.rows() - dimGrid - 1 + j ) = 0; - } - - // Calculate coefficients gamma and beta - - matrix B; - - matrix valuesLU( H.rows(), values.cols() ); - valuesLU.setZero(); - valuesLU.topLeftCorner( values.rows(), values.cols() ) = values; - - lu.compute( H.selfadjointView() ); - B = lu.solve( valuesLU ); - - // Evaluate Phi_BA which contains the evaluation of the radial basis function - // This method is only used by the greedy algorithm, and the matrix Phi - // is therefore enlarged at every greedy step. - - buildPhi( positions, positionsInterpolation ); - - if ( polynomialTerm ) - { - // Include polynomial contributions in matrix Phi - - for ( int i = 0; i < Phi.rows(); i++ ) - Phi( i, n_A ) = 1; - - Phi.topRightCorner( n_B, dimGrid ) = positionsInterpolation.block( 0, 0, n_B, dimGrid ); - } - - valuesInterpolation.noalias() = Phi * B; - - computed = true; - } - - void RBFInterpolation::buildPhi( - const matrix & positions, - const matrix & positionsInterpolation - ) - { - n_A = positions.rows(); - n_B = positionsInterpolation.rows(); - int phiColsOld = Phi.cols(); - - if ( polynomialTerm ) - Phi.conservativeResize( n_B, n_A + dimGrid + 1 ); - else - Phi.conservativeResize( n_B, n_A ); - - int nNewPoints = Phi.cols() - phiColsOld; - - if ( nNewPoints == Phi.cols() ) - nNewPoints = n_A; - - scalar r = 0; - - for ( int i = 0; i < nNewPoints; i++ ) - { - int index = Phi.cols() - (i + 1); - - if ( polynomialTerm ) - index = Phi.cols() - 1 - dimGrid - (i + 1); - - for ( int j = 0; j < n_B; j++ ) - { - r = ( positions.row( index ) - positionsInterpolation.row( j ) ).norm(); - Phi( j, index ) = rbfFunction->evaluate( r ); - } - } - } - - /* - * This function is only called by the RBFCoarsening class. - * It is assumed that the polynomial term is included in the - * interpolation, and that the fullPivLu decomposition is - * used to solve for the coefficients B. - */ - void RBFInterpolation::interpolate2( - const matrix & values, - matrix & valuesInterpolation - ) - { - assert( computed ); - - matrix valuesLU( values.rows(), values.cols() ); - - // resize valuesLU if polynomial is used - if ( polynomialTerm ) - { - valuesLU.conservativeResize( values.rows() + values.cols() + 1, values.cols() ); - } - - valuesLU.setZero(); // initialize all values zero - - // Set correct part of valuesLU equal to values - if ( polynomialTerm ) - { - valuesLU.topLeftCorner( values.rows(), values.cols() ) = values; - } - else - { - valuesLU = values; - } - - valuesInterpolation.noalias() = Phi * lu.solve( valuesLU ); - - assert( valuesInterpolation.rows() == n_B ); - assert( values.cols() == valuesInterpolation.cols() ); - } -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.H deleted file mode 100644 index 9f580d58e..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFInterpolation.H +++ /dev/null @@ -1,87 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef RBFInterpolation_H -#define RBFInterpolation_H - -#include -#include -#include "RBFFunctionInterface.H" -#include "fvCFD.H" - -namespace rbf -{ - typedef Eigen::Matrix matrix; - typedef Eigen::Matrix vector; - - class RBFInterpolation - { - public: - RBFInterpolation(); - - explicit RBFInterpolation( std::shared_ptr rbfFunction ); - - RBFInterpolation( - std::shared_ptr rbfFunction, - bool polynomialTerm, - bool cpu - ); - - void compute( - const matrix & positions, - const matrix & positionsInterpolation - ); - - void interpolate( - const matrix & values, - matrix & valuesInterpolation - ); - - void interpolate( - const matrix & positions, - const matrix & positionsInterpolation, - const matrix & values, - matrix & valuesInterpolation - ); - - void interpolate2( - const matrix & values, - matrix & valuesInterpolation - ); - - void buildPhi( - const matrix & positions, - const matrix & positionsInterpolation - ); - - std::shared_ptr rbfFunction; - bool polynomialTerm; - bool cpu; - bool computed; - int n_A; - int n_B; - int dimGrid; - matrix Hhat; - matrix Phi; - Eigen::FullPivLU lu; - matrix positions; - matrix positionsInterpolation; - - private: - void evaluateH( - const matrix & positions, - matrix & H - ); - - void evaluatePhi( - const matrix & positions, - const matrix & positionsInterpolation, - matrix & Phi - ); - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.C deleted file mode 100644 index 6a60e70ce..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.C +++ /dev/null @@ -1,833 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "RBFMeshMotionSolver.H" -#include - -using namespace Foam; - -defineTypeNameAndDebug( RBFMeshMotionSolver, 0 ); - -addToRunTimeSelectionTable -( - motionSolver, - RBFMeshMotionSolver, - dictionary -); - -RBFMeshMotionSolver::RBFMeshMotionSolver( - const polyMesh & mesh, - Istream & msData - ) - : - motionSolver( mesh ), - motionCenters( mesh.boundaryMesh().size(), vectorField( 0 ) ), - staticPatches( lookup( "staticPatches" ) ), - staticPatchIDs( staticPatches.size() ), - movingPatches( lookup( "movingPatches" ) ), - movingPatchIDs( movingPatches.size() ), - fixedPatches( lookup( "fixedPatches" ) ), - fixedPatchIDs( fixedPatches.size() ), - newPoints( mesh.points().size(), vector::zero ), - rbf( false ), - nbGlobalFaceCenters( Pstream::nProcs(), 0 ), - nbGlobalMovingFaceCenters( Pstream::nProcs(), 0 ), - nbGlobalStaticFaceCenters( Pstream::nProcs(), 0 ), - nbGlobalFixedFaceCenters( Pstream::nProcs(), 0 ), - globalMovingPointsLabelList( mesh.boundaryMesh().size(), labelList( 0 ) ), - twoDCorrector( mesh ), - nbPoints( 0 ), - faceCellCenters( true ), - cpu( false ), - // timeIntegrationScheme( false ), - corrector( false ), - k( 0 ) -{ - // Find IDs of staticPatches - forAll( staticPatches, patchI ) - { - label patchIndex = mesh.boundaryMesh().findPatchID( staticPatches[patchI] ); - - assert( patchIndex >= 0 ); - - staticPatchIDs[patchI] = patchIndex; - } - - // Find IDs of movingPatches - forAll( movingPatches, patchI ) - { - label patchIndex = mesh.boundaryMesh().findPatchID( movingPatches[patchI] ); - - assert( patchIndex >= 0 ); - - movingPatchIDs[patchI] = patchIndex; - } - - // Find IDs of fixedPatches - forAll( fixedPatches, patchI ) - { - label patchIndex = mesh.boundaryMesh().findPatchID( fixedPatches[patchI] ); - - assert( patchIndex >= 0 ); - - fixedPatchIDs[patchI] = patchIndex; - } - - // Verify that a patch is not defined as a static and a moving patch - - forAll( staticPatchIDs, staticPatchI ) - { - // Search the moving patches for static patchI - forAll( movingPatchIDs, movingPatchI ) - { - assert( movingPatchIDs[movingPatchI] != staticPatchIDs[staticPatchI] ); - } - - // Search the fixed patches for static patchI - forAll( fixedPatchIDs, fixedPatchI ) - { - assert( fixedPatchIDs[fixedPatchI] != staticPatchIDs[staticPatchI] ); - } - } - - forAll( fixedPatchIDs, fixedPatchI ) - { - // Search the moving patches for fixed patchI - forAll( movingPatchIDs, movingPatchI ) - { - assert( movingPatchIDs[movingPatchI] != fixedPatchIDs[fixedPatchI] ); - } - } - - // Initialize RBF interpolator - - dictionary & dict = subDict( "interpolation" ); - - word function = dict.lookup( "function" ); - - assert( function == "TPS" || function == "WendlandC0" || function == "WendlandC2" || function == "WendlandC4" || function == "WendlandC6" ); - - std::shared_ptr rbfFunction; - - Info << "Radial Basis Function interpolation: Selecting RBF function: " << function << endl; - - if ( function == "TPS" ) - rbfFunction = std::shared_ptr ( new rbf::TPSFunction() ); - - if ( function == "WendlandC0" ) - { - scalar radius = readScalar( dict.lookup( "radius" ) ); - rbfFunction = std::shared_ptr ( new rbf::WendlandC0Function( radius ) ); - } - - if ( function == "WendlandC2" ) - { - scalar radius = readScalar( dict.lookup( "radius" ) ); - rbfFunction = std::shared_ptr ( new rbf::WendlandC2Function( radius ) ); - } - - if ( function == "WendlandC4" ) - { - scalar radius = readScalar( dict.lookup( "radius" ) ); - rbfFunction = std::shared_ptr ( new rbf::WendlandC4Function( radius ) ); - } - - if ( function == "WendlandC6" ) - { - scalar radius = readScalar( dict.lookup( "radius" ) ); - rbfFunction = std::shared_ptr ( new rbf::WendlandC6Function( radius ) ); - } - - assert( rbfFunction ); - - bool polynomialTerm = dict.lookupOrDefault( "polynomial", false ); - bool cpu = dict.lookupOrDefault( "cpu", false ); - this->cpu = dict.lookupOrDefault( "fullCPU", false ); - std::shared_ptr rbfInterpolator( new rbf::RBFInterpolation( rbfFunction, polynomialTerm, cpu ) ); - - if ( this->cpu == true ) - assert( cpu == true ); - - bool coarsening = readBool( subDict( "coarsening" ).lookup( "enabled" ) ); - scalar tol = 0.1; - scalar tolLivePointSelection = 0.1; - bool livePointSelection = false; - bool exportSelectedPoints = false; - int coarseningMinPoints = 1; - int coarseningMaxPoints = 2; - bool twoPointSelection = false; - bool surfaceCorrection = false; - scalar ratioRadiusError = 10.0; - - if ( coarsening ) - { - tol = readScalar( subDict( "coarsening" ).lookup( "tol" ) ); - coarseningMinPoints = readLabel( subDict( "coarsening" ).lookup( "minPoints" ) ); - coarseningMaxPoints = readLabel( subDict( "coarsening" ).lookup( "maxPoints" ) ); - livePointSelection = readBool( subDict( "coarsening" ).lookup( "livePointSelection" ) ); - exportSelectedPoints = readBool( subDict( "coarsening" ).lookup( "exportSelectedPoints" ) ); - twoPointSelection = subDict( "coarsening" ).lookupOrDefault( "twoPointSelection", false ); - } - - if ( livePointSelection ) - { - tolLivePointSelection = readScalar( subDict( "coarsening" ).lookup( "tolLivePointSelection" ) ); - surfaceCorrection = subDict( "coarsening" ).lookupOrDefault( "surfaceCorrection", false ); - - if ( surfaceCorrection ) - { - ratioRadiusError = subDict( "coarsening" ).lookupOrDefault( "ratioRadiusError", 10.0 ); - } - } - - rbf = std::shared_ptr ( new rbf::RBFCoarsening( rbfInterpolator, coarsening, livePointSelection, true, tol, tolLivePointSelection, coarseningMinPoints, coarseningMaxPoints, twoPointSelection, surfaceCorrection, ratioRadiusError, exportSelectedPoints ) ); - - faceCellCenters = lookupOrDefault( "faceCellCenters", true ); - - Info << "RBF mesh deformation settings:" << endl; - Info << " interpolation function = " << function << endl; - Info << " interpolation polynomial term = " << polynomialTerm << endl; - Info << " interpolation cpu formulation = " << cpu << endl; - Info << " coarsening = " << coarsening << endl; - Info << " coarsening tolerance = " << tol << endl; - Info << " coarsening reselection tolerance = " << tolLivePointSelection << endl; - Info << " coarsening two-point selection = " << twoPointSelection << endl; -} - -RBFMeshMotionSolver::~RBFMeshMotionSolver() -{} - -tmp RBFMeshMotionSolver::curPoints() const -{ - // Prepare new points: same as old point - tmp tnewPoints - ( - new vectorField( mesh().nPoints(), vector::zero ) - ); - - pointField & newPoints = tnewPoints(); - - newPoints = this->newPoints; - - // Add old point positions - newPoints += mesh().points(); - - return tnewPoints; -} - -// As a first step, the motion is defined in the -void RBFMeshMotionSolver::setMotion( const Field & motion ) -{ - // Input checking - - assert( motion.size() == mesh().boundaryMesh().size() ); - - forAll( motion, ipatch ) - { - const vectorField & mpatch = motion[ipatch]; - - // Check whether the size of patch motion is equal to number of face centers in patch - if ( faceCellCenters && mpatch.size() > 0 ) - assert( mpatch.size() == mesh().boundaryMesh()[ipatch].faceCentres().size() ); - - if ( not faceCellCenters && mpatch.size() > 0 ) - assert( mpatch.size() == mesh().boundaryMesh()[ipatch].meshPoints().size() ); - - // Check whether the size of a moving patch is equal to the number of face centers in the patch - // First check if patchid is a moving patch - bool movingPatch = false; - forAll( movingPatchIDs, movingPatchI ) - { - if ( movingPatchIDs[movingPatchI] == ipatch ) - movingPatch = true; - } - - if ( faceCellCenters && movingPatch ) - assert( mpatch.size() == mesh().boundaryMesh()[ipatch].faceCentres().size() ); - - if ( not faceCellCenters && movingPatch ) - assert( mpatch.size() == mesh().boundaryMesh()[ipatch].meshPoints().size() ); - } - - motionCenters = motion; -} - -void RBFMeshMotionSolver::updateMesh( const mapPolyMesh & ) -{ - assert( false ); -} - -void RBFMeshMotionSolver::solve() -{ - assert( motionCenters.size() == mesh().boundaryMesh().size() ); - - /* - * RBF interpolator from face centers to local complete mesh vertices - * The interpolation consists of the following steps: - * 1. Build a matrix with the face center positions of the static patches and the moving patches - * 2. Build a matrix with the positions of every vertex in the local mesh - * 3. Build a matrix with the displacement/motion of the face center positions of the static patches and the moving patches - * 4. Perform the interpolation from the face centers to the complete mesh - * 5. Correct the mesh vertices of the static patches. Set these displacement to zero. - * 6. Set the motion of the mesh vertices - */ - - /* - * Step 1: Build a matrix with the face center positions of the static patches and the moving patches - * The order of the matrix is defined as first a list of the moving patch face centers, - * thereafter the static patch face centers. These are the control points used by the - * radial basis function interpolation. - * The control points should be exactly the same at each processor, and are therefore communicated - * to each process. As only an absolute RBF interpolation is implemented, this communication step is only - * performed once per simulation. - * The global ordering of the data is first the information of the moving patches, - * thereafter the static patches. - */ - - unsigned int nbFaceCenters = 0; - unsigned int nbMovingFaceCenters = 0; - unsigned int nbStaticFaceCenters = 0; - unsigned int nbFixedFaceCenters = 0; - - std::unordered_map staticControlPointLabels; - std::unordered_map fixedControlPointLabels; - std::vector movingControlPointLabelsVector; - std::unordered_map movingControlPointLabelsMap; - - std::vector movingControlPointPatchIds; - std::vector movingControlPointIndices; - - std::unordered_map staticControlGlobalPointLabels; - std::unordered_map fixedControlGlobalPointLabels; - std::unordered_map movingControlGlobalPointLabelsMap; - - labelList globalStaticPointsListEnabled( nbStaticFaceCenters, 0 ); - labelList globalFixedPointsListEnabled( nbFixedFaceCenters, 0 ); - labelList globalMovingPointsListEnabled( nbMovingFaceCenters, 0 ); - unsigned int globalStaticOffsetNonUnique = 0; - unsigned int globalFixedOffsetNonUnique = 0; - unsigned int globalMovingOffsetNonUnique = 0; - - if ( sum( nbGlobalFaceCenters ) == 0 ) - { - // Determine the number of face centers - // The total number of face centers is simply the sum of the face centers - // on each processor. - nbFaceCenters = 0; - - // First add the static patches, thereafter the fixed patches, and - // the moving patches as last. - - forAll( staticPatchIDs, i ) - { - const labelList & meshPoints = mesh().boundaryMesh()[staticPatchIDs[i]].meshPoints(); - - forAll( meshPoints, j ) - { - if ( twoDCorrector.marker()[meshPoints[j]] != 0 ) - continue; - - if ( staticControlPointLabels.find( meshPoints[j] ) == staticControlPointLabels.end() ) - staticControlPointLabels[meshPoints[j]] = staticControlPointLabels.size() - 1; - } - } - - nbStaticFaceCenters = staticControlPointLabels.size(); - - forAll( fixedPatchIDs, i ) - { - const labelList & meshPoints = mesh().boundaryMesh()[fixedPatchIDs[i]].meshPoints(); - - forAll( meshPoints, j ) - { - if ( twoDCorrector.marker()[meshPoints[j]] != 0 ) - continue; - - if ( staticControlPointLabels.find( meshPoints[j] ) == staticControlPointLabels.end() - && fixedControlPointLabels.find( meshPoints[j] ) == fixedControlPointLabels.end() ) - fixedControlPointLabels[meshPoints[j]] = fixedControlPointLabels.size() - 1; - } - } - - nbFixedFaceCenters = fixedControlPointLabels.size(); - - if ( faceCellCenters ) - { - forAll( movingPatchIDs, i ) - { - nbMovingFaceCenters += mesh().boundaryMesh()[movingPatchIDs[i]].faceCentres().size(); - } - } - - if ( not faceCellCenters ) - { - forAll( movingPatchIDs, patchI ) - { - const labelList & meshPoints = mesh().boundaryMesh()[movingPatchIDs[patchI]].meshPoints(); - globalMovingPointsLabelList[movingPatchIDs[patchI]] = labelList( meshPoints.size(), 0 ); - - forAll( meshPoints, j ) - { - if ( twoDCorrector.marker()[meshPoints[j]] != 0 ) - continue; - - if ( staticControlPointLabels.find( meshPoints[j] ) == staticControlPointLabels.end() - && fixedControlPointLabels.find( meshPoints[j] ) == fixedControlPointLabels.end() - && movingControlPointLabelsMap.find( meshPoints[j] ) == movingControlPointLabelsMap.end() ) - { - movingControlPointLabelsMap[meshPoints[j]] = movingControlPointLabelsMap.size() - 1; - movingControlPointLabelsVector.push_back( meshPoints[j] ); - movingControlPointPatchIds.push_back( movingPatchIDs[patchI] ); - movingControlPointIndices.push_back( j ); - globalMovingPointsLabelList[movingPatchIDs[patchI]][j] = 1; - } - } - } - - nbMovingFaceCenters = movingControlPointLabelsVector.size(); - } - - if ( Pstream::nProcs() == 1 ) - { - globalStaticPointsListEnabled.resize( nbStaticFaceCenters ); - globalStaticPointsListEnabled = 1; - globalFixedPointsListEnabled.resize( nbFixedFaceCenters ); - globalFixedPointsListEnabled = 1; - globalMovingPointsListEnabled.resize( nbMovingFaceCenters ); - globalMovingPointsListEnabled = 1; - } - - if ( Pstream::nProcs() > 1 ) - { - IOobject addrHeader - ( - "pointProcAddressing", - mesh().facesInstance(), - mesh().meshSubDir, - mesh(), - IOobject::MUST_READ - ); - - assert( addrHeader.headerOk() ); - labelIOList pointProcAddressing( addrHeader ); - - assert( pointProcAddressing.size() == mesh().points().size() ); - - // Count the number of global static points including scalar points - nbGlobalStaticFaceCenters[Pstream::myProcNo()] = nbStaticFaceCenters; - nbGlobalFixedFaceCenters[Pstream::myProcNo()] = nbFixedFaceCenters; - nbGlobalMovingFaceCenters[Pstream::myProcNo()] = nbMovingFaceCenters; - reduce( nbGlobalStaticFaceCenters, sumOp() ); - reduce( nbGlobalFixedFaceCenters, sumOp() ); - reduce( nbGlobalMovingFaceCenters, sumOp() ); - nbStaticFaceCenters = sum( nbGlobalStaticFaceCenters ); - nbFixedFaceCenters = sum( nbGlobalFixedFaceCenters ); - - if ( not faceCellCenters ) - nbMovingFaceCenters = sum( nbGlobalMovingFaceCenters ); - - // Construct a list with all the global point labels, thus including - // also scalar points. Thereafter, construct a list of static control - // list which indicates whether the point is already included or not. - // Use this later to build a list of the unique static control points. - labelList globalStaticPointsList( nbStaticFaceCenters, 0 ); - labelList globalFixedPointsList( nbFixedFaceCenters, 0 ); - labelList globalMovingPointsList( nbMovingFaceCenters, 0 ); - labelList globalMovingPointsPatchIds( nbMovingFaceCenters, 0 ); - labelList globalMovingPointsIndices( nbMovingFaceCenters, 0 ); - - globalStaticOffsetNonUnique = 0; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalStaticOffsetNonUnique += nbGlobalStaticFaceCenters[i]; - - globalFixedOffsetNonUnique = 0; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalFixedOffsetNonUnique += nbGlobalFixedFaceCenters[i]; - - globalMovingOffsetNonUnique = 0; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalMovingOffsetNonUnique += nbGlobalMovingFaceCenters[i]; - - for ( auto label : staticControlPointLabels ) - { - globalStaticPointsList[label.second + globalStaticOffsetNonUnique] = pointProcAddressing[label.first]; - } - - for ( auto label : fixedControlPointLabels ) - { - globalFixedPointsList[label.second + globalFixedOffsetNonUnique] = pointProcAddressing[label.first]; - } - - for ( unsigned int i = 0; i < movingControlPointLabelsVector.size(); ++i ) - { - globalMovingPointsList[i + globalMovingOffsetNonUnique] = pointProcAddressing[movingControlPointLabelsVector[i]]; - globalMovingPointsPatchIds[i + globalMovingOffsetNonUnique] = movingControlPointPatchIds[i]; - globalMovingPointsIndices[i + globalMovingOffsetNonUnique] = movingControlPointIndices[i]; - } - - reduce( globalStaticPointsList, sumOp() ); - reduce( globalFixedPointsList, sumOp() ); - - if ( not faceCellCenters ) - { - reduce( globalMovingPointsList, sumOp() ); - reduce( globalMovingPointsPatchIds, sumOp() ); - reduce( globalMovingPointsIndices, sumOp() ); - } - - // Construct a list of static control points which indicate whether - // should be included or not. - - globalStaticPointsListEnabled.resize( nbStaticFaceCenters ); - globalStaticPointsListEnabled = 0; - globalFixedPointsListEnabled.resize( nbFixedFaceCenters ); - globalFixedPointsListEnabled = 0; - globalMovingPointsListEnabled.resize( nbMovingFaceCenters ); - globalMovingPointsListEnabled = 0; - forAll( globalStaticPointsList, i ) - { - if ( staticControlGlobalPointLabels.find( globalStaticPointsList[i] ) == staticControlGlobalPointLabels.end() ) - { - staticControlGlobalPointLabels[globalStaticPointsList[i]] = staticControlGlobalPointLabels.size() - 1; - globalStaticPointsListEnabled[i] = 1; - } - } - - forAll( globalFixedPointsList, i ) - { - if ( staticControlGlobalPointLabels.find( globalFixedPointsList[i] ) == staticControlGlobalPointLabels.end() - && fixedControlGlobalPointLabels.find( globalFixedPointsList[i] ) == fixedControlGlobalPointLabels.end() ) - { - fixedControlGlobalPointLabels[globalFixedPointsList[i]] = fixedControlGlobalPointLabels.size() - 1; - globalFixedPointsListEnabled[i] = 1; - } - } - - if ( not faceCellCenters ) - { - forAll( movingPatchIDs, patchI ) - { - const labelList & meshPoints = mesh().boundaryMesh()[movingPatchIDs[patchI]].meshPoints(); - globalMovingPointsLabelList[movingPatchIDs[patchI]] = labelList( meshPoints.size(), 0 ); - } - - forAll( globalMovingPointsList, i ) - { - if ( staticControlGlobalPointLabels.find( globalMovingPointsList[i] ) == staticControlGlobalPointLabels.end() - && fixedControlGlobalPointLabels.find( globalMovingPointsList[i] ) == fixedControlGlobalPointLabels.end() - && movingControlGlobalPointLabelsMap.find( globalMovingPointsList[i] ) == movingControlGlobalPointLabelsMap.end() ) - { - movingControlGlobalPointLabelsMap[globalMovingPointsList[i]] = movingControlGlobalPointLabelsMap.size() - 1; - globalMovingPointsListEnabled[i] = 1; - - if ( static_cast(i) < movingControlPointLabelsVector.size() + globalMovingOffsetNonUnique - && static_cast(i) >= globalMovingOffsetNonUnique ) - { - label patchId = globalMovingPointsPatchIds[i]; - label index = globalMovingPointsIndices[i]; - globalMovingPointsLabelList[patchId][index] = 1; - } - } - } - } - - // Count the number of local unique static points - nbStaticFaceCenters = 0; - - for ( auto label : staticControlPointLabels ) - { - if ( globalStaticPointsListEnabled[label.second + globalStaticOffsetNonUnique] == 1 ) - nbStaticFaceCenters++; - } - - nbFixedFaceCenters = 0; - - for ( auto label : fixedControlPointLabels ) - { - if ( globalFixedPointsListEnabled[label.second + globalFixedOffsetNonUnique] == 1 ) - nbFixedFaceCenters++; - } - - if ( not faceCellCenters ) - { - nbMovingFaceCenters = 0; - - for ( unsigned int i = 0; i < movingControlPointLabelsVector.size(); ++i ) - { - if ( globalMovingPointsListEnabled[i + globalMovingOffsetNonUnique] == 1 ) - nbMovingFaceCenters++; - } - } - } - - // Calculate sum of all faces on each processor - nbGlobalStaticFaceCenters = 0; - nbGlobalFixedFaceCenters = 0; - nbGlobalMovingFaceCenters = 0; - nbGlobalMovingFaceCenters[Pstream::myProcNo()] = nbMovingFaceCenters; - nbGlobalStaticFaceCenters[Pstream::myProcNo()] = nbStaticFaceCenters; - nbGlobalFixedFaceCenters[Pstream::myProcNo()] = nbFixedFaceCenters; - nbGlobalFaceCenters[Pstream::myProcNo()] = nbMovingFaceCenters + nbStaticFaceCenters + nbFixedFaceCenters; - - reduce( nbGlobalMovingFaceCenters, sumOp() ); - reduce( nbGlobalStaticFaceCenters, sumOp() ); - reduce( nbGlobalFixedFaceCenters, sumOp() ); - reduce( nbGlobalFaceCenters, sumOp() ); - } - - nbMovingFaceCenters = sum( nbGlobalMovingFaceCenters ); - nbStaticFaceCenters = sum( nbGlobalStaticFaceCenters ); - nbFixedFaceCenters = sum( nbGlobalFixedFaceCenters ); - nbFaceCenters = sum( nbGlobalFaceCenters ); - - // Determine the offset taking into account multiple processors - - int globalMovingOffset = 0; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalMovingOffset += nbGlobalMovingFaceCenters[i]; - - int globalStaticOffset = nbMovingFaceCenters; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalStaticOffset += nbGlobalStaticFaceCenters[i]; - - int globalFixedOffset = nbMovingFaceCenters + nbStaticFaceCenters; - - for ( int i = 0; i < Pstream::myProcNo(); i++ ) - globalFixedOffset += nbGlobalFixedFaceCenters[i]; - - if ( !rbf->rbf->computed ) - { - rbf::matrix positions( nbFaceCenters, mesh().nGeometricD() ); - positions.setZero(); - - const Foam::pointField & points = mesh().points(); - - vectorField positionsField( positions.rows(), vector::zero ); - - if ( faceCellCenters ) - { - int offset = 0; - - forAll( movingPatchIDs, i ) - { - const Foam::vectorField::subField faceCentres = mesh().boundaryMesh()[movingPatchIDs[i]].faceCentres(); - - // Set the positions for patch i - forAll( faceCentres, j ) - { - positionsField[j + offset + globalMovingOffset] = faceCentres[j]; - } - - offset += faceCentres.size(); - } - } - - int index = 0; - - if ( not faceCellCenters ) - { - for ( unsigned int i = 0; i < movingControlPointLabelsVector.size(); ++i ) - { - if ( globalMovingPointsListEnabled[i + globalMovingOffsetNonUnique] == 1 ) - { - assert( index + globalMovingOffset < positionsField.size() ); - positionsField[index + globalMovingOffset] = points[movingControlPointLabelsVector[i]]; - index++; - } - } - - assert( index == nbGlobalMovingFaceCenters[Pstream::myProcNo()] ); - } - - index = 0; - - for ( auto label : staticControlPointLabels ) - { - if ( globalStaticPointsListEnabled[label.second + globalStaticOffsetNonUnique] == 1 ) - { - positionsField[index + globalStaticOffset] = points[label.first]; - index++; - } - } - - assert( index == nbGlobalStaticFaceCenters[Pstream::myProcNo()] ); - - index = 0; - - for ( auto label : fixedControlPointLabels ) - { - if ( globalFixedPointsListEnabled[label.second + globalFixedOffsetNonUnique] == 1 ) - { - assert( index + globalFixedOffset < positionsField.size() ); - positionsField[index + globalFixedOffset] = points[label.first]; - index++; - } - } - - assert( index == nbGlobalFixedFaceCenters[Pstream::myProcNo()] ); - - reduce( positionsField, sumOp() ); - - // Copy the FOAM vector field to an Eigen matrix - for ( int i = 0; i < positions.rows(); i++ ) - for ( int j = 0; j < positions.cols(); j++ ) - positions( i, j ) = positionsField[i][j]; - - /* - * Step 2: Build a matrix with the positions of every vertex in the local mesh. - * This is only local information and does not need to be communicated to other - * processors. - */ - - // Determine the number of points by using the 2d corrector - nbPoints = 0; - forAll( points, i ) - { - if ( twoDCorrector.marker()[i] == 0 ) - nbPoints++; - } - - rbf::matrix positionsInterpolation( nbPoints, positions.cols() ); - - index = 0; - forAll( points, i ) - { - if ( twoDCorrector.marker()[i] == 0 ) - { - for ( int j = 0; j < positionsInterpolation.cols(); j++ ) - positionsInterpolation( index, j ) = points[i][j]; - - index++; - } - } - - rbf->compute( positions, positionsInterpolation ); - - rbf->setNbMovingAndStaticFaceCenters( nbMovingFaceCenters, nbStaticFaceCenters + nbFixedFaceCenters ); - } - - /* - * Step 3: Build a matrix with the displacement/motion of the face center - * positions of the static patches and the moving patches. - * The motion needs to be communicated to every process at every mesh deformation. - * This is considered to be the most expensive step with regards to parallel - * scalability of the overall algorithm. - */ - - rbf::matrix values( nbFaceCenters, mesh().nGeometricD() ); - values.setZero(); - - vectorField valuesField( values.rows(), vector::zero ); - - if ( faceCellCenters ) - { - int offset = 0; - - forAll( movingPatchIDs, i ) - { - const Foam::vectorField::subField faceCentres = mesh().boundaryMesh()[movingPatchIDs[i]].faceCentres(); - - forAll( motionCenters[movingPatchIDs[i]], j ) - { - valuesField[j + offset + globalMovingOffset] = motionCenters[movingPatchIDs[i]][j]; - } - - offset += faceCentres.size(); - } - } - - if ( not faceCellCenters ) - { - int index = 0; - - forAll( movingPatchIDs, patchI ) - { - forAll( globalMovingPointsLabelList[movingPatchIDs[patchI]], j ) - { - if ( globalMovingPointsLabelList[movingPatchIDs[patchI]][j] == 1 ) - { - valuesField[index + globalMovingOffset] = motionCenters[movingPatchIDs[patchI]][j]; - index++; - } - } - } - - assert( index == nbGlobalMovingFaceCenters[Pstream::myProcNo()] ); - } - - reduce( valuesField, sumOp() ); - - // Copy the FOAM vector field to an Eigen matrix - for ( int i = 0; i < values.rows(); i++ ) - for ( int j = 0; j < values.cols(); j++ ) - values( i, j ) = valuesField[i][j]; - - /* - * Step 4: Perform the interpolation from the face centers to the complete mesh - */ - - rbf::matrix valuesInterpolation( nbPoints, values.cols() ); - valuesInterpolation.setZero(); - - if ( cpu ) - rbf->rbf->computed = false; - - rbf->interpolate( values, valuesInterpolation ); - - // Apply the 2d correction - - vectorField valuesInterpolationField( mesh().points().size(), Foam::vector::zero ); - int index = 0; - forAll( valuesInterpolationField, i ) - { - if ( twoDCorrector.marker()[i] == 0 ) - { - for ( int j = 0; j < valuesInterpolation.cols(); j++ ) - valuesInterpolationField[i][j] = valuesInterpolation( index, j ); - - index++; - } - } - - twoDCorrector.setShadowSide( valuesInterpolationField ); - - /* - * Step 5: Correct the mesh vertices of the fixed patches. Set these displacements to zero. - */ - - // Loop over all the patches, and set the fixed patches to zero. - - forAll( mesh().boundaryMesh(), i ) - { - const labelList & meshPoints = mesh().boundaryMesh()[i].meshPoints(); - - bool isFixedPatch = false; - forAll( fixedPatchIDs, j ) - { - if ( i == fixedPatchIDs[j] ) - isFixedPatch = true; - } - - if ( isFixedPatch ) - { - for ( int j = 0; j < meshPoints.size(); j++ ) - valuesInterpolationField[meshPoints[j]] = Foam::vector::zero; - } - } - - /* - * Step 6: Set the motion of the mesh vertices - */ - - assert( newPoints.size() == valuesInterpolationField.size() ); - - newPoints = valuesInterpolationField; -} diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.H deleted file mode 100644 index 58cfbc6e3..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/RBFMeshMotionSolver.H +++ /dev/null @@ -1,97 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef RBFMeshMotionSolver_H -#define RBFMeshMotionSolver_H - -#include "motionSolver.H" -#include "polyMesh.H" -#include "addToRunTimeSelectionTable.H" -#include "RBFInterpolation.H" -#include "RBFCoarsening.H" -#include "TPSFunction.H" -#include "WendlandC0Function.H" -#include "WendlandC2Function.H" -#include "WendlandC4Function.H" -#include "WendlandC6Function.H" -#include "twoDPointCorrectorRBF.H" -#include -#include -// #include "TimeIntegrationScheme.H" - -namespace Foam -{ - class RBFMeshMotionSolver : public motionSolver - { - protected: - // Disallow default bitwise copy construct - RBFMeshMotionSolver( const RBFMeshMotionSolver & ); - - // Disallow default bitwise assignment - void operator=( const RBFMeshMotionSolver & ); - - Field motionCenters; - - wordList staticPatches; - labelList staticPatchIDs; - wordList movingPatches; - labelList movingPatchIDs; - wordList fixedPatches; - labelList fixedPatchIDs; - - pointField newPoints; - - std::shared_ptr rbf; - - labelList nbGlobalFaceCenters; - labelList nbGlobalMovingFaceCenters; - labelList nbGlobalStaticFaceCenters; - labelList nbGlobalFixedFaceCenters; - - List globalMovingPointsLabelList; - - // 2d corrector - twoDPointCorrectorRBF twoDCorrector; - int nbPoints; - bool faceCellCenters; - bool cpu; - - public: - // Runtime type information - TypeName( "RBFMeshMotionSolver" ); - - // Constructors - - // Construct from polyMesh - RBFMeshMotionSolver( - const polyMesh & mesh, - Istream & msData - ); - - // Destructor - virtual ~RBFMeshMotionSolver(); - - // Set the motion of every boundary patch, where m is equal to number of patches and with non-empty vectorFields for moving patches. - // The motion is defined at the face centers of the boundary patch. - void setMotion( const Field & motion ); - - // Return point location obtained from the current motion field - virtual tmp curPoints() const; - - // Solve for motion - virtual void solve(); - - // Update the mesh corresponding to given map - virtual void updateMesh( const mapPolyMesh & ); - - // std::shared_ptr timeIntegrationScheme; - - bool corrector; - int k; - }; -} - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.C deleted file mode 100644 index 209583b51..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.C +++ /dev/null @@ -1,98 +0,0 @@ - -#include "twoDPointCorrectorRBF.H" -#include "polyMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - - // Set marker to -1 if is on the "back" side - void twoDPointCorrectorRBF::setMarker() - { - if ( !required() ) - return; - - // Change size of useablePointIDs_ and shadowPointIDs_ - useablePointIDs_.setSize( mesh_.nPoints() / 2 ); - shadowPointIDs_.setSize( mesh_.nPoints() / 2 ); - - // Get reference to edges - const edgeList & meshEdges = mesh_.edges(); - const pointField & points( mesh_.points() ); - - const labelList & neIndices = normalEdgeIndices(); - const vector & pn = planeNormal(); - - forAll( neIndices, edgeI ) - { - const label & pStartInd = meshEdges[neIndices[edgeI]].start(); - const label & pEndInd = meshEdges[neIndices[edgeI]].end(); - const point & pStart = points[pStartInd]; - const point & pEnd = points[pEndInd]; - - // calculate average point position - const point A = 0.5 * (pStart + pEnd); - - // Calculate inner product with plane normal - scalar pStartInner = ( pn & (pStart - A) ); - scalar pEndInner = ( pn & (pEnd - A) ); - - if ( pStartInner > 0 && pEndInner < 0 ) - { - pointMarker_[pEndInd] = -1; - useablePointIDs_[edgeI] = pStartInd; - shadowPointIDs_[edgeI] = pEndInd; - } - else - if ( pEndInner > 0 && pStartInner < 0 ) - { - pointMarker_[pStartInd] = -1; - useablePointIDs_[edgeI] = pEndInd; - shadowPointIDs_[edgeI] = pStartInd; - } - else - { - FatalErrorIn( "void twoDPointCorrectorRBF::setMarker()" ) - << "Both points give back a negative value with the inner product. Programming error?" - << abort( FatalError ); - } - } - } - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - - twoDPointCorrectorRBF::twoDPointCorrectorRBF( const polyMesh & mesh ) - : - twoDPointCorrector( mesh ), - mesh_( mesh ), - pointMarker_( mesh.nPoints(), 0 ), - useablePointIDs_( mesh.nPoints(), 0 ), - shadowPointIDs_( mesh.nPoints(), 0 ) - { - setMarker(); - } - - // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - - twoDPointCorrectorRBF::~twoDPointCorrectorRBF() - {} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - const labelList & twoDPointCorrectorRBF::marker() const - { - return pointMarker_; - } - - void twoDPointCorrectorRBF::setShadowSide( vectorField & newpoints ) const - { - forAll( useablePointIDs_, ipoint ) - { - newpoints[shadowPointIDs_[ipoint]] = newpoints[useablePointIDs_[ipoint]]; - } - } -} // End namespace Foam diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.H b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.H deleted file mode 100644 index 541562955..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/RBFMeshMotionSolver/twoDPointCorrectorRBF.H +++ /dev/null @@ -1,71 +0,0 @@ - -#ifndef twoDPointCorrectorRBF_H -#define twoDPointCorrectorRBF_H - -#include "pointField.H" -#include "labelList.H" -#include "vector.H" -#include "twoDPointCorrector.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // Forward class declarations - class polyMesh; - - /*---------------------------------------------------------------------------*\ - * Class twoDPointCorrectorRBF Declaration - \*---------------------------------------------------------------------------*/ - - class twoDPointCorrectorRBF - : - public twoDPointCorrector - { - // Private data - - // - Reference to moving mesh - const polyMesh & mesh_; - - // - Holder for marker value: -1 if on "back" side and 0 otherwise - labelList pointMarker_; - - // - Point IDs marked with 0 - labelList useablePointIDs_; - - // - Point IDs marked with -1 - labelList shadowPointIDs_; - - - // Private Member Functions - - // - Disallow default bitwise copy construct - twoDPointCorrectorRBF( const twoDPointCorrectorRBF & ); - - // - Disallow default bitwise assignment - void operator=( const twoDPointCorrectorRBF & ); - - void setMarker(); - - public: - // Constructors - - // - Construct from components - twoDPointCorrectorRBF( const polyMesh & mesh ); - - - // Destructor - - virtual ~twoDPointCorrectorRBF(); - - - // Member Functions - - const labelList & marker() const; - - void setShadowSide( vectorField & newpoints ) const; - }; -} // End namespace Foam - - -#endif diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/EulerDdtScheme.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/EulerDdtScheme.C deleted file mode 100644 index aee60beda..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/EulerDdtScheme.C +++ /dev/null @@ -1,281 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "EulerDdtScheme.H" -#include "surfaceInterpolate.H" -#include "fvcDiv.H" -#include "fvMatrices.H" -#include "slipFvPatchFields.H" -#include "symmetryFvPatchFields.H" -#include "wedgeFvPatchFields.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace fv -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template<> -tmp > -EulerDdtScheme::fvcDdt -( - const GeometricField& vf -) -{ - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - - IOobject ddtIOobject - ( - "ddt("+vf.name()+')', - mesh().time().timeName(), - mesh() - ); - - if (mesh().moving()) - { - if - ( - mesh().objectRegistry::found("grad(" + vf.name() + ")") - && mesh().objectRegistry::found("meshU") - ) - { - const volTensorField& gradVf = - mesh().objectRegistry::lookupObject - ( - "grad(" + vf.name() + ")" - ); - - const volVectorField& meshU = - mesh().objectRegistry::lookupObject - ( - "meshU" - ); - - return tmp > - ( - new GeometricField - ( - ddtIOobject, - rDeltaT*(vf - vf.oldTime()) - (meshU&gradVf.oldTime()) - ) - ); - } - else - { - return tmp > - ( - new GeometricField - ( - ddtIOobject, - rDeltaT*(vf - vf.oldTime()) - ) - ); - } - -// return tmp > -// ( -// new GeometricField -// ( -// ddtIOobject, -// mesh(), -// rDeltaT.dimensions()*vf.dimensions(), -// rDeltaT.value()* -// ( -// vf.internalField() -// - vf.oldTime().internalField() -// ), -// rDeltaT.value()* -// ( -// vf.boundaryField() - vf.oldTime().boundaryField() -// ) -// ) -// ); - } - else - { - return tmp > - ( - new GeometricField - ( - ddtIOobject, - rDeltaT*(vf - vf.oldTime()) - ) - ); - } -} - - -template<> -tmp EulerDdtScheme::fvcDdtPhiCorr -( - const volScalarField& rA, - const GeometricField& U, - const surfaceScalarField& phi -) -{ - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - - IOobject ddtIOobject - ( - "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')', - mesh().time().timeName(), - mesh() - ); - - surfaceScalarField ddtPhiCoeff - ( - IOobject - ( - "ddtPhiCoeff", - mesh().time().timeName(), - mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh(), - dimensioned("1", dimless, 1.0) - ); - - forAll(U.boundaryField(), patchI) - { - if - ( - U.boundaryField()[patchI].fixesValue() - || isA(U.boundaryField()[patchI]) - || isA(U.boundaryField()[patchI]) - || isA(U.boundaryField()[patchI]) - ) - { - ddtPhiCoeff.boundaryField()[patchI] = 0.0; - } - } - - if (mesh().moving()) - { - Info << "Moving mesh ddtPhiCorr " << U.name() << endl; - - volScalarField V0oV - ( - IOobject - ( - "V0oV", - mesh().time().timeName(), - mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh(), - dimless, - zeroGradientFvPatchScalarField::typeName - ); - - V0oV.internalField() = mesh().V0()/mesh().V(); - V0oV.correctBoundaryConditions(); - - const surfaceVectorField& Sf = - mesh().objectRegistry::lookupObject("Sf"); - - // Non-conservative cell-face velocity - surfaceVectorField U0 = - fvc::interpolate(V0oV*U.oldTime(), "interpolate(U)"); - forAll(U0.boundaryField(), patchI) - { - if (!U.boundaryField()[patchI].coupled()) - { - U0.boundaryField()[patchI] = - U.oldTime().boundaryField()[patchI] - .patchInternalField() - *V0oV.boundaryField()[patchI]; - } - } - - // Conservataive cell-face velocity - surfaceVectorField U0c = - fvc::interpolate(U.oldTime(), "interpolate(U)"); - U0c -= (Sf.oldTime()&U0c)*Sf.oldTime()/magSqr(Sf.oldTime()); - U0c += phi.oldTime()*Sf.oldTime()/magSqr(Sf.oldTime()); - - return tmp - ( - new surfaceScalarField - ( - ddtIOobject, - rDeltaT*ddtPhiCoeff - *((fvc::interpolate(V0oV)*U0c - U0) & mesh().Sf()) - /fvc::interpolate(1/rA, "interpolate(U)") - ) - ); - } - else - { - Info << "Fixed mesh ddtPhiCorr " << U.name() << endl; - - // Non-conservative cell-face velocity - surfaceVectorField U0 = - fvc::interpolate(U.oldTime(), "interpolate(U)"); - forAll(U0.boundaryField(), patchI) - { - if (!U.boundaryField()[patchI].coupled()) - { - U0.boundaryField()[patchI] = - U.oldTime().boundaryField()[patchI] - .patchInternalField(); - } - } - - // Conservataive cell-face velocity - surfaceVectorField U0c = U0; -// fvc::interpolate(U.oldTime(), "interpolate(U)"); - U0c -= (mesh().Sf() & U0c)*mesh().Sf()/magSqr(mesh().Sf()); - U0c += phi.oldTime()*mesh().Sf()/magSqr(mesh().Sf()); - - return tmp - ( - new surfaceScalarField - ( - ddtIOobject, - rDeltaT*ddtPhiCoeff - *((U0c - U0) & mesh().Sf()) - /fvc::interpolate(1/rA, "interpolate(U)") - ) - ); - } -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace fv - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/backwardDdtScheme.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/backwardDdtScheme.C deleted file mode 100644 index 741643380..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/finiteVolume/ddtSchemes/backwardDdtScheme.C +++ /dev/null @@ -1,420 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "backwardDdtScheme.H" -#include "surfaceInterpolate.H" -#include "fvcDiv.H" -#include "fvMatrices.H" -#include "symmetryFvPatchFields.H" -#include "slipFvPatchFields.H" -#include "wedgeFvPatchFields.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace fv -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// template<> -// tmp > -// backwardDdtScheme::fvmDdt -// ( -// GeometricField& vf -// ) -// { -// tmp > tfvm -// ( -// new fvMatrix -// ( -// vf, -// vf.dimensions()*dimVol/dimTime -// ) -// ); - -// fvMatrix& fvm = tfvm(); - -// scalar rDeltaT = 1.0/deltaT_(); - -// scalar deltaT = deltaT_(); -// scalar deltaT0 = deltaT0_(vf); - -// scalar coefft = 1 + deltaT/(deltaT + deltaT0); -// scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); -// scalar coefft0 = coefft + coefft00; - -// fvm.diag() = rDeltaT*mesh().V(); -// // fvm.diag() = (coefft*rDeltaT)*mesh().V(); - -// if (mesh().moving()) -// { -// Info << "Corrected backward ddt" << endl; - -// fvm.source() = rDeltaT* -// ( -// vf.oldTime().internalField()*mesh().V0() -// - ( -// // backward -// ( -// coefft*vf.internalField()*mesh().V() -// - coefft0*vf.oldTime().internalField()*mesh().V0() -// + coefft00*vf.oldTime().oldTime().internalField() -// *mesh().V00() -// ) -// // Euler -// - ( -// vf.internalField()*mesh().V() -// - vf.oldTime().internalField()*mesh().V0() -// ) -// ) -// ); -// } -// else -// { -// fvm.source() = rDeltaT*mesh().V()* -// ( -// coefft0*vf.oldTime().internalField() -// - coefft00*vf.oldTime().oldTime().internalField() -// ); -// } - -// return tfvm; -// } - - - -template<> -tmp -backwardDdtScheme::fvcDdtPhiCorr -( - const volScalarField& rA, - const GeometricField& U, - const surfaceScalarField& phi -) -{ - Info << "Consistent backwardDdtPhiCorr" << endl; - - dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); - - IOobject ddtIOobject - ( - "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')', - mesh().time().timeName(), - mesh() - ); - - scalar deltaT = deltaT_(); - scalar deltaT0 = deltaT0_(U); - - scalar coefft = 1 + deltaT/(deltaT + deltaT0); - scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0)); - scalar coefft0 = coefft + coefft00; - - surfaceScalarField ddtPhiCoeff - ( - IOobject - ( - "ddtPhiCoeff", - mesh().time().timeName(), - mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh(), - dimensioned("1", dimless, 1.0) - ); - - forAll(U.boundaryField(), patchI) - { -// if (!U.boundaryField()[patchI].coupled()) -// { -// ddtPhiCoeff.boundaryField()[patchI] = 0.0; -// } - if - ( - U.boundaryField()[patchI].fixesValue() - || isA(U.boundaryField()[patchI]) - || isA(U.boundaryField()[patchI]) - || isA(U.boundaryField()[patchI]) - ) - { - ddtPhiCoeff.boundaryField()[patchI] = 0.0; - } - } - - if (mesh().moving()) - { - volScalarField V0oV - ( - IOobject - ( - "V0oV", - mesh().time().timeName(), - mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh(), - dimless, - zeroGradientFvPatchScalarField::typeName - ); - - V0oV.internalField() = mesh().V0()/mesh().V(); - V0oV.correctBoundaryConditions(); - - volScalarField V00oV - ( - IOobject - ( - "V00oV", - mesh().time().timeName(), - mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh(), - dimless, - zeroGradientFvPatchScalarField::typeName - ); - - V00oV.internalField() = mesh().V00()/mesh().V(); - V00oV.correctBoundaryConditions(); - - if - ( - U.dimensions() == dimVelocity - && phi.dimensions() == dimVelocity*dimArea - ) - { - surfaceVectorField U0 = - fvc::interpolate(U.oldTime(), "interpolate(U)"); - - surfaceVectorField U00 = - fvc::interpolate(U.oldTime().oldTime(), "interpolate(U)"); - - -// surfaceVectorField dU0 = -// fvc::interpolate(U.oldTime()); -// forAll(dU0.boundaryField(), patchI) -// { -// if (!U.boundaryField()[patchI].coupled()) -// { -// dU0.boundaryField()[patchI] = -// U.oldTime().boundaryField()[patchI] -// .patchInternalField(); -// } -// } - -// surfaceVectorField dU00 = -// fvc::interpolate(U.oldTime().oldTime()); -// forAll(dU00.boundaryField(), patchI) -// { -// if (!U.boundaryField()[patchI].coupled()) -// { -// dU00.boundaryField()[patchI] = -// U.oldTime().oldTime().boundaryField()[patchI] -// .patchInternalField(); -// } -// } - - const surfaceVectorField& Sf = - mesh().objectRegistry::lookupObject - ( - "Sf" - ); - - - U0 += Sf.oldTime() - *(phi.oldTime() - (Sf.oldTime()&U0)) - /( - magSqr(Sf.oldTime()) - + dimensionedScalar("SMALL", sqr(dimArea), SMALL) - ); - - U00 += Sf.oldTime().oldTime() - *(phi.oldTime().oldTime() - (Sf.oldTime().oldTime()&U00)) - /( - magSqr(Sf.oldTime().oldTime()) - + dimensionedScalar("SMALL", sqr(dimArea), SMALL) - ); - -// dU0 = Sf.oldTime() -// *(phi.oldTime() - (Sf.oldTime()&dU0)) -// /( -// magSqr(Sf.oldTime()) -// + dimensionedScalar("SMALL", sqr(dimArea), SMALL) -// ); - -// dU00 = Sf.oldTime().oldTime() -// *(phi.oldTime().oldTime() - (Sf.oldTime().oldTime()&dU00)) -// /( -// magSqr(Sf.oldTime().oldTime()) -// + dimensionedScalar("SMALL", sqr(dimArea), SMALL) -// ); - - return tmp - ( - new surfaceScalarField - ( - ddtIOobject, - rDeltaT*ddtPhiCoeff - *( - ( - ( - coefft0*fvc::interpolate(V0oV)*U0 - - coefft00*fvc::interpolate(V00oV)*U00 - ) & mesh().Sf() - ) - - ( - fvc::interpolate - ( - coefft0*V0oV*U.oldTime() - - coefft00*V00oV*U.oldTime().oldTime(), - "interpolate(U)" - ) & mesh().Sf() - ) - )/fvc::interpolate(1/rA, "interpolate(U)") - ) - ); - -// return tmp -// ( -// new surfaceScalarField -// ( -// ddtIOobject, -// rDeltaT*ddtPhiCoeff -// *( -// coefft0*fvc::interpolate(V0oV) -// *(mesh().Sf()&dU0) -// - coefft00 -// *fvc::interpolate(V00oV) -// *(mesh().Sf()&dU00) -// ) -// /fvc::interpolate(1.0/rA) -// ) -// ); - } - else - { - FatalErrorIn - ( - "backwardDdtScheme::fvcDdtPhiCorr" - ) << "dimensions of phi are not correct" - << abort(FatalError); - - return surfaceScalarField::null(); - } - } - else - { - if - ( - U.dimensions() == dimVelocity - && phi.dimensions() == dimVelocity*dimArea - ) - { - surfaceVectorField dU0 = - fvc::interpolate(U.oldTime(), "interpolate(U)"); - forAll(dU0.boundaryField(), patchI) - { - if (!U.boundaryField()[patchI].coupled()) - { - dU0.boundaryField()[patchI] = - U.oldTime().boundaryField()[patchI] - .patchInternalField(); - } - } - - surfaceVectorField dU00 = - fvc::interpolate(U.oldTime().oldTime(), "interpolate(U)"); - forAll(dU00.boundaryField(), patchI) - { - if (!U.boundaryField()[patchI].coupled()) - { - dU00.boundaryField()[patchI] = - U.oldTime().oldTime().boundaryField()[patchI] - .patchInternalField(); - } - } - - dU0 = mesh().Sf() - *(phi.oldTime() - (mesh().Sf()&dU0)) - /( - magSqr(mesh().Sf()) - + dimensionedScalar("SMALL", sqr(dimArea), SMALL) - ); - - dU00 = mesh().Sf() - *(phi.oldTime().oldTime() - (mesh().Sf()&dU00)) - /( - magSqr(mesh().Sf()) - + dimensionedScalar("SMALL", sqr(dimArea), SMALL) - ); - - return tmp - ( - new surfaceScalarField - ( - ddtIOobject, - rDeltaT*ddtPhiCoeff - *( - coefft0*(mesh().Sf()&dU0) - - coefft00*(mesh().Sf()&dU00) - ) - /fvc::interpolate(1.0/rA, "interpolate(U)") - ) - ); - } - else - { - FatalErrorIn - ( - "backwardDdtScheme::fvcDdtPhiCorr" - ) << "dimensions of phi are not correct" - << abort(FatalError); - - return surfaceScalarField::null(); - } - } -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace fv - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/fluidSolver/fluidSolver.C b/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/fluidSolver/fluidSolver.C deleted file mode 100644 index f3a73b5fd..000000000 --- a/applications/solvers/FSI/fluidSolidInteraction/fluidSolvers/fluidSolver/fluidSolver.C +++ /dev/null @@ -1,409 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "fluidSolver.H" -#include "volFields.H" -#include "wedgePolyPatch.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(fluidSolver, 0); - defineRunTimeSelectionTable(fluidSolver, dictionary); -} - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::fluidSolver::calcGlobalFaceZones() const -{ - // Find global face zones - if (globalFaceZonesPtr_) - { - FatalErrorIn - ( - "void fluidSolver::calcGlobalFaceZones() const" - ) - << "Global face zones already fonud" - << abort(FatalError); - } - - SLList