diff --git a/applications/solvers/combustion/dieselEngineFoam/createSpray.H b/applications/solvers/combustion/dieselEngineFoam/createSpray.H index ea4d0121b6..6b78689f09 100644 --- a/applications/solvers/combustion/dieselEngineFoam/createSpray.H +++ b/applications/solvers/combustion/dieselEngineFoam/createSpray.H @@ -30,7 +30,7 @@ scalar gasMass0 = fvc::domainIntegrate(rho).value(); if (dieselSpray.twoD()) { - gasMass0 *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge(); + gasMass0 *= constant::math::twoPi/dieselSpray.angleOfWedge(); } gasMass0 -= diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C index e3f9c0b118..709b70f7a6 100644 --- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C +++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C @@ -43,6 +43,7 @@ Description #include "OFstream.H" #include "volPointInterpolation.H" #include "thermoPhysicsTypes.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H index 31337d0a77..2bacaa91ee 100644 --- a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H +++ b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H @@ -20,7 +20,7 @@ if (dieselSpray.twoD()) { - gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge(); + gasMass *= constant::math::twoPi/dieselSpray.angleOfWedge(); } scalar addedMass = gasMass - gasMass0; diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C index 86460c1e9a..8d5dd0d9fa 100644 --- a/applications/solvers/combustion/dieselFoam/dieselFoam.C +++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C @@ -41,6 +41,7 @@ Description #include "IFstream.H" #include "OFstream.H" #include "Switch.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C index 9adf73cc67..8ff998217d 100644 --- a/applications/solvers/combustion/engineFoam/engineFoam.C +++ b/applications/solvers/combustion/engineFoam/engineFoam.C @@ -58,6 +58,7 @@ Description #include "ignition.H" #include "Switch.H" #include "OFstream.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C index b2d007094e..276cac6ae0 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C +++ b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C @@ -28,7 +28,7 @@ License #include "addToRunTimeSelectionTable.H" #include "fvPatchFieldMapper.H" #include "volFields.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -182,7 +182,7 @@ void smoluchowskiJumpTFvPatchScalarField::updateCoeffs() } Field C2 = pmu/prho - *sqrt(ppsi*mathematicalConstant::pi/2.0) + *sqrt(ppsi*constant::math::piByTwo) *2.0*gamma_/Pr.value()/(gamma_ + 1.0) *(2.0 - accommodationCoeff_)/accommodationCoeff_; diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C index a087d0c964..7cf2ab1f25 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C +++ b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C @@ -28,7 +28,7 @@ Description #include "maxwellSlipUFvPatchVectorField.H" #include "addToRunTimeSelectionTable.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "fvPatchFieldMapper.H" #include "volFields.H" #include "surfaceFields.H" @@ -147,7 +147,7 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs() const fvPatchField& ppsi = patch().lookupPatchField("psi"); - Field C1 = sqrt(ppsi*mathematicalConstant::pi/2.0) + Field C1 = sqrt(ppsi*constant::math::piByTwo) *(2.0 - accommodationCoeff_)/accommodationCoeff_; Field pnu = pmu/prho; diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C index 8a7fb290b8..f6794e92b2 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C +++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C @@ -38,7 +38,6 @@ Description #include "CoalCloud.H" #include "psiChemistryModel.H" #include "chemistrySolver.H" -#include "thermoPhysicsTypes.H" #include "timeActivatedExplicitCellSource.H" #include "radiationModel.H" diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H index 0a5fd91b75..dbd7cf9659 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H @@ -1,5 +1,5 @@ Info<< "\nConstructing coal cloud" << endl; -CoalCloud coalParcels +thermoCoalCloud coalParcels ( "coalCloud1", rho, diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H index 728e605e57..2accb8e1c9 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H @@ -1,5 +1,5 @@ Info<< "\nConstructing reacting cloud" << endl; -BasicReactingCloud parcels +icoPoly8ThermoReactingCloud parcels ( "reactingCloud1", rho, diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C index 4242452440..bf0a16af28 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C @@ -46,7 +46,6 @@ Description #include "BasicReactingCloud.H" #include "rhoChemistryModel.H" #include "chemistrySolver.H" -#include "thermoPhysicsTypes.H" #include "radiationModel.H" #include "porousZones.H" #include "timeActivatedExplicitMulticomponentPointSource.H" diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H index 4ae0633ab7..f3e043143c 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H @@ -1,5 +1,5 @@ Info<< "\nConstructing reacting cloud" << endl; -BasicReactingCloud parcels +thermoReactingCloud parcels ( "reactingCloud1", rho, diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C index 62f8f4834a..4dfd367d09 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C @@ -37,7 +37,6 @@ Description #include "BasicReactingCloud.H" #include "psiChemistryModel.H" #include "chemistrySolver.H" -#include "thermoPhysicsTypes.H" #include "radiationModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C index 02d220138b..03e46a5c63 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "SchnerrSauer.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -77,7 +77,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb { return pow ( - ((4*mathematicalConstant::pi*n_)/3) + ((4*constant::math::pi*n_)/3) *limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1), 1.0/3.0 ); @@ -87,7 +87,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb Foam::dimensionedScalar Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::alphaNuc() const { - dimensionedScalar Vnuc = n_*mathematicalConstant::pi*pow3(dNuc_)/6; + dimensionedScalar Vnuc = n_*constant::math::pi*pow3(dNuc_)/6; return Vnuc/(1 + Vnuc); } diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index 613dbfbaef..240207c8ab 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C @@ -29,11 +29,12 @@ License #include "Time.H" #include "subCycle.H" #include "fvCFD.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * // const scalar Foam::multiphaseMixture::convertToRad = - Foam::mathematicalConstant::pi/180.0; + Foam::constant::math::pi/180.0; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C index 5eb399afbf..fe24603713 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "GidaspowConductivity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -69,13 +69,13 @@ Foam::tmp Foam::GidaspowConductivity::kappa const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); return rhoa*da*sqrt(Theta)* ( 2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi + (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha) - + (15.0/16.0)*sqrtPi*alpha + + (15.0/16.0)*sqrtPi*alpha + (25.0/64.0)*sqrtPi/((1.0 + e)*g0) ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C index 3fdd409f5b..2d0acca9fa 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "HrenyaSinclairConductivity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -74,9 +74,9 @@ Foam::tmp Foam::HrenyaSinclairConductivity::kappa const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); - volScalarField lamda = + volScalarField lamda = scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_; return rhoa*da*sqrt(Theta)* diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C index 8e555dd487..40c296f7e0 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "SyamlalConductivity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -69,7 +69,7 @@ Foam::tmp Foam::SyamlalConductivity::kappa const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); return rhoa*da*sqrt(Theta)* ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index 5a6d4aa52d..dd3f345179 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -26,7 +26,7 @@ License #include "kineticTheoryModel.H" #include "surfaceInterpolate.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "fvCFD.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -204,7 +204,7 @@ void Foam::kineticTheoryModel::solve() volScalarField alpha = alpha_; alpha.max(1.0e-6); - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C index b1a127e2e5..2a1c00b188 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "GidaspowViscosity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -70,7 +70,7 @@ Foam::kineticTheoryModels::GidaspowViscosity::mua const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); return rhoa*da*sqrt(Theta)* ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C index 8b43147895..b5e60e94c2 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "HrenyaSinclairViscosity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -77,7 +77,7 @@ Foam::kineticTheoryModels::HrenyaSinclairViscosity::mua const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); volScalarField lamda = scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C index 28c2c56b07..892319df6c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "SyamlalViscosity.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -69,7 +69,7 @@ Foam::tmp Foam::kineticTheoryModels::SyamlalViscosity::mua const dimensionedScalar& e ) const { - const scalar sqrtPi = sqrt(mathematicalConstant::pi); + const scalar sqrtPi = sqrt(constant::math::pi); return rhoa*da*sqrt(Theta)* ( diff --git a/applications/test/graphXi/graphXi.C b/applications/test/graphXi/graphXi.C index 2c3a5750f3..c717f6cc1f 100644 --- a/applications/test/graphXi/graphXi.C +++ b/applications/test/graphXi/graphXi.C @@ -32,7 +32,7 @@ Description #include "graph.H" #include "OFstream.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" using namespace Foam; @@ -50,7 +50,7 @@ int main() scalarField b = 0.5*(1.0 + erf(x)); scalarField c = 1.0 - b; - scalarField gradb = (1/::sqrt(mathematicalConstant::pi))*exp(-sqr(x)); + scalarField gradb = (1/::sqrt(constant::math::pi))*exp(-sqr(x)); scalarField lapb = -2*x*gradb; r = lapb*b*c/(gradb*gradb); diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index d4d6d11539..0905771d19 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -47,7 +47,7 @@ Description #include "polyTopoChanger.H" #include "polyMesh.H" #include "mapPolyMesh.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "PackedBoolList.H" #include "SortableList.H" @@ -467,7 +467,7 @@ int main(int argc, char *argv[]) scalar angle(readScalar(IStringStream(args.additionalArgs()[1])())); bool overwrite = args.optionFound("overwrite"); - scalar maxCos = Foam::cos(angle*180/mathematicalConstant::pi); + scalar maxCos = Foam::cos(angle*180/constant::math::pi); Info<< "Merging:" << nl << " edges with length less than " << minLen << " meters" << nl diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C index 60c9439561..de7af2b077 100644 --- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C +++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C @@ -53,7 +53,7 @@ Description #include "removePoints.H" #include "polyMesh.H" #include "mapPolyMesh.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" using namespace Foam; @@ -445,12 +445,12 @@ int main(int argc, char *argv[]) scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); - scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.0); + scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0); scalar concaveAngle = defaultConcaveAngle; args.optionReadIfPresent("concaveAngle", concaveAngle); - scalar concaveSin = Foam::sin(concaveAngle*mathematicalConstant::pi/180.0); + scalar concaveSin = Foam::sin(concaveAngle*constant::math::pi/180.0); bool snapMeshDict = args.optionFound("snapMesh"); bool overwrite = args.optionFound("overwrite"); diff --git a/applications/utilities/mesh/advanced/splitCells/splitCells.C b/applications/utilities/mesh/advanced/splitCells/splitCells.C index 85afb3a771..1d9c0e5b54 100644 --- a/applications/utilities/mesh/advanced/splitCells/splitCells.C +++ b/applications/utilities/mesh/advanced/splitCells/splitCells.C @@ -49,7 +49,7 @@ Description #include "cellSet.H" #include "cellModeller.H" #include "meshCutter.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "geomCellLooper.H" #include "plane.H" #include "edgeVertex.H" @@ -539,7 +539,7 @@ int main(int argc, char *argv[]) scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); - scalar radAngle = featureAngle * mathematicalConstant::pi/180.0; + scalar radAngle = featureAngle*constant::math::pi/180.0; scalar minCos = Foam::cos(radAngle); scalar minSin = Foam::sin(radAngle); diff --git a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C index 1da7ba049e..996cd67d4c 100644 --- a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C +++ b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C @@ -43,7 +43,7 @@ Description #include "symmetryPolyPatch.H" #include "wedgePolyPatch.H" #include "cyclicPolyPatch.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" using namespace Foam; diff --git a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H index d315dce4a0..f1ec25a59c 100644 --- a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H +++ b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H @@ -223,7 +223,7 @@ const char* kivaPatchNames[nBCs] = "cylinderHead", "axis", "wedge", - "inflow", + "inflow", "outflow", "presin", "presout", @@ -434,7 +434,7 @@ if (pFaces[WEDGE].size() && pFaces[WEDGE][0].size()) { // Distribute the points to be +/- 2.5deg from the x-z plane - scalar tanTheta = Foam::tan(2.5*mathematicalConstant::pi/180.0); + scalar tanTheta = Foam::tan(2.5*constant::math::pi/180.0); SLList::iterator iterf = pFaces[WEDGE][0].begin(); SLList::iterator iterb = pFaces[WEDGE][1].begin(); diff --git a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C index 2707469a6c..65dd5c952d 100644 --- a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C +++ b/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C @@ -59,7 +59,7 @@ Usage #include "Time.H" #include "timeSelector.H" #include "fvMesh.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "polyTopoChange.H" #include "mapPolyMesh.H" #include "PackedBoolList.H" @@ -91,7 +91,7 @@ void simpleMarkFeatures labelList& multiCellFeaturePoints ) { - scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0); + scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0); const polyBoundaryMesh& patches = mesh.boundaryMesh(); @@ -387,7 +387,7 @@ int main(int argc, char *argv[]) scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); - scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0); + scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0); Info<< "Feature:" << featureAngle << endl << "minCos :" << minCos << endl diff --git a/applications/utilities/mesh/conversion/starToFoam/createCoupleMatches.C b/applications/utilities/mesh/conversion/starToFoam/createCoupleMatches.C index 910df915d8..0b135188e0 100644 --- a/applications/utilities/mesh/conversion/starToFoam/createCoupleMatches.C +++ b/applications/utilities/mesh/conversion/starToFoam/createCoupleMatches.C @@ -33,7 +33,7 @@ Description #include "IOmanip.H" #include "boundBox.H" #include "Map.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -99,7 +99,7 @@ void starMesh::createCoupleMatches() << coupleI << ". STAR couple ID: " << couples_[coupleI].coupleID() << endl << "The angle between face normals is " - << Foam::acos(faceAreaAngle)/mathematicalConstant::pi*180 + << Foam::acos(faceAreaAngle)/constant::math::pi*180 << " deg." << endl << "master cell: " << fp.masterCell() << " STAR number: " << starCellID_[fp.masterCell()] diff --git a/applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C b/applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C index 886281274b..2039ffb0bb 100644 --- a/applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C +++ b/applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C @@ -29,7 +29,7 @@ Description \*---------------------------------------------------------------------------*/ #include "arcEdge.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -78,7 +78,7 @@ Foam::cylindricalCS Foam::arcEdge::calcAngle() // find angles scalar tmp = (r3&r1)/(mag(r3)*mag(r1)); - angle_ = acos(tmp)*180.0/mathematicalConstant::pi; + angle_ = acos(tmp)*180.0/constant::math::pi; // check if the vectors define an exterior or an interior arcEdge if (((r1 ^ r2)&(r1 ^ r3)) < 0.0) angle_ = 360 - angle_; @@ -162,7 +162,7 @@ Foam::vector Foam::arcEdge::position(const scalar lambda) const //- Return the length of the curve Foam::scalar Foam::arcEdge::length() const { - return angle_*radius_*mathematicalConstant::pi/180.0; + return angle_*radius_*constant::math::pi/180.0; } diff --git a/applications/utilities/mesh/generation/extrudeMesh/Make/options b/applications/utilities/mesh/generation/extrudeMesh/Make/options index ce0e27f401..1cdc68a2c7 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/Make/options +++ b/applications/utilities/mesh/generation/extrudeMesh/Make/options @@ -1,10 +1,14 @@ EXE_INC = \ -IextrudedMesh \ -IextrudeModel/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude EXE_LIBS = \ + -lfiniteVolume \ + -lsurfMesh \ -lmeshTools \ -ldynamicMesh \ -lextrudeModel diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C index d149f85923..4d7d6bc81d 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C @@ -36,13 +36,14 @@ Description #include "Time.H" #include "dimensionedTypes.H" #include "IFstream.H" -#include "faceMesh.H" #include "polyTopoChange.H" #include "polyTopoChanger.H" #include "edgeCollapser.H" -#include "mathematicalConstants.H" #include "globalMeshData.H" #include "perfectInterface.H" +#include "addPatchCellLayer.H" +#include "fvMesh.H" +#include "MeshedSurfaces.H" #include "extrudedMesh.H" #include "extrudeModel.H" @@ -50,14 +51,148 @@ Description using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +enum ExtrudeMode +{ + MESH, + PATCH, + SURFACE +}; + +template<> +const char* NamedEnum::names[] = +{ + "mesh", + "patch", + "surface" +}; +static const NamedEnum ExtrudeModeNames; + + +void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName) +{ + // Create dummy system/fv* + { + IOobject io + ( + "fvSchemes", + mesh.time().system(), + regionName, + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + Info<< "Testing:" << io.objectPath() << endl; + + if (!io.headerOk()) + { + Info<< "Writing dummy " << regionName/io.name() << endl; + dictionary dummyDict; + dictionary divDict; + dummyDict.add("divSchemes", divDict); + dictionary gradDict; + dummyDict.add("gradSchemes", gradDict); + dictionary laplDict; + dummyDict.add("laplacianSchemes", laplDict); + + IOdictionary(io, dummyDict).regIOobject::write(); + } + } + { + IOobject io + ( + "fvSolution", + mesh.time().system(), + regionName, + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + if (!io.headerOk()) + { + Info<< "Writing dummy " << regionName/io.name() << endl; + dictionary dummyDict; + IOdictionary(io, dummyDict).regIOobject::write(); + } + } +} + + +label findPatchID(const polyBoundaryMesh& patches, const word& name) +{ + label patchID = patches.findPatchID(name); + + if (patchID == -1) + { + FatalErrorIn("findPatchID(const polyBoundaryMesh&, const word&)") + << "Cannot find patch " << name + << " in the source mesh.\n" + << "Valid patch names are " << patches.names() + << exit(FatalError); + } + return patchID; +} + + +labelList patchFaces(const polyBoundaryMesh& patches, const word& name) +{ + label patchID = findPatchID(patches, name); + const polyPatch& pp = patches[patchID]; + + return identity(pp.size()) + pp.start(); +} + + +void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels) +{ + const labelList& reverseMap = map.reverseFaceMap(); + + labelList newFaceLabels(faceLabels.size()); + label newI = 0; + + forAll(faceLabels, i) + { + label oldFaceI = faceLabels[i]; + + if (reverseMap[oldFaceI] >= 0) + { + newFaceLabels[newI++] = reverseMap[oldFaceI]; + } + } + newFaceLabels.setSize(newI); + faceLabels.transfer(newFaceLabels); +} + + + // Main program: int main(int argc, char *argv[]) { + #include "addRegionOption.H" #include "setRootCase.H" #include "createTimeExtruded.H" - autoPtr meshPtr(NULL); + // Get optional regionName + word regionName; + word regionDir; + if (args.optionReadIfPresent("region", regionName)) + { + regionDir = regionName; + Info<< "Create mesh " << regionName << " for time = " + << runTimeExtruded.timeName() << nl << endl; + } + else + { + regionName = fvMesh::defaultRegion; + Info<< "Create mesh for time = " + << runTimeExtruded.timeName() << nl << endl; + } + IOdictionary dict ( @@ -65,26 +200,44 @@ int main(int argc, char *argv[]) ( "extrudeProperties", runTimeExtruded.constant(), + regionDir, runTimeExtruded, IOobject::MUST_READ ) ); + // Point generator autoPtr model(extrudeModel::New(dict)); - const word sourceType(dict.lookup("constructFrom")); + // Whether to flip normals + const Switch flipNormals(dict.lookup("flipNormals")); - autoPtr fMesh; + // What to extrude + const ExtrudeMode mode = ExtrudeModeNames.read + ( + dict.lookup("constructFrom") + ); - if (sourceType == "patch") + + // Generated mesh (one of either) + autoPtr meshFromMesh; + autoPtr meshFromSurface; + + // Faces on front and back for stitching (in case of mergeFaces) + word frontPatchName; + labelList frontPatchFaces; + word backPatchName; + labelList backPatchFaces; + + if (mode == PATCH || mode == MESH) { fileName sourceCasePath(dict.lookup("sourceCase")); sourceCasePath.expand(); fileName sourceRootDir = sourceCasePath.path(); fileName sourceCaseDir = sourceCasePath.name(); - word patchName(dict.lookup("sourcePatch")); + dict.lookup("sourcePatch") >> frontPatchName; - Info<< "Extruding patch " << patchName + Info<< "Extruding patch " << frontPatchName << " on mesh " << sourceCasePath << nl << endl; @@ -94,31 +247,183 @@ int main(int argc, char *argv[]) sourceRootDir, sourceCaseDir ); - #include "createPolyMesh.H" + #include "createMesh.H" - label patchID = mesh.boundaryMesh().findPatchID(patchName); + const polyBoundaryMesh& patches = mesh.boundaryMesh(); - if (patchID == -1) + + // Topo change container. Either copy an existing mesh or start + // with empty storage (number of patches only needed for checking) + autoPtr meshMod + ( + ( + mode == MESH + ? new polyTopoChange(mesh) + : new polyTopoChange(patches.size()) + ) + ); + + // Extrusion engine. Either adding to existing mesh or + // creating separate mesh. + addPatchCellLayer layerExtrude(mesh, (mode == MESH)); + + indirectPrimitivePatch extrudePatch + ( + IndirectList + ( + mesh.faces(), + patchFaces(patches, frontPatchName) + ), + mesh.points() + ); + + // Only used for addPatchCellLayer into new mesh + labelList exposedPatchIDs; + if (mode == PATCH) { - FatalErrorIn(args.executable()) - << "Cannot find patch " << patchName - << " in the source mesh.\n" - << "Valid patch names are " << mesh.boundaryMesh().names() - << exit(FatalError); + dict.lookup("exposedPatchName") >> backPatchName; + exposedPatchIDs.setSize + ( + extrudePatch.size(), + findPatchID(patches, backPatchName) + ); } - const polyPatch& pp = mesh.boundaryMesh()[patchID]; - fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints())); + pointField layer0Points(extrudePatch.nPoints()); + pointField displacement(extrudePatch.nPoints()); + forAll(displacement, pointI) { - fileName surfName(runTime.path()/patchName + ".sMesh"); - Info<< "Writing patch as surfaceMesh to " - << surfName << nl << endl; - OFstream os(surfName); - os << fMesh() << nl; + const vector& patchNormal = extrudePatch.pointNormals()[pointI]; + + // layer0 point + layer0Points[pointI] = model() + ( + extrudePatch.localPoints()[pointI], + patchNormal, + 0 + ); + // layerN point + point extrudePt = model() + ( + extrudePatch.localPoints()[pointI], + patchNormal, + model().nLayers() + ); + displacement[pointI] = extrudePt - layer0Points[pointI]; } + + if (flipNormals) + { + Info<< "Flipping faces." << nl << endl; + displacement = -displacement; + } + + // Check if wedge (has layer0 different from original patch points) + // If so move the mesh to starting position. + if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL) + { + Info<< "Moving mesh to layer0 points since differ from original" + << " points - this can happen for wedge extrusions." << nl + << endl; + + pointField newPoints(mesh.points()); + forAll(extrudePatch.meshPoints(), i) + { + newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i]; + } + mesh.movePoints(newPoints); + } + + + // Layers per face + labelList nFaceLayers(extrudePatch.size(), model().nLayers()); + // Layers per point + labelList nPointLayers(extrudePatch.nPoints(), model().nLayers()); + // Displacement for first layer + vectorField firstLayerDisp = displacement*model().sumThickness(1); + // Expansion ratio not used. + scalarField ratio(extrudePatch.nPoints(), 1.0); + + layerExtrude.setRefinement + ( + ratio, // expansion ratio + extrudePatch, // patch faces to extrude + exposedPatchIDs, // if new mesh: patches for exposed faces + nFaceLayers, + nPointLayers, + firstLayerDisp, + meshMod() + ); + + // Reset points according to extrusion model + forAll(layerExtrude.addedPoints(), pointI) + { + const labelList& pPoints = layerExtrude.addedPoints()[pointI]; + forAll(pPoints, pPointI) + { + label meshPointI = pPoints[pPointI]; + + point modelPt + ( + model() + ( + extrudePatch.localPoints()[pointI], + extrudePatch.pointNormals()[pointI], + pPointI+1 // layer + ) + ); + + const_cast&> + ( + meshMod().points() + )[meshPointI] = modelPt; + } + } + + // Store faces on front and exposed patch (if mode=patch there are + // only added faces so cannot used map to old faces) + const labelListList& layerFaces = layerExtrude.layerFaces(); + backPatchFaces.setSize(layerFaces.size()); + frontPatchFaces.setSize(layerFaces.size()); + forAll(backPatchFaces, i) + { + backPatchFaces[i] = layerFaces[i][0]; + frontPatchFaces[i] = layerFaces[i][layerFaces[i].size()-1]; + } + + // Create dummy fvSchemes, fvSolution + createDummyFvMeshFiles(mesh, regionDir); + + // Create actual mesh from polyTopoChange container + autoPtr map = meshMod().makeMesh + ( + meshFromMesh, + IOobject + ( + regionName, + runTimeExtruded.constant(), + runTimeExtruded, + IOobject::NO_READ, + IOobject::AUTO_WRITE, + false + ), + mesh + ); + + // Calculate face labels for front and back. + frontPatchFaces = renumber + ( + map().reverseFaceMap(), + frontPatchFaces + ); + backPatchFaces = renumber + ( + map().reverseFaceMap(), + backPatchFaces + ); } - else if (sourceType == "surface") + else { // Read from surface fileName surfName(dict.lookup("surface")); @@ -126,52 +431,62 @@ int main(int argc, char *argv[]) Info<< "Extruding surfaceMesh read from file " << surfName << nl << endl; - IFstream is(surfName); + MeshedSurface fMesh(surfName); - fMesh.reset(new faceMesh(is)); - - Info<< "Read patch from file " << surfName << nl - << endl; - } - else - { - FatalErrorIn(args.executable()) - << "Illegal 'constructFrom' specification. Should either be " - << "patch or surface." << exit(FatalError); - } - - Switch flipNormals(dict.lookup("flipNormals")); - - if (flipNormals) - { - Info<< "Flipping faces." << nl << endl; - - faceList faces(fMesh().size()); - forAll(faces, i) + if (flipNormals) { - faces[i] = fMesh()[i].reverseFace(); + Info<< "Flipping faces." << nl << endl; + faceList& faces = const_cast(fMesh.faces()); + forAll(faces, i) + { + faces[i] = fMesh[i].reverseFace(); + } } - fMesh.reset(new faceMesh(faces, fMesh().localPoints())); + + Info<< "Extruding surface with :" << nl + << " points : " << fMesh.points().size() << nl + << " faces : " << fMesh.size() << nl + << " normals[0] : " << fMesh.faceNormals()[0] + << nl + << endl; + + meshFromSurface.reset + ( + new extrudedMesh + ( + IOobject + ( + extrudedMesh::defaultRegion, + runTimeExtruded.constant(), + runTimeExtruded + ), + fMesh, + model() + ) + ); + + + // Get the faces on front and back + frontPatchName = "originalPatch"; + frontPatchFaces = patchFaces + ( + meshFromSurface().boundaryMesh(), + frontPatchName + ); + backPatchName = "otherSide"; + backPatchFaces = patchFaces + ( + meshFromSurface().boundaryMesh(), + backPatchName + ); } - Info<< "Extruding patch with :" << nl - << " points : " << fMesh().points().size() << nl - << " faces : " << fMesh().size() << nl - << " normals[0] : " << fMesh().faceNormals()[0] - << nl - << endl; - - extrudedMesh mesh + polyMesh& mesh = ( - IOobject - ( - extrudedMesh::defaultRegion, - runTimeExtruded.constant(), - runTimeExtruded - ), - fMesh(), - model() + meshFromMesh.valid() + ? meshFromMesh() + : meshFromSurface() ); @@ -184,17 +499,6 @@ int main(int argc, char *argv[]) << "Merge distance : " << mergeDim << nl << endl; - const polyBoundaryMesh& patches = mesh.boundaryMesh(); - - const label origPatchID = patches.findPatchID("originalPatch"); - const label otherPatchID = patches.findPatchID("otherSide"); - - if (origPatchID == -1 || otherPatchID == -1) - { - FatalErrorIn(args.executable()) - << "Cannot find patch originalPatch or otherSide." << nl - << "Valid patches are " << patches.names() << exit(FatalError); - } // Collapse edges // ~~~~~~~~~~~~~~ @@ -237,6 +541,10 @@ int main(int argc, char *argv[]) // Update fields mesh.updateMesh(map); + // Update stored data + updateFaceLabels(map(), frontPatchFaces); + updateFaceLabels(map(), backPatchFaces); + // Move mesh (if inflation used) if (map().hasMotionPoints()) { @@ -252,22 +560,33 @@ int main(int argc, char *argv[]) Switch mergeFaces(dict.lookup("mergeFaces")); if (mergeFaces) { + if (mode == MESH) + { + FatalErrorIn(args.executable()) + << "Cannot stitch front and back of extrusion since" + << " in 'mesh' mode (extrusion appended to mesh)." + << exit(FatalError); + } + Info<< "Assuming full 360 degree axisymmetric case;" << " stitching faces on patches " - << patches[origPatchID].name() << " and " - << patches[otherPatchID].name() << " together ..." << nl << endl; + << frontPatchName << " and " + << backPatchName << " together ..." << nl << endl; + + if (frontPatchFaces.size() != backPatchFaces.size()) + { + FatalErrorIn(args.executable()) + << "Differing number of faces on front (" + << frontPatchFaces.size() << ") and back (" + << backPatchFaces.size() << ")" + << exit(FatalError); + } + + polyTopoChanger stitcher(mesh); stitcher.setSize(1); - // Make list of masterPatch faces - labelList isf(patches[origPatchID].size()); - - forAll (isf, i) - { - isf[i] = patches[origPatchID].start() + i; - } - const word cutZoneName("originalCutFaceZone"); List fz @@ -276,8 +595,8 @@ int main(int argc, char *argv[]) new faceZone ( cutZoneName, - isf, - boolList(isf.size(), false), + frontPatchFaces, + boolList(frontPatchFaces.size(), false), 0, mesh.faceZones() ) @@ -286,26 +605,58 @@ int main(int argc, char *argv[]) mesh.addZones(List(0), fz, List(0)); // Add the perfect interface mesh modifier - stitcher.set + perfectInterface perfectStitcher ( + "couple", 0, - new perfectInterface - ( - "couple", - 0, - stitcher, - cutZoneName, - patches[origPatchID].name(), - patches[otherPatchID].name() - ) + stitcher, + cutZoneName, + word::null, // dummy patch name + word::null // dummy patch name ); - // Execute all polyMeshModifiers - autoPtr morphMap = stitcher.changeMesh(true); + // Topo change container + polyTopoChange meshMod(mesh); - mesh.movePoints(morphMap->preMotionPoints()); + perfectStitcher.setRefinement + ( + indirectPrimitivePatch + ( + IndirectList + ( + mesh.faces(), + frontPatchFaces + ), + mesh.points() + ), + indirectPrimitivePatch + ( + IndirectList + ( + mesh.faces(), + backPatchFaces + ), + mesh.points() + ), + meshMod + ); + + // Construct new mesh from polyTopoChange. + autoPtr map = meshMod.changeMesh(mesh, false); + + // Update fields + mesh.updateMesh(map); + + // Move mesh (if inflation used) + if (map().hasMotionPoints()) + { + mesh.movePoints(map().preMotionPoints()); + } } + mesh.setInstance(runTimeExtruded.constant()); + Info<< "Writing mesh to " << mesh.objectPath() << nl << endl; + if (!mesh.write()) { FatalErrorIn(args.executable()) << "Failed writing mesh" diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C index 7305431b96..5550c56b15 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C @@ -43,6 +43,7 @@ Foam::extrudeModel::extrudeModel ) : nLayers_(readLabel(dict.lookup("nLayers"))), + expansionRatio_(readScalar(dict.lookup("expansionRatio"))), dict_(dict), coeffDict_(dict.subDict(modelType + "Coeffs")) {} @@ -62,5 +63,28 @@ Foam::label Foam::extrudeModel::nLayers() const } +Foam::scalar Foam::extrudeModel::expansionRatio() const +{ + return expansionRatio_; +} + + +Foam::scalar Foam::extrudeModel::sumThickness(const label layer) const +{ + // 1+r+r^2+ .. +r^(n-1) = (1-r^n)/(1-r) + + if (mag(1.0-expansionRatio_) < SMALL) + { + return scalar(layer)/nLayers_; + } + else + { + return + (1.0-pow(expansionRatio_, layer)) + / (1.0-pow(expansionRatio_, nLayers_)); + } +} + + // ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H index d8f3699fd5..ac96bbeea6 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H @@ -58,6 +58,8 @@ protected: const label nLayers_; + const scalar expansionRatio_; + const dictionary& dict_; const dictionary& coeffDict_; @@ -113,9 +115,15 @@ public: label nLayers() const; + scalar expansionRatio() const; + // Member Operators + //- Helper: calculate cumulative relative thickness for layer. + // (layer=0 -> 0; layer=nLayers -> 1) + scalar sumThickness(const label layer) const; + virtual point operator() ( const point& surfacePoint, diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C index c36098ab1e..cd42e91816 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C @@ -72,7 +72,8 @@ point linearNormal::operator() const label layer ) const { - scalar d = thickness_*layer/nLayers_; + //scalar d = thickness_*layer/nLayers_; + scalar d = thickness_*sumThickness(layer); return surfacePoint + d*surfaceNormal; } diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H index 949547311b..2dbcdccb16 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H @@ -64,7 +64,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary linearNormal(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C index e8d406b110..7d4fcff9ba 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C @@ -67,8 +67,7 @@ point linearRadial::operator() scalar rs = mag(surfacePoint); vector rsHat = surfacePoint/rs; - scalar delta = (R_ - rs)/nLayers_; - scalar r = rs + layer*delta; + scalar r = rs + (R_ - rs)*sumThickness(layer); return r*rsHat; } diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H index 0fda12e1b1..66758c1b12 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H @@ -61,7 +61,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary linearRadial(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C index fc6f040305..d35ad4e539 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C @@ -49,7 +49,13 @@ sigmaRadial::sigmaRadial(const dictionary& dict) RTbyg_(readScalar(coeffDict_.lookup("RTbyg"))), pRef_(readScalar(coeffDict_.lookup("pRef"))), pStrat_(readScalar(coeffDict_.lookup("pStrat"))) -{} +{ + if (mag(expansionRatio() - 1.0) > SMALL) + { + WarningIn("sigmaRadial::sigmaRadial(const dictionary&)") + << "Ignoring expansionRatio setting." << endl; + } +} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H index 525e89d11c..c74ad02457 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H @@ -63,7 +63,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary sigmaRadial(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C index 3d2c883ea1..07ddf34eaf 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C @@ -26,7 +26,7 @@ License #include "wedge.H" #include "addToRunTimeSelectionTable.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,7 +52,7 @@ wedge::wedge(const dictionary& dict) angle_ ( readScalar(coeffDict_.lookup("angle")) - *mathematicalConstant::pi/180.0 + *constant::math::pi/180.0 ) {} @@ -88,8 +88,8 @@ point wedge::operator() } else { - //sliceAngle = angle_*(layer + 1)/nLayers_; - sliceAngle = angle_*layer/nLayers_; + //sliceAngle = angle_*layer/nLayers_; + sliceAngle = angle_*sumThickness(layer); } // Find projection onto axis (or rather decompose surfacePoint diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H index 39c4b4a2d3..18707c32d1 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H @@ -77,7 +77,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary wedge(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties index 45618034de..e1283ab8b5 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties @@ -14,24 +14,28 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Where to get surface from: either from surface ('surface') or -// from (flipped) patch of existing case ('patch') -constructFrom patch; //surface; +// What to extrude: +// patch : from patch of another case ('sourceCase') +// mesh : as above but with original case included +// surface : from externally read surface -// If construct from (flipped) patch -sourceCase "$FOAM_RUN/icoFoam/cavity"; +//constructFrom mesh; +constructFrom patch; +//constructFrom surface; + +// If construct from patch/mesh: +sourceCase "../cavity"; sourcePatch movingWall; +// If construct from patch: patch to use for back (can be same as sourcePatch) +exposedPatchName movingWall; +// If construct from surface: +surface "movingWall.stl"; + + // Flip surface normals before usage. flipNormals false; -// If construct from surface -surface "movingWall.sMesh"; - - -// Do front and back need to be merged? Usually only makes sense for 360 -// degree wedges. -mergeFaces true; //- Linear extrusion in point-normal direction @@ -46,11 +50,13 @@ extrudeModel wedge; //- Extrudes into sphere with grading according to pressure (atmospherics) //extrudeModel sigmaRadial; -nLayers 20; +nLayers 10; + +expansionRatio 1.0; //0.9; wedgeCoeffs { - axisPt (0 0.1 0); + axisPt (0 0.1 -0.05); axis (-1 0 0); angle 360; // For nLayers=1 assume symmetry so angle/2 on each side } @@ -73,5 +79,10 @@ sigmaRadialCoeffs } +// Do front and back need to be merged? Usually only makes sense for 360 +// degree wedges. +mergeFaces true; + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index cd7c13d783..92c92d72e5 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -280,11 +280,14 @@ addLayersControls nBufferCellsNoExtrude 0; - // Overall max number of layer addition iterations + // Overall max number of layer addition iterations. The mesher will exit + // if it reaches this number of iterations; possibly with an illegal + // mesh. nLayerIter 50; // Max number of iterations after which relaxed meshQuality controls - // get used. + // get used. Up to nRelaxIter it uses the settings in meshQualityControls, + // after nRelaxIter it uses the values in meshQualityControls::relaxed. nRelaxedIter 20; } diff --git a/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C b/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C index 758592590c..3c48f6d7d7 100644 --- a/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C +++ b/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C @@ -33,7 +33,7 @@ Description #include "Time.H" #include "boundaryMesh.H" #include "repatchPolyTopoChanger.H" -#include "mathematicalConstants.H" +#include "mathConstants.H" #include "OFstream.H" #include "ListOps.H" @@ -93,7 +93,7 @@ int main(int argc, char *argv[]) scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); bool overwrite = args.optionFound("overwrite"); - scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0); + scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0); Info<< "Feature:" << featureAngle << endl << "minCos :" << minCos << endl diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C index 0ba1af5bea..c6fa201218 100644 --- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C +++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C @@ -125,6 +125,7 @@ label findPatchID(const polyMesh& mesh, const word& name) int main(int argc, char *argv[]) { +# include "addRegionOption.H" argList::validArgs.append("faceZone"); argList::validArgs.append("patch"); argList::validOptions.insert("additionalPatches", "(patch2 .. patchN)"); @@ -134,7 +135,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" runTime.functionObjects().off(); -# include "createMesh.H" +# include "createNamedMesh.H" const word oldInstance = mesh.pointsInstance(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); @@ -146,6 +147,14 @@ int main(int argc, char *argv[]) Info<< "Converting faces on zone " << zoneID.name() << " into baffles." << nl << endl; + if (zoneID.index() == -1) + { + FatalErrorIn(args.executable()) << "Cannot find faceZone " + << zoneID.name() << endl + << "Valid zones are " << faceZones.names() + << exit(FatalError); + } + const faceZone& fZone = faceZones[zoneID.index()]; Info<< "Found " << returnReduce(fZone.size(), sumOp