Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/chemFoam/solveChemistry.H

@@ -1,7 +1,3 @@
-    dtChem = chemistry.solve
-    (
-        runTime.value() - runTime.deltaT().value(),
-        runTime.deltaT().value()
-    );
+    dtChem = chemistry.solve(runTime.deltaT().value());
     scalar Sh = chemistry.Sh()()[0]/rho[0];
     integratedHeat += Sh*runTime.deltaT().value();

applications/solvers/combustion/reactingFoam/Allwmake

@@ -5,5 +5,6 @@ set -x
 wmake
 wmake rhoReactingFoam
 wmake rhoReactingBuoyantFoam
+wmake LTSReactingFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -16,7 +16,6 @@
           : -dpdt
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
-//    - fvm::laplacian(turbulence->muEff(), he)  // unit lewis no.
      ==
         reaction->Sh()
       + fvOptions(rho, he)

applications/solvers/combustion/reactingFoam/LTSReactingFoam/LTSReactingFoam.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    LTSReactingFoam
+
+Description
+    Local time stepping (LTS) solver for steady, compressible, laminar or
+    turbulent reacting and non-reacting flow.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulenceModel.H"
+#include "psiCombustionModel.H"
+#include "multivariateScheme.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialrDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "setrDeltaT.H"
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/files

@@ -0,0 +1,3 @@
+LTSReactingFoam.C
+
+EXE = $(FOAM_APPBIN)/LTSReactingFoam

applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/options

@@ -0,0 +1,28 @@
+EXE_INC = -ggdb3 \
+    -I.. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lreactionThermophysicalModels \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lcombustionModels

applications/solvers/combustion/reactingFoam/LTSReactingFoam/readTimeControls.H

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// Maximum flow Courant number
+scalar maxCo(readScalar(pimple.dict().lookup("maxCo")));
+
+// Maximum time scale
+scalar maxDeltaT(pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT));
+
+// Smoothing parameter (0-1) when smoothing iterations > 0
+scalar rDeltaTSmoothingCoeff
+(
+    pimple.dict().lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 1)
+);
+
+// Damping coefficient (1-0)
+scalar rDeltaTDampingCoeff
+(
+    pimple.dict().lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1)
+);
+
+// Maximum change in cell temperature per iteration (relative to previous value)
+scalar alphaTemp(pimple.dict().lookupOrDefault("alphaTemp", 0.05));
+
+
+// ************************************************************************* //

applications/solvers/combustion/reactingFoam/LTSReactingFoam/setInitialrDeltaT.H

@@ -0,0 +1,14 @@
+    volScalarField rDeltaT
+    (
+        IOobject
+        (
+            "rDeltaT",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless/dimTime, 1),
+        zeroGradientFvPatchScalarField::typeName
+    );

applications/solvers/combustion/reactingFoam/LTSReactingFoam/setrDeltaT.H

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    Info<< "Time scales min/max:" << endl;
+
+    // Cache old reciprocal time scale field
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Flow time scale
+    {
+        rDeltaT.dimensionedInternalField() =
+        (
+            fvc::surfaceSum(mag(phi))().dimensionedInternalField()
+           /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
+        );
+
+        // Limit the largest time scale
+        rDeltaT.max(1/maxDeltaT);
+
+        Info<< "    Flow        = "
+            << gMin(1/rDeltaT.internalField()) << ", "
+            << gMax(1/rDeltaT.internalField()) << endl;
+    }
+
+    // Reaction source time scale
+    if (alphaTemp < 1.0)
+    {
+        volScalarField::DimensionedInternalField rDeltaTT
+        (
+            mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T)
+        );
+
+        Info<< "    Temperature = "
+            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+
+        rDeltaT.dimensionedInternalField() = max
+        (
+            rDeltaT.dimensionedInternalField(),
+            rDeltaTT
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    // Spatially smooth the time scale field
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+    }
+
+    Info<< "    Overall     = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+}
+
+// ************************************************************************* //

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -115,5 +115,5 @@ K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
 {
-    dpdt = fvc::ddt(p);
+    dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);
 }

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso multiphaseMixtureThermo
+wclean
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso multiphaseMixtureThermo
+wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files

@@ -0,0 +1,3 @@
+compressibleMultiphaseInterFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleMultiphaseInterFoam

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -ImultiphaseMixtureThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmultiphaseMixtureThermo \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -linterfaceProperties \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T)
+      + fvm::div(multiphaseProperties.rhoPhi(), T)
+      - fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T)
+      + (
+            fvc::div(fvc::absolute(phi, U), p)
+          + fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K)
+        )*multiphaseProperties.rCv()
+    );
+
+    TEqn.relax();
+    TEqn.solve();
+
+    multiphaseProperties.correct();
+}

applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H

@@ -0,0 +1,27 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(multiphaseProperties.rhoPhi(), U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    multiphaseProperties.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        K = 0.5*magSqr(U);
+    }

applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        multiphaseProperties.nearInterface()().internalField()
+      * fvc::surfaceSum(mag(phi))().internalField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    compressibleMultiphaseInterFoam
+
+Description
+    Solver for n compressible, non-isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixtureThermo.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            multiphaseProperties.solve();
+
+            solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi()));
+
+            #include "UEqn.H"
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H

@@ -0,0 +1,71 @@
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+    Info<< "Constructing multiphaseMixtureThermo\n" << endl;
+    multiphaseMixtureThermo multiphaseProperties(U, phi);
+
+    volScalarField& p = multiphaseProperties.p();
+    volScalarField& T = multiphaseProperties.T();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        multiphaseProperties.rho()
+    );
+
+    Info<< max(rho) << min(rho);
+
+    dimensionedScalar pMin(multiphaseProperties.lookup("pMin"));
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    // Construct compressible turbulence model
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            multiphaseProperties.rhoPhi(),
+            multiphaseProperties
+        )
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files

@@ -0,0 +1,5 @@
+phaseModel/phaseModel.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixtureThermo.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseMixtureThermo

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixtureThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixtureThermo::interfacePair,
+        multiphaseMixtureThermo::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H

@@ -0,0 +1,434 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multiphaseMixtureThermo
+
+Description
+
+SourceFiles
+    multiphaseMixtureThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixtureThermo_H
+#define multiphaseMixtureThermo_H
+
+#include "phaseModel.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "rhoThermo.H"
+#include "psiThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class multiphaseMixtureThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixtureThermo
+:
+    public psiThermo
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phaseModel& alpha1, const phaseModel& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phaseModel> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        surfaceScalarField rhoPhi_;
+
+        volScalarField alphas_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2,
+            surfaceVectorField::GeometricBoundaryField& nHatb
+        ) const;
+
+        tmp<volScalarField> K
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseMixtureThermo");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixtureThermo
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixtureThermo()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phaseModel>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return non-const access to the phases
+        PtrDictionary<phaseModel>& phases()
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+        //- Update densities for given pressure change
+        void correctRho(const volScalarField& dp);
+
+        //- Return true if the equation of state is incompressible
+        //  i.e. rho != f(p)
+        virtual bool incompressible() const;
+
+        //- Return true if the equation of state is isochoric
+        //  i.e. rho = const
+        virtual bool isochoric() const;
+
+
+        // Access to thermodynamic state variables
+
+            //- Enthalpy/Internal energy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& he()
+            {
+                notImplemented("multiphaseMixtureThermo::he()");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy [J/kg]
+            virtual const volScalarField& he() const
+            {
+                notImplemented("multiphaseMixtureThermo::he() const");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy
+            //  for given pressure and temperature [J/kg]
+            virtual tmp<volScalarField> he
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) const;
+
+            //- Enthalpy/Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy/Internal energy for patch [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
+            //- Temperature from enthalpy/internal energy for cell-set
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const labelList& cells
+            ) const;
+
+            //- Temperature from enthalpy/internal energy for patch
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const label patchi
+            ) const;
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Density [kg/m^3]
+            virtual tmp<volScalarField> rho() const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- gamma = Cp/Cv []
+            virtual tmp<volScalarField> gamma() const;
+
+            //- gamma = Cp/Cv for patch []
+            virtual tmp<scalarField> gamma
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure/volume [J/kg/K]
+            virtual tmp<volScalarField> Cpv() const;
+
+            //- Heat capacity at constant pressure/volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity ratio []
+            virtual tmp<volScalarField> CpByCpv() const;
+
+            //- Heat capacity ratio for patch []
+            virtual tmp<scalarField> CpByCpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+
+        // Fields derived from transport state variables
+
+            //- Thermal diffusivity for temperature of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappa() const;
+
+            //- Thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappa
+            (
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> alphaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> alphaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+
+        //- Return the phase-averaged reciprocal Cv
+        tmp<volScalarField> rCv() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseModel::phaseModel
+(
+    const word& phaseName,
+    const volScalarField& p,
+    const volScalarField& T
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh()
+    ),
+    name_(phaseName),
+    p_(p),
+    T_(T),
+    thermo_(NULL),
+    dgdt_
+    (
+        IOobject
+        (
+            IOobject::groupName("dgdt", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh(),
+        dimensionedScalar("0", dimless/dimTime, 0)
+    )
+{
+    {
+        volScalarField Tp(IOobject::groupName("T", phaseName), T);
+        Tp.write();
+    }
+
+    thermo_ = rhoThermo::New(p.mesh(), phaseName);
+    thermo_->validate(phaseName, "e");
+
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::clone() const
+{
+    notImplemented("phaseModel::clone() const");
+    return autoPtr<phaseModel>(NULL);
+}
+
+
+void Foam::phaseModel::correct()
+{
+    thermo_->he() = thermo_->he(p_, T_);
+    thermo_->correct();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseModel
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phaseModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseModel_H
+#define phaseModel_H
+
+#include "rhoThermo.H"
+#include "volFields.H"
+#include "dictionaryEntry.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phaseModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseModel
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        const volScalarField& p_;
+        const volScalarField& T_;
+        autoPtr<rhoThermo> thermo_;
+        volScalarField dgdt_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        phaseModel
+        (
+            const word& phaseName,
+            const volScalarField& p,
+            const volScalarField& T
+        );
+
+        //- Return clone
+        autoPtr<phaseModel> clone() const;
+
+        //- Return a pointer to a new phaseModel created on freestore
+        //  from Istream
+        class iNew
+        {
+            const volScalarField& p_;
+            const volScalarField& T_;
+
+        public:
+
+            iNew
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            )
+            :
+                p_(p),
+                T_(T)
+            {}
+
+            autoPtr<phaseModel> operator()(Istream& is) const
+            {
+                return autoPtr<phaseModel>(new phaseModel(is, p_, T_));
+            }
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+
+        //- Return const-access to phase rhoThermo
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        //- Return access to phase rhoThermo
+        rhoThermo& thermo()
+        {
+            return thermo_();
+        }
+
+        //- Return const-access to phase divergence
+        const volScalarField& dgdt() const
+        {
+            return dgdt_;
+        }
+
+        //- Return access to phase divergence
+        volScalarField& dgdt()
+        {
+            return dgdt_;
+        }
+
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -0,0 +1,142 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+
+    surfaceScalarField phig
+    (
+        (
+            multiphaseProperties.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    PtrList<fvScalarMatrix> p_rghEqnComps(multiphaseProperties.phases().size());
+
+    label phasei = 0;
+    forAllConstIter
+    (
+        PtrDictionary<phaseModel>,
+        multiphaseProperties.phases(),
+        phase
+    )
+    {
+        const rhoThermo& thermo = phase().thermo();
+        const volScalarField& rho = thermo.rho()();
+
+        p_rghEqnComps.set
+        (
+            phasei,
+            (
+                fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
+              + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
+            ).ptr()
+        );
+
+        phasei++;
+    }
+
+    // Cache p_rgh prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        tmp<fvScalarMatrix> p_rghEqnComp;
+
+        phasei = 0;
+        forAllConstIter
+        (
+            PtrDictionary<phaseModel>,
+            multiphaseProperties.phases(),
+            phase
+        )
+        {
+            tmp<fvScalarMatrix> hmm
+            (
+                (max(phase(), scalar(0))/phase().thermo().rho())
+               *p_rghEqnComps[phasei]
+            );
+
+            if (phasei == 0)
+            {
+                p_rghEqnComp = hmm;
+            }
+            else
+            {
+                p_rghEqnComp() += hmm;
+            }
+
+            phasei++;
+        }
+
+        solve
+        (
+            p_rghEqnComp
+          + p_rghEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+            // p_rgh = p - multiphaseProperties.rho()*gh;
+
+            phasei = 0;
+            forAllIter
+            (
+                PtrDictionary<phaseModel>,
+                multiphaseProperties.phases(),
+                phase
+            )
+            {
+                phase().dgdt() =
+                    pos(phase())
+                  *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
+            }
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+
+    // Update densities from change in p_rgh
+    multiphaseProperties.correctRho(p_rgh - p_rgh_0);
+    rho = multiphaseProperties.rho();
+
+    K = 0.5*magSqr(U);
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+}

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -121,9 +121,11 @@
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
-        rDeltaT =
-             rDeltaT0
-            *max(rDeltaT/rDeltaT0, scalar(1.0) - rDeltaTDampingCoeff);
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
 
         Info<< "Damped flow time scale min/max = "
             << gMin(1/rDeltaT.internalField())

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -126,11 +126,25 @@ void Foam::multiphaseSystem::solveAlphas()
         }
 
         // Ensure that the flux at inflow BCs is preserved
-        phiAlphaCorr.boundaryField() = min
-        (
-            phase1.phi().boundaryField()*alpha1.boundaryField(),
-            phiAlphaCorr.boundaryField()
-        );
+        forAll(phiAlphaCorr.boundaryField(), patchi)
+        {
+            fvsPatchScalarField& phiAlphaCorrp =
+                phiAlphaCorr.boundaryField()[patchi];
+
+            if (!phiAlphaCorrp.coupled())
+            {
+                const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
+                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                forAll(phiAlphaCorrp, facei)
+                {
+                    if (phi1p[facei] < 0)
+                    {
+                        phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
+                    }
+                }
+            }
+        }
 
         MULES::limit
         (

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -2,6 +2,8 @@
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
+    alpha1.correctBoundaryConditions();
+
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
 
@@ -93,11 +95,25 @@
             );
 
             // Ensure that the flux at inflow BCs is preserved
-            alphaPhic1.boundaryField() = min
-            (
-                phi1.boundaryField()*alpha1.boundaryField(),
-                alphaPhic1.boundaryField()
-            );
+            forAll(alphaPhic1.boundaryField(), patchi)
+            {
+                fvsPatchScalarField& alphaPhic1p =
+                    alphaPhic1.boundaryField()[patchi];
+
+                if (!alphaPhic1p.coupled())
+                {
+                    const scalarField& phi1p = phi1.boundaryField()[patchi];
+                    const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                    forAll(alphaPhic1p, facei)
+                    {
+                        if (phi1p[facei] < 0)
+                        {
+                            alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
+                        }
+                    }
+                }
+            }
 
             MULES::explicitSolve
             (

applications/test/PtrList/Test-PtrList.C

@@ -32,6 +32,8 @@ Description
 #include "scalar.H"
 #include "IOstreams.H"
 #include "PtrList.H"
+#include "plane.H"
+#include "DynamicList.H"
 
 using namespace Foam;
 
@@ -76,6 +78,7 @@ int main(int argc, char *argv[])
 {
     PtrList<Scalar> list1(10);
     PtrList<Scalar> list2(15);
+    PtrList<Scalar> listApp;
 
     forAll(list1, i)
     {
@@ -87,9 +90,14 @@ int main(int argc, char *argv[])
         list2.set(i, new Scalar(10 + 1.3*i));
     }
 
+    for (label i = 0; i < 5; ++i)
+    {
+        listApp.append(new Scalar(1.3*i));
+    }
 
     Info<<"list1: " << list1 << endl;
     Info<<"list2: " << list2 << endl;
+    Info<<"listApp: " << listApp << endl;
 
     Info<<"indirectly delete some items via set(.., 0) :" << endl;
     for (label i = 0; i < 3; i++)
@@ -115,6 +123,13 @@ int main(int argc, char *argv[])
     Info<<"list2: " << list2 << endl;
     Info<<"list3: " << list3 << endl;
 
+    PtrList<plane> planes;
+    planes.append(new plane(vector::one, vector::one));
+    planes.append(new plane(vector(1,2,3), vector::one));
+
+    forAll(planes, p)
+        Info<< "plane " << planes[p] << endl;
+
     Info<< nl << "Done." << endl;
     return 0;
 }

applications/test/UniformField/Make/files

@@ -0,0 +1,3 @@
+Test-UniformField.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-UniformField

applications/test/UniformField/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */
+/* EXE_LIBS = -lfiniteVolume */

applications/test/UniformField/Test-UniformField.C

@@ -0,0 +1,15 @@
+#include "UniformField.H"
+#include "vector.H"
+#include "IOstreams.H"
+
+using namespace Foam;
+
+int main()
+{
+    UniformField<scalar> uf1(13.1);
+    UniformField<vector> uf2(vector(1, 2, 3));
+
+    Info<< "uf1 = " << uf1[22] << "; uf2 = " << uf2[1002] << endl;
+
+    return 0;
+}

applications/test/fieldMapping/Make/files

@@ -0,0 +1,7 @@
+/*
+fvMeshGeometry.C
+fvMesh.C
+*/
+Test-fieldMapping.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-fieldMapping

applications/test/fieldMapping/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -DFULLDEBUG -g -O0 \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -ldynamicMesh

applications/test/fieldMapping/README.txt

@@ -0,0 +1,8 @@
+Test application for volField and surfaceField mapping with topology
+changes.
+
+Run
+
+    pipe1D/Allrun
+
+to compile and map a few fields

applications/test/fieldMapping/Test-fieldMapping.C

@@ -0,0 +1,335 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-fieldMapping
+
+Description
+    Test app for mapping of fields.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "Time.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "removeFaces.H"
+#include "mapPolyMesh.H"
+#include "polyTopoChange.H"
+#include "fvcDiv.H"
+#include "Random.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool notEqual(const scalar s1, const scalar s2, const scalar tol)
+{
+    return mag(s1-s2) > tol;
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+    argList::validArgs.append("inflate (true|false)");
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    const Switch inflate(args.args()[1]);
+
+    if (inflate)
+    {
+        Info<< "Deleting cells using inflation/deflation" << nl << endl;
+    }
+    else
+    {
+        Info<< "Deleting cells, introducing points at new position" << nl
+            << endl;
+    }
+
+
+    Random rndGen(0);
+
+
+
+    // Test mapping
+    // ------------
+    // Mapping is volume averaged
+
+
+    // 1. uniform field stays uniform
+    volScalarField one
+    (
+        IOobject
+        (
+            "one",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+    Info<< "Writing one field "
+        << one.name() << " in " << runTime.timeName() << endl;
+    one.write();
+
+
+    // 2. linear profile gets preserved
+    volScalarField ccX
+    (
+        IOobject
+        (
+            "ccX",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh.C().component(0)
+    );
+    Info<< "Writing x component of cell centres to "
+        << ccX.name()
+        << " in " << runTime.timeName() << endl;
+    ccX.write();
+
+
+    // Uniform surface field
+    surfaceScalarField surfaceOne
+    (
+        IOobject
+        (
+            "surfaceOne",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 1.0),
+        calculatedFvsPatchScalarField::typeName
+    );
+    Info<< "Writing surface one field "
+        << surfaceOne.name() << " in " << runTime.timeName() << endl;
+    surfaceOne.write();
+
+
+
+
+    // Force allocation of V. Important for any mesh changes since otherwise
+    // old time volumes are not stored
+    const scalar totalVol = gSum(mesh.V());
+
+    // Face removal engine. No checking for not merging boundary faces.
+    removeFaces faceRemover(mesh, GREAT);
+
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (!mesh.nInternalFaces())
+        {
+            break;
+        }
+
+        // Remove face
+        label candidateFaceI = rndGen.integer(0, mesh.nInternalFaces()-1);
+        Info<< "Wanting to delete face " << mesh.faceCentres()[candidateFaceI]
+            << nl << endl;
+
+        labelList candidates(1, candidateFaceI);
+
+
+        // Get compatible set of faces and connected sets of cells.
+        labelList cellRegion;
+        labelList cellRegionMaster;
+        labelList facesToRemove;
+
+        faceRemover.compatibleRemoves
+        (
+            candidates,
+            cellRegion,
+            cellRegionMaster,
+            facesToRemove
+        );
+
+        // Topo changes container
+        polyTopoChange meshMod(mesh);
+
+        // Insert mesh refinement into polyTopoChange.
+        faceRemover.setRefinement
+        (
+            facesToRemove,
+            cellRegion,
+            cellRegionMaster,
+            meshMod
+        );
+
+        // Change mesh and inflate
+        Info<< "Actually changing mesh" << nl << endl;
+        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, inflate);
+
+        Info<< "Mapping fields" << nl << endl;
+        mesh.updateMesh(morphMap);
+
+        // Move mesh (since morphing does not do this)
+        if (morphMap().hasMotionPoints())
+        {
+            Info<< "Moving mesh" << nl << endl;
+            mesh.movePoints(morphMap().preMotionPoints());
+        }
+
+        // Update numbering of cells/vertices.
+        faceRemover.updateMesh(morphMap);
+
+
+        Info<< "Writing fields" << nl << endl;
+        runTime.write();
+
+
+        // Check mesh volume conservation
+        if (mesh.moving())
+        {
+            #include "volContinuity.H"
+        }
+        else
+        {
+            if (mesh.V().size() != mesh.nCells())
+            {
+                FatalErrorIn(args.executable())
+                    << "Volume not mapped. V:" << mesh.V().size()
+                    << " nCells:" << mesh.nCells()
+                    << exit(FatalError);
+            }
+
+            const scalar newVol = gSum(mesh.V());
+            Info<< "Initial volume = " << totalVol
+                << "  New volume = " << newVol
+                << endl;
+
+            if (mag(newVol-totalVol)/totalVol > 1e-10)
+            {
+                FatalErrorIn(args.executable())
+                    << "Volume loss: old volume:" << totalVol
+                    << "  new volume:" << newVol
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Volume check OK" << nl << endl;
+            }
+        }
+
+
+        // Check constant profile
+        {
+            const scalar max = gMax(one);
+            const scalar min = gMin(one);
+
+            Info<< "Uniform one field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10))
+            {
+                FatalErrorIn(args.executable())
+                    << "Uniform volVectorField not preserved."
+                    << " Min and max should both be 1.0. min:" << min
+                    << " max:" << max
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Uniform field mapping check OK" << nl << endl;
+            }
+        }
+
+        // Check linear profile
+        {
+            const scalarField diff = ccX-mesh.C().component(0);
+
+            const scalar max = gMax(diff);
+            const scalar min = gMin(diff);
+
+            Info<< "Linear profile field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 0.0, 1e-10) || notEqual(min, 0.0, 1e-10))
+            {
+                FatalErrorIn(args.executable())
+                    << "Linear profile not preserved."
+                    << " Min and max should both be 0.0. min:" << min
+                    << " max:" << max
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Linear profile mapping check OK" << nl << endl;
+            }
+        }
+
+        // Check face field mapping
+        if (surfaceOne.size())
+        {
+            const scalar max = gMax(surfaceOne.internalField());
+            const scalar min = gMin(surfaceOne.internalField());
+
+            Info<< "Uniform surface field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10))
+            {
+                FatalErrorIn(args.executable())
+                    << "Uniform surfaceScalarField not preserved."
+                    << " Min and max should both be 1.0. min:" << min
+                    << " max:" << max
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Uniform surfaceScalarField mapping check OK" << nl
+                    << endl;
+            }
+        }
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/Allrun

@@ -0,0 +1,23 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Get application name
+application=Test-fieldMapping
+
+# Compile
+wmake ..
+
+runApplication blockMesh
+
+# Run with inflation
+runApplication $application true
+mv "log.$application" "log.$application-inflate"
+
+# Run without inflation
+runApplication $application false
+
+
+# ----------------------------------------------------------------- end-of-file

applications/test/fieldMapping/pipe1D/constant/manualDecomposition

@@ -0,0 +1,424 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       labelList;
+    location    "constant";
+    object      cellDecomposition;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+400
+(
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+)
+
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/constant/polyMesh/blockMeshDict

@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 0.1 0)
+    (0 0.1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 0.1 0.1)
+    (0 0.1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 1 1) simpleGrading (10 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+);
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/constant/polyMesh/boundary

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+3
+(
+    inlet
+    {
+        type            patch;
+        nFaces          1;
+        startFace       19;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          1;
+        startFace       20;
+    }
+    defaultFaces
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          80;
+        startFace       21;
+    }
+)
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;
+
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/system/controlDict

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+//    primitiveMesh   1;
+//    polyMesh        1;
+//    fvMesh          1;
+}
+
+application     Test-fieldMapping;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/system/fvSchemes

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

applications/test/fieldMapping/pipe1D/system/fvSolution

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //

applications/test/hexRef8/Make/files

@@ -0,0 +1,3 @@
+Test-hexRef8.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-hexRef8

applications/test/hexRef8/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -DFULLDEBUG -g -O0 \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -ldynamicMesh

applications/test/hexRef8/README.txt

@@ -0,0 +1,11 @@
+Test application for volField and surfaceField mapping with
+refinement/unrefinement
+
+Run
+
+    block/Allrun
+
+to compile and map a few fields. Note that hexRef8 cannot be
+run in inflation mode - there is the problem of getting
+a set of faces to sweep that is consistent for a cell and
+all its neighbours.

applications/test/hexRef8/Test-hexRef8.C

@@ -0,0 +1,371 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-hexRef8
+
+Description
+    Test app for refinement and unrefinement. Runs a few iterations refining
+    and unrefining.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "hexRef8.H"
+#include "mapPolyMesh.H"
+#include "polyTopoChange.H"
+#include "Random.H"
+#include "zeroGradientFvPatchFields.H"
+#include "fvcDiv.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool notEqual(const scalar s1, const scalar s2, const scalar tol)
+{
+    return mag(s1-s2) > tol;
+}
+
+
+// Main program:
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+    argList::validArgs.append("inflate (true|false)");
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+
+    const Switch inflate(args.args()[1]);
+
+    if (inflate)
+    {
+        Info<< "Splitting/deleting cells using inflation/deflation" << nl
+            << endl;
+    }
+    else
+    {
+        Info<< "Splitting/deleting cells, introducing points at new position"
+            << nl << endl;
+    }
+
+
+    Random rndGen(0);
+
+
+    // Force generation of V()
+    (void)mesh.V();
+
+
+    // Test mapping
+    // ------------
+
+    // 1. uniform field stays uniform
+    volScalarField one
+    (
+        IOobject
+        (
+            "one",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+    Info<< "Writing one field "
+        << one.name() << " in " << runTime.timeName() << endl;
+    one.write();
+
+
+    // 2. linear profile gets preserved
+    volScalarField ccX
+    (
+        IOobject
+        (
+            "ccX",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh.C().component(0)
+    );
+    Info<< "Writing x component of cell centres to "
+        << ccX.name()
+        << " in " << runTime.timeName() << endl;
+    ccX.write();
+
+
+    // Uniform surface field
+    surfaceScalarField surfaceOne
+    (
+        IOobject
+        (
+            "surfaceOne",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless, 1.0),
+        calculatedFvsPatchScalarField::typeName
+    );
+    Info<< "Writing surface one field "
+        << surfaceOne.name() << " in " << runTime.timeName() << endl;
+    surfaceOne.write();
+
+
+    // Force allocation of V. Important for any mesh changes since otherwise
+    // old time volumes are not stored
+    const scalar totalVol = gSum(mesh.V());
+
+
+    // Construct refiner. Read initial cell and point levels.
+    hexRef8 meshCutter(mesh);
+
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (mesh.globalData().nTotalCells() == 0)
+        {
+            break;
+        }
+
+
+
+        // Mesh changing engine.
+        polyTopoChange meshMod(mesh);
+
+        if (rndGen.bit())
+        {
+            // Refine
+            label nRefine = mesh.nCells()/20;
+            DynamicList<label> refineCandidates(nRefine);
+
+            for (label i=0; i<nRefine; i++)
+            {
+                refineCandidates.append(rndGen.integer(0, mesh.nCells()-1));
+            }
+
+            labelList cellsToRefine
+            (
+                meshCutter.consistentRefinement
+                (
+                    refineCandidates,
+                    true                  // buffer layer
+                )
+            );
+            Info<< nl << "-- selected "
+                << returnReduce(cellsToRefine.size(), sumOp<label>())
+                << " cells out of " << mesh.globalData().nTotalCells()
+                << " for refinement" << endl;
+
+            // Play refinement commands into mesh changer.
+            meshCutter.setRefinement(cellsToRefine, meshMod);
+        }
+        else
+        {
+            // Unrefine
+            labelList allSplitPoints(meshCutter.getSplitPoints());
+
+            label nUnrefine = allSplitPoints.size()/20;
+            labelHashSet candidates(2*nUnrefine);
+
+            for (label i=0; i<nUnrefine; i++)
+            {
+                candidates.insert
+                (
+                    allSplitPoints[rndGen.integer(0, allSplitPoints.size()-1)]
+                );
+            }
+
+            labelList splitPoints = meshCutter.consistentUnrefinement
+            (
+                candidates.toc(),
+                false
+            );
+            Info<< nl << "-- selected "
+                << returnReduce(splitPoints.size(), sumOp<label>())
+                << " points out of "
+                << returnReduce(allSplitPoints.size(), sumOp<label>())
+                << " for unrefinement" << endl;
+
+            // Play refinement commands into mesh changer.
+            meshCutter.setUnrefinement(splitPoints, meshMod);
+        }
+
+
+
+
+        // Create mesh, return map from old to new mesh.
+        Info<< nl << "-- actually changing mesh" << endl;
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, inflate);
+
+        // Update fields
+        Info<< nl << "-- mapping mesh data" << endl;
+        mesh.updateMesh(map);
+
+        // Inflate mesh
+        if (map().hasMotionPoints())
+        {
+            Info<< nl << "-- moving mesh" << endl;
+            mesh.movePoints(map().preMotionPoints());
+        }
+
+        // Update numbering of cells/vertices.
+        Info<< nl << "-- mapping hexRef8 data" << endl;
+        meshCutter.updateMesh(map);
+
+
+
+        Info<< "Writing fields" << nl << endl;
+        runTime.write();
+
+
+
+        // Check mesh volume conservation
+        if (mesh.moving())
+        {
+            #include "volContinuity.H"
+        }
+        else
+        {
+            if (mesh.V().size() != mesh.nCells())
+            {
+                FatalErrorIn(args.executable())
+                    << "Volume not mapped. V:" << mesh.V().size()
+                    << " nCells:" << mesh.nCells()
+                    << exit(FatalError);
+            }
+
+            const scalar newVol = gSum(mesh.V());
+            Info<< "Initial volume = " << totalVol
+                << "  New volume = " << newVol
+                << endl;
+
+            if (mag(newVol-totalVol)/totalVol > 1e-10)
+            {
+                FatalErrorIn(args.executable())
+                    << "Volume loss: old volume:" << totalVol
+                    << "  new volume:" << newVol
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Volume check OK" << nl << endl;
+            }
+        }
+
+
+        // Check constant profile
+        {
+            const scalar max = gMax(one);
+            const scalar min = gMin(one);
+
+            Info<< "Uniform one field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10))
+            {
+                FatalErrorIn(args.executable())
+                    << "Uniform volVectorField not preserved."
+                    << " Min and max should both be 1.0. min:" << min
+                    << " max:" << max
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Uniform field mapping check OK" << nl << endl;
+            }
+        }
+
+        // Check linear profile
+        {
+            const scalarField diff = ccX-mesh.C().component(0);
+
+            const scalar max = gMax(diff);
+            const scalar min = gMin(diff);
+
+            Info<< "Linear profile field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 0.0, 1e-10) || notEqual(min, 0.0, 1e-10))
+            {
+                Info<< "Linear profile not preserved."
+                    << " Min and max should both be 0.0. min:" << min
+                    << " max:" << max << nl << endl;
+            }
+            else
+            {
+                Info<< "Linear profile mapping check OK" << nl << endl;
+            }
+        }
+
+        // Check face field mapping
+        if (surfaceOne.size())
+        {
+            const scalar max = gMax(surfaceOne.internalField());
+            const scalar min = gMin(surfaceOne.internalField());
+
+            Info<< "Uniform surface field min = " << min
+                << "  max = " << max << endl;
+
+            if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10))
+            {
+                FatalErrorIn(args.executable())
+                    << "Uniform surfaceScalarField not preserved."
+                    << " Min and max should both be 1.0. min:" << min
+                    << " max:" << max
+                    << exit(FatalError);
+            }
+            else
+            {
+                Info<< "Uniform surfaceScalarField mapping check OK" << nl
+                    << endl;
+            }
+        }
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/hexRef8/block/Allrun

@@ -0,0 +1,23 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Get application name
+application=Test-hexRef8
+
+# Compile
+runApplication wmake ..
+
+runApplication blockMesh
+
+# Run with inflation
+runApplication $application true
+mv "log.$application" "log.$application-inflate"
+
+# Run without inflation
+runApplication $application false
+
+
+# ----------------------------------------------------------------- end-of-file

applications/test/hexRef8/block/constant/manualDecomposition

@@ -0,0 +1,424 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       labelList;
+    location    "constant";
+    object      cellDecomposition;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+400
+(
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+)
+
+
+// ************************************************************************* //

applications/test/hexRef8/block/constant/polyMesh/blockMeshDict

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 1)
+    (1 0 1)
+    (1 1 1)
+    (0 1 1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (5 5 5) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    allWalls
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //

applications/test/hexRef8/block/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;
+
+
+// ************************************************************************* //

applications/test/hexRef8/block/system/controlDict

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+    primitiveMesh   1;
+    polyMesh        1;
+    fvMesh          1;
+}
+
+application     icoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         30;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

applications/test/hexRef8/block/system/decomposeParDict

@@ -0,0 +1,139 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  2;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+//- Keep owner and neighbour on same processor for faces in patches:
+//  (makes sense only for cyclic patches)
+//preservePatches (cyclic_half0 cyclic_half1);
+
+//- Keep all of faceSet on a single processor. This puts all cells
+//  connected with a point, edge or face on the same processor.
+//  (just having face connected cells might not guarantee a balanced
+//  decomposition)
+// The processor can be -1 (the decompositionMethod chooses the processor
+// for a good load balance) or explicitly provided (upsets balance).
+//singleProcessorFaceSets ((f0 -1));
+
+
+//- Use the volScalarField named here as a weight for each cell in the
+//  decomposition.  For example, use a particle population field to decompose
+//  for a balanced number of particles in a lagrangian simulation.
+// weightField dsmcRhoNMean;
+
+//method          scotch;
+//method          hierarchical;
+// method          simple;
+//method          metis;
+method          manual;
+// method          multiLevel;
+// method          structured;  // does 2D decomposition of structured mesh
+
+multiLevelCoeffs
+{
+    // Decomposition methods to apply in turn. This is like hierarchical but
+    // fully general - every method can be used at every level.
+
+    level0
+    {
+        numberOfSubdomains  64;
+        //method simple;
+        //simpleCoeffs
+        //{
+        //    n           (2 1 1);
+        //    delta       0.001;
+        //}
+        method scotch;
+    }
+    level1
+    {
+        numberOfSubdomains  4;
+        method scotch;
+    }
+}
+
+// Desired output
+
+simpleCoeffs
+{
+    n           (2 1 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (1 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+metisCoeffs
+{
+    method recursive;
+ /*
+    processorWeights
+    (
+        1
+        1
+        1
+        1
+    );
+  */
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "manualDecomposition";
+}
+
+structuredCoeffs
+{
+    // Patches to do 2D decomposition on. Structured mesh only; cells have
+    // to be in 'columns' on top of patches.
+    patches     (movingWall);
+
+    // Method to use on the 2D subset
+    method      scotch;
+}
+
+//// Is the case distributed? Note: command-line argument -roots takes
+//// precedence
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+// ************************************************************************* //

applications/test/hexRef8/block/system/fvSchemes

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

applications/test/hexRef8/block/system/fvSolution

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //

applications/test/sort/Test-sortList.C

@@ -50,11 +50,18 @@ int main(int argc, char *argv[])
     labelList a(orig);
     sortedOrder(a, order);
 
+    SortableList<label> aReverse(a.size());
+    aReverse = a;
+
     Info<< "unsorted: " << a << endl;
     sort(a);
     Info<< "sorted:   " << a << endl;
     Info<< "indices:  " << order << endl;
 
+    aReverse.reverseSort();
+    Info<< "reverse sorted:   " << aReverse << endl;
+    Info<< "reverse indices:  " << aReverse.indices() << endl;
+
     SortableList<label> b(orig);
     Info<< "unsorted: " << orig << endl;
     Info<< "sorted:   " << b << endl;

applications/test/volField/README.txt

@@ -1 +1,2 @@
-Test field reading and manipulation. Run on cavity subdirectory.
+Test field reading and manipulation.
+See cavity/Allrun in the subdirectory.

applications/test/volField/Test-volField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    volFieldTest
+    Test-volField
 
 \*---------------------------------------------------------------------------*/
 
@@ -52,8 +52,6 @@ int main(int argc, char *argv[])
         mesh
     );
 
-    //Info<< min(p, p);
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -85,9 +83,12 @@ int main(int argc, char *argv[])
         zeroGradientFvPatchSymmTensorField::typeName
     );
 
-    //Info<< fvc::div(st) << endl;
-
-    solve(fvm::ddt(st) + fvm::div(phi, st) - fvm::laplacian(st));
+    solve
+    (
+        fvm::ddt(st)
+      + fvm::div(phi, st)
+      - fvm::laplacian(dimensionedScalar("D", sqr(dimLength)/dimTime, 1), st)
+    );
 
     return 0;
 }

applications/test/volField/cavity/Allrun

@@ -0,0 +1,17 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Get application name
+application=Test-volField
+
+# Compile
+runApplication wmake ..
+
+runApplication blockMesh
+
+runApplication $application
+
+# ----------------------------------------------------------------- end-of-file

applications/test/volField/cavity/system/fvSchemes

@@ -23,26 +23,22 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
 }
 
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linear;
+    div(phi,st)     Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(nu,U) Gauss linear orthogonal;
-    laplacian((1|A(U)),p) Gauss linear orthogonal;
+    default         Gauss linear orthogonal;
 }
 
 interpolationSchemes
 {
     default         linear;
-    interpolate(HbyA) linear;
 }
 
 snGradSchemes

applications/test/volField/cavity/system/fvSolution

@@ -17,15 +17,7 @@ FoamFile
 
 solvers
 {
-    p
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-06;
-        relTol          0;
-    }
-
-    U
+    st
     {
         solver          PBiCG;
         preconditioner  DILU;

applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C

@@ -52,6 +52,7 @@ Usage
 #include "polyTopoChange.H"
 #include "fvMesh.H"
 #include "polyMeshFilter.H"
+#include "faceSet.H"
 
 using namespace Foam;
 
@@ -74,9 +75,9 @@ int main(int argc, char *argv[])
 
     argList::addOption
     (
-        "collapseFaceZone",
-        "zoneName",
-        "Collapse faces that are in the supplied face zone"
+        "collapseFaceSet",
+        "faceSet",
+        "Collapse faces that are in the supplied face set"
     );
 
 #   include "addOverwriteOption.H"
@@ -93,46 +94,115 @@ int main(int argc, char *argv[])
     const bool overwrite = args.optionFound("overwrite");
 
     const bool collapseFaces = args.optionFound("collapseFaces");
-    const bool collapseFaceZone = args.optionFound("collapseFaceZone");
+    const bool collapseFaceSet = args.optionFound("collapseFaceSet");
 
+    if (collapseFaces && collapseFaceSet)
+    {
+        FatalErrorIn("main(int, char*[])")
+            << "Both face zone collapsing and face collapsing have been"
+            << "selected. Choose only one of:" << nl
+            << "    -collapseFaces" << nl
+            << "    -collapseFaceSet <faceSet>"
+            << abort(FatalError);
+    }
+
+
+    // maintain indirectPatchFaces if it is there (default) or force
+    // (if collapseFaceSet option provided)
+    word faceSetName("indirectPatchFaces");
+    IOobject::readOption readFlag = IOobject::READ_IF_PRESENT;
+
+    if (args.optionReadIfPresent("collapseFaceSet", faceSetName))
+    {
+        readFlag = IOobject::MUST_READ;
+    }
+
+
+
+    labelIOList pointPriority
+    (
+        IOobject
+        (
+            "pointPriority",
+            runTime.timeName(),
+            runTime,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        labelList(mesh.nPoints(), labelMin)
+    );
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
 
         Info<< "Time = " << runTime.timeName() << endl;
 
-        polyMeshFilter meshFilter(mesh);
+        autoPtr<polyMeshFilter> meshFilterPtr;
+
+        label nBadFaces = 0;
+
+        faceSet indirectPatchFaces
+        (
+            mesh,
+            faceSetName,
+            readFlag,
+            IOobject::AUTO_WRITE
+        );
+        Info<< "Read faceSet " << indirectPatchFaces.name()
+            << " with "
+            << returnReduce(indirectPatchFaces.size(), sumOp<label>())
+            << " faces" << endl;
+
+
+        {
+            meshFilterPtr.set(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
 
             // newMesh will be empty until it is filtered
             const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
 
             // Filter small edges only. This reduces the number of faces so that
             // the face filtering is sped up.
-        label nBadFaces = meshFilter.filterEdges(0);
+            nBadFaces = meshFilter.filterEdges(0);
             {
-            polyTopoChange meshMod(newMesh);
+                polyTopoChange meshMod(newMesh());
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
 
-        if (collapseFaceZone)
-        {
-            const word faceZoneName = args.optionRead<word>("collapseFaceZone");
+            pointPriority = meshFilter.pointPriority();
+        }
 
-            const faceZone& fZone = mesh.faceZones()[faceZoneName];
+        if (collapseFaceSet)
+        {
+            meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
+
+            const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
 
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
-            meshFilter.filterFaceZone(fZone);
+            meshFilter.filter(indirectPatchFaces);
             {
                 polyTopoChange meshMod(newMesh);
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
+
+            pointPriority = meshFilter.pointPriority();
         }
 
         if (collapseFaces)
         {
+            meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
+
+            const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
+
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
             meshFilter.filter(nBadFaces);
@@ -140,7 +210,11 @@ int main(int argc, char *argv[])
                 polyTopoChange meshMod(newMesh);
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
+
+            pointPriority = meshFilter.pointPriority();
         }
 
         // Write resulting mesh
@@ -157,6 +231,7 @@ int main(int argc, char *argv[])
             << runTime.timeName() << nl << endl;
 
         mesh.write();
+        pointPriority.write();
     }
 
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/utilities/mesh/generation/Allwmake

@@ -0,0 +1,18 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake all blockMesh
+wmake all extrude
+wmake all extrude2DMesh
+wmake all snappyHexMesh
+
+if [ -d "$CGAL_ARCH_PATH" ]
+then
+    wmake libso foamyHexMesh/conformalVoronoiMesh
+    wmake all foamyHexMesh
+    wmake all foamyQuadMesh
+fi
+
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/foamyHexMesh/Make/options

@@ -10,8 +10,10 @@ include $(GENERAL_RULES)/CGAL
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_EXACT} \
     ${CGAL_INEXACT} \
     ${CGAL_INC} \
+    ${c++CGALWARN} \
     -IconformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C

@@ -190,20 +190,25 @@ void Foam::DelaunayMesh<Triangulation>::reset()
     resetVertexCount();
     resetCellCount();
 
-    insertPoints(vertices);
+    insertPoints(vertices, false);
 
     Info<< "Inserted " << vertexCount() << " fixed points" << endl;
 }
 
 
 template<class Triangulation>
-void Foam::DelaunayMesh<Triangulation>::insertPoints(const List<Vb>& vertices)
+Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::insertPoints
+(
+    const List<Vb>& vertices,
+    const bool reIndex
+)
 {
-    rangeInsertWithInfo
+    return rangeInsertWithInfo
     (
         vertices.begin(),
         vertices.end(),
-        false
+        false,
+        reIndex
     );
 }
 
@@ -268,7 +273,7 @@ const
 
 template<class Triangulation>
 template<class PointIterator>
-void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
+Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 (
     PointIterator begin,
     PointIterator end,
@@ -307,6 +312,10 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 
     Vertex_handle hint;
 
+    Map<label> oldToNewIndex(points.size());
+
+    label maxIndex = -1;
+
     for
     (
         typename vectorPairPointIndex::const_iterator p = points.begin();
@@ -335,11 +344,14 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
         {
             if (reIndex)
             {
+                const label oldIndex = vert.index();
                 hint->index() = getNewVertexIndex();
+                oldToNewIndex.insert(oldIndex, hint->index());
             }
             else
             {
                 hint->index() = vert.index();
+                maxIndex = max(maxIndex, vert.index());
             }
             hint->type() = vert.type();
             hint->procIndex() = vert.procIndex();
@@ -347,6 +359,13 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
             hint->alignment() = vert.alignment();
         }
     }
+
+    if (!reIndex)
+    {
+        vertexCount_ = maxIndex + 1;
+    }
+
+    return oldToNewIndex;
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H

@@ -228,7 +228,11 @@ public:
             void reset();
 
             //- Insert the list of vertices (calls rangeInsertWithInfo)
-            void insertPoints(const List<Vb>& vertices);
+            Map<label> insertPoints
+            (
+                const List<Vb>& vertices,
+                const bool reIndex
+            );
 
             //- Function inserting points into a triangulation and setting the
             //  index and type data of the point in the correct order. This is
@@ -237,7 +241,7 @@ public:
             //  Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
             //  Sebastien Loriot (Geometry Factory).
             template<class PointIterator>
-            void rangeInsertWithInfo
+            Map<label> rangeInsertWithInfo
             (
                 PointIterator begin,
                 PointIterator end,

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C

@@ -536,7 +536,8 @@ Foam::label Foam::DistributedDelaunayMesh<Triangulation>::referVertices
     labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
     (
         referredVertices.begin(),
-        referredVertices.end()
+        referredVertices.end(),
+        true
     );
 
     if (!pointsNotInserted.empty())
@@ -966,6 +967,15 @@ Foam::DistributedDelaunayMesh<Triangulation>::rangeInsertReferredWithInfo
             )   << "Point is outside affine hull! pt = " << pointToInsert
                 << endl;
         }
+        else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
+        {
+            // @todo Can this be optimised?
+            //
+            // Only want to insert if a connection is formed between
+            // pointToInsert and an internal or internal boundary point.
+            hint = Triangulation::insert(pointToInsert, c);
+            inserted = true;
+        }
         else
         {
             // Get the cells that conflict with p in a vector V,

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files

@@ -76,6 +76,7 @@ faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
 
 searchableSurfaceFeatures/searchableSurfaceFeatures.C
 searchableSurfaceFeatures/searchableBoxFeatures.C
+searchableSurfaceFeatures/searchablePlateFeatures.C
 searchableSurfaceFeatures/triSurfaceMeshFeatures.C
 
 LIB = $(FOAM_LIBBIN)/libconformalVoronoiMesh

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options

@@ -11,8 +11,10 @@ FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_EXACT} \
     ${CGAL_INEXACT} \
     ${CGAL_INC} \
+    ${c++CGALWARN} \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C

@@ -628,49 +628,6 @@ Foam::tensorField Foam::cellShapeControlMesh::dumpAlignments() const
 }
 
 
-void Foam::cellShapeControlMesh::insertBoundingPoints
-(
-    const boundBox& bb,
-    const cellSizeAndAlignmentControls& sizeControls
-)
-{
-    // Loop over bound box points and get cell size and alignment
-    const pointField bbPoints(bb.points());
-
-    forAll(bbPoints, pI)
-    {
-        const Foam::point& pt = bbPoints[pI];
-
-        // Cell size here will return default cell size
-        const scalar cellSize = sizeControls.cellSize(pt);
-
-        if (debug)
-        {
-            Info<< "Insert Bounding Point: " << pt << " " << cellSize << endl;
-        }
-
-        // Get the cell size of the nearest surface.
-//        geometryToConformTo_.findSurfaceNearest
-//        (
-//            pt,
-//            GREAT,
-//            surfHit,
-//            hitSurface
-//        );
-
-        const tensor alignment = tensor::I;
-
-        insert
-        (
-            pt,
-            cellSize,
-            alignment,
-            Vb::vtInternalNearBoundary
-        );
-    }
-}
-
-
 void Foam::cellShapeControlMesh::write() const
 {
     Info<< "Writing " << meshSubDir << endl;
@@ -781,4 +738,78 @@ void Foam::cellShapeControlMesh::write() const
 }
 
 
+Foam::label Foam::cellShapeControlMesh::estimateCellCount
+(
+    const autoPtr<backgroundMeshDecomposition>& decomposition
+) const
+{
+    // Loop over all the tets and estimate the cell count in each one
+
+    scalar cellCount = 0;
+
+    for
+    (
+        Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if (!cit->hasFarPoint() && !is_infinite(cit))
+        {
+            // @todo Check if tet centre is on the processor..
+            CGAL::Tetrahedron_3<baseK> tet
+            (
+                cit->vertex(0)->point(),
+                cit->vertex(1)->point(),
+                cit->vertex(2)->point(),
+                cit->vertex(3)->point()
+            );
+
+            pointFromPoint centre = topoint(CGAL::centroid(tet));
+
+            if
+            (
+                Pstream::parRun()
+             && !decomposition().positionOnThisProcessor(centre)
+            )
+            {
+                continue;
+            }
+
+            scalar volume = CGAL::to_double(tet.volume());
+
+            scalar averagedPointCellSize = 0;
+            //scalar averagedPointCellSize = 1;
+
+            // Get an average volume by averaging the cell size of the vertices
+            for (label vI = 0; vI < 4; ++vI)
+            {
+                averagedPointCellSize += cit->vertex(vI)->targetCellSize();
+                //averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
+            }
+
+            averagedPointCellSize /= 4;
+            //averagedPointCellSize = ::sqrt(averagedPointCellSize);
+
+//            if (averagedPointCellSize < SMALL)
+//            {
+//                Pout<< "Volume = " << volume << endl;
+//
+//                for (label vI = 0; vI < 4; ++vI)
+//                {
+//                    Pout<< "Point " << vI
+//                        << ", point = " << topoint(cit->vertex(vI)->point())
+//                        << ", size = " << cit->vertex(vI)->targetCellSize()
+//                        << endl;
+//                }
+//            }
+
+            cellCount += volume/pow(averagedPointCellSize, 3);
+        }
+    }
+
+    return cellCount;
+}
+
+
 // ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H

@@ -153,15 +153,14 @@ public:
 
             tensorField dumpAlignments() const;
 
-            void insertBoundingPoints
-            (
-                const boundBox& bb,
-                const cellSizeAndAlignmentControls& sizeControls
-            );
-
             void writeTriangulation();
 
             void write() const;
+
+            label estimateCellCount
+            (
+                const autoPtr<backgroundMeshDecomposition>& decomposition
+            ) const;
 };
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -429,6 +429,11 @@ void Foam::searchableSurfaceControl::initialVertices
                 vectorField normals(1);
                 searchableSurface_.getNormal(infoList, normals);
 
+                if (mag(normals[0]) < SMALL)
+                {
+                    normals[0] = vector(1, 1, 1);
+                }
+
                 pointAlignment.set(new triad(normals[0]));
 
                 if (infoList[0].hit())

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C

@@ -768,7 +768,7 @@ Foam::label Foam::controlMeshRefinement::refineMesh
         }
     }
 
-    mesh_.insertPoints(verts);
+    mesh_.insertPoints(verts, false);
 
     return verts.size();
 }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -42,8 +42,25 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(conformalVoronoiMesh, 0);
+
+    template<>
+    const char* NamedEnum
+    <
+        conformalVoronoiMesh::dualMeshPointType,
+        5
+    >::names[] =
+    {
+        "internal",
+        "surface",
+        "featureEdge",
+        "featurePoint",
+        "constrained"
+    };
 }
 
+const Foam::NamedEnum<Foam::conformalVoronoiMesh::dualMeshPointType, 5>
+    Foam::conformalVoronoiMesh::dualMeshPointTypeNames_;
+
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -162,16 +179,24 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
 }
 
 
-void Foam::conformalVoronoiMesh::insertPoints
+Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
 (
     List<Vb>& vertices,
-    bool distribute
+    bool distribute,
+    bool reIndex
 )
 {
     if (Pstream::parRun() && distribute)
     {
+        autoPtr<mapDistribute> mapDist =
             decomposition_().distributePoints(vertices);
 
+        // Re-index the point pairs if one or both have been distributed.
+        // If both, remove
+
+        // If added a point, then need to know its point pair
+        // If one moved, then update procIndex locally
+
         forAll(vertices, vI)
         {
             vertices[vI].procIndex() = Pstream::myProcNo();
@@ -180,7 +205,8 @@ void Foam::conformalVoronoiMesh::insertPoints
 
     label preReinsertionSize(number_of_vertices());
 
-    this->DelaunayMesh<Delaunay>::insertPoints(vertices);
+    Map<label> oldToNewIndices =
+        this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
 
     const label nReinserted = returnReduce
     (
@@ -191,17 +217,18 @@ void Foam::conformalVoronoiMesh::insertPoints
     Info<< "    Reinserted " << nReinserted << " vertices out of "
         << returnReduce(vertices.size(), sumOp<label>())
         << endl;
+
+    return oldToNewIndices;
 }
 
 
 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 (
     const pointIndexHitAndFeatureList& surfaceHits,
-    const fileName fName
+    const fileName fName,
+    DynamicList<Vb>& pts
 )
 {
-    DynamicList<Vb> pts(2.0*surfaceHits.size());
-
     forAll(surfaceHits, i)
     {
         vectorField norm(1);
@@ -229,6 +256,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 normal,
+                true,
                 pts
             );
         }
@@ -239,6 +267,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 normal,
+                true,
                 pts
             );
         }
@@ -249,6 +278,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 -normal,
+                true,
                 pts
             );
         }
@@ -263,8 +293,6 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
         }
     }
 
-    insertPoints(pts, true);
-
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
         DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@@ -275,11 +303,10 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 void Foam::conformalVoronoiMesh::insertEdgePointGroups
 (
     const pointIndexHitAndFeatureList& edgeHits,
-    const fileName fName
+    const fileName fName,
+    DynamicList<Vb>& pts
 )
 {
-    DynamicList<Vb> pts(3.0*edgeHits.size());
-
     forAll(edgeHits, i)
     {
         if (edgeHits[i].first().hit())
@@ -301,10 +328,6 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
         }
     }
 
-    pts.shrink();
-
-    insertPoints(pts, true);
-
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
         DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@@ -331,6 +354,28 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
 }
 
 
+bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
+(
+    const Foam::point& pt
+) const
+{
+    scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
+
+    pointIndexHit info;
+    label featureHit;
+
+    geometryToConformTo_.findEdgeNearest
+    (
+        pt,
+        exclusionRangeSqr,
+        info,
+        featureHit
+    );
+
+    return info.hit();
+}
+
+
 void Foam::conformalVoronoiMesh::insertInitialPoints()
 {
     Info<< nl << "Inserting initial points" << endl;
@@ -842,6 +887,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         geometryToConformTo_
     ),
     limitBounds_(),
+    ptPairs_(*this),
     ftPtConformer_(*this),
     edgeLocationTreePtr_(),
     surfacePtLocationTreePtr_(),
@@ -1091,7 +1137,7 @@ void Foam::conformalVoronoiMesh::move()
                 vB->alignment().T() & cartesianDirections
             );
 
-            Field<vector> alignmentDirs(3);
+            Field<vector> alignmentDirs(alignmentDirsA);
 
             forAll(alignmentDirsA, aA)
             {
@@ -1249,7 +1295,7 @@ void Foam::conformalVoronoiMesh::move()
                     if
                     (
                         (
-                            (vA->internalPoint() || vB->internalPoint())
+                            (vA->internalPoint() && vB->internalPoint())
                          && (!vA->referred() || !vB->referred())
 //                         ||
 //                            (
@@ -1315,14 +1361,19 @@ void Foam::conformalVoronoiMesh::move()
                     {
                         // Point removal
                         if
+                        (
                             (
                                 vA->internalPoint()
                              && !vA->referred()
                              && !vA->fixed()
-                         && vB->internalPoint()
+                            )
+                         &&
+                            (
+                                vB->internalPoint()
                              && !vB->referred()
                              && !vB->fixed()
                             )
+                        )
                         {
                             // Only insert a point at the midpoint of
                             // the short edge if neither attached
@@ -1490,9 +1541,7 @@ void Foam::conformalVoronoiMesh::move()
         Info<< "Writing point displacement vectors to file." << endl;
         OFstream str
         (
-            time().path()
-          + "displacements_" + runTime_.timeName()
-          + ".obj"
+            time().path()/"displacements_" + runTime_.timeName() + ".obj"
         );
 
         for
@@ -1590,6 +1639,74 @@ void Foam::conformalVoronoiMesh::move()
                 time().path()/"boundaryPoints_" + time().timeName() + ".obj",
                 *this
             );
+
+            DelaunayMeshTools::writeOBJ
+            (
+                time().path()/"internalBoundaryPoints_" + time().timeName()
+              + ".obj",
+                *this,
+                Foam::indexedVertexEnum::vtInternalSurface,
+                Foam::indexedVertexEnum::vtInternalFeaturePoint
+            );
+
+            DelaunayMeshTools::writeOBJ
+            (
+                time().path()/"externalBoundaryPoints_" + time().timeName()
+              + ".obj",
+                *this,
+                Foam::indexedVertexEnum::vtExternalSurface,
+                Foam::indexedVertexEnum::vtExternalFeaturePoint
+            );
+
+            OBJstream multipleIntersections
+            (
+                "multipleIntersections_"
+              + time().timeName()
+              + ".obj"
+            );
+
+            for
+            (
+                Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+                eit != finite_edges_end();
+                ++eit
+            )
+            {
+                Cell_handle c = eit->first;
+                Vertex_handle vA = c->vertex(eit->second);
+                Vertex_handle vB = c->vertex(eit->third);
+
+                Foam::point ptA(topoint(vA->point()));
+                Foam::point ptB(topoint(vB->point()));
+
+                List<pointIndexHit> surfHits;
+                labelList hitSurfaces;
+
+                geometryToConformTo().findSurfaceAllIntersections
+                (
+                    ptA,
+                    ptB,
+                    surfHits,
+                    hitSurfaces
+                );
+
+                if
+                (
+                    surfHits.size() >= 2
+                 || (
+                     surfHits.size() == 0
+                  && (
+                          (vA->externalBoundaryPoint()
+                       && vB->internalBoundaryPoint())
+                       || (vB->externalBoundaryPoint()
+                       && vA->internalBoundaryPoint())
+                     )
+                    )
+                )
+                {
+                    multipleIntersections.write(linePointRef(ptA, ptB));
+                }
+            }
         }
     }
 
@@ -1600,7 +1717,6 @@ void Foam::conformalVoronoiMesh::move()
         printVertexInfo(Info);
     }
 
-    // Write the intermediate mesh, do not filter the dual faces.
     if (time().outputTime())
     {
         writeMesh(time().timeName());
@@ -1613,49 +1729,6 @@ void Foam::conformalVoronoiMesh::move()
 }
 
 
-//Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
-//(
-//    const List<Foam::point>& centres,
-//    const List<scalar>& radiusSqrs
-//) const
-//{
-//    if (!Pstream::parRun())
-//    {
-//        return labelListList(centres.size(), labelList(0));
-//    }
-//
-////    DynamicList<Foam::point> pts(number_of_vertices());
-//
-////    for
-////    (
-////        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-////        vit != finite_vertices_end();
-////        vit++
-////    )
-////    {
-////        pts.append(topoint(vit->point()));
-////    }
-////
-////    dynamicIndexedOctree<dynamicTreeDataPoint> vertexOctree
-////    (
-////        dynamicTreeDataPoint(pts),
-////        treeBoundBox(min(pts), max(pts)),
-////        10, // maxLevel
-////        10, // leafSize
-////        3.0 // duplicity
-////    );
-//
-//    return decomposition_().overlapsProcessors
-//    (
-//        centres,
-//        radiusSqrs,
-//        *this,
-//        false//,
-////        vertexOctree
-//    );
-//}
-
-
 void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
 {
     typedef CGAL::Exact_predicates_exact_constructions_kernel   Kexact;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -77,7 +77,7 @@ SourceFiles
 #include "Pair.H"
 #include "DistributedDelaunayMesh.H"
 #include "featurePointConformer.H"
-
+#include "pointPairs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,12 +86,10 @@ namespace Foam
 
 // Forward declaration of classes
 class initialPointsMethod;
-
 class relaxationModel;
-
 class faceAreaWeightModel;
-
 class backgroundMeshDecomposition;
+class OBJstream;
 
 /*---------------------------------------------------------------------------*\
                     Class conformalVoronoiMesh Declaration
@@ -115,6 +113,19 @@ public:
     typedef List<pointIndexHitAndFeature>        pointIndexHitAndFeatureList;
     typedef DynamicList<pointIndexHitAndFeature> pointIndexHitAndFeatureDynList;
 
+    // Static data
+
+        enum dualMeshPointType
+        {
+            internal     = 0,
+            surface      = 1,
+            featureEdge  = 2,
+            featurePoint = 3,
+            constrained  = 4
+        };
+
+        static const NamedEnum<dualMeshPointType, 5> dualMeshPointTypeNames_;
+
 
 private:
 
@@ -152,6 +163,9 @@ private:
         //- Limiting bound box before infinity begins
         treeBoundBox limitBounds_;
 
+        //-
+        mutable pointPairs<Delaunay> ptPairs_;
+
         //-
         featurePointConformer ftPtConformer_;
 
@@ -216,10 +230,11 @@ private:
             const bool distribute = false
         );
 
-        void insertPoints
+        Map<label> insertPointPairs
         (
             List<Vb>& vertices,
-            bool distribute
+            bool distribute,
+            bool reIndex
         );
 
         //- Create a point-pair at a ppDist distance either side of
@@ -229,6 +244,7 @@ private:
             const scalar ppDist,
             const Foam::point& surfPt,
             const vector& n,
+            const bool ptPair,
             DynamicList<Vb>& pts
         ) const;
 
@@ -241,24 +257,20 @@ private:
             const scalar ppDist,
             const Foam::point& surfPt,
             const vector& n,
+            const bool ptPair,
             DynamicList<Vb>& pts
         ) const;
 
         //- Check internal point is completely inside the meshable region
         inline bool internalPointIsInside(const Foam::point& pt) const;
 
-        inline bool isPointPair
-        (
-            const Vertex_handle& vA,
-            const Vertex_handle& vB
-        ) const;
-
         //- Insert pairs of points on the surface with the given normals, at the
         //  specified spacing
         void insertSurfacePointPairs
         (
             const pointIndexHitAndFeatureList& surfaceHits,
-            const fileName fName = fileName::null
+            const fileName fName,
+            DynamicList<Vb>& pts
         );
 
         //- Insert groups of points to conform to an edge given a list of
@@ -267,7 +279,8 @@ private:
         void insertEdgePointGroups
         (
             const pointIndexHitAndFeatureList& edgeHits,
-            const fileName fName = fileName::null
+            const fileName fName,
+            DynamicList<Vb>& pts
         );
 
         void createEdgePointGroupByCirculating
@@ -338,6 +351,9 @@ private:
         //- Check if a location is in exclusion range around a feature point
         bool nearFeaturePt(const Foam::point& pt) const;
 
+        //- Check if a surface point is in exclusion range around a feature edge
+        bool surfacePtNearFeatureEdge(const Foam::point& pt) const;
+
         //- Insert the initial points into the triangulation, based on the
         //  initialPointsMethod
         void insertInitialPoints();
@@ -461,14 +477,11 @@ private:
             label& hitSurface
         ) const;
 
-        //- Find the "worst" incursion of the dual cell of a non-internal or
-        //  boundary point through the surface, subject to the
-        //  maxSurfaceProtrusion tolerance
         void dualCellLargestSurfaceIncursion
         (
             const Delaunay::Finite_vertices_iterator& vit,
-            pointIndexHit& surfHitLargest,
-            label& hitSurfaceLargest
+            pointIndexHit& surfHit,
+            label& hitSurface
         ) const;
 
         //- Write out vertex-processor occupancy information for debugging
@@ -683,13 +696,12 @@ private:
         void mergeIdenticalDualVertices
         (
             const pointField& pts,
-            const labelList& boundaryPts
+            labelList& boundaryPts
         );
 
         label mergeIdenticalDualVertices
         (
             const pointField& pts,
-            const labelList& boundaryPts,
             Map<label>& dualPtIndexMap
         ) const;
 
@@ -716,6 +728,8 @@ private:
         //  elements damage the mesh quality, allowing backtracking.
         labelHashSet checkPolyMeshQuality(const pointField& pts) const;
 
+        label classifyBoundaryPoint(Cell_handle cit) const;
+
         //- Index all of the the Delaunay cells and calculate their
         //- dual points
         void indexDualVertices
@@ -725,7 +739,11 @@ private:
         );
 
         //- Re-index all of the Delaunay cells
-        void reindexDualVertices(const Map<label>& dualPtIndexMap);
+        void reindexDualVertices
+        (
+            const Map<label>& dualPtIndexMap,
+            labelList& boundaryPts
+        );
 
         label createPatchInfo
         (
@@ -835,6 +853,8 @@ private:
             const PtrList<dictionary>& patchDicts
         ) const;
 
+        void writePointPairs(const fileName& fName) const;
+
         //- Disallow default bitwise copy construct
         conformalVoronoiMesh(const conformalVoronoiMesh&);
 
@@ -981,7 +1001,7 @@ public:
                 const wordList& patchNames,
                 const PtrList<dictionary>& patchDicts,
                 const pointField& cellCentres,
-                const PackedBoolList& boundaryFacesToRemove
+                PackedBoolList& boundaryFacesToRemove
             ) const;
 
             //- Calculate and write a field of the target cell size,

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -30,480 +30,11 @@ License
 #include "indexedCellChecks.H"
 #include "OBJstream.H"
 #include "indexedCellOps.H"
+#include "ListOps.H"
 #include "DelaunayMeshTools.H"
 
-#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
-#include "CGAL/Gmpq.h"
-
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-void Foam::conformalVoronoiMesh::checkCells()
-{
-    List<List<FixedList<Foam::point, 4> > > cellListList(Pstream::nProcs());
-
-    List<FixedList<Foam::point, 4> > cells(number_of_finite_cells());
-
-    globalIndex gIndex(number_of_vertices());
-
-    label count = 0;
-    for
-    (
-        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (tetrahedron(cit).volume() == 0)
-        {
-            Pout<< "ZERO VOLUME TET" << endl;
-            Pout<< cit->info();
-            Pout<< cit->dual();
-        }
-
-        if (cit->hasFarPoint())
-        {
-            continue;
-        }
-
-        List<labelPair> cellVerticesPair(4);
-        List<Foam::point> cellVertices(4);
-
-        for (label vI = 0; vI < 4; ++vI)
-        {
-            cellVerticesPair[vI] = labelPair
-            (
-                cit->vertex(vI)->procIndex(),
-                cit->vertex(vI)->index()
-            );
-            cellVertices[vI] = topoint(cit->vertex(vI)->point());
-        }
-
-        List<Foam::point> cellVerticesOld(cellVertices);
-        labelList oldToNew;
-        sortedOrder(cellVerticesPair, oldToNew);
-        oldToNew = invert(oldToNew.size(), oldToNew);
-        inplaceReorder(oldToNew, cellVerticesPair);
-        inplaceReorder(oldToNew, cellVertices);
-
-//        FixedList<label, 4> globalTetCell
-//        (
-//            cit->globallyOrderedCellVertices(gIndex)
-//        );
-//
-//        FixedList<Point, 4> cellVertices(Point(0,0,0));
-//
-//        forAll(globalTetCell, gvI)
-//        {
-//            label gI = globalTetCell[gvI];
-//
-//            cellVertices[gvI] = cit->vertex(gI)->point();
-//        }
-
-//        if (cit->hasFarPoint())
-//        {
-//            continue;
-//        }
-
-        for (label i = 0; i < 4; ++i)
-        {
-            //cells[count][i] = topoint(cit->vertex(i)->point());
-            cells[count][i] = cellVertices[i];
-        }
-
-        count++;
-    }
-
-    cells.setSize(count);
-
-    cellListList[Pstream::myProcNo()] = cells;
-
-    Pstream::gatherList(cellListList);
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the cells of each " << nl
-            << "processor point list against the master cell list." << nl
-            << "There are " << cellListList.size() << " processors" << nl
-            << "The size of each processor's cell list is:" << endl;
-
-        forAll(cellListList, cfI)
-        {
-            Info<< "    Proc " << cfI << " has " << cellListList[cfI].size()
-                << " cells" << endl;
-        }
-
-        label nMatches = 0, nMatchFoundDiffOrder = 0;
-
-        forAll(cellListList[0], cmI)
-        {
-            const FixedList<Foam::point, 4>& masterCell = cellListList[0][cmI];
-
-            bool matchFound = false;
-            bool matchFoundDiffOrder = false;
-
-            forAll(cellListList, cpI)
-            {
-                if (cpI == 0)
-                {
-                    continue;
-                }
-
-                forAll(cellListList[cpI], csI)
-                {
-                    const FixedList<Foam::point, 4>& slaveCell
-                        = cellListList[cpI][csI];
-
-                    if (masterCell == slaveCell)
-                    {
-                        matchFound = true;
-                        break;
-                    }
-                    else
-                    {
-                        label samePt = 0;
-
-                        forAll(masterCell, mI)
-                        {
-                            const Foam::point& mPt = masterCell[mI];
-
-                            forAll(slaveCell, sI)
-                            {
-                                const Foam::point& sPt = slaveCell[sI];
-
-                                if (mPt == sPt)
-                                {
-                                    samePt++;
-                                }
-                            }
-                        }
-
-                        if (samePt == 4)
-                        {
-                            matchFoundDiffOrder = true;
-
-                            Pout<< masterCell << nl << slaveCell << endl;
-
-                            break;
-                        }
-                    }
-                }
-            }
-
-            if (matchFound)
-            {
-                nMatches++;
-            }
-
-            if (matchFoundDiffOrder)
-            {
-                nMatchFoundDiffOrder++;
-            }
-        }
-
-        Info<< "Found " << nMatches << " matching cells and "
-            << nMatchFoundDiffOrder << " matching cells with different "
-            << "vertex ordering"<< endl;
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::checkDuals()
-{
-    List<List<Point> > pointFieldList(Pstream::nProcs());
-
-    List<Point> duals(number_of_finite_cells());
-
-//    PackedBoolList bPoints(number_of_finite_cells());
-
-//    indexDualVertices(duals, bPoints);
-
-    label count = 0;//duals.size();
-
-    duals.setSize(number_of_finite_cells());
-
-    globalIndex gIndex(number_of_vertices());
-
-    for
-    (
-        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (cit->hasFarPoint())
-        {
-            continue;
-        }
-
-        duals[count++] = cit->circumcenter();
-
-//        List<labelPair> cellVerticesPair(4);
-//        List<Point> cellVertices(4);
-//
-//        for (label vI = 0; vI < 4; ++vI)
-//        {
-//            cellVerticesPair[vI] = labelPair
-//            (
-//                cit->vertex(vI)->procIndex(),
-//                cit->vertex(vI)->index()
-//            );
-//            cellVertices[vI] = cit->vertex(vI)->point();
-//        }
-//
-//        labelList oldToNew;
-//        sortedOrder(cellVerticesPair, oldToNew);
-//        oldToNew = invert(oldToNew.size(), oldToNew);
-//        inplaceReorder(oldToNew, cellVerticesPair);
-//        inplaceReorder(oldToNew, cellVertices);
-//
-//        duals[count++] = CGAL::circumcenter
-//        (
-//            cellVertices[0],
-//            cellVertices[1],
-//            cellVertices[2],
-//            cellVertices[3]
-//        );
-
-//        To_exact to_exact;
-//        Back_from_exact back_from_exact;
-//        EK::Construct_circumcenter_3 exact_circumcenter =
-//            EK().construct_circumcenter_3_object();
-//
-//        duals[count++] = topoint
-//        (
-//            back_from_exact
-//            (
-//                exact_circumcenter
-//                (
-//                    to_exact(cit->vertex(0)->point()),
-//                    to_exact(cit->vertex(1)->point()),
-//                    to_exact(cit->vertex(2)->point()),
-//                    to_exact(cit->vertex(3)->point())
-//                )
-//            )
-//        );
-    }
-
-    Pout<< "Duals Calculated " << count << endl;
-
-    duals.setSize(count);
-
-    pointFieldList[Pstream::myProcNo()] = duals;
-
-    Pstream::gatherList(pointFieldList);
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the dual locations of each " << nl
-            << "processor point list against the master dual list." << nl
-            << "There are " << pointFieldList.size() << " processors" << nl
-            << "The size of each processor's dual list is:" << endl;
-
-        forAll(pointFieldList, pfI)
-        {
-            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
-                << " duals" << endl;
-        }
-
-        label nNonMatches = 0;
-        label nNearMatches = 0;
-        label nExactMatches = 0;
-
-        forAll(pointFieldList[0], pI)
-        {
-            const Point& masterPoint = pointFieldList[0][pI];
-
-            bool foundMatch = false;
-            bool foundNearMatch = false;
-
-            scalar minCloseness = GREAT;
-            Point closestPoint(0, 0, 0);
-
-            forAll(pointFieldList, pfI)
-            {
-                if (pfI == 0)
-                {
-                    continue;
-                }
-
-//                label pfI = 1;
-
-                forAll(pointFieldList[pfI], pISlave)
-                {
-                    const Point& slavePoint
-                        = pointFieldList[pfI][pISlave];
-
-                    if (masterPoint == slavePoint)
-                    {
-                        foundMatch = true;
-                        break;
-                    }
-
-                    const scalar closeness = mag
-                    (
-                        topoint(masterPoint) - topoint(slavePoint)
-                    );
-
-                    if (closeness < 1e-12)
-                    {
-                        foundNearMatch = true;
-                    }
-                    else
-                    {
-                        if (closeness < minCloseness)
-                        {
-                            minCloseness = closeness;
-                            closestPoint = slavePoint;
-                        }
-                    }
-                }
-
-                if (!foundMatch)
-                {
-                    if (foundNearMatch)
-                    {
-                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-
-                        std::cout<< "master = " << x << " " << y << " " << z
-                            << std::endl;
-
-                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-                        std::cout<< "slave  = " << xs << " " << ys << " " << zs
-                            << std::endl;
-
-                        nNearMatches++;
-                    }
-                    else
-                    {
-                        nNonMatches++;
-                        Info<< "    Closest point to " << masterPoint << " is "
-                            << closestPoint << nl
-                            << "    Separation is " << minCloseness << endl;
-
-                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-
-                        std::cout<< "master = " << x << " " << y << " " << z
-                                 << std::endl;
-
-                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-                        std::cout<< "slave  = " << xs << " " << ys << " " << zs
-                                 << std::endl;
-                    }
-                }
-                else
-                {
-                    nExactMatches++;
-                }
-            }
-        }
-
-        Info<< "Found " << nNonMatches << " non-matching duals" << nl
-            << " and " << nNearMatches << " near matches"
-            << " and " << nExactMatches << " exact matches" << endl;
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::checkVertices()
-{
-    List<pointField> pointFieldList(Pstream::nProcs());
-
-    pointField points(number_of_vertices());
-
-    labelPairHashSet duplicateVertices;
-
-    label count = 0;
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (duplicateVertices.found(labelPair(vit->procIndex(), vit->index())))
-        {
-            Pout<< "DUPLICATE " << vit->procIndex() << vit->index() << endl;
-        }
-        else
-        {
-            duplicateVertices.insert(labelPair(vit->procIndex(), vit->index()));
-        }
-
-        points[count++] = topoint(vit->point());
-    }
-
-    pointFieldList[Pstream::myProcNo()] = points;
-
-    Pstream::gatherList(pointFieldList);
-
-    OFstream str("missingPoints.obj");
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the point locations of each " << nl
-            << "processor point list against the master point list." << nl
-            << "There are " << pointFieldList.size() << " processors" << nl
-            << "The size of each processor's point list is:" << endl;
-
-        forAll(pointFieldList, pfI)
-        {
-            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
-                << " points" << endl;
-        }
-
-        label nNonMatches = 0;
-
-        forAll(pointFieldList[0], pI)
-        {
-            const Foam::point& masterPoint = pointFieldList[0][pI];
-
-            forAll(pointFieldList, pfI)
-            {
-                if (pI == 0)
-                {
-                    continue;
-                }
-
-                bool foundMatch = false;
-
-                forAll(pointFieldList[pfI], pISlave)
-                {
-                    const Foam::point& slavePoint
-                        = pointFieldList[pfI][pISlave];
-
-                    if (masterPoint == slavePoint)
-                    {
-                        foundMatch = true;
-                        break;
-                    }
-                }
-
-                if (!foundMatch)
-                {
-                    Info<< "    Proc " << pfI << " Master != Slave -> "
-                        << masterPoint << endl;
-
-                    meshTools::writeOBJ(str, masterPoint);
-
-                    nNonMatches++;
-                }
-            }
-        }
-
-        Info<< "Found a total of " << nNonMatches << " non-matching points"
-            << endl;
-    }
-}
-
-
 void Foam::conformalVoronoiMesh::calcDualMesh
 (
     pointField& points,
@@ -521,53 +52,6 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 {
     timeCheck("Start calcDualMesh");
 
-//    if (debug)
-//    {
-//        Pout<< nl << "Perfoming some checks . . ." << nl << nl
-//            << "Total number of vertices = " << number_of_vertices() << nl
-//            << "Total number of cells    = " << number_of_finite_cells()
-//            << endl;
-//
-//        checkVertices();
-//        checkCells();
-//        checkDuals();
-//
-//        Info<< nl << "Finished checks" << nl << endl;
-//    }
-
-//    OFstream str("attachedToFeature.obj");
-//    label offset = 0;
-//
-//    for
-//    (
-//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-//        vit != finite_vertices_end();
-//        ++vit
-//    )
-//    {
-//        if (vit->featurePoint())
-//        {
-//            std::list<Cell_handle> adjacentCells;
-//
-//            finite_incident_cells(vit, std::back_inserter(adjacentCells));
-//
-//            for
-//            (
-//                std::list<Cell_handle>::iterator acit = adjacentCells.begin();
-//                acit != adjacentCells.end();
-//                ++acit
-//            )
-//            {
-//                if ((*acit)->real())
-//                {
-//                    drawDelaunayCell(str, (*acit), offset);
-//                    offset++;
-////                    meshTools::writeOBJ(str, topoint((*acit)->dual()));
-//                }
-//            }
-//        }
-//    }
-
     setVertexSizeAndAlignment();
 
     timeCheck("After setVertexSizeAndAlignment");
@@ -814,7 +298,7 @@ void Foam::conformalVoronoiMesh::calcTetMesh
 void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 (
     const pointField& pts,
-    const labelList& boundaryPts
+    labelList& boundaryPts
 )
 {
     // Assess close points to be merged
@@ -829,11 +313,10 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
         nPtsMerged = mergeIdenticalDualVertices
         (
             pts,
-            boundaryPts,
             dualPtIndexMap
         );
 
-        reindexDualVertices(dualPtIndexMap);
+        reindexDualVertices(dualPtIndexMap, boundaryPts);
 
         reduce(nPtsMerged, sumOp<label>());
 
@@ -851,7 +334,6 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 (
     const pointField& pts,
-    const labelList& boundaryPts,
     Map<label>& dualPtIndexMap
 ) const
 {
@@ -883,6 +365,19 @@ Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 
             if (p1 == p2)
             {
+//                if (c1->parallelDualVertex() || c2->parallelDualVertex())
+//                {
+//                    if (c1->vertexLowestProc() < c2->vertexLowestProc())
+//                    {
+//                        dualPtIndexMap.insert(c1I, c1I);
+//                        dualPtIndexMap.insert(c2I, c1I);
+//                    }
+//                    else
+//                    {
+//                        dualPtIndexMap.insert(c1I, c2I);
+//                        dualPtIndexMap.insert(c2I, c2I);
+//                    }
+//                }
                 if (c1I < c2I)
                 {
                     dualPtIndexMap.insert(c1I, c1I);
@@ -1228,7 +723,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
         false
     );
 
-    //createCellCentres(cellCentres);
     cellCentres = DelaunayMeshTools::allPoints(*this);
 
     labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
@@ -1310,24 +804,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
 
     pMesh.addPatches(patches);
 
-    // Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
-
-    // forAll(patches, p)
-    // {
-    //     patches[p] = new polyPatch
-    //     (
-    //         patchNames[p],
-    //         patchSizes[p],
-    //         patchStarts[p],
-    //         p,
-    //         pMesh.boundaryMesh()
-    //     );
-    // }
-
-    // pMesh.addPatches(patches, false);
-
-    // pMesh.overrideCellCentres(cellCentres);
-
     return meshPtr;
 }
 
@@ -1338,7 +814,7 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
 
     timeCheck("Start of Cell Sizing");
 
-    labelList boundaryPts(number_of_finite_cells(), -1);
+    labelList boundaryPts(number_of_finite_cells(), internal);
     pointField ptsField;
 
     indexDualVertices(ptsField, boundaryPts);
@@ -1723,6 +1199,41 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
+(
+    Cell_handle cit
+) const
+{
+    if (cit->boundaryDualVertex())
+    {
+        if (cit->featurePointDualVertex())
+        {
+            return featurePoint;
+        }
+        else if (cit->featureEdgeDualVertex())
+        {
+            return featureEdge;
+        }
+        else
+        {
+            return surface;
+        }
+    }
+    else if (cit->baffleSurfaceDualVertex())
+    {
+        return surface;
+    }
+    else if (cit->baffleEdgeDualVertex())
+    {
+        return featureEdge;
+    }
+    else
+    {
+        return internal;
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::indexDualVertices
 (
     pointField& pts,
@@ -1755,7 +1266,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
     boundaryPts.setSize
     (
         number_of_finite_cells() + nConstrainedVertices,
-        -1
+        internal
     );
 
     if (foamyHexMeshControls().guardFeaturePoints())
@@ -1774,7 +1285,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                     topoint(vit->point());
 
                 boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
-                    1;
+                    constrained;
 
                 nConstrainedVertices++;
             }
@@ -1972,17 +1483,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 //                }
 //            }
 
-            if (cit->boundaryDualVertex())
-            {
-                if (cit->featureEdgeDualVertex())
-                {
-                    boundaryPts[cit->cellIndex()] = 1;
-                }
-                else
-                {
-                    boundaryPts[cit->cellIndex()] = 0;
-                }
-            }
+            boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
         }
         else
         {
@@ -1998,7 +1499,8 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
 void Foam::conformalVoronoiMesh::reindexDualVertices
 (
-    const Map<label>& dualPtIndexMap
+    const Map<label>& dualPtIndexMap,
+    labelList& boundaryPts
 )
 {
     for
@@ -2011,6 +1513,12 @@ void Foam::conformalVoronoiMesh::reindexDualVertices
         if (dualPtIndexMap.found(cit->cellIndex()))
         {
             cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
+            boundaryPts[cit->cellIndex()] =
+                max
+                (
+                    boundaryPts[cit->cellIndex()],
+                    boundaryPts[dualPtIndexMap[cit->cellIndex()]]
+                );
         }
     }
 }
@@ -2225,11 +1733,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
     bool includeEmptyPatches
 ) const
 {
-    const label defaultPatchIndex = createPatchInfo
-    (
-        patchNames,
-        patchDicts
-    );
+    const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);
 
     const label nPatches = patchNames.size();
 
@@ -2254,6 +1758,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
     List<DynamicList<bool> > indirectPatchFace(nPatches, DynamicList<bool>(0));
 
+
     faces.setSize(number_of_finite_edges());
     owner.setSize(number_of_finite_edges());
     neighbour.setSize(number_of_finite_edges());
@@ -2723,7 +2228,12 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
                     // If the two vertices are a pair, then the patch face is
                     // a desired one.
-                    if (!isPointPair(vA, vB))
+                    if
+                    (
+                        vA->boundaryPoint() && vB->boundaryPoint()
+                     && !ptPairs_.isPointPair(vA, vB)
+                     && !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
+                    )
                     {
                         indirectPatchFace[patchIndex].append(true);
                     }
@@ -2755,6 +2265,18 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
                     if (patchIndex != -1)
                     {
+//                        if
+//                        (
+//                            vA->boundaryPoint() && vB->boundaryPoint()
+//                         && !ptPairs_.isPointPair(vA, vB)
+//                        )
+//                        {
+//                            indirectPatchFace[patchIndex].append(true);
+//                        }
+//                        else
+//                        {
+//                            indirectPatchFace[patchIndex].append(false);
+//                        }
 //                        patchFaces[patchIndex].append(newDualFace);
 //                        patchOwners[patchIndex].append(own);
 //                        indirectPatchFace[patchIndex].append(false);
@@ -2771,8 +2293,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 //                        {
 //                            patchPPSlaves[patchIndex].append(vA->index());
 //                        }
-
-//                        baffleFaces[dualFaceI] = patchIndex;
                     }
 //                    else
                     {
@@ -2947,7 +2467,6 @@ void Foam::conformalVoronoiMesh::sortProcPatches
         faceList& faces = patchFaces[patchI];
         labelList& owner = patchOwners[patchI];
         DynamicList<label>& slaves = patchPointPairSlaves[patchI];
-
         DynamicList<Pair<labelPair> >& sortingIndices
             = patchSortingIndices[patchI];
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -62,11 +62,13 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
         if (Pstream::parRun())
         {
-            sync(decomposition_().procBounds());
+            sync(decomposition().procBounds());
         }
     }
     else
     {
+        ptPairs_.clear();
+
         // Rebuild, insert and store new surface conformation
         buildSurfaceConformation();
 
@@ -74,7 +76,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
         {
             if (Pstream::parRun())
             {
-                sync(decomposition_().procBounds());
+                sync(decomposition().procBounds());
             }
         }
 
@@ -358,10 +360,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
         }
 
+        DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
+
         insertSurfacePointPairs
         (
             surfaceHits,
-            "surfaceConformationLocations_initial.obj"
+            "surfaceConformationLocations_initial.obj",
+            pts
         );
 
         // In parallel, synchronise the edge trees
@@ -373,9 +378,21 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
         insertEdgePointGroups
         (
             featureEdgeHits,
-            "edgeConformationLocations_initial.obj"
+            "edgeConformationLocations_initial.obj",
+            pts
         );
 
+        pts.shrink();
+
+        Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
+
+        // Re-index the point pairs
+        ptPairs_.reIndex(oldToNewIndices);
+
+        //writePointPairs("pointPairs_initial.obj");
+
+        // Remove location from surface/edge tree
+
         timeCheck("After initial conformation");
 
         initialTotalHits = nSurfHits + nFeatEdHits;
@@ -550,215 +567,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             }
         }
 
-//        for
-//        (
-//            Delaunay::Finite_edges_iterator eit = finite_edges_begin();
-//            eit != finite_edges_end();
-//            ++eit
-//        )
-//        {
-//            Cell_handle c = eit->first;
-//            Vertex_handle vA = c->vertex(eit->second);
-//            Vertex_handle vB = c->vertex(eit->third);
-//
-//            if
-//            (
-//                vA->referred()
-//             || vB->referred()
-//            )
-//            {
-//                continue;
-//            }
-//
-//            if
-//            (
-//                (vA->internalPoint() && vA->referred())
-//             || (vB->internalPoint() && vB->referred())
-//            )
-//            {
-//                continue;
-//            }
-//
-//            if
-//            (
-//                (vA->internalPoint() && vB->externalBoundaryPoint())
-//             || (vB->internalPoint() && vA->externalBoundaryPoint())
-//             || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
-//            )
-//            {
-//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-//                pointIndexHit surfHit;
-//                label hitSurface;
-//
-//                geometryToConformTo_.findSurfaceNearest
-//                (
-//                    (
-//                        vA->internalPoint()
-//                      ? topoint(vA->point())
-//                      : topoint(vB->point())
-//                    ),
-//                    magSqr(topoint(vA->point()) - topoint(vB->point())),
-//                    surfHit,
-//                    hitSurface
-//                );
-//
-//                if (surfHit.hit())
-//                {
-//                    surfaceIntersections.append
-//                    (
-//                        pointIndexHitAndFeature(surfHit, hitSurface)
-//                    );
-//
-//                    addSurfaceAndEdgeHits
-//                    (
-//                        (
-//                            vA->internalPoint()
-//                          ? topoint(vA->point())
-//                          : topoint(vB->point())
-//                        ),
-//                        surfaceIntersections,
-//                        surfacePtReplaceDistCoeffSqr,
-//                        edgeSearchDistCoeffSqr,
-//                        surfaceHits,
-//                        featureEdgeHits,
-//                        surfaceToTreeShape,
-//                        edgeToTreeShape,
-//                        surfacePtToEdgePtDist,
-//                        false
-//                    );
-//                }
-//            }
-//        }
-
-        for
-        (
-            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-            cit != finite_cells_end();
-            ++cit
-        )
-        {
-            label nInternal = 0;
-            for (label vI = 0; vI < 4; vI++)
-            {
-                if (cit->vertex(vI)->internalPoint())
-                {
-                    nInternal++;
-                }
-            }
-
-            if (nInternal == 1 && cit->hasBoundaryPoint())
-            //if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
-            {
-                const Foam::point& pt = cit->dual();
-
-                const Foam::point cellCentre =
-                    topoint
-                    (
-                        CGAL::centroid
-                        (
-                            CGAL::Tetrahedron_3<baseK>
-                            (
-                                cit->vertex(0)->point(),
-                                cit->vertex(1)->point(),
-                                cit->vertex(2)->point(),
-                                cit->vertex(3)->point()
-                            )
-                        )
-                    );
-
-                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-                pointIndexHit surfHit;
-                label hitSurface;
-
-                geometryToConformTo_.findSurfaceNearestIntersection
-                (
-                    cellCentre,
-                    pt,
-                    surfHit,
-                    hitSurface
-                );
-
-                if (surfHit.hit())
-                {
-                    surfaceIntersections.append
-                    (
-                        pointIndexHitAndFeature(surfHit, hitSurface)
-                    );
-
-                    addSurfaceAndEdgeHits
-                    (
-                        pt,
-                        surfaceIntersections,
-                        surfacePtReplaceDistCoeffSqr,
-                        edgeSearchDistCoeffSqr,
-                        surfaceHits,
-                        featureEdgeHits,
-                        surfaceToTreeShape,
-                        edgeToTreeShape,
-                        surfacePtToEdgePtDist,
-                        false
-                    );
-                }
-            }
-        }
-
-//        for
-//        (
-//            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-//            cit != finite_cells_end();
-//            ++cit
-//        )
-//        {
-//            if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
-//            {
-//                const Foam::point& pt = cit->dual();
-//
-//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-//                pointIndexHit surfHit;
-//                label hitSurface;
-//
-//                geometryToConformTo_.findSurfaceNearest
-//                (
-//                    pt,
-//                    sqr(0.1*targetCellSize(pt)),
-//                    surfHit,
-//                    hitSurface
-//                );
-//
-//                if (!surfHit.hit())
-//                {
-//                    geometryToConformTo_.findSurfaceNearest
-//                    (
-//                        pt,
-//                        GREAT,
-//                        surfHit,
-//                        hitSurface
-//                    );
-//
-//                    if (surfHit.hit())
-//                    {
-//                        surfaceIntersections.append
-//                        (
-//                            pointIndexHitAndFeature(surfHit, hitSurface)
-//                        );
-//
-//                        addSurfaceAndEdgeHits
-//                        (
-//                            pt,
-//                            surfaceIntersections,
-//                            surfacePtReplaceDistCoeffSqr,
-//                            edgeSearchDistCoeffSqr,
-//                            surfaceHits,
-//                            featureEdgeHits,
-//                            surfaceToTreeShape,
-//                            edgeToTreeShape,
-//                            false
-//                        );
-//                    }
-//                }
-//            }
-//        }
-
         label nVerts = number_of_vertices();
         label nSurfHits = surfaceHits.size();
         label nFeatEdHits = featureEdgeHits.size();
@@ -789,10 +597,16 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                     synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
             }
 
+            DynamicList<Vb> pts
+            (
+                2*surfaceHits.size() + 3*featureEdgeHits.size()
+            );
+
             insertSurfacePointPairs
             (
                 surfaceHits,
-                "surfaceConformationLocations_" + name(iterationNo) + ".obj"
+                "surfaceConformationLocations_" + name(iterationNo) + ".obj",
+                pts
             );
 
             // In parallel, synchronise the edge trees
@@ -805,9 +619,19 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             insertEdgePointGroups
             (
                 featureEdgeHits,
-                "edgeConformationLocations_" + name(iterationNo) + ".obj"
+                "edgeConformationLocations_" + name(iterationNo) + ".obj",
+                pts
             );
 
+            pts.shrink();
+
+            Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
+
+            // Reindex the point pairs
+            ptPairs_.reIndex(oldToNewIndices);
+
+            //writePointPairs("pointPairs_" + name(iterationNo) + ".obj");
+
             if (Pstream::parRun())
             {
                 sync
@@ -1155,8 +979,6 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
          || is_infinite(fit->first->neighbor(fit->second))
          || !fit->first->hasInternalPoint()
          || !fit->first->neighbor(fit->second)->hasInternalPoint()
-//         || fit->first->hasFarPoint()
-//         || fit->first->neighbor(fit->second)->hasFarPoint()
         )
         {
             continue;
@@ -1209,14 +1031,13 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
             const scalar d = p.normalIntersect(r);
 
-            const Foam::point newPoint = vertex + d*n;
+            Foam::point newPoint = vertex + d*n;
 
             pointIndexHitAndFeature info;
-
             geometryToConformTo_.findSurfaceNearest
             (
                 newPoint,
-                surfaceSearchDistanceSqr(newPoint),
+                4.0*magSqr(newPoint - vertex),
                 info.first(),
                 info.second()
             );
@@ -1359,7 +1180,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
+         || (
+                !c1->internalOrBoundaryDualVertex()
+             || !c2->internalOrBoundaryDualVertex()
+            )
          || !c1->real() || !c2->real()
         )
         {
@@ -1369,19 +1193,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
 //        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
         Foam::point endPt = c1->dual();
 
-        if
-        (
-            magSqr(vert - c1->dual())
-          < magSqr(vert - c2->dual())
-        )
+        if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
         {
             endPt = c2->dual();
         }
 
+        if
+        (
+            magSqr(vert - endPt)
+          > magSqr(geometryToConformTo().globalBounds().mag())
+        )
+        {
+            continue;
+        }
+
         pointIndexHit surfHit;
         label hitSurface;
 
-        geometryToConformTo_.findSurfaceAnyIntersection
+        geometryToConformTo_.findSurfaceNearestIntersection
         (
             vert,
             endPt,
@@ -1401,18 +1230,49 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
 
             const vector& n = norm[0];
 
-            const scalar normalProtrusionDistance =
-                (endPt - surfHit.hitPoint()) & n;
+            const scalar normalProtrusionDistance
+            (
+                (endPt - surfHit.hitPoint()) & n
+            );
 
             if (normalProtrusionDistance > maxProtrusionDistance)
             {
-                surfHitLargest = surfHit;
-                hitSurfaceLargest = hitSurface;
+                const plane p(surfHit.hitPoint(), n);
+
+                const plane::ray r(endPt, -n);
+
+                const scalar d = p.normalIntersect(r);
+
+                Foam::point newPoint = endPt - d*n;
+
+                pointIndexHitAndFeature info;
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    newPoint,
+                    4.0*magSqr(newPoint - endPt),
+                    info.first(),
+                    info.second()
+                );
+
+                if (info.first().hit())
+                {
+                    if
+                    (
+                        surfaceLocationConformsToInside
+                        (
+                            pointIndexHitAndFeature(info.first(), info.second())
+                        )
+                    )
+                    {
+                        surfHitLargest = info.first();
+                        hitSurfaceLargest = info.second();
 
                         maxProtrusionDistance = normalProtrusionDistance;
                     }
                 }
             }
+        }
+    }
 
     // Relying on short-circuit evaluation to not call for hitPoint when this
     // is a miss
@@ -1465,7 +1325,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
+         || (
+                !c1->internalOrBoundaryDualVertex()
+             || !c2->internalOrBoundaryDualVertex()
+            )
          || !c1->real() || !c2->real()
         )
         {
@@ -1475,19 +1338,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 //        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
         Foam::point endPt = c1->dual();
 
-        if
-        (
-            magSqr(vert - c1->dual())
-          < magSqr(vert - c2->dual())
-        )
+        if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
         {
             endPt = c2->dual();
         }
 
+        if
+        (
+            magSqr(vert - endPt)
+          > magSqr(geometryToConformTo().globalBounds().mag())
+        )
+        {
+            continue;
+        }
+
         pointIndexHit surfHit;
         label hitSurface;
 
-        geometryToConformTo_.findSurfaceAnyIntersection
+        geometryToConformTo_.findSurfaceNearestIntersection
         (
             vert,
             endPt,
@@ -1507,19 +1375,46 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 
             const vector& n = norm[0];
 
-            scalar normalIncursionDistance = (endPt - surfHit.hitPoint()) & n;
+            scalar normalIncursionDistance
+            (
+                (endPt - surfHit.hitPoint()) & n
+            );
 
             if (normalIncursionDistance < minIncursionDistance)
             {
-                surfHitLargest = surfHit;
-                hitSurfaceLargest = hitSurface;
+                const plane p(surfHit.hitPoint(), n);
+
+                const plane::ray r(endPt, n);
+
+                const scalar d = p.normalIntersect(r);
+
+                Foam::point newPoint = endPt + d*n;
+
+                pointIndexHitAndFeature info;
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    newPoint,
+                    4.0*magSqr(newPoint - endPt),
+                    info.first(),
+                    info.second()
+                );
+
+                if (info.first().hit())
+                {
+                    if
+                    (
+                        surfaceLocationConformsToInside
+                        (
+                            pointIndexHitAndFeature(info.first(), info.second())
+                        )
+                    )
+                    {
+                        surfHitLargest = info.first();
+                        hitSurfaceLargest = info.second();
 
                         minIncursionDistance = normalIncursionDistance;
-
-                // Info<< nl << "# Incursion: " << endl;
-                // meshTools::writeOBJ(Info, vert);
-                // meshTools::writeOBJ(Info, edgeMid);
-                // Info<< "l Na Nb" << endl;
+                    }
+                }
             }
         }
     }
@@ -2162,9 +2057,11 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
 
         bool isNearFeaturePt = nearFeaturePt(surfPt);
 
+        bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
+
         bool isNearSurfacePt = nearSurfacePoint(surfHitI);
 
-        if (isNearFeaturePt || isNearSurfacePt)
+        if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
         {
             keepSurfacePoint = false;
         }
@@ -2219,6 +2116,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                             // feature edge, give control to edge control points
                             // instead, this will prevent "pits" forming.
 
+                            // Allow if different surfaces
+
+
                             keepSurfacePoint = false;
                         }
 
@@ -2298,11 +2198,13 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             surfaceHits.append(surfHitI);
             appendToSurfacePtTree(surfPt);
             surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
+
+//            addedPoints.write(surfPt);
+        }
+        else
+        {
+//            removedPoints.write(surfPt);
         }
-//        else
-//        {
-//            surfaceToTreeShape.remove();
-//        }
     }
 }
 
@@ -2338,7 +2240,9 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
                 Vb
                 (
                     vit->point(),
-                    vit->type()
+                    vit->index(),
+                    vit->type(),
+                    Pstream::myProcNo()
                 )
             );
         }
@@ -2362,24 +2266,40 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
 {
     Info<< nl << "Reinserting stored surface conformation" << endl;
 
-    const label preReinsertionSize(number_of_vertices());
+    Map<label> oldToNewIndices =
+        insertPointPairs(surfaceConformationVertices_, true, true);
 
-    // It is assumed that the stored surface conformation is on the correct
-    // processor and does not need distributed
-    rangeInsertWithInfo
+    ptPairs_.reIndex(oldToNewIndices);
+
+    PackedBoolList selectedElems(surfaceConformationVertices_.size(), true);
+
+    forAll(surfaceConformationVertices_, vI)
+    {
+        Vb& v = surfaceConformationVertices_[vI];
+        label& vIndex = v.index();
+
+        Map<label>::const_iterator iter = oldToNewIndices.find(vIndex);
+
+        if (iter != oldToNewIndices.end())
+        {
+            const label newIndex = iter();
+
+            if (newIndex != -1)
+            {
+                vIndex = newIndex;
+            }
+            else
+            {
+                selectedElems[vI] = false;
+            }
+        }
+    }
+
+    inplaceSubset<PackedBoolList, List<Vb> >
     (
-        surfaceConformationVertices_.begin(),
-        surfaceConformationVertices_.end(),
-        true
+        selectedElems,
+        surfaceConformationVertices_
     );
-
-    const label nInserted = label(number_of_vertices()) - preReinsertionSize;
-    const label nFailed = surfaceConformationVertices_.size() - nInserted;
-
-    Info<< "    " << returnReduce(nInserted, sumOp<label>())
-        << " points reinserted, failed to insert "
-        << returnReduce(nFailed, sumOp<label>())
-        << endl;
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C

@@ -29,6 +29,7 @@ License
 #include "tetrahedron.H"
 #include "const_circulator.H"
 #include "DelaunayMeshTools.H"
+#include "OBJstream.H"
 
 using namespace Foam::vectorTools;
 
@@ -225,54 +226,53 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
         vector masterPtVec(normalDir + nextNormalDir);
         masterPtVec /= mag(masterPtVec) + SMALL;
 
-        if (((normalDir ^ nextNormalDir) & edDir) < SMALL)
+        if
+        (
+            ((normalDir ^ nextNormalDir) & edDir) < SMALL
+         || mag(masterPtVec) < SMALL
+        )
         {
 //            Info<< "    IGNORE REGION" << endl;
             addedMasterPreviously = false;
 
-            if (circ.size() == 2 && mag((normal & nextNormal) - 1) < SMALL)
+            if
+            (
+                circ.size() == 2
+             && mag((normal & nextNormal) - 1) < SMALL
+            )
             {
-                // Add an extra point
-
                 const vector n = 0.5*(normal + nextNormal);
 
-                vector s = ppDist*(edDir ^ n);
+                const vector s = ppDist*(edDir ^ n);
 
-                createPointPair(ppDist, edgePt + s, n, pts);
-                createPointPair(ppDist, edgePt - s, n, pts);
-
-                break;
-            }
-
-            continue;
-        }
-
-        if (mag(masterPtVec) < SMALL)
+                if (volType == extendedFeatureEdgeMesh::BOTH)
                 {
-            if (circ.size() == 2)
-            {
-                // Add an extra point
-
-                // Average normal to remove any bias to one side, although as
-                // it is a flat edge, the normals should be essentially the same
-                const vector n = 0.5*(normal + nextNormal);
-
-                // Direction along the surface to the control point, sense of
-                // direction not important, as +s and -s can be used because it
-                // is a flat edge
-                vector s = ppDist*(edDir ^ n);
-
-                createPointPair(ppDist, edgePt + s, n, pts);
-                createPointPair(ppDist, edgePt - s, n, pts);
-
-                break;
+                    createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+                    createBafflePointPair(ppDist, edgePt - s, n, true, pts);
                 }
                 else
                 {
-//                Info<< "    IGNORE REGION" << endl;
-                addedMasterPreviously = false;
-                continue;
+                    WarningIn
+                    (
+                        "Foam::conformalVoronoiMesh::"
+                        "createEdgePointGroupByCirculating"
+                        "("
+                        "    const extendedFeatureEdgeMesh&,"
+                        "    const pointIndexHit&,"
+                        "    DynamicList<Vb>&"
+                        ")"
+                    )   << "Failed to insert flat/open edge as volType is "
+                        << extendedFeatureEdgeMesh::sideVolumeTypeNames_
+                           [
+                               volType
+                           ]
+                        << endl;
                 }
+
+                break;
+            }
+
+            continue;
         }
 
         const Foam::point masterPt = edgePt + ppDist*masterPtVec;
@@ -490,11 +490,19 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is an external edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
     const vector& nB = feNormals[edNormalIs[1]];
 
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
+        normalVolumeTypes[edNormalIs[0]];
+
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
+        normalVolumeTypes[edNormalIs[1]];
+
     if (areParallel(nA, nB))
     {
         // The normals are nearly parallel, so this is too sharp a feature to
@@ -524,13 +532,6 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
     // Convex. So refPt will be inside domain and hence a master point
     Foam::point refPt = edgePt - ppDist*refVec;
 
-    // Result when the points are eventually inserted.
-    // Add number_of_vertices() at insertion of first vertex to all numbers:
-    // pt           index type
-    // refPt        0     1
-    // reflectedA   1     0
-    // reflectedB   2     0
-
     // Insert the master point pairing the the first slave
 
     if (!geometryToConformTo_.inside(refPt))
@@ -559,7 +560,11 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
         (
             reflectedA,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeA == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
@@ -571,10 +576,26 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
         (
             reflectedB,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeB == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 1].index()
+    );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 2].index()
+    );
 }
 
 
@@ -591,11 +612,19 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is an external edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
     const vector& nB = feNormals[edNormalIs[1]];
 
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
+        normalVolumeTypes[edNormalIs[0]];
+
+//    const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
+//        normalVolumeTypes[edNormalIs[1]];
+
     if (areParallel(nA, nB))
     {
         // The normals are nearly parallel, so this is too sharp a feature to
@@ -705,14 +734,30 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
         (
             reflMasterPt,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeA == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
 
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 2].index(),
+        pts[pts.size() - 1].index()
+    );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 1].index()
+    );
+
     if (nAddPoints == 1)
     {
-        // One additinal point is the reflection of the slave point,
+        // One additional point is the reflection of the slave point,
         // i.e. the original reference point
         pts.append
         (
@@ -767,6 +812,8 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is a flat edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
@@ -781,8 +828,21 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
     // is a flat edge
     vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
 
-    createPointPair(ppDist, edgePt + s, n, pts);
-    createPointPair(ppDist, edgePt - s, n, pts);
+    if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::OUTSIDE)
+    {
+        createPointPair(ppDist, edgePt + s, -n, true, pts);
+        createPointPair(ppDist, edgePt - s, -n, true, pts);
+    }
+    else if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::BOTH)
+    {
+        createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+        createBafflePointPair(ppDist, edgePt - s, n, true, pts);
+    }
+    else
+    {
+        createPointPair(ppDist, edgePt + s, n, true, pts);
+        createPointPair(ppDist, edgePt - s, n, true, pts);
+    }
 }
 
 
@@ -814,16 +874,23 @@ void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
 
         plane facePlane(edgePt, n);
 
-        Foam::point pt1 = edgePt + s + ppDist*n;
-        Foam::point pt2 = edgePt - s + ppDist*n;
+        const label initialPtsSize = pts.size();
 
-        Foam::point pt3 = facePlane.mirror(pt1);
-        Foam::point pt4 = facePlane.mirror(pt2);
+        if
+        (
+            !geometryToConformTo_.inside(edgePt)
+        )
+        {
+            return;
+        }
 
-        pts.append(Vb(pt1, Vb::vtInternalSurface));
-        pts.append(Vb(pt2, Vb::vtInternalSurface));
-        pts.append(Vb(pt3, Vb::vtInternalSurface));
-        pts.append(Vb(pt4, Vb::vtInternalSurface));
+        createBafflePointPair(ppDist, edgePt - s, n, true, pts);
+        createBafflePointPair(ppDist, edgePt + s, n, false, pts);
+
+        for (label ptI = initialPtsSize; ptI < pts.size(); ++ptI)
+        {
+            pts[ptI].type() = Vb::vtInternalFeatureEdge;
+        }
     }
     else
     {
@@ -846,26 +913,245 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
 
     const scalar ppDist = pointPairDistance(edgePt);
 
+    const vector edDir = feMesh.edgeDirections()[edHit.index()];
+
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const labelList& normalDirs = feMesh.normalDirections()[edHit.index()];
 
-    Info<< edNormalIs.size() << endl;
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
-    // As this is a flat edge, there are two normals by definition
-    const vector& nA = feNormals[edNormalIs[0]];
-    const vector& nB = feNormals[edNormalIs[1]];
+    labelList nNormalTypes(4, 0);
 
-    // Average normal to remove any bias to one side, although as this
-    // is a flat edge, the normals should be essentially the same
-    const vector n = 0.5*(nA + nB);
+    forAll(edNormalIs, edgeNormalI)
+    {
+        const extendedFeatureEdgeMesh::sideVolumeType sType =
+            normalVolumeTypes[edNormalIs[edgeNormalI]];
 
-    // Direction along the surface to the control point, sense of edge
-    // direction not important, as +s and -s can be used because this
-    // is a flat edge
-    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+        nNormalTypes[sType]++;
+    }
 
-    createPointPair(ppDist, edgePt + s, n, pts);
-    createPointPair(ppDist, edgePt - s, n, pts);
+    if (nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 4)
+    {
+        label masterEdgeNormalIndex = -1;
+
+        forAll(edNormalIs, edgeNormalI)
+        {
+            const extendedFeatureEdgeMesh::sideVolumeType sType =
+                normalVolumeTypes[edNormalIs[edgeNormalI]];
+
+            if (sType == extendedFeatureEdgeMesh::BOTH)
+            {
+                masterEdgeNormalIndex = edgeNormalI;
+                break;
+            }
+        }
+
+        const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
+
+        label nDir = normalDirs[masterEdgeNormalIndex];
+
+        vector normalDir =
+            (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
+        normalDir *= nDir/mag(normalDir);
+
+        Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
+        Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
+
+        plane plane3(edgePt, normalDir);
+
+        Foam::point pt3 = plane3.mirror(pt1);
+        Foam::point pt4 = plane3.mirror(pt2);
+
+        pts.append
+        (
+            Vb
+            (
+                pt1,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+        pts.append
+        (
+            Vb
+            (
+                pt2,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt3,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt4,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
+    else if
+    (
+        nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 1
+     && nNormalTypes[extendedFeatureEdgeMesh::INSIDE] == 2
+    )
+    {
+        label masterEdgeNormalIndex = -1;
+
+        forAll(edNormalIs, edgeNormalI)
+        {
+            const extendedFeatureEdgeMesh::sideVolumeType sType =
+                normalVolumeTypes[edNormalIs[edgeNormalI]];
+
+            if (sType == extendedFeatureEdgeMesh::BOTH)
+            {
+                masterEdgeNormalIndex = edgeNormalI;
+                break;
+            }
+        }
+
+        const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
+
+        label nDir = normalDirs[masterEdgeNormalIndex];
+
+        vector normalDir =
+            (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
+        normalDir *= nDir/mag(normalDir);
+
+        const label nextNormalI =
+            (masterEdgeNormalIndex + 1) % edNormalIs.size();
+        if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
+        {
+            normalDir *= -1;
+        }
+
+        Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
+        Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
+
+        plane plane3(edgePt, normalDir);
+
+        Foam::point pt3 = plane3.mirror(pt1);
+        Foam::point pt4 = plane3.mirror(pt2);
+
+        pts.append
+        (
+            Vb
+            (
+                pt1,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+        pts.append
+        (
+            Vb
+            (
+                pt2,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt3,
+                vertexCount() + pts.size(),
+                Vb::vtExternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt4,
+                vertexCount() + pts.size(),
+                Vb::vtExternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
+
+
+//    // As this is a flat edge, there are two normals by definition
+//    const vector& nA = feNormals[edNormalIs[0]];
+//    const vector& nB = feNormals[edNormalIs[1]];
+//
+//    // Average normal to remove any bias to one side, although as this
+//    // is a flat edge, the normals should be essentially the same
+//    const vector n = 0.5*(nA + nB);
+//
+//    // Direction along the surface to the control point, sense of edge
+//    // direction not important, as +s and -s can be used because this
+//    // is a flat edge
+//    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+//
+//    createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+//    createBafflePointPair(ppDist, edgePt - s, n, true, pts);
 }
 
 
@@ -884,7 +1170,10 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
     const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
 
     // Insert the created points directly as already distributed.
-    this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices);
+    Map<label> oldToNewIndices =
+        this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices, true);
+
+    ftPtConformer_.reIndexPointPairs(oldToNewIndices);
 
     label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
     reduce(nFeatureVertices, sumOp<label>());

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -225,6 +225,7 @@ inline void Foam::conformalVoronoiMesh::createPointPair
     const scalar ppDist,
     const Foam::point& surfPt,
     const vector& n,
+    const bool ptPair,
     DynamicList<Vb>& pts
 ) const
 {
@@ -260,14 +261,14 @@ inline void Foam::conformalVoronoiMesh::createPointPair
             )
         );
 
-//        if (ptPair)
-//        {
-//            ptPairs_.addPointPair
-//            (
-//                pts[pts.size() - 2].index(),
-//                pts[pts.size() - 1].index() // external 0 -> slave
-//            );
-//        }
+        if (ptPair)
+        {
+            ptPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts[pts.size() - 1].index() // external 0 -> slave
+            );
+        }
     }
 //    else
 //    {
@@ -305,6 +306,7 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
     const scalar ppDist,
     const Foam::point& surfPt,
     const vector& n,
+    const bool ptPair,
     DynamicList<Vb>& pts
 ) const
 {
@@ -315,8 +317,8 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
         Vb
         (
             surfPt - ppDistn,
-            vertexCount() + 1 + pts.size(),
-            Vb::vtInternalSurface,
+            vertexCount() + pts.size(),
+            Vb::vtInternalSurfaceBaffle,
             Pstream::myProcNo()
         )
     );
@@ -326,11 +328,20 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
         Vb
         (
             surfPt + ppDistn,
-            vertexCount() + 1 + pts.size(),
-            Vb::vtInternalSurface,
+            vertexCount() + pts.size(),
+            Vb::vtExternalSurfaceBaffle,
             Pstream::myProcNo()
         )
     );
+
+    if (ptPair)
+    {
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
 }
 
 
@@ -341,8 +352,8 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
 {
     if
     (
-        !geometryToConformTo_.inside(pt)
-     || !geometryToConformTo_.globalBounds().contains(pt)
+        !geometryToConformTo_.globalBounds().contains(pt)
+     || !geometryToConformTo_.inside(pt)
     )
     {
         return false;
@@ -352,76 +363,6 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
 }
 
 
-inline bool Foam::conformalVoronoiMesh::isPointPair
-(
-    const Vertex_handle& vA,
-    const Vertex_handle& vB
-) const
-{
-    // Want to do this topologically, but problem if vertices are redistributed
-    // so that one of the point pair is one processor and the other is on
-    // another.
-
-    const Foam::point& ptA = topoint(vA->point());
-    const Foam::point& ptB = topoint(vB->point());
-
-    if
-    (
-        (
-            vA->type() == Vb::vtInternalSurface
-         && vB->type() == Vb::vtExternalSurface
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalSurface
-         && vA->type() == Vb::vtExternalSurface
-        )
-     ||
-        (
-            vA->type() == Vb::vtInternalFeatureEdge
-         && vB->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalFeatureEdge
-         && vA->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vA->type() == Vb::vtInternalSurface
-         && vB->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalSurface
-         && vA->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vA->type() == Vb::vtExternalSurface
-         && vB->type() == Vb::vtInternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtExternalSurface
-         && vA->type() == Vb::vtInternalFeatureEdge
-        )
-    )
-    {
-        const scalar distSqr = magSqr(ptA - ptB);
-
-        const scalar ppDistSqr = sqr(2*pointPairDistance(0.5*(ptA + ptB)));
-
-        if (distSqr > 1.001*ppDistSqr)
-        {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-
 inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 (
     const Delaunay::Finite_edges_iterator& eit

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -36,6 +36,8 @@ License
 #include "indexedVertexOps.H"
 #include "DelaunayMeshTools.H"
 #include "syncTools.H"
+#include "faceSet.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -114,7 +116,6 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
         wordList patchNames;
         PtrList<dictionary> patchDicts;
         pointField cellCentres;
-
         PackedBoolList boundaryFacesToRemove;
 
         calcDualMesh
@@ -789,7 +790,7 @@ void Foam::conformalVoronoiMesh::writeMesh
     const wordList& patchNames,
     const PtrList<dictionary>& patchDicts,
     const pointField& cellCentres,
-    const PackedBoolList& boundaryFacesToRemove
+    PackedBoolList& boundaryFacesToRemove
 ) const
 {
     if (foamyHexMeshControls().objOutput())
@@ -912,18 +913,64 @@ void Foam::conformalVoronoiMesh::writeMesh
 
 
 
-    // Add indirectPatchFaces to a face zone
+    Info<< indent << "Add pointZones" << endl;
     {
+        label sz = mesh.pointZones().size();
+
+        DynamicList<label> bPts(boundaryPts.size());
+
+        forAll(dualMeshPointTypeNames_, typeI)
+        {
+            forAll(boundaryPts, ptI)
+            {
+                const label& bPtType = boundaryPts[ptI];
+
+                if (bPtType == typeI)
+                {
+                    bPts.append(ptI);
+                }
+            }
+
+//            syncTools::syncPointList(mesh, bPts, maxEqOp<label>(), -1);
+
+            Info<< incrIndent << indent
+                << "Adding " << bPts.size()
+                << " points of type " << dualMeshPointTypeNames_.words()[typeI]
+                << decrIndent << endl;
+
+            mesh.pointZones().append
+            (
+                new pointZone
+                (
+                    dualMeshPointTypeNames_.words()[typeI],
+                    bPts,
+                    sz + typeI,
+                    mesh.pointZones()
+                )
+            );
+
+            bPts.clear();
+        }
+    }
+
+
+
+    // Add indirectPatchFaces to a face zone
+    Info<< indent << "Adding indirect patch faces set" << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh,
+        boundaryFacesToRemove,
+        orEqOp<unsigned int>()
+    );
+
     labelList addr(boundaryFacesToRemove.count());
     label count = 0;
 
     forAll(boundaryFacesToRemove, faceI)
     {
-            if
-            (
-                boundaryFacesToRemove[faceI]
-             && mesh.faceZones().whichZone(faceI) == -1
-            )
+        if (boundaryFacesToRemove[faceI])
         {
             addr[count++] = faceI;
         }
@@ -931,22 +978,18 @@ void Foam::conformalVoronoiMesh::writeMesh
 
     addr.setSize(count);
 
-        label sz = mesh.faceZones().size();
-        boolList flip(addr.size(), false);
-        mesh.faceZones().setSize(sz + 1);
-        mesh.faceZones().set
-        (
-            sz,
-            new faceZone
+    faceSet indirectPatchFaces
     (
+        mesh,
         "indirectPatchFaces",
         addr,
-                flip,
-                sz,
-                mesh.faceZones()
-            )
+        IOobject::AUTO_WRITE
     );
-    }
+
+    indirectPatchFaces.sync(mesh);
+
+
+    Info<< decrIndent;
 
     timeCheck("Before fvMesh filtering");
 
@@ -958,7 +1001,7 @@ void Foam::conformalVoronoiMesh::writeMesh
     {
         Info<< nl << "Filtering edges on polyMesh" << nl << endl;
 
-        meshFilter.reset(new polyMeshFilter(mesh));
+        meshFilter.reset(new polyMeshFilter(mesh, boundaryPts));
 
         // Filter small edges only. This reduces the number of faces so that
         // the face filtering is sped up.
@@ -966,25 +1009,48 @@ void Foam::conformalVoronoiMesh::writeMesh
         {
             const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
 
-            polyTopoChange meshMod(newMesh);
+            polyTopoChange meshMod(newMesh());
 
-            meshMod.changeMesh(mesh, false);
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+            polyMeshFilter::copySets(newMesh(), mesh);
         }
     }
 
     if (foamyHexMeshControls().filterFaces())
     {
+        labelIOList boundaryPtsIO
+        (
+            IOobject
+            (
+                "pointPriority",
+                instance,
+                time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(mesh.nPoints(), labelMin)
+        );
+
+        forAll(mesh.points(), ptI)
+        {
+            boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
+        }
+
+
         Info<< nl << "Filtering faces on polyMesh" << nl << endl;
 
-        meshFilter.reset(new polyMeshFilter(mesh));
+        meshFilter.reset(new polyMeshFilter(mesh, boundaryPtsIO));
 
         meshFilter().filter(nInitialBadFaces);
         {
             const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
 
-            polyTopoChange meshMod(newMesh);
+            polyTopoChange meshMod(newMesh());
 
-            meshMod.changeMesh(mesh, false);
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+            polyMeshFilter::copySets(newMesh(), mesh);
         }
     }
 
@@ -1005,158 +1071,43 @@ void Foam::conformalVoronoiMesh::writeMesh
             << endl;
     }
 
-
-//    volTensorField alignments
-//    (
-//        IOobject
-//        (
-//            "alignmentsField",
-//            runTime_.timeName(),
-//            runTime_,
-//            IOobject::NO_READ,
-//            IOobject::AUTO_WRITE
-//        ),
-//        mesh,
-//        tensor::zero
-//    );
-//
-//    forAll(mesh.cellCentres(), pI)
-//    {
-//        Vertex_handle nearV =
-//            nearest_vertex
-//            (
-//                toPoint<Point>(mesh.cellCentres()[pI])
-//            );
-//        alignments[pI] = nearV->alignment();
-//    }
-//    alignments.write();
-//
-//    {
-//        volVectorField alignmentx
-//        (
-//            IOobject
-//            (
-//                "alignmentsx",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmentx, aI)
-//        {
-//            alignmentx[aI] = alignments[aI].x();
-//        }
-//        alignmentx.write();
-//    }
-//    {
-//        volVectorField alignmenty
-//        (
-//            IOobject
-//            (
-//                "alignmentsy",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmenty, aI)
-//        {
-//            alignmenty[aI] = alignments[aI].y();
-//        }
-//        alignmenty.write();
-//    }
-//    {
-//        volVectorField alignmentz
-//        (
-//            IOobject
-//            (
-//                "alignmentsz",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmentz, aI)
-//        {
-//            alignmentz[aI] = alignments[aI].z();
-//        }
-//        alignmentz.write();
-//    }
-
-
-    labelIOList boundaryIOPts
+    {
+        pointScalarField boundaryPtsScalarField
         (
             IOobject
             (
-            "boundaryPoints",
+                "boundaryPoints_collapsed",
                 instance,
-            runTime_,
+                time(),
                 IOobject::NO_READ,
                 IOobject::AUTO_WRITE
             ),
-        boundaryPts
+            pointMesh::New(mesh),
+            scalar(labelMin)
         );
 
-    // Dump list of boundary points
-    forAll(mesh.boundaryMesh(), patchI)
-    {
-        const polyPatch& pp = mesh.boundaryMesh()[patchI];
+        labelIOList boundaryPtsIO
+        (
+            IOobject
+            (
+                "pointPriority",
+                instance,
+                time(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            labelList(mesh.nPoints(), labelMin)
+        );
 
-        if (!isA<coupledPolyPatch>(pp))
+        forAll(mesh.points(), ptI)
         {
-            forAll(pp, fI)
-            {
-                const face& boundaryFace = pp[fI];
-
-                forAll(boundaryFace, pI)
-                {
-                    const label boundaryPointI = boundaryFace[pI];
-
-                    boundaryIOPts[boundaryPointI] = boundaryPts[boundaryPointI];
-                }
-            }
-        }
+            boundaryPtsScalarField[ptI] = mesh.pointZones().whichZone(ptI);
+            boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
         }
 
-    boundaryIOPts.write();
-
-//    forAllConstIter(labelHashSet, pointsInPatch, pI)
-//    {
-//        const Foam::point& ptMaster = mesh.points()[pI.key()];
-//
-//        forAllConstIter(labelHashSet, pointsInPatch, ptI)
-//        {
-//            if (ptI.key() != pI.key())
-//            {
-//                const Foam::point& ptSlave = mesh.points()[ptI.key()];
-//
-//                const scalar dist = mag(ptMaster - ptSlave);
-//                if (ptMaster == ptSlave)
-//                {
-//                    Pout<< "Point(" << pI.key() << ") " << ptMaster
-//                        << " == "
-//                        << "(" << ptI.key() << ") " << ptSlave
-//                        << endl;
-//                }
-//                else if (dist == 0)
-//                {
-//                    Pout<< "Point(" << pI.key() << ") " << ptMaster
-//                        << " ~= "
-//                        << "(" << ptI.key() << ") " << ptSlave
-//                        << endl;
-//                }
-//            }
-//        }
-//    }
+        boundaryPtsScalarField.write();
+        boundaryPtsIO.write();
+    }
 
 //    writeCellSizes(mesh);
 
@@ -1164,7 +1115,7 @@ void Foam::conformalVoronoiMesh::writeMesh
 
 //    writeCellCentres(mesh);
 
-//    findRemainingProtrusionSet(mesh);
+    findRemainingProtrusionSet(mesh);
 }
 
 
@@ -1471,4 +1422,33 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findRemainingProtrusionSet
 }
 
 
+void Foam::conformalVoronoiMesh::writePointPairs
+(
+    const fileName& fName
+) const
+{
+    OBJstream os(fName);
+
+    for
+    (
+        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if (ptPairs_.isPointPair(vA, vB))
+        {
+            os.write
+            (
+                linePointRef(topoint(vA->point()), topoint(vB->point()))
+            );
+        }
+    }
+}
+
+
 // ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C

@@ -146,7 +146,7 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
             )
         );
 
-//        const label masterIndex = pts[pts.size() - 1].index();
+        const label masterIndex = pts.last().index();
 
         //meshTools::writeOBJ(strMasters, masterPt);
 
@@ -180,6 +180,11 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
                 )
             );
 
+            ftPtPairs_.addPointPair
+            (
+                masterIndex,
+                pts.last().index()
+            );
 
             //meshTools::writeOBJ(strSlaves, slavePt);
         }
@@ -525,7 +530,8 @@ Foam::featurePointConformer::featurePointConformer
     foamyHexMesh_(foamyHexMesh),
     foamyHexMeshControls_(foamyHexMesh.foamyHexMeshControls()),
     geometryToConformTo_(foamyHexMesh.geometryToConformTo()),
-    featurePointVertices_()
+    featurePointVertices_(),
+    ftPtPairs_(foamyHexMesh)
 {
     Info<< nl
         << "Conforming to feature points" << endl;
@@ -598,6 +604,30 @@ void Foam::featurePointConformer::distribute
     {
         featurePointVertices_[vI].procIndex() = Pstream::myProcNo();
     }
+
+    // Also update the feature point pairs
+}
+
+
+void Foam::featurePointConformer::reIndexPointPairs
+(
+    const Map<label>& oldToNewIndices
+)
+{
+    forAll(featurePointVertices_, vI)
+    {
+        const label currentIndex = featurePointVertices_[vI].index();
+
+        Map<label>::const_iterator newIndexIter =
+            oldToNewIndices.find(currentIndex);
+
+        if (newIndexIter != oldToNewIndices.end())
+        {
+            featurePointVertices_[vI].index() = newIndexIter();
+        }
+    }
+
+    ftPtPairs_.reIndex(oldToNewIndices);
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H

@@ -44,6 +44,7 @@ SourceFiles
 #include "DynamicList.H"
 #include "List.H"
 #include "extendedFeatureEdgeMesh.H"
+#include "pointPairs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,6 +84,9 @@ class featurePointConformer
         //  triangulation clear.
         List<Vb> featurePointVertices_;
 
+        //-
+        mutable pointPairs<Delaunay> ftPtPairs_;
+
 
     // Private Member Functions
 
@@ -163,12 +167,18 @@ public:
             //  triangulation.
             inline const List<Vb>& featurePointVertices() const;
 
+            //- Return the feature point pair table
+            inline const pointPairs<Delaunay>& featurePointPairs() const;
+
 
         // Edit
 
             //- Distribute the feature point vertices according to the
             //  supplied background mesh
             void distribute(const backgroundMeshDecomposition& decomposition);
+
+            //- Reindex the feature point pairs using the map.
+            void reIndexPointPairs(const Map<label>& oldToNewIndices);
 };
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H

@@ -28,5 +28,11 @@ const Foam::List<Vb>& Foam::featurePointConformer::featurePointVertices() const
     return featurePointVertices_;
 }
 
+const Foam::pointPairs<Delaunay>&
+Foam::featurePointConformer::featurePointPairs() const
+{
+    return ftPtPairs_;
+}
+
 
 // ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C

@@ -168,9 +168,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtA, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                internalPtA,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtAIndex(pts.last().index());
+
         const Foam::point internalPtB =
             concaveEdgeExternalPt
           - 2.0*planeB.distance(concaveEdgeExternalPt)
@@ -178,9 +186,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtB, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                internalPtB,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtBIndex(pts.last().index());
+
         // Add the external points
 
         Foam::point externalPtD;
@@ -288,7 +304,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtD, Vb::vtExternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtD,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtExternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtAIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -319,7 +347,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtE, Vb::vtExternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtE,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtExternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtBIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -328,9 +368,29 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(concaveEdgeExternalPt, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                concaveEdgeExternalPt,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        ftPtPairs_.addPointPair
+        (
+            internalPtBIndex,
+            pts.last().index()
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtAIndex,
+            pts.last().index()
+        );
+
+        const label concaveEdgeExternalPtIndex(pts.last().index());
+
         const scalar totalAngle = radToDeg
         (
             constant::mathematical::pi
@@ -377,7 +437,21 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(internalPtF, Vb::vtInternalFeaturePoint)
+                Vb
+                (
+                    internalPtF,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtInternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            const label internalPtFIndex(pts.last().index());
+
+            ftPtPairs_.addPointPair
+            (
+                concaveEdgeExternalPtIndex,
+                pts.last().index()
             );
 
             const Foam::point externalPtG =
@@ -386,7 +460,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(externalPtG, Vb::vtExternalFeaturePoint)
+                Vb
+                (
+                    externalPtG,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtExternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                internalPtFIndex,
+                pts.last().index()
             );
         }
 
@@ -520,9 +606,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtA, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                internalPtA,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtAIndex(pts.last().index());
+
         const Foam::point internalPtB =
             convexEdgeExternalPt
           + 2.0*planeB.distance(convexEdgeExternalPt)
@@ -530,9 +624,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtB, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                internalPtB,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtBIndex(pts.last().index());
+
         // Add the internal points
 
         Foam::point externalPtD;
@@ -643,7 +745,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtD, Vb::vtInternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtD,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtInternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtAIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -674,7 +788,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtE, Vb::vtInternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtE,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtInternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtBIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -683,7 +809,25 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(convexEdgeExternalPt, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                convexEdgeExternalPt,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtBIndex,
+            pts.last().index()
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtAIndex,
+            pts.last().index()
         );
 
         const scalar totalAngle = radToDeg
@@ -732,7 +876,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(internalPtF, Vb::vtExternalFeaturePoint)
+                Vb
+                (
+                    internalPtF,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtExternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts.last().index()
             );
 
             const Foam::point externalPtG =
@@ -741,7 +897,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(externalPtG, Vb::vtInternalFeaturePoint)
+                Vb
+                (
+                    externalPtG,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtInternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts.last().index()
             );
         }
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H

@@ -176,6 +176,9 @@ public:
         //- Does the Delaunay cell have a far point
         inline bool hasFarPoint() const;
 
+        //- Does the Delaunay cell have a referred point
+        inline bool hasReferredPoint() const;
+
         //- Does the Delaunay cell have a feature point
         inline bool hasFeaturePoint() const;
 
@@ -190,6 +193,8 @@ public:
         //- Does the Dual vertex form part of a processor patch
         inline bool parallelDualVertex() const;
 
+        inline Foam::label vertexLowestProc() const;
+
         //- Using the globalIndex object, return a list of four (sorted) global
         //  Delaunay vertex indices that uniquely identify this tet in parallel
         inline Foam::tetCell vertexGlobalIndices
@@ -216,6 +221,10 @@ public:
         //  least one Delaunay vertex outside and at least one inside
         inline bool boundaryDualVertex() const;
 
+        inline bool baffleSurfaceDualVertex() const;
+
+        inline bool baffleEdgeDualVertex() const;
+
         //- A dual vertex on a feature edge will result from this Delaunay cell
         inline bool featureEdgeDualVertex() const;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H

@@ -189,6 +189,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasFarPoint() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasReferredPoint() const
+{
+    return
+    (
+        this->vertex(0)->referred()
+     || this->vertex(1)->referred()
+     || this->vertex(2)->referred()
+     || this->vertex(3)->referred()
+    );
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::hasFeaturePoint() const
 {
@@ -278,6 +291,23 @@ inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
 }
 
 
+template<class Gt, class Cb>
+inline Foam::label CGAL::indexedCell<Gt, Cb>::vertexLowestProc() const
+{
+    Foam::label lowestProc = -1;
+
+    for (int i = 0; i < 4; ++i)
+    {
+        if (this->vertex(i)->referred())
+        {
+            lowestProc = min(lowestProc, this->vertex(i)->procIndex());
+        }
+    }
+
+    return lowestProc;
+}
+
+
 template<class Gt, class Cb>
 inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
 (
@@ -372,6 +402,14 @@ inline bool CGAL::indexedCell<Gt, Cb>::anyInternalOrBoundaryDualVertex() const
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
 {
+//    return
+//    (
+//        this->vertex(0)->boundaryPoint()
+//     && this->vertex(1)->boundaryPoint()
+//     && this->vertex(2)->boundaryPoint()
+//     && this->vertex(3)->boundaryPoint()
+//    );
+
     return
     (
         (
@@ -387,6 +425,62 @@ inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
         || this->vertex(3)->externalBoundaryPoint()
         )
     );
+
+//    Foam::label nBoundaryPoints = 0;
+//
+//    for (Foam::label i = 0; i < 4; ++i)
+//    {
+//        Vertex_handle v = this->vertex(i);
+//
+//        if (v->boundaryPoint())
+//        {
+//            nBoundaryPoints++;
+//        }
+//    }
+//
+//    return (nBoundaryPoints > 1);
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::baffleSurfaceDualVertex() const
+{
+    return
+    (
+        (
+           this->vertex(0)->internalBaffleSurfacePoint()
+        || this->vertex(1)->internalBaffleSurfacePoint()
+        || this->vertex(2)->internalBaffleSurfacePoint()
+        || this->vertex(3)->internalBaffleSurfacePoint()
+        )
+     && (
+           this->vertex(0)->externalBaffleSurfacePoint()
+        || this->vertex(1)->externalBaffleSurfacePoint()
+        || this->vertex(2)->externalBaffleSurfacePoint()
+        || this->vertex(3)->externalBaffleSurfacePoint()
+        )
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::baffleEdgeDualVertex() const
+{
+    return
+    (
+        (
+           this->vertex(0)->internalBaffleEdgePoint()
+        || this->vertex(1)->internalBaffleEdgePoint()
+        || this->vertex(2)->internalBaffleEdgePoint()
+        || this->vertex(3)->internalBaffleEdgePoint()
+        )
+     && (
+           this->vertex(0)->externalBaffleEdgePoint()
+        || this->vertex(1)->externalBaffleEdgePoint()
+        || this->vertex(2)->externalBaffleEdgePoint()
+        || this->vertex(3)->externalBaffleEdgePoint()
+        )
+    );
 }
 
 
@@ -400,6 +494,31 @@ inline bool CGAL::indexedCell<Gt, Cb>::featureEdgeDualVertex() const
      && this->vertex(2)->featureEdgePoint()
      && this->vertex(3)->featureEdgePoint()
     );
+//        (
+//            this->vertex(0)->featureEdgePoint()
+//         || this->vertex(1)->featureEdgePoint()
+//         || this->vertex(2)->featureEdgePoint()
+//         || this->vertex(3)->featureEdgePoint()
+//        )
+//     && (
+//            (
+//                this->vertex(0)->featureEdgePoint()
+//             || this->vertex(0)->featurePoint()
+//            )
+//         && (
+//                this->vertex(1)->featureEdgePoint()
+//             || this->vertex(1)->featurePoint()
+//            )
+//         && (
+//                this->vertex(2)->featureEdgePoint()
+//             || this->vertex(2)->featurePoint()
+//            )
+//         && (
+//                this->vertex(3)->featureEdgePoint()
+//             || this->vertex(3)->featurePoint()
+//            )
+//        )
+//    );
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H

@@ -243,8 +243,12 @@ public:
         inline bool surfacePoint() const;
 
         inline bool internalBoundaryPoint() const;
+        inline bool internalBaffleSurfacePoint() const;
+        inline bool internalBaffleEdgePoint() const;
 
         inline bool externalBoundaryPoint() const;
+        inline bool externalBaffleSurfacePoint() const;
+        inline bool externalBaffleEdgePoint() const;
 
         inline bool constrained() const;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C

@@ -30,13 +30,17 @@ License
 
 template<>
 const char*
-Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>::names[] =
+Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 15>::names[] =
 {
     "Unassigned",
     "Internal",
     "InternalNearBoundary",
     "InternalSurface",
+    "InternalSurfaceBaffle",
+    "ExternalSurfaceBaffle",
     "InternalFeatureEdge",
+    "InternalFeatureEdgeBaffle",
+    "ExternalFeatureEdgeBaffle",
     "InternalFeaturePoint",
     "ExternalSurface",
     "ExternalFeatureEdge",
@@ -45,7 +49,7 @@ Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>::names[] =
     "Constrained"
 };
 
-const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>
+const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 15>
 Foam::indexedVertexEnum::vertexTypeNames_;
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.H

@@ -53,13 +53,17 @@ public:
         vtInternal                  = 1,
         vtInternalNearBoundary      = 2,
         vtInternalSurface           = 3,
-        vtInternalFeatureEdge   = 4,
-        vtInternalFeaturePoint  = 5,
-        vtExternalSurface       = 6,
-        vtExternalFeatureEdge   = 7,
-        vtExternalFeaturePoint  = 8,
-        vtFar                   = 9,
-        vtConstrained           = 10
+        vtInternalSurfaceBaffle     = 4,
+        vtExternalSurfaceBaffle     = 5,
+        vtInternalFeatureEdge       = 6,
+        vtInternalFeatureEdgeBaffle = 7,
+        vtExternalFeatureEdgeBaffle = 8,
+        vtInternalFeaturePoint      = 9,
+        vtExternalSurface           = 10,
+        vtExternalFeatureEdge       = 11,
+        vtExternalFeaturePoint      = 12,
+        vtFar                       = 13,
+        vtConstrained               = 14
     };
 
     enum vertexMotion
@@ -68,7 +72,7 @@ public:
         movable = 1
     };
 
-    static const Foam::NamedEnum<vertexType, 11> vertexTypeNames_;
+    static const Foam::NamedEnum<vertexType, 15> vertexTypeNames_;
 
     static const Foam::NamedEnum<vertexMotion, 2> vertexMotionNames_;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H

@@ -307,6 +307,17 @@ inline bool CGAL::indexedVertex<Gt, Vb>::internalBoundaryPoint() const
     return type_ >= vtInternalSurface && type_ <= vtInternalFeaturePoint;
 }
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalBaffleSurfacePoint() const
+{
+    return (type_ == vtInternalSurfaceBaffle);
+}
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalBaffleEdgePoint() const
+{
+    return (type_ == vtInternalFeatureEdgeBaffle);
+}
 
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
@@ -314,6 +325,18 @@ inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
     return type_ >= vtExternalSurface && type_ <= vtExternalFeaturePoint;
 }
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::externalBaffleSurfacePoint() const
+{
+    return (type_ == vtExternalSurfaceBaffle);
+}
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::externalBaffleEdgePoint() const
+{
+    return (type_ == vtExternalFeatureEdgeBaffle);
+}
+
 
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::constrained() const

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -83,6 +83,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
 
             const pointField& surfPts = triSurf.points();
 
+            Info<< "    Checking " << surface.name() << endl;
+
+            label nBaffles = 0;
+
+            Info<< "        Index = " << surfaces_[s] << endl;
+            Info<< "        Offset = " << regionOffset_[s] << endl;
+
             forAll(triSurf, sI)
             {
                 const label patchID =
@@ -98,7 +105,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
                 {
                     sum += triSurf[sI].normal(surfPts);
                 }
+                else
+                {
+                    nBaffles++;
                 }
+            }
+            Info<< "        has " << nBaffles << " baffles out of "
+                << triSurf.size() << " triangles" << endl;
 
             totalTriangles += triSurf.size();
         }
@@ -153,8 +166,9 @@ void Foam::conformationSurfaces::readFeatures
     {
         const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
-        Info<< "    features: " << surface.name() << " of type "
-            << surface.type() << endl;
+        Info<< "    features: " << surface.name()
+            << " of type " << surface.type()
+            << ", id: " << featureIndex << endl;
 
         autoPtr<searchableSurfaceFeatures> ssFeatures
         (
@@ -210,7 +224,8 @@ void Foam::conformationSurfaces::readFeatures
     {
         fileName feMeshName(featureDict.lookup("extendedFeatureEdgeMesh"));
 
-        Info<< "    features: " << feMeshName << endl;
+        Info<< "    features: " << feMeshName << ", id: " << featureIndex
+            << endl;
 
         features_.set
         (
@@ -716,7 +731,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
 
             const searchableSurface& surface(allGeometry_[surfaces_[s]]);
 
-            if (!surface.hasVolumeType())
+            if
+            (
+                !surface.hasVolumeType()
+             //&& surfaceVolumeTests[s][i] == volumeType::UNKNOWN
+            )
             {
                 pointField sample(1, samplePts[i]);
                 scalarField nearestDistSqr(1, GREAT);
@@ -733,7 +752,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 findSurfaceNearestIntersection
                 (
                     samplePts[i],
-                    info[0].rawPoint(),
+                    info[0].rawPoint() - 1e-3*mag(hitDir)*hitDir,
                     surfHit,
                     hitSurface
                 );

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C

@@ -187,7 +187,7 @@ Foam::label Foam::autoDensity::recurseAndFill
     word recursionName
 ) const
 {
-    label maxDepth = 0;
+    label treeDepth = 0;
 
     for (direction i = 0; i < 8; i++)
     {
@@ -201,10 +201,10 @@ Foam::label Foam::autoDensity::recurseAndFill
         {
             if (levelLimit > 0)
             {
-                maxDepth =
+                treeDepth =
                     max
                     (
-                        maxDepth,
+                        treeDepth,
                         recurseAndFill
                         (
                             initialPoints,
@@ -229,10 +229,10 @@ Foam::label Foam::autoDensity::recurseAndFill
 
                 if (!fillBox(initialPoints, subBB, true))
                 {
-                    maxDepth =
+                    treeDepth =
                         max
                         (
-                            maxDepth,
+                            treeDepth,
                             recurseAndFill
                             (
                                 initialPoints,
@@ -259,10 +259,10 @@ Foam::label Foam::autoDensity::recurseAndFill
 
             if (!fillBox(initialPoints, subBB, false))
             {
-                maxDepth =
+                treeDepth =
                     max
                     (
-                        maxDepth,
+                        treeDepth,
                         recurseAndFill
                         (
                             initialPoints,
@@ -286,7 +286,7 @@ Foam::label Foam::autoDensity::recurseAndFill
         }
     }
 
-    return maxDepth + 1;
+    return treeDepth + 1;
 }
 
 
@@ -941,7 +941,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         Pout<< "    Filling box " << hierBB << endl;
     }
 
-    label maxDepth = recurseAndFill
+    label treeDepth = recurseAndFill
     (
         initialPoints,
         hierBB,
@@ -966,7 +966,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         << scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
         << " success rate" << nl
         << indent
-        << returnReduce(maxDepth, maxOp<label>())
+        << returnReduce(treeDepth, maxOp<label>())
         << " levels of recursion (maximum)"
         << decrIndent << decrIndent
         << endl;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C

@@ -49,7 +49,6 @@ void rayShooting::splitLine
 {
     Foam::point midPoint(l.centre());
     const scalar localCellSize(cellShapeControls().cellSize(midPoint));
-    const scalar lineLength(l.mag());
 
     const scalar minDistFromSurfaceSqr
     (
@@ -64,6 +63,8 @@ void rayShooting::splitLine
     )
     {
         // Add extra points if line length is much bigger than local cell size
+//        const scalar lineLength(l.mag());
+//
 //        if (lineLength > 4.0*localCellSize)
 //        {
 //            splitLine

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointPairs.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class Triangulation>
+inline Foam::Pair<Foam::labelPair>
+Foam::pointPairs<Triangulation>::orderPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    return
+    (
+        (vA < vB)
+      ? Pair<labelPair>(vA, vB)
+      : Pair<labelPair>(vB, vA)
+    );
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::findPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    if (vA == vB)
+    {
+        return false;
+    }
+    else if (find(orderPointPair(vA, vB)) == end())
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::insertPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+)
+{
+    if (vA == vB)
+    {
+        return false;
+    }
+
+    return insert(orderPointPair(vA, vB));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+Foam::pointPairs<Triangulation>::pointPairs(const Triangulation& triangulation)
+:
+    ptPairTable(),
+    triangulation_(triangulation)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+Foam::pointPairs<Triangulation>::~pointPairs()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+)
+{
+    return insertPointPair(vA, vB);
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const labelPair& master,
+    const DynamicList<labelPair>& slaves
+)
+{
+    forAll(slaves, sI)
+    {
+        const labelPair& slave = slaves[sI];
+
+        addPointPair(master, slave);
+    }
+
+    return true;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const label vA,
+    const label vB
+)
+{
+    const label procNo = Pstream::myProcNo();
+
+    labelPair a(vA, procNo);
+    labelPair b(vB, procNo);
+
+    return addPointPair(a, b);
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::isPointPair
+(
+    const Vertex_handle& vA,
+    const Vertex_handle& vB
+) const
+{
+    if (vA->boundaryPoint() && vB->boundaryPoint())
+    {
+        labelPair a(vA->index(), vA->procIndex());
+        labelPair b(vB->index(), vB->procIndex());
+
+        return findPointPair(a, b);
+    }
+
+    return false;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::isPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    return findPointPair(vA, vB);
+}
+
+
+template<class Triangulation>
+void Foam::pointPairs<Triangulation>::reIndex(const Map<label>& oldToNewIndices)
+{
+    pointPairs<Triangulation> newTable(triangulation_);
+
+    forAllConstIter(pointPairs, *this, iter)
+    {
+        Pair<labelPair> e = iter.key();
+
+        labelPair& start = e.first();
+        labelPair& end = e.second();
+
+        bool insert = true;
+
+        if (start.second() == Pstream::myProcNo())
+        {
+            Map<label>::const_iterator iter2 =
+                oldToNewIndices.find(start.first());
+
+            if (iter2 != oldToNewIndices.end())
+            {
+                if (iter2() != -1)
+                {
+                    start.first() = iter2();
+                }
+                else
+                {
+                    insert = false;
+                }
+            }
+        }
+
+        if (end.second() == Pstream::myProcNo())
+        {
+            Map<label>::const_iterator iter2 =
+                oldToNewIndices.find(end.first());
+
+            if (iter2 != oldToNewIndices.end())
+            {
+                if (iter2() != -1)
+                {
+                    end.first() = iter2();
+                }
+                else
+                {
+                    insert = false;
+                }
+            }
+        }
+
+        if (insert)
+        {
+            if (e.first() < e.second())
+            {
+                newTable.insert(e);
+            }
+            else if (e.first() > e.second())
+            {
+                newTable.insert(reverse(e));
+            }
+        }
+    }
+
+    this->transfer(newTable);
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    pointPairs
+
+Description
+    HashSet of unique edges. The edges are stored as a pair of pairs:
+
+        ( (local index, processor index) (local index, processor index) )
+
+    e.g.,
+
+        ( (0 1) (3 1) )
+        ( (0 2) (5 1) )
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointPairs_H
+#define pointPairs_H
+
+#include "labelPair.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef HashSet
+<
+    Pair<labelPair>,
+    FixedList<labelPair, 2>::Hash<>
+> ptPairTable;
+
+/*---------------------------------------------------------------------------*\
+                         Class pointPairs Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Triangulation>
+class pointPairs
+:
+    public ptPairTable
+{
+    // Private typedefs
+
+        typedef typename Triangulation::Vertex_handle Vertex_handle;
+
+
+    // Private data
+
+        const Triangulation& triangulation_;
+
+
+    // Private Member Functions
+
+        inline Pair<labelPair> orderPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        ) const;
+
+        inline bool insertPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        );
+
+        inline bool findPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from triangulation
+        pointPairs(const Triangulation& triangulation);
+
+
+    //- Destructor
+    ~pointPairs();
+
+
+    // Member Functions
+
+        // Access
+
+            inline bool isPointPair
+            (
+                const Vertex_handle& vA,
+                const Vertex_handle& vB
+            ) const;
+
+            inline bool isPointPair
+            (
+                const labelPair& vA,
+                const labelPair& vB
+            ) const;
+
+
+        // Edit
+
+            inline bool addPointPair
+            (
+                const labelPair& vA,
+                const labelPair& vB
+            );
+
+            inline bool addPointPair
+            (
+                const labelPair& master,
+                const DynamicList<labelPair>& slaves
+            );
+
+            inline bool addPointPair
+            (
+                const label vA,
+                const label vB
+            );
+
+            void reIndex(const Map<label>& oldToNewIndices);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "pointPairs.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C

@@ -104,11 +104,14 @@ Foam::searchableBoxFeatures::features() const
     edgeNormals[10][0] = 1; edgeNormals[10][1] = 3;
     edgeNormals[11][0] = 3; edgeNormals[11][1] = 0;
 
+    tmp<pointField> surfacePointsTmp(surface().points());
+    pointField& surfacePoints = surfacePointsTmp();
+
     forAll(edgeDirections, eI)
     {
         edgeDirections[eI] =
-            surface().points()()[treeBoundBox::edges[eI].end()]
-          - surface().points()()[treeBoundBox::edges[eI].start()];
+            surfacePoints[treeBoundBox::edges[eI].end()]
+          - surfacePoints[treeBoundBox::edges[eI].start()];
 
         normalDirections[eI] = labelList(2, 0);
         for (label j = 0; j < 2; ++j)
@@ -116,8 +119,8 @@ Foam::searchableBoxFeatures::features() const
             const vector cross =
                 (faceNormals[edgeNormals[eI][j]] ^ edgeDirections[eI]);
             const vector fC0tofE0 =
-                0.5*(max(surface().points()() + min(surface().points()())))
-              - surface().points()()[treeBoundBox::edges[eI].start()];
+                0.5*(max(surfacePoints + min(surfacePoints)))
+              - surfacePoints[treeBoundBox::edges[eI].start()];
 
             normalDirections[eI][j] =
                 (