diff --git a/applications/solvers/lagrangian/coalChemistryFoam/Make/options b/applications/solvers/lagrangian/coalChemistryFoam/Make/options index b4758dace7..c6daa2b570 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options +++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options @@ -14,6 +14,7 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/surfaceFilmModels/lnInclude \ @@ -36,6 +37,7 @@ EXE_LIBS = \ -lsolidMixture \ -lthermophysicalFunctions \ -lreactionThermophysicalModels \ + -lSLGThermo \ -lchemistryModel \ -lradiation \ -lsurfaceFilmModels \ diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C index 65fb0b616b..5f373cc224 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C +++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C @@ -39,11 +39,12 @@ Description #include "hCombustionThermo.H" #include "turbulenceModel.H" #include "basicThermoCloud.H" -#include "CoalCloud.H" +#include "coalCloud.H" #include "psiChemistryModel.H" #include "chemistrySolver.H" #include "timeActivatedExplicitSource.H" #include "radiationModel.H" +#include "SLGThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H index dbd7cf9659..5dcfe1df4f 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H @@ -1,11 +1,11 @@ Info<< "\nConstructing coal cloud" << endl; -thermoCoalCloud coalParcels +coalCloud coalParcels ( "coalCloud1", rho, U, g, - thermo + slgThermo ); Info<< "\nConstructing limestone cloud" << endl; @@ -15,5 +15,5 @@ basicThermoCloud limestoneParcels rho, U, g, - thermo + slgThermo ); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H index ea325dabe0..97e409215c 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H @@ -8,6 +8,8 @@ hsCombustionThermo& thermo = chemistry.thermo(); + SLGThermo slgThermo(mesh, thermo); + basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options index 3639d884a2..6e117fc63b 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options @@ -14,6 +14,7 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/ODE/lnInclude \ @@ -35,6 +36,7 @@ EXE_LIBS = \ -lsolidMixture \ -lthermophysicalFunctions \ -lreactionThermophysicalModels \ + -lSLGThermo \ -lchemistryModel \ -lradiation \ -lODE \ diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H index 74a66b63e8..954b74e069 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H @@ -1,9 +1,9 @@ Info<< "\nConstructing reacting cloud" << endl; -icoPoly8ThermoReactingMultiphaseCloud parcels +basicReactingMultiphaseCloud parcels ( "reactingCloud1", rho, U, g, - thermo + slgThermo ); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H index 3ad06c22b3..22d7c1f219 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H @@ -8,6 +8,8 @@ hsReactionThermo& thermo = chemistry.thermo(); + SLGThermo slgThermo(mesh, thermo); + basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C index c16ba1ec1a..35dfc7e616 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C @@ -33,7 +33,6 @@ Description - reacting multiphase parcel cloud - porous media - mass, momentum and energy sources - - polynomial based, incompressible thermodynamics (f(T)) Note: ddtPhiCorr not used here when porous zones are active - not well defined for porous calculations @@ -43,12 +42,13 @@ Description #include "fvCFD.H" #include "hReactionThermo.H" #include "turbulenceModel.H" -#include "BasicReactingMultiphaseCloud.H" +#include "basicReactingMultiphaseCloud.H" #include "rhoChemistryModel.H" #include "chemistrySolver.H" #include "radiationModel.H" #include "porousZones.H" #include "timeActivatedExplicitSource.H" +#include "SLGThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options index 8b15c3e9ae..230a08dfad 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options @@ -1,5 +1,3 @@ -DEV_PATH=./../.. - EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ @@ -13,6 +11,7 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/surfaceFilmModels/lnInclude \ @@ -31,6 +30,7 @@ EXE_LIBS = \ -lsolidMixture \ -lthermophysicalFunctions \ -lreactionThermophysicalModels \ + -lSLGThermo \ -lchemistryModel \ -lradiation \ -lsurfaceFilmModels \ diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H index f3e043143c..c568be12a1 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H @@ -1,9 +1,9 @@ Info<< "\nConstructing reacting cloud" << endl; -thermoReactingCloud parcels +basicReactingCloud parcels ( "reactingCloud1", rho, U, g, - thermo + slgThermo ); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H index 909d0351f7..2e3208b0ee 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H @@ -8,6 +8,8 @@ hsCombustionThermo& thermo = chemistry.thermo(); + SLGThermo slgThermo(mesh, thermo); + basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C index f39d8f070a..f788da57e4 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C @@ -33,11 +33,12 @@ Description #include "fvCFD.H" #include "hCombustionThermo.H" #include "turbulenceModel.H" -#include "BasicReactingCloud.H" +#include "basicReactingCloud.H" #include "surfaceFilmModel.H" #include "psiChemistryModel.H" #include "chemistrySolver.H" #include "radiationModel.H" +#include "SLGThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options index 33d575609b..3954a34ea3 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options +++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options @@ -13,6 +13,7 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/ODE/lnInclude \ @@ -34,6 +35,7 @@ EXE_LIBS = \ -lsolidMixture \ -lthermophysicalFunctions \ -lreactionThermophysicalModels \ + -lSLGThermo \ -lchemistryModel \ -lradiation \ -lODE \ diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H index f3e043143c..c568be12a1 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H @@ -1,9 +1,9 @@ Info<< "\nConstructing reacting cloud" << endl; -thermoReactingCloud parcels +basicReactingCloud parcels ( "reactingCloud1", rho, U, g, - thermo + slgThermo ); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H index 0fd631fb7f..089489a014 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H @@ -8,6 +8,8 @@ hsCombustionThermo& thermo = chemistry.thermo(); + SLGThermo slgThermo(mesh, thermo); + basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C index e9f33f6229..38a1bf908e 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C @@ -33,10 +33,11 @@ Description #include "fvCFD.H" #include "hCombustionThermo.H" #include "turbulenceModel.H" -#include "BasicReactingCloud.H" +#include "basicReactingCloud.H" #include "psiChemistryModel.H" #include "chemistrySolver.H" #include "radiationModel.H" +#include "SLGThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H index 9f0417045b..c32855904a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H @@ -4,12 +4,24 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime); { { volTensorField gradUaT = fvc::grad(Ua)().T(); + + if (kineticTheory.on()) + { + kineticTheory.solve(gradUaT); + nuEffa = kineticTheory.mua()/rhoa; + } + else // If not using kinetic theory is using Ct model + { + nuEffa = sqr(Ct)*nutb + nua; + } + volTensorField Rca ( "Rca", ((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT ); + if (kineticTheory.on()) { Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H new file mode 100644 index 0000000000..dab75ee07d --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H @@ -0,0 +1,62 @@ +if(turbulence) +{ + if (mesh.changing()) + { + y.correct(); + } + + tmp tgradUb = fvc::grad(Ub); + volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb()))); + tgradUb.clear(); + + #include "wallFunctions.H" + + // Dissipation equation + fvScalarMatrix epsEqn + ( + fvm::ddt(beta, epsilon) + + fvm::div(phib, epsilon) + - fvm::laplacian + ( + alphaEps*nuEffb, epsilon, + "laplacian(DepsilonEff,epsilon)" + ) + == + C1*beta*G*epsilon/k + - fvm::Sp(C2*beta*epsilon/k, epsilon) + ); + + #include "wallDissipation.H" + + epsEqn.relax(); + epsEqn.solve(); + + epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15)); + + + // Turbulent kinetic energy equation + fvScalarMatrix kEqn + ( + fvm::ddt(beta, k) + + fvm::div(phib, k) + - fvm::laplacian + ( + alphak*nuEffb, k, + "laplacian(DkEff,k)" + ) + == + beta*G + - fvm::Sp(beta*epsilon/k, k) + ); + kEqn.relax(); + kEqn.solve(); + + k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8)); + + //- Re-calculate turbulence viscosity + nutb = Cmu*sqr(k)/epsilon; + + #include "wallViscosity.H" +} + +nuEffb = nutb + nub; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C index af326ffd5b..bcd41ac268 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C @@ -102,7 +102,7 @@ Foam::tmp Foam::JohnsonJacksonFrictionalStress::muf const volScalarField& alpha, const dimensionedScalar& alphaMax, const volScalarField& pf, - const volTensorField& D, + const volSymmTensorField& D, const dimensionedScalar& phi ) const { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H index 08733247c5..d758f2f85c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H @@ -93,7 +93,7 @@ public: const volScalarField& alpha, const dimensionedScalar& alphaMax, const volScalarField& pf, - const volTensorField& D, + const volSymmTensorField& D, const dimensionedScalar& phi ) const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C index dbc658f1d4..0f39146345 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C @@ -99,7 +99,7 @@ Foam::tmp Foam::SchaefferFrictionalStress::muf const volScalarField& alpha, const dimensionedScalar& alphaMax, const volScalarField& pf, - const volTensorField& D, + const volSymmTensorField& D, const dimensionedScalar& phi ) const { @@ -124,9 +124,9 @@ Foam::tmp Foam::SchaefferFrictionalStress::muf volScalarField& muff = tmuf(); - forAll(D, celli) + forAll (D, celli) { - if (alpha[celli] > alphaMax.value()-5e-2) + if (alpha[celli] > alphaMax.value() - 5e-2) { muff[celli] = 0.5*pf[celli]*sin(phi.value()) diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H index e9b0c4efd7..89b0e25ded 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H @@ -93,7 +93,7 @@ public: const volScalarField& alpha, const dimensionedScalar& alphaMax, const volScalarField& pf, - const volTensorField& D, + const volSymmTensorField& D, const dimensionedScalar& phi ) const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H index 8176e6aac2..b496ee9cf6 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H @@ -48,7 +48,7 @@ namespace Foam class frictionalStressModel { - // Private Member Functions + // Private member functions //- Disallow default bitwise copy construct frictionalStressModel(const frictionalStressModel&); @@ -127,7 +127,7 @@ public: const volScalarField& alpha, const dimensionedScalar& alphaMax, const volScalarField& pf, - const volTensorField& D, + const volSymmTensorField& D, const dimensionedScalar& phi ) const = 0; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index abdd6080be..ee2300d91e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -56,12 +56,13 @@ Foam::kineticTheoryModel::kineticTheoryModel "kineticTheoryProperties", Ua_.time().constant(), Ua_.mesh(), - IOobject::MUST_READ_IF_MODIFIED, + IOobject::MUST_READ, IOobject::NO_WRITE ) ), kineticTheory_(kineticTheoryProperties_.lookup("kineticTheory")), equilibrium_(kineticTheoryProperties_.lookup("equilibrium")), + viscosityModel_ ( kineticTheoryModels::viscosityModel::New @@ -192,24 +193,19 @@ Foam::kineticTheoryModel::~kineticTheoryModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::kineticTheoryModel::solve() +void Foam::kineticTheoryModel::solve(const volTensorField& gradUat) { if (!kineticTheory_) { return; } - word scheme("div(phi,Theta)"); - - volScalarField alpha = alpha_; - alpha.max(1.0e-6); const scalar sqrtPi = sqrt(constant::mathematical::pi); surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_); - volTensorField dU = fvc::grad(Ua_); - volTensorField dUT = dU.T(); - volTensorField D = 0.5*(dU + dUT); + volTensorField dU = gradUat.T();//fvc::grad(Ua_); + volSymmTensorField D = symm(dU); // NB, drag = K*alpha*beta, // (the alpha and beta has been extracted from the drag function for @@ -220,45 +216,52 @@ void Foam::kineticTheoryModel::solve() // Calculating the radial distribution function (solid volume fraction is // limited close to the packing limit, but this needs improvements) // The solution is higly unstable close to the packing limit. - gs0_ = radialModel_->g0(min(alpha, alphaMax_-1.0e-2), alphaMax_); + gs0_ = radialModel_->g0 + ( + min(max(alpha_, 1e-6), alphaMax_ - 0.01), + alphaMax_ + ); // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45) - volScalarField PsCoeff = - granularPressureModel_->granularPressureCoeff(alpha_,gs0_,rhoa_,e_ ); + volScalarField PsCoeff = granularPressureModel_->granularPressureCoeff + ( + alpha_, + gs0_, + rhoa_, + e_ + ); + + // 'thermal' conductivity (Table 3.3, p. 49) + kappa_ = conductivityModel_->kappa(alpha_, Theta_, gs0_, rhoa_, da_, e_); + + // particle viscosity (Table 3.2, p.47) + mua_ = viscosityModel_->mua(alpha_, Theta_, gs0_, rhoa_, da_, e_); dimensionedScalar Tsmall ( "small", - dimensionSet(0,2,-2,0,0,0,0), + dimensionSet(0 , 2 ,-2 ,0 , 0, 0, 0), 1.0e-6 ); dimensionedScalar TsmallSqrt = sqrt(Tsmall); volScalarField ThetaSqrt = sqrt(Theta_); - // 'thermal' conductivity (Table 3.3, p. 49) - kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rhoa_, da_, e_); - - // particle viscosity (Table 3.2, p.47) - mua_ = viscosityModel_->mua(alpha, Theta_, gs0_, rhoa_, da_, e_); - // dissipation (Eq. 3.24, p.50) volScalarField gammaCoeff = - 12.0*(1.0 - e_*e_)*sqr(alpha)*rhoa_*gs0_*(1.0/da_) - *ThetaSqrt/sqrtPi; + 12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi; // Eq. 3.25, p. 50 Js = J1 - J2 volScalarField J1 = 3.0*betaPrim; volScalarField J2 = 0.25*sqr(betaPrim)*da_*sqr(Ur) - /(alpha*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt)); + /(max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt)); // bulk viscosity p. 45 (Lun et al. 1984). lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi; - // stress tensor, Definitions, Table 3.1, p. 43 - volTensorField tau = 2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I; + volSymmTensorField tau = 2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I; if (!equilibrium_) { @@ -268,8 +271,8 @@ void Foam::kineticTheoryModel::solve() // wrong sign infront of laplacian fvScalarMatrix ThetaEqn ( - fvm::ddt(1.5*alpha*rhoa_, Theta_) - + fvm::div(phi, Theta_, scheme) + fvm::ddt(1.5*alpha_*rhoa_, Theta_) + + fvm::div(phi, Theta_, "div(phi,Theta)") == fvm::SuSp(-((PsCoeff*I) && dU), Theta_) + (tau && dU) @@ -290,33 +293,31 @@ void Foam::kineticTheoryModel::solve() volScalarField K3 = 0.5*da_*rhoa_* ( (sqrtPi/(3.0*(3.0-e_))) - *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_) - + 1.6*alpha*gs0_*(1.0 + e_)/sqrtPi + *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_) + +1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi ); volScalarField K2 = - 4.0*da_*rhoa_*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0; + 4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0; - volScalarField K4 = 12.0*(1.0 - e_*e_)*rhoa_*gs0_/(da_*sqrtPi); + volScalarField K4 = 12.0*(1.0 - sqr(e_))*rhoa_*gs0_/(da_*sqrtPi); volScalarField trD = tr(D); - volTensorField D2 = D & D; - volScalarField tr2D = trD*trD; - volScalarField trD2 = tr(D2); + volScalarField tr2D = sqr(trD); + volScalarField trD2 = tr(D & D); - volScalarField t1 = K1*alpha + rhoa_; + volScalarField t1 = K1*alpha_ + rhoa_; volScalarField l1 = -t1*trD; volScalarField l2 = sqr(t1)*tr2D; - volScalarField l3 = 4.0*K4*alpha*(2.0*K3*trD2 + K2*tr2D); + volScalarField l3 = 4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D); - Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha + 1.0e-4)*K4)); + Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4)); } Theta_.max(1.0e-15); Theta_.min(1.0e+3); - volScalarField pf = - frictionalStressModel_->frictionalPressure + volScalarField pf = frictionalStressModel_->frictionalPressure ( alpha_, alphaMinFriction_, @@ -344,13 +345,11 @@ void Foam::kineticTheoryModel::solve() phi_ ); - // add frictional stress for alpha > alphaMinFriction - mua_ = viscosityModel_->mua(alpha, Theta_, gs0_, rhoa_, da_, e_) + muf; + // add frictional stress + mua_ += muf; mua_.min(1.0e+2); mua_.max(0.0); - lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi; - Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl; volScalarField ktn = mua_/rhoa_; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H index 2659d8c1d6..2863c2d038 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H @@ -156,7 +156,7 @@ public: // Member Functions - void solve(); + void solve(const volTensorField& gradUat); bool on() const { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index b2bc0d817d..39a133d39b 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -5,6 +5,9 @@ volScalarField rUaA = 1.0/UaEqn.A(); volScalarField rUbA = 1.0/UbEqn.A(); + phia == (fvc::interpolate(Ua) & mesh.Sf()); + phib == (fvc::interpolate(Ub) & mesh.Sf()); + rUaAf = fvc::interpolate(rUaA); surfaceScalarField rUbAf = fvc::interpolate(rUbA); @@ -47,11 +50,10 @@ surfaceScalarField Dp ( - "(rho*(1|A(U)))", - alphaf*rUaAf/rhoa + betaf*rUbAf/rhob + "(rho*(1|A(U)))", alphaf*rUaAf/rhoa + betaf*rUbAf/rhob ); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index b3e0eb0009..bd0310765d 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C @@ -92,11 +92,6 @@ int main(int argc, char *argv[]) #include "kEpsilon.H" - if (kineticTheory.on()) - { - kineticTheory.solve(); - nuEffa += kineticTheory.mua()/rhoa; - } #include "write.H" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L index 6c3a0c315a..5fd94ef38d 100644 --- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L +++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L @@ -321,7 +321,7 @@ endOfSection {space}")"{space} pointGroupEndIndex.append(strtol(endPtr, &endPtr, 16) - 1); // point group type skipped - strtol(endPtr, &endPtr, 16); + (void)strtol(endPtr, &endPtr, 16); pointi = pointGroupStartIndex.last(); @@ -435,7 +435,7 @@ endOfSection {space}")"{space} faceGroupEndIndex.append(strtol(endPtr, &endPtr, 16) - 1); // face group type - strtol(endPtr, &endPtr, 16); + (void)strtol(endPtr, &endPtr, 16); faceGroupElementType = strtol(endPtr, &endPtr, 16); @@ -583,7 +583,7 @@ endOfSection {space}")"{space} cellGroupType.append(strtol(endPtr, &endPtr, 16)); // Note. Potentially skip cell set if type is zero. - strtol(endPtr, &endPtr, 16); + (void)strtol(endPtr, &endPtr, 16); Info<< "CellGroup: " << cellGroupZoneID.last() diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options index c00ca0534b..7c1ab5a64a 100644 --- a/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options +++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options @@ -1,11 +1,9 @@ EXE_INC = \ - -InonuniformTransformCyclic/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude EXE_LIBS = \ - -lnonuniformTransformCyclic \ -lfiniteVolume \ -lmeshTools \ -ldynamicMesh diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C index 49fd63f8e1..fb019694da 100644 --- a/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C +++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C @@ -955,6 +955,12 @@ int main(int argc, char *argv[]) argList::validArgs.append("faceZones"); argList::validArgs.append("thickness"); + Foam::argList::addBoolOption + ( + "oneD", + "generate columns of 1D cells" + ); + #include "addRegionOption.H" #include "addOverwriteOption.H" #include "setRootCase.H" @@ -966,6 +972,7 @@ int main(int argc, char *argv[]) const wordList zoneNames(IStringStream(args.additionalArgs()[1])()); scalar thickness = readScalar(IStringStream(args.additionalArgs()[2])()); bool overwrite = args.optionFound("overwrite"); + bool oneD = args.optionFound("oneD"); Info<< "Extruding zones " << zoneNames @@ -1225,9 +1232,25 @@ int main(int argc, char *argv[]) label nSide = 0; forAll(zoneSidePatch, zoneI) { - if (zoneSidePatch[zoneI] > 0) + if (oneD) + { + // Always add empty patches, one per zone. + word patchName = faceZones[zoneI].name() + "_" + "side"; + + zoneSidePatch[zoneI] = addPatch + ( + mesh, + patchName + ); + + Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl; + + nSide++; + } + else if (zoneSidePatch[zoneI] > 0) { word patchName = faceZones[zoneI].name() + "_" + "side"; + zoneSidePatch[zoneI] = addPatch ( mesh, @@ -1257,47 +1280,53 @@ int main(int argc, char *argv[]) ); label nInter = 0; - forAll(zoneZonePatch_min, minZone) + if (!oneD) { - for (label maxZone = minZone; maxZone < faceZones.size(); maxZone++) + forAll(zoneZonePatch_min, minZone) { - label index = minZone*faceZones.size()+maxZone; - - if (zoneZonePatch_min[index] > 0) + for (label maxZone = minZone; maxZone < faceZones.size(); maxZone++) { - word minToMax = - faceZones[minZone].name() - + "_to_" - + faceZones[maxZone].name(); - word maxToMin = - faceZones[maxZone].name() - + "_to_" - + faceZones[minZone].name(); - { - transformDict.set("neighbourPatch", maxToMin); - zoneZonePatch_min[index] = - addPatch - ( - mesh, - minToMax, - transformDict - ); - Info<< zoneZonePatch_min[index] << '\t' << minToMax << nl; - nInter++; - } - { - transformDict.set("neighbourPatch", minToMax); - zoneZonePatch_max[index] = - addPatch - ( - mesh, - maxToMin, - transformDict - ); - Info<< zoneZonePatch_max[index] << '\t' << maxToMin << nl; - nInter++; - } + label index = minZone*faceZones.size()+maxZone; + if (zoneZonePatch_min[index] > 0) + { + word minToMax = + faceZones[minZone].name() + + "_to_" + + faceZones[maxZone].name(); + word maxToMin = + faceZones[maxZone].name() + + "_to_" + + faceZones[minZone].name(); + + { + transformDict.set("neighbourPatch", maxToMin); + zoneZonePatch_min[index] = + addPatch + ( + mesh, + minToMax, + transformDict + ); + Info<< zoneZonePatch_min[index] << '\t' << minToMax + << nl; + nInter++; + } + { + transformDict.set("neighbourPatch", minToMax); + zoneZonePatch_max[index] = + addPatch + ( + mesh, + maxToMin, + transformDict + ); + Info<< zoneZonePatch_max[index] << '\t' << maxToMin + << nl; + nInter++; + } + + } } } } @@ -1323,7 +1352,16 @@ int main(int argc, char *argv[]) labelList& ePatches = extrudeEdgePatches[edgeI]; - if (eFaces.size() == 2) + if (oneD) + { + nonManifoldEdge[edgeI] = 1; + ePatches.setSize(eFaces.size()); + forAll(eFaces, i) + { + ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]]; + } + } + else if (eFaces.size() == 2) { label zone0 = zoneID[eFaces[0]]; label zone1 = zoneID[eFaces[1]]; diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/files b/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/files deleted file mode 100644 index ffacacb4da..0000000000 --- a/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/files +++ /dev/null @@ -1,8 +0,0 @@ -fvPatches/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicFvPatch.C -pointPatches/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicPointPatch.C -polyPatches/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicPolyPatch.C -pointPatchFields/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicPointPatchFields.C -fvsPatchFields/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicFvsPatchFields.C -fvPatchFields/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicFvPatchFields.C - -LIB = $(FOAM_LIBBIN)/libnonuniformTransformCyclic diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/options b/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/options deleted file mode 100644 index 71b7873964..0000000000 --- a/applications/utilities/mesh/generation/extrudeToRegionMesh/nonuniformTransformCyclic/Make/options +++ /dev/null @@ -1,5 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude - -LIB_LIBS = \ - -lfiniteVolume diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index 00913b8817..4c411fb150 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -381,6 +381,7 @@ int main(int argc, char *argv[]) "write cellMap, faceMap, pointMap in polyMesh/" ); +# include "addRegionOption.H" # include "addOverwriteOption.H" # include "addTimeOptions.H" @@ -396,7 +397,7 @@ int main(int argc, char *argv[]) runTime.setTime(Times[startTime], startTime); -# include "createMesh.H" +# include "createNamedMesh.H" const word oldInstance = mesh.pointsInstance(); const bool blockOrder = args.optionFound("blockOrder"); diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files index 23c2c6738c..2c10cfe303 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files +++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files @@ -1,9 +1,10 @@ +surfaceMeshWriter.C + foamToVTK.C internalWriter.C lagrangianWriter.C patchWriter.C writeFuns.C -writePatchGeom.C writeFaceSet.C writePointSet.C writeSurfFields.C diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C index ad67095ad9..272307588a 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C +++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C @@ -157,7 +157,7 @@ Note #include "writeFaceSet.H" #include "writePointSet.H" -#include "writePatchGeom.H" +#include "surfaceMeshWriter.H" #include "writeSurfFields.H" @@ -963,20 +963,42 @@ int main(int argc, char *argv[]) if (doFaceZones) { + PtrList ssf; + readFields + ( + vMesh, + vMesh.baseMesh(), + objects, + selectedFields, + ssf + ); + print(" surfScalarFields :", Info, ssf); + + PtrList svf; + readFields + ( + vMesh, + vMesh.baseMesh(), + objects, + selectedFields, + svf + ); + print(" surfVectorFields :", Info, svf); + const faceZoneMesh& zones = mesh.faceZones(); forAll(zones, zoneI) { - const faceZone& pp = zones[zoneI]; + const faceZone& fz = zones[zoneI]; - mkDir(fvPath/pp.name()); + mkDir(fvPath/fz.name()); fileName patchFileName; if (vMesh.useSubMesh()) { patchFileName = - fvPath/pp.name()/cellSetName + fvPath/fz.name()/cellSetName + "_" + timeDesc + ".vtk"; @@ -984,7 +1006,7 @@ int main(int argc, char *argv[]) else { patchFileName = - fvPath/pp.name()/pp.name() + fvPath/fz.name()/fz.name() + "_" + timeDesc + ".vtk"; @@ -992,18 +1014,31 @@ int main(int argc, char *argv[]) Info<< " FaceZone : " << patchFileName << endl; - std::ofstream ostr(patchFileName.c_str()); - - writeFuns::writeHeader(ostr, binary, pp.name()); - ostr<< "DATASET POLYDATA" << std::endl; - - writePatchGeom + indirectPrimitivePatch pp ( - binary, - pp().localFaces(), - pp().localPoints(), - ostr + IndirectList(mesh.faces(), fz), + mesh.points() ); + + surfaceMeshWriter writer + ( + vMesh, + binary, + pp, + fz.name(), + patchFileName + ); + + // Number of fields + writeFuns::writeCellDataHeader + ( + writer.os(), + pp.size(), + ssf.size()+svf.size() + ); + + writer.write(ssf); + writer.write(svf); } } @@ -1171,7 +1206,11 @@ int main(int argc, char *argv[]) + "_" + procFile.name() ); - system(cmd.c_str()); + if (system(cmd.c_str()) == -1) + { + WarningIn(args.executable()) + << "Could not execute command " << cmd << endl; + } } } } diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriter.C similarity index 59% rename from applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C rename to applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriter.C index 0dd2bd9200..d947eb3de3 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C +++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriter.C @@ -23,60 +23,63 @@ License \*---------------------------------------------------------------------------*/ -#include "writePatchGeom.H" -#include "OFstream.H" -#include "floatScalar.H" +#include "surfaceMeshWriter.H" #include "writeFuns.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -namespace Foam -{ - -// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // - -void writePatchGeom +Foam::surfaceMeshWriter::surfaceMeshWriter ( + const vtkMesh& vMesh, const bool binary, - const faceList& faces, - const pointField& points, - std::ofstream& ostr + const indirectPrimitivePatch& pp, + const word& name, + const fileName& fName ) +: + vMesh_(vMesh), + binary_(binary), + pp_(pp), + fName_(fName), + os_(fName.c_str()) { - ostr<< "POINTS " << points.size() << " float" << std::endl; - - DynamicList ptField(3*points.size()); - - writeFuns::insert(points, ptField); - - writeFuns::write(ostr, binary, ptField); + // Write header + writeFuns::writeHeader(os_, binary_, name); + os_ << "DATASET POLYDATA" << std::endl; + // Write topology label nFaceVerts = 0; - forAll(faces, faceI) + forAll(pp, faceI) { - nFaceVerts += faces[faceI].size() + 1; + nFaceVerts += pp[faceI].size() + 1; } - ostr<< "POLYGONS " << faces.size() << ' ' << nFaceVerts << std::endl; + os_ << "POINTS " << pp.nPoints() << " float" << std::endl; + + DynamicList ptField(3*pp.nPoints()); + writeFuns::insert(pp.localPoints(), ptField); + writeFuns::write(os_, binary, ptField); + + + os_ << "POLYGONS " << pp.size() << ' ' << nFaceVerts << std::endl; DynamicList