Merge branch 'particleInteractions'

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -64,11 +64,17 @@ Usage
 #include "IOobjectList.H"
 #include "domainDecomposition.H"
 #include "labelIOField.H"
+#include "labelFieldIOField.H"
 #include "scalarIOField.H"
+#include "scalarFieldIOField.H"
 #include "vectorIOField.H"
+#include "vectorFieldIOField.H"
 #include "sphericalTensorIOField.H"
+#include "sphericalTensorFieldIOField.H"
 #include "symmTensorIOField.H"
+#include "symmTensorFieldIOField.H"
 #include "tensorIOField.H"
+#include "tensorFieldIOField.H"
 #include "pointFields.H"
 
 #include "readFields.H"
@@ -374,20 +380,42 @@ int main(int argc, char *argv[])
     PtrList< List<SLList<indexedParticle*>*> > cellParticles(cloudDirs.size());
 
     PtrList<PtrList<labelIOField> > lagrangianLabelFields(cloudDirs.size());
+    PtrList<PtrList<labelIOFieldField> > lagrangianLabelFieldFields
+    (
+        cloudDirs.size()
+    );
     PtrList<PtrList<scalarIOField> > lagrangianScalarFields(cloudDirs.size());
+    PtrList<PtrList<scalarIOFieldField> > lagrangianScalarFieldFields
+    (
+        cloudDirs.size()
+    );
     PtrList<PtrList<vectorIOField> > lagrangianVectorFields(cloudDirs.size());
+    PtrList<PtrList<vectorIOFieldField> > lagrangianVectorFieldFields
+    (
+        cloudDirs.size()
+    );
     PtrList<PtrList<sphericalTensorIOField> > lagrangianSphericalTensorFields
     (
         cloudDirs.size()
     );
+    PtrList<PtrList<sphericalTensorIOFieldField> >
+        lagrangianSphericalTensorFieldFields(cloudDirs.size());
     PtrList<PtrList<symmTensorIOField> > lagrangianSymmTensorFields
     (
         cloudDirs.size()
     );
+    PtrList<PtrList<symmTensorIOFieldField> > lagrangianSymmTensorFieldFields
+    (
+        cloudDirs.size()
+    );
     PtrList<PtrList<tensorIOField> > lagrangianTensorFields
     (
         cloudDirs.size()
     );
+    PtrList<PtrList<tensorIOFieldField> > lagrangianTensorFieldFields
+    (
+        cloudDirs.size()
+    );
 
     label cloudI = 0;
 
@@ -487,6 +515,13 @@ int main(int argc, char *argv[])
                 lagrangianLabelFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianLabelFieldFields
+            );
+
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -494,6 +529,14 @@ int main(int argc, char *argv[])
                 lagrangianScalarFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianScalarFieldFields
+            );
+
+
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -501,6 +544,13 @@ int main(int argc, char *argv[])
                 lagrangianVectorFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianVectorFieldFields
+            );
+
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -508,6 +558,13 @@ int main(int argc, char *argv[])
                 lagrangianSphericalTensorFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianSphericalTensorFieldFields
+            );
+
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -515,6 +572,13 @@ int main(int argc, char *argv[])
                 lagrangianSymmTensorFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianSymmTensorFieldFields
+            );
+
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -522,6 +586,13 @@ int main(int argc, char *argv[])
                 lagrangianTensorFields
             );
 
+            lagrangianFieldDecomposer::readFieldFields
+            (
+                cloudI,
+                lagrangianObjects,
+                lagrangianTensorFieldFields
+            );
+
             cloudI++;
         }
     }
@@ -529,11 +600,17 @@ int main(int argc, char *argv[])
     lagrangianPositions.setSize(cloudI);
     cellParticles.setSize(cloudI);
     lagrangianLabelFields.setSize(cloudI);
+    lagrangianLabelFieldFields.setSize(cloudI);
     lagrangianScalarFields.setSize(cloudI);
+    lagrangianScalarFieldFields.setSize(cloudI);
     lagrangianVectorFields.setSize(cloudI);
+    lagrangianVectorFieldFields.setSize(cloudI);
     lagrangianSphericalTensorFields.setSize(cloudI);
+    lagrangianSphericalTensorFieldFields.setSize(cloudI);
     lagrangianSymmTensorFields.setSize(cloudI);
+    lagrangianSymmTensorFieldFields.setSize(cloudI);
     lagrangianTensorFields.setSize(cloudI);
+    lagrangianTensorFieldFields.setSize(cloudI);
 
 
     // Any uniform data to copy/link?
@@ -725,11 +802,17 @@ int main(int argc, char *argv[])
                 if
                 (
                     lagrangianLabelFields[cloudI].size()
+                 || lagrangianLabelFieldFields[cloudI].size()
                  || lagrangianScalarFields[cloudI].size()
+                 || lagrangianScalarFieldFields[cloudI].size()
                  || lagrangianVectorFields[cloudI].size()
+                 || lagrangianVectorFieldFields[cloudI].size()
                  || lagrangianSphericalTensorFields[cloudI].size()
+                 || lagrangianSphericalTensorFieldFields[cloudI].size()
                  || lagrangianSymmTensorFields[cloudI].size()
+                 || lagrangianSymmTensorFieldFields[cloudI].size()
                  || lagrangianTensorFields[cloudI].size()
+                 || lagrangianTensorFieldFields[cloudI].size()
                 )
                 {
                     fieldDecomposer.decomposeFields
@@ -737,31 +820,61 @@ int main(int argc, char *argv[])
                         cloudDirs[cloudI],
                         lagrangianLabelFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianLabelFieldFields[cloudI]
+                    );
                     fieldDecomposer.decomposeFields
                     (
                         cloudDirs[cloudI],
                         lagrangianScalarFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianScalarFieldFields[cloudI]
+                    );
                     fieldDecomposer.decomposeFields
                     (
                         cloudDirs[cloudI],
                         lagrangianVectorFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianVectorFieldFields[cloudI]
+                    );
                     fieldDecomposer.decomposeFields
                     (
                         cloudDirs[cloudI],
                         lagrangianSphericalTensorFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianSphericalTensorFieldFields[cloudI]
+                    );
                     fieldDecomposer.decomposeFields
                     (
                         cloudDirs[cloudI],
                         lagrangianSymmTensorFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianSymmTensorFieldFields[cloudI]
+                    );
                     fieldDecomposer.decomposeFields
                     (
                         cloudDirs[cloudI],
                         lagrangianTensorFields[cloudI]
                     );
+                    fieldDecomposer.decomposeFieldFields
+                    (
+                        cloudDirs[cloudI],
+                        lagrangianTensorFieldFields[cloudI]
+                    );
                 }
             }
         }

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H

@@ -37,6 +37,7 @@ SourceFiles
 #define lagrangianFieldDecomposer_H
 
 #include "Cloud.H"
+#include "IOFieldField.H"
 #include "indexedParticle.H"
 #include "passiveParticle.H"
 
@@ -102,6 +103,19 @@ public:
 //            PtrList<IOField<Type> >& lagrangianFields
         );
 
+        template<class Type>
+        static void readFieldFields
+        (
+            const label cloudI,
+            const IOobjectList& lagrangianObjects,
+            PtrList
+            <
+                PtrList<IOFieldField<Field<Type>, Type> >
+            >& lagrangianFields
+//            PtrList<IOFieldField<Field<Type>, Type > >& lagrangianFields
+        );
+
+
         //- Decompose volume field
         template<class Type>
         tmp<IOField<Type> > decomposeField
@@ -110,12 +124,27 @@ public:
             const IOField<Type>& field
         ) const;
 
+        template<class Type>
+        tmp<IOFieldField<Field<Type>, Type> > decomposeFieldField
+        (
+            const word& cloudName,
+            const IOFieldField<Field<Type>, Type>& field
+        ) const;
+
+
         template<class GeoField>
         void decomposeFields
         (
             const word& cloudName,
             const PtrList<GeoField>& fields
         ) const;
+
+        template<class GeoField>
+        void decomposeFieldFields
+        (
+            const word& cloudName,
+            const PtrList<GeoField>& fields
+        ) const;
 };
 
 

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C

@@ -63,6 +63,56 @@ void Foam::lagrangianFieldDecomposer::readFields
 }
 
 
+template<class Type>
+void Foam::lagrangianFieldDecomposer::readFieldFields
+(
+    const label cloudI,
+    const IOobjectList& lagrangianObjects,
+    PtrList<PtrList<IOFieldField<Field<Type>, Type> > >& lagrangianFields
+)
+{
+    // Search list of objects for lagrangian fields
+    IOobjectList lagrangianTypeObjectsA
+    (
+        lagrangianObjects.lookupClass(IOField<Field<Type> >::typeName)
+    );
+
+    IOobjectList lagrangianTypeObjectsB
+    (
+        lagrangianObjects.lookupClass(IOFieldField<Field<Type>, Type>::typeName)
+    );
+
+    lagrangianFields.set
+    (
+        cloudI,
+        new PtrList<IOFieldField<Field<Type>, Type> >
+        (
+            lagrangianTypeObjectsA.size() + lagrangianTypeObjectsB.size()
+        )
+    );
+
+    label lagrangianFieldi=0;
+
+    forAllIter(IOobjectList, lagrangianTypeObjectsA, iter)
+    {
+        lagrangianFields[cloudI].set
+        (
+            lagrangianFieldi++,
+            new IOFieldField<Field<Type>, Type>(*iter())
+        );
+    }
+
+    forAllIter(IOobjectList, lagrangianTypeObjectsB, iter)
+    {
+        lagrangianFields[cloudI].set
+        (
+            lagrangianFieldi++,
+            new IOFieldField<Field<Type>, Type>(*iter())
+        );
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::IOField<Type> >
 Foam::lagrangianFieldDecomposer::decomposeField
@@ -94,6 +144,37 @@ Foam::lagrangianFieldDecomposer::decomposeField
 }
 
 
+template<class Type>
+Foam::tmp<Foam::IOFieldField<Foam::Field<Type>, Type> >
+Foam::lagrangianFieldDecomposer::decomposeFieldField
+(
+    const word& cloudName,
+    const IOFieldField<Field<Type>, Type>& field
+) const
+{
+    // Create and map the internal field values
+    Field<Field<Type> > procField(field, particleIndices_);
+
+    // Create the field for the processor
+    return tmp<IOFieldField<Field<Type>, Type> >
+    (
+        new IOFieldField<Field<Type>, Type>
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                cloud::prefix/cloudName,
+                procMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            procField
+        )
+    );
+}
+
+
 template<class GeoField>
 void Foam::lagrangianFieldDecomposer::decomposeFields
 (
@@ -111,4 +192,21 @@ void Foam::lagrangianFieldDecomposer::decomposeFields
 }
 
 
+template<class GeoField>
+void Foam::lagrangianFieldDecomposer::decomposeFieldFields
+(
+    const word& cloudName,
+    const PtrList<GeoField>& fields
+) const
+{
+    if (particleIndices_.size())
+    {
+        forAll(fields, fieldI)
+        {
+            decomposeFieldField(cloudName, fields[fieldI])().write();
+        }
+    }
+}
+
+
 // ************************************************************************* //

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -402,6 +402,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<label>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                     reconstructLagrangianFields<scalar>
                     (
                         cloudName,
@@ -409,6 +416,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<scalar>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                     reconstructLagrangianFields<vector>
                     (
                         cloudName,
@@ -416,6 +430,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<vector>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                     reconstructLagrangianFields<sphericalTensor>
                     (
                         cloudName,
@@ -423,6 +444,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<sphericalTensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                     reconstructLagrangianFields<symmTensor>
                     (
                         cloudName,
@@ -430,6 +458,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<symmTensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                     reconstructLagrangianFields<tensor>
                     (
                         cloudName,
@@ -437,6 +472,13 @@ int main(int argc, char *argv[])
                         procMeshes.meshes(),
                         sprayObjs
                     );
+                    reconstructLagrangianFieldFields<tensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
                 }
             }
             else

src/OpenFOAM/Make/files

@@ -497,10 +497,15 @@ $(Fields)/scalarField/scalarFieldIOField.C
 $(Fields)/vectorField/vectorIOField.C
 $(Fields)/vectorField/vectorFieldIOField.C
 $(Fields)/vector2DField/vector2DIOField.C
+$(Fields)/vector2DField/vector2DFieldIOField.C
 $(Fields)/sphericalTensorField/sphericalTensorIOField.C
+$(Fields)/sphericalTensorField/sphericalTensorFieldIOField.C
 $(Fields)/diagTensorField/diagTensorIOField.C
+$(Fields)/diagTensorField/diagTensorFieldIOField.C
 $(Fields)/symmTensorField/symmTensorIOField.C
+$(Fields)/symmTensorField/symmTensorFieldIOField.C
 $(Fields)/tensorField/tensorIOField.C
+$(Fields)/tensorField/tensorFieldIOField.C
 $(Fields)/transformField/transformField.C
 
 pointPatchFields = fields/pointPatchFields

src/OpenFOAM/fields/Fields/diagTensorField/diagTensorFieldIOField.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Description
+    diagTensorField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "diagTensorFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        diagTensorFieldIOField,
+        "diagTensorFieldField",
+        0
+    );
+
+    defineTemplateTypeNameAndDebugWithName
+    (
+        diagTensorIOFieldField,
+        "diagTensorCompactFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/diagTensorField/diagTensorFieldIOField.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Typedef
+    Foam::diagTensorFieldIOField
+
+Description
+    diagTensorFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diagTensorFieldIOField_H
+#define diagTensorFieldIOField_H
+
+#include "diagTensorField.H"
+#include "IOField.H"
+#include "IOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<diagTensorField> diagTensorFieldIOField;
+    typedef IOFieldField<diagTensorField, diagTensor> diagTensorIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/sphericalTensorField/sphericalTensorFieldIOField.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Description
+    sphericalTensorField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sphericalTensorFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        sphericalTensorFieldIOField,
+        "sphericalTensorFieldField",
+        0
+    );
+
+    defineTemplateTypeNameAndDebugWithName
+    (
+        sphericalTensorIOFieldField,
+        "sphericalTensorCompactFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/sphericalTensorField/sphericalTensorFieldIOField.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Typedef
+    Foam::sphericalTensorFieldIOField
+
+Description
+    sphericalTensorFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sphericalTensorFieldIOField_H
+#define sphericalTensorFieldIOField_H
+
+#include "sphericalTensorField.H"
+#include "IOField.H"
+#include "IOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<sphericalTensorField> sphericalTensorFieldIOField;
+
+    typedef IOFieldField<sphericalTensorField, sphericalTensor>
+    sphericalTensorIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/sphericalTensorField/sphericalTensorIOField.C

@@ -32,7 +32,12 @@ Description
 
 namespace Foam
 {
-    defineTemplateTypeNameAndDebugWithName(sphericalTensorIOField, "sphericalTensorField", 0);
+    defineTemplateTypeNameAndDebugWithName
+    (
+        sphericalTensorIOField,
+        "sphericalTensorField",
+        0
+    );
 }
 
 // ************************************************************************* //

src/OpenFOAM/fields/Fields/symmTensorField/symmTensorFieldIOField.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Description
+    symmTensorField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "symmTensorFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        symmTensorFieldIOField,
+        "symmTensorFieldField",
+        0
+    );
+
+    defineTemplateTypeNameAndDebugWithName
+    (
+        symmTensorIOFieldField,
+        "symmTensorCompactFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/symmTensorField/symmTensorFieldIOField.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Typedef
+    Foam::symmTensorFieldIOField
+
+Description
+    symmTensorFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef symmTensorFieldIOField_H
+#define symmTensorFieldIOField_H
+
+#include "symmTensorField.H"
+#include "IOField.H"
+#include "IOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<symmTensorField> symmTensorFieldIOField;
+    typedef IOFieldField<symmTensorField, symmTensor> symmTensorIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/tensorField/tensorFieldIOField.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Description
+    tensorField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tensorFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        tensorFieldIOField,
+        "tensorFieldField",
+        0
+    );
+
+    defineTemplateTypeNameAndDebugWithName
+    (
+        tensorIOFieldField,
+        "tensorCompactFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/tensorField/tensorFieldIOField.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Typedef
+    Foam::tensorFieldIOField
+
+Description
+    tensorFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tensorFieldIOField_H
+#define tensorFieldIOField_H
+
+#include "tensorField.H"
+#include "IOField.H"
+#include "IOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<tensorField> tensorFieldIOField;
+    typedef IOFieldField<tensorField, tensor> tensorIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/vector2DField/vector2DFieldIOField.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Description
+    vector2DField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vector2DFieldIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        vector2DFieldIOField,
+        "vector2DFieldField",
+        0
+    );
+
+    defineTemplateTypeNameAndDebugWithName
+    (
+        vector2DIOFieldField,
+        "vector2DCompactFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/vector2DField/vector2DFieldIOField.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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/>.
+
+Typedef
+    Foam::vector2DFieldIOField
+
+Description
+    vector2DFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vector2DFieldIOField_H
+#define vector2DFieldIOField_H
+
+#include "vector2DField.H"
+#include "IOField.H"
+#include "IOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<vector2DField> vector2DFieldIOField;
+    typedef IOFieldField<vector2DField, vector2D> vector2DIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/vector2DField/vector2DIOField.C

@@ -32,7 +32,12 @@ Description
 
 namespace Foam
 {
-    defineTemplateTypeNameAndDebugWithName(vector2DIOField, "vector2DField", 0);
+    defineTemplateTypeNameAndDebugWithName
+    (
+        vector2DIOField,
+        "vector2DField",
+        0
+    );
 }
 
 // ************************************************************************* //

src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.C

@@ -27,6 +27,7 @@ License
 #include "ListOps.H"
 #include "Pstream.H"
 #include "commSchedule.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -42,6 +43,8 @@ Foam::labelList Foam::ProcessorTopology<Patch, ProcPatch>::procNeighbours
 
     label maxNb = 0;
 
+    boolList isNeighbourProc(Pstream::nProcs(), false);
+
     forAll(patches, patchi)
     {
         const Patch& patch = patches[patchi];
@@ -51,19 +54,34 @@ Foam::labelList Foam::ProcessorTopology<Patch, ProcPatch>::procNeighbours
             const ProcPatch& procPatch =
                 refCast<const ProcPatch>(patch);
 
+            label pNeighbProcNo = procPatch.neighbProcNo();
+
+            if (!isNeighbourProc[pNeighbProcNo])
+            {
                 nNeighbours++;
 
                 maxNb = max(maxNb, procPatch.neighbProcNo());
+
+                isNeighbourProc[pNeighbProcNo] = true;
+            }
         }
     }
 
-    labelList neighbours(nNeighbours);
+    labelList neighbours(nNeighbours, -1);
+
+    nNeighbours = 0;
+
+    forAll(isNeighbourProc, procI)
+    {
+        if (isNeighbourProc[procI])
+        {
+            neighbours[nNeighbours++] = procI;
+        }
+    }
 
     procPatchMap_.setSize(maxNb + 1);
     procPatchMap_ = -1;
 
-    nNeighbours = 0;
-
     forAll(patches, patchi)
     {
         const Patch& patch = patches[patchi];
@@ -73,8 +91,6 @@ Foam::labelList Foam::ProcessorTopology<Patch, ProcPatch>::procNeighbours
             const ProcPatch& procPatch =
                 refCast<const ProcPatch>(patch);
 
-            neighbours[nNeighbours++] = procPatch.neighbProcNo();
-
             // Construct reverse map
             procPatchMap_[procPatch.neighbProcNo()] = patchi;
         }

src/OpenFOAM/meshes/ProcessorTopology/ProcessorTopology.H

@@ -30,6 +30,9 @@ Description
 
     *this[procI] gives the list of neighbouring processors.
 
+    TODO: This does not currently correctly support multiple processor
+    patches connecting two processors.
+
 SourceFiles
     ProcessorTopology.C
 

src/lagrangian/basic/Cloud/Cloud.C

@@ -101,9 +101,29 @@ template<class ParticleType>
 template<class TrackingData>
 void Foam::Cloud<ParticleType>::move(TrackingData& td)
 {
+    const polyBoundaryMesh& pbm = pMesh().boundaryMesh();
     const globalMeshData& pData = polyMesh_.globalData();
-    const labelList& processorPatches = pData.processorPatches();
+
-    const labelList& processorPatchIndices = pData.processorPatchIndices();
+    // Which patches are processor patches
+    const labelList& procPatches = pData.processorPatches();
+
+    // Indexing of patches into the procPatches list
+    const labelList& procPatchIndices = pData.processorPatchIndices();
+
+    // Indexing of equivalent patch on neighbour processor into the
+    // procPatches list on the neighbour
+    const labelList& procPatchNeighbours = pData.processorPatchNeighbours();
+
+    // Which processors this processor is connected to
+    const labelList& neighbourProcs = pData[Pstream::myProcNo()];
+
+    // Indexing from the processor number into the neighbourProcs list
+    labelList neighbourProcIndices(Pstream::nProcs(), -1);
+
+    forAll(neighbourProcs, i)
+    {
+        neighbourProcIndices[neighbourProcs[i]] = i;
+    }
 
     // Initialise the stepFraction moved for the particles
     forAllIter(typename Cloud<ParticleType>, *this, pIter)
@@ -114,9 +134,19 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
     // While there are particles to transfer
     while (true)
     {
-        // List of lists of particles to be transfered for all the processor
+        // List of lists of particles to be transfered for all of the
-        // patches
+        // neighbour processors
-        List<IDLList<ParticleType> > transferList(processorPatches.size());
+        List<IDLList<ParticleType> > particleTransferLists
+        (
+            neighbourProcs.size()
+        );
+
+        // List of destination processorPatches indices for all of the
+        // neighbour processors
+        List<DynamicList<label> > patchIndexTransferLists
+        (
+            neighbourProcs.size()
+        );
 
         // Loop over all particles
         forAllIter(typename Cloud<ParticleType>, *this, pIter)
@@ -134,15 +164,28 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
                 // boundary face
                 if (Pstream::parRun() && p.facei_ >= pMesh().nInternalFaces())
                 {
-                    label patchi = pMesh().boundaryMesh().whichPatch(p.facei_);
+                    label patchi = pbm.whichPatch(p.facei_);
-                    label n = processorPatchIndices[patchi];
 
                     // ... and the face is on a processor patch
                     // prepare it for transfer
-                    if (n != -1)
+                    if (procPatchIndices[patchi] != -1)
                     {
+                        label n = neighbourProcIndices
+                        [
+                            refCast<const processorPolyPatch>
+                            (
+                                pbm[patchi]
+                            ).neighbProcNo()
+                        ];
+
                         p.prepareForParallelTransfer(patchi, td);
-                        transferList[n].append(this->remove(&p));
+
+                        particleTransferLists[n].append(this->remove(&p));
+
+                        patchIndexTransferLists[n].append
+                        (
+                            procPatchNeighbours[patchi]
+                        );
                     }
                 }
             }
@@ -157,31 +200,30 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
             break;
         }
 
-
         // Allocate transfer buffers
         PstreamBuffers pBufs(Pstream::nonBlocking);
 
         // Stream into send buffers
-        forAll(transferList, i)
+        forAll(particleTransferLists, i)
         {
-            if (transferList[i].size())
+            if (particleTransferLists[i].size())
             {
                 UOPstream particleStream
                 (
-                    refCast<const processorPolyPatch>
+                    neighbourProcs[i],
-                    (
-                        pMesh().boundaryMesh()[processorPatches[i]]
-                    ).neighbProcNo(),
                     pBufs
                 );
 
-                particleStream << transferList[i];
+                particleStream
+                    << patchIndexTransferLists[i]
+                    << particleTransferLists[i];
             }
         }
 
         // Set up transfers when in non-blocking mode. Returns sizes (in bytes)
         // to be sent/received.
         labelListList allNTrans(Pstream::nProcs());
+
         pBufs.finishedSends(allNTrans);
 
         bool transfered = false;
@@ -203,34 +245,35 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
             break;
         }
 
-
         // Retrieve from receive buffers
-        forAll(processorPatches, i)
+        forAll(neighbourProcs, i)
         {
-            label patchi = processorPatches[i];
+            label neighbProci = neighbourProcs[i];
 
-            const processorPolyPatch& procPatch =
+            label nRec = allNTrans[neighbProci][Pstream::myProcNo()];
-                refCast<const processorPolyPatch>
-                (pMesh().boundaryMesh()[patchi]);
 
-            label neighbProci = procPatch.neighbProcNo();
+            if (nRec)
-
-            label nRecPs = allNTrans[neighbProci][Pstream::myProcNo()];
-
-            if (nRecPs)
             {
                 UIPstream particleStream(neighbProci, pBufs);
 
+                labelList receivePatchIndex(particleStream);
+
                 IDLList<ParticleType> newParticles
                 (
                     particleStream,
                     typename ParticleType::iNew(*this)
                 );
 
+                label pI = 0;
+
                 forAllIter(typename Cloud<ParticleType>, newParticles, newpIter)
                 {
                     ParticleType& newp = newpIter();
+
+                    label patchi = procPatches[receivePatchIndex[pI++]];
+
                     newp.correctAfterParallelTransfer(patchi, td);
+
                     addParticle(newParticles.remove(&newp));
                 }
             }

src/lagrangian/basic/Cloud/Cloud.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "cloud.H"
 #include "IDLList.H"
 #include "IOField.H"
+#include "IOFieldField.H"
 #include "polyMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -267,6 +268,14 @@ public:
                 const IOField<DataType>& data
             ) const;
 
+            //- Check lagrangian data fieldfield
+            template<class DataType>
+            void checkFieldFieldIOobject
+            (
+                const Cloud<ParticleType>& c,
+                const IOFieldField<Field<DataType>, DataType>& data
+            ) const;
+
             //- Read the field data for the cloud of particles. Dummy at
             //  this level.
             virtual void readFields();

src/lagrangian/basic/Cloud/CloudIO.C

@@ -205,6 +205,31 @@ void Foam::Cloud<ParticleType>::checkFieldIOobject
 }
 
 
+template<class ParticleType>
+template<class DataType>
+void Foam::Cloud<ParticleType>::checkFieldFieldIOobject
+(
+    const Cloud<ParticleType>& c,
+    const IOFieldField<Field<DataType>, DataType>& data
+) const
+{
+    if (data.size() != c.size())
+    {
+        FatalErrorIn
+        (
+            "void Cloud<ParticleType>::checkFieldFieldIOobject"
+            "("
+                "const Cloud<ParticleType>&, "
+                "const IOFieldField<Field<DataType>, DataType>&"
+            ") const"
+        )   << "Size of " << data.name()
+            << " field " << data.size()
+            << " does not match the number of particles " << c.size()
+            << abort(FatalError);
+    }
+}
+
+
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::readFields()
 {}

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/CollisionRecordList.C

@@ -51,6 +51,104 @@ Foam::CollisionRecordList<PairType, WallType>::CollisionRecordList(Istream& is)
 }
 
 
+template<class PairType, class WallType>
+Foam::CollisionRecordList<PairType, WallType>::CollisionRecordList
+(
+    const labelField& pairAccessed,
+    const labelField& pairOrigProcOfOther,
+    const labelField& pairOrigIdOfOther,
+    const Field<PairType>& pairData,
+    const labelField& wallAccessed,
+    const vectorField& wallPRel,
+    const Field<WallType>& wallData
+)
+:
+    pairRecords_(),
+    wallRecords_()
+{
+    label nPair = pairAccessed.size();
+
+    if
+    (
+        pairOrigProcOfOther.size() != nPair
+     || pairOrigIdOfOther.size() != nPair
+     || pairData.size() != nPair
+    )
+    {
+        FatalErrorIn
+        (
+            "Foam::CollisionRecordList<PairType, WallType>::CollisionRecordList"
+            "("
+                "const labelField& pairAccessed,"
+                "const labelField& pairOrigProcOfOther,"
+                "const labelField& pairOrigIdOfOther,"
+                "const Field<PairType>& pairData,"
+                "const labelField& wallAccessed,"
+                "const vectorField& wallPRel,"
+                "const Field<WallType>& wallData"
+            ")"
+        )
+            << "Pair field size mismatch." << nl
+            << pairAccessed << nl
+            << pairOrigProcOfOther << nl
+            << pairOrigIdOfOther << nl
+            << pairData << nl
+            << abort(FatalError);
+    }
+
+    forAll(pairAccessed, i)
+    {
+        pairRecords_.append
+        (
+            PairCollisionRecord<PairType>
+            (
+                pairAccessed[i],
+                pairOrigProcOfOther[i],
+                pairOrigIdOfOther[i],
+                pairData[i]
+            )
+        );
+    }
+
+    label nWall = wallAccessed.size();
+
+    if (wallPRel.size() != nWall || wallData.size() != nWall)
+    {
+        FatalErrorIn
+        (
+            "Foam::CollisionRecordList<PairType, WallType>::CollisionRecordList"
+            "("
+                "const labelField& pairAccessed,"
+                "const labelField& pairOrigProcOfOther,"
+                "const labelField& pairOrigIdOfOther,"
+                "const Field<PairType>& pairData,"
+                "const labelField& wallAccessed,"
+                "const vectorField& wallPRel,"
+                "const Field<WallType>& wallData"
+            ")"
+        )
+            << "Wall field size mismatch." << nl
+            << wallAccessed << nl
+            << wallPRel << nl
+            << wallData << nl
+            << abort(FatalError);
+    }
+
+    forAll(wallAccessed, i)
+    {
+        wallRecords_.append
+        (
+            WallCollisionRecord<WallType>
+            (
+                wallAccessed[i],
+                wallPRel[i],
+                wallData[i]
+            )
+        );
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * /
 
 template<class PairType, class WallType>
@@ -60,6 +158,111 @@ Foam::CollisionRecordList<PairType, WallType>::~CollisionRecordList()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class PairType, class WallType>
+Foam::labelField
+Foam::CollisionRecordList<PairType, WallType>::pairAccessed() const
+{
+    labelField f(pairRecords_.size());
+
+    forAll(pairRecords_, i)
+    {
+        f[i] = pairRecords_[i].accessed();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::labelField
+Foam::CollisionRecordList<PairType, WallType>::pairOrigProcOfOther() const
+{
+    labelField f(pairRecords_.size());
+
+    forAll(pairRecords_, i)
+    {
+        f[i] = pairRecords_[i].origProcOfOther();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::labelField
+Foam::CollisionRecordList<PairType, WallType>::pairOrigIdOfOther() const
+{
+    labelField f(pairRecords_.size());
+
+    forAll(pairRecords_, i)
+    {
+        f[i] = pairRecords_[i].origIdOfOther();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::Field<PairType>
+Foam::CollisionRecordList<PairType, WallType>::pairData() const
+{
+    Field<PairType> f(pairRecords_.size());
+
+    forAll(pairRecords_, i)
+    {
+        f[i] = pairRecords_[i].collisionData();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::labelField
+Foam::CollisionRecordList<PairType, WallType>::wallAccessed() const
+{
+    labelField f(wallRecords_.size());
+
+    forAll(wallRecords_, i)
+    {
+        f[i] = wallRecords_[i].accessed();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::vectorField
+Foam::CollisionRecordList<PairType, WallType>::wallPRel() const
+{
+    vectorField f(wallRecords_.size());
+
+    forAll(wallRecords_, i)
+    {
+        f[i] = wallRecords_[i].pRel();
+    }
+
+    return f;
+}
+
+
+template<class PairType, class WallType>
+Foam::Field<WallType>
+Foam::CollisionRecordList<PairType, WallType>::wallData() const
+{
+    Field<WallType> f(wallRecords_.size());
+
+    forAll(wallRecords_, i)
+    {
+        f[i] = wallRecords_[i].collisionData();
+    }
+
+    return f;
+}
+
+
 template<class PairType, class WallType>
 Foam::PairCollisionRecord<PairType>&
 Foam::CollisionRecordList<PairType, WallType>::matchPairRecord
@@ -85,12 +288,12 @@ Foam::CollisionRecordList<PairType, WallType>::matchPairRecord
     }
 
     // Record not found, create a new one and return it as the last
-    // member of the list.  The status of the record will be accessed
+    // member of the list.  Setting the status of the record to be accessed
-    // by construction.
+    // on construction.
 
     pairRecords_.append
     (
-        PairCollisionRecord<PairType>(origProcOfOther, origIdOfOther)
+        PairCollisionRecord<PairType>(true, origProcOfOther, origIdOfOther)
     );
 
     return pairRecords_.last();
@@ -121,10 +324,10 @@ Foam::CollisionRecordList<PairType, WallType>::matchWallRecord
     }
 
     // Record not found, create a new one and return it as the last
-    // member of the list.  The status of the record will be accessed
+    // member of the list.  Setting the status of the record to be accessed
-    // by construction.
+    // on construction.
 
-    wallRecords_.append(WallCollisionRecord<WallType>(pRel));
+    wallRecords_.append(WallCollisionRecord<WallType>(true, pRel));
 
     return wallRecords_.last();
 }

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/CollisionRecordList.H

@@ -95,6 +95,17 @@ public:
         //- Construct from Istream
         CollisionRecordList(Istream&);
 
+        //- Construct from component fields (for IO)
+        CollisionRecordList
+        (
+            const labelField& pairAccessed,
+            const labelField& pairOrigProcOfOther,
+            const labelField& pairOrigIdOfOther,
+            const Field<PairType>& pairData,
+            const labelField& wallAccessed,
+            const vectorField& wallPRel,
+            const Field<WallType>& wallData
+        );
 
     //- Destructor
     ~CollisionRecordList();
@@ -102,6 +113,50 @@ public:
 
     // Member Functions
 
+        //- Return the active pair collisions
+        inline const DynamicList<PairCollisionRecord<PairType> >&
+        pairRecords() const;
+
+        //- Return the active wall collisions
+        inline const DynamicList<WallCollisionRecord<WallType> >&
+        wallRecords() const;
+
+        // Fields representing the data from each record, i.e if the
+        // records 0-N containing each data members {a, b, c, d...}
+        // are organised:
+        //
+        // a0 b0 c0 d0 ...
+        // a1 b1 c1 d1 ...
+        // a2 b2 c2 d2 ...
+        // ...
+        // aN bN cN dN ...
+        //
+        // Then these field return, for example, (c0, c1, c2,... cN)
+
+        //- Return field of pair accessed from each record, used for
+        //  field IO
+        labelField pairAccessed() const;
+
+        //- Return field of pair origProcOfOther from each record,
+        //  used for field IO
+        labelField pairOrigProcOfOther() const;
+
+        //- Return field of pair origIdOfOther from each record, used
+        //  for field IO
+        labelField pairOrigIdOfOther() const;
+
+        //- Return field of pair data from each record, used for field IO
+        Field<PairType> pairData() const;
+
+        //- Return field of wall accessed from each record, used for field IO
+        labelField wallAccessed() const;
+
+        //- Return field of wall pRel from each record, used for field IO
+        vectorField wallPRel() const;
+
+        //- Return field of wall data from each record, used for field IO
+        Field<WallType> wallData() const;
+
         //- Enquires if the proc and id pair of the other particle are
         //  present in the records.  If so, return non-const access to
         //  the PairCollisionRecord (hence the data) and mark the
@@ -172,6 +227,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "CollisionRecordListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #ifdef NoRepository
 #   include "CollisionRecordList.C"
 #endif

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/CollisionRecordListI.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class PairType, class WallType>
+const Foam::DynamicList<Foam::PairCollisionRecord<PairType> >&
+Foam::CollisionRecordList<PairType, WallType>::pairRecords() const
+{
+    return pairRecords_;
+}
+
+
+template<class PairType, class WallType>
+const Foam::DynamicList<Foam::WallCollisionRecord<WallType> >&
+Foam::CollisionRecordList<PairType, WallType>::wallRecords() const
+{
+    return wallRecords_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/PairCollisionRecord/PairCollisionRecord.C

@@ -30,8 +30,8 @@ License
 template<class Type>
 Foam::PairCollisionRecord<Type>::PairCollisionRecord()
 :
-    origProcOfOther_(-VGREAT),
+    origProcOfOther_(0),
-    origIdOfOther_(-VGREAT),
+    origIdOfOther_(-1),
     data_(pTraits<Type>::zero)
 {}
 
@@ -39,6 +39,7 @@ Foam::PairCollisionRecord<Type>::PairCollisionRecord()
 template<class Type>
 Foam::PairCollisionRecord<Type>::PairCollisionRecord
 (
+    bool accessed,
     label origProcOfOther,
     label origIdOfOther,
     const Type& data
@@ -47,7 +48,14 @@ Foam::PairCollisionRecord<Type>::PairCollisionRecord
     origProcOfOther_(origProcOfOther + 1),
     origIdOfOther_(origIdOfOther),
     data_(data)
-{}
+{
+    // Default assignment to origProcOfOther_ assumes accessed is true
+
+    if (!accessed)
+    {
+        setUnaccessed();
+    }
+}
 
 
 template<class Type>
@@ -56,7 +64,7 @@ Foam::PairCollisionRecord<Type>::PairCollisionRecord
     const PairCollisionRecord<Type>& pCR
 )
 :
-    origProcOfOther_(pCR.origProcOfOther() + 1),
+    origProcOfOther_(pCR.origProcOfOther_),
     origIdOfOther_(pCR.origIdOfOther_),
     data_(pCR.data_)
 {}

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/PairCollisionRecord/PairCollisionRecord.H

@@ -108,6 +108,7 @@ public:
         //- Construct from components
         PairCollisionRecord
         (
+            bool accessed,
             label origProcOfOther,
             label origIdOfOther,
             const Type& data = pTraits<Type>::zero

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/WallCollisionRecord/WallCollisionRecord.C

@@ -25,6 +25,12 @@ License
 
 #include "WallCollisionRecord.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class Type>
+const Foam::scalar Foam::WallCollisionRecord<Type>::errorCosAngle(1.0 + 1e-6);
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -39,11 +45,12 @@ Foam::WallCollisionRecord<Type>::WallCollisionRecord()
 template<class Type>
 Foam::WallCollisionRecord<Type>::WallCollisionRecord
 (
+    bool accessed,
     const vector& pRel,
     const Type& data
 )
 :
-    accessed_(true),
+    accessed_(accessed),
     pRel_(pRel),
     data_(data)
 {}

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/WallCollisionRecord/WallCollisionRecord.H

@@ -82,7 +82,7 @@ class WallCollisionRecord
         // //- Recording whether or not this record has been accessed
         bool accessed_;
 
-        //- The position of wall impact relative to the cell centre
+        //- The position of wall impact relative to the particle centre
         vector pRel_;
 
         //- Collision data, stored as if the storing particle was the
@@ -92,6 +92,12 @@ class WallCollisionRecord
 
 public:
 
+    // Static data members
+
+        //- Tolerance for detecting seriously erroneous wall matches
+        static const scalar errorCosAngle;
+
+
     // Constructors
 
         //- Construct null
@@ -100,6 +106,7 @@ public:
         //- Construct from components
         WallCollisionRecord
         (
+            bool accessed,
             const vector& pRel,
             const Type& data = pTraits<Type>::zero
         );

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/CollisionRecordList/WallCollisionRecord/WallCollisionRecordI.H

@@ -39,7 +39,7 @@ inline bool Foam::WallCollisionRecord<Type>::match
     // Using the new data as the acceptance criterion
     scalar cosAcceptanceAngle = magpRel/radius;
 
-    if (cosAcceptanceAngle > 1.0)
+    if (cosAcceptanceAngle > errorCosAngle)
     {
         Info<< "pRel_ " << pRel_ << " " << magpRel_ << nl
             << "pRel " << pRel << " " << magpRel << nl
@@ -75,6 +75,14 @@ inline bool Foam::WallCollisionRecord<Type>::match
 }
 
 
+template<class Type>
+inline const Foam::vector&
+Foam::WallCollisionRecord<Type>::pRel() const
+{
+    return pRel_;
+}
+
+
 template<class Type>
 inline const Type&
 Foam::WallCollisionRecord<Type>::collisionData() const

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -33,7 +33,7 @@ Description
     - drag
     - turbulent dispersion
     - wall interactions
-    - many-body collisions
+    - many-body collisions, including memory of data from previous collision
 
 SourceFiles
     KinematicParcelI.H
@@ -53,12 +53,18 @@ SourceFiles
 
 #include "KinematicCloud.H"
 #include "CollisionRecordList.H"
+#include "labelFieldIOField.H"
+#include "vectorFieldIOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef CollisionRecordList<vector, vector> collisionRecordList;
+typedef vectorIOFieldField pairDataIOFieldField;
+typedef vectorIOFieldField wallDataIOFieldField;
+
 template<class ParcelType>
 class KinematicParcel;
 
@@ -268,7 +274,7 @@ protected:
             vector UTurb_;
 
             //- Particle collision records
-            CollisionRecordList<vector, vector> collisionRecords_;
+            collisionRecordList collisionRecords_;
 
 
         // Cell-based quantities
@@ -400,8 +406,7 @@ public:
             inline const vector& UTurb() const;
 
             //- Return const access to the collision records
-            inline const CollisionRecordList<vector, vector>&
+            inline const collisionRecordList& collisionRecords() const;
-            collisionRecords() const;
 
 
         // Edit
@@ -440,7 +445,7 @@ public:
             inline vector& UTurb();
 
             //- Return access to collision records
-            inline CollisionRecordList<vector, vector>& collisionRecords();
+            inline collisionRecordList& collisionRecords();
 
 
         // Helper functions

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -348,7 +348,7 @@ inline bool& Foam::KinematicParcel<ParcelType>::active()
 
 
 template <class ParcelType>
-inline const Foam::CollisionRecordList<Foam::vector, Foam::vector>&
+inline const Foam::collisionRecordList&
 Foam::KinematicParcel<ParcelType>::collisionRecords() const
 {
     return collisionRecords_;
@@ -426,7 +426,7 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
 
 
 template <class ParcelType>
-inline Foam::CollisionRecordList<Foam::vector, Foam::vector>&
+inline Foam::collisionRecordList&
 Foam::KinematicParcel<ParcelType>::collisionRecords()
 {
     return collisionRecords_;

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C

@@ -43,7 +43,14 @@ Foam::string Foam::KinematicParcel<ParcelType>::propHeader =
   + " (torquex torquey torquez)"
   + " rho"
   + " tTurb"
-  + " (UTurbx UTurby UTurbz)";
+  + " (UTurbx UTurby UTurbz)"
+  + " collisionRecordsPairAccessed"
+  + " collisionRecordsPairOrigProcOfOther"
+  + " collisionRecordsPairOrigIdOfOther"
+  + " (collisionRecordsPairData)"
+  + " collisionRecordsWallAccessed"
+  + " collisionRecordsWallPRel"
+  + " (collisionRecordsWallData)";
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -167,7 +174,58 @@ void Foam::KinematicParcel<ParcelType>::readFields(Cloud<ParcelType>& c)
     IOField<vector> UTurb(c.fieldIOobject("UTurb", IOobject::MUST_READ));
     c.checkFieldIOobject(c, UTurb);
 
+    labelIOFieldField collisionRecordsPairAccessed
+    (
+        c.fieldIOobject("collisionRecordsPairAccessed", IOobject::MUST_READ)
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsPairAccessed);
+
+    labelIOFieldField collisionRecordsPairOrigProcOfOther
+    (
+        c.fieldIOobject
+        (
+            "collisionRecordsPairOrigProcOfOther",
+            IOobject::MUST_READ
+        )
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsPairOrigProcOfOther);
+
+    labelIOFieldField collisionRecordsPairOrigIdOfOther
+    (
+        c.fieldIOobject
+        (
+            "collisionRecordsPairOrigIdOfOther",
+            IOobject::MUST_READ
+        )
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsPairOrigProcOfOther);
+
+    pairDataIOFieldField collisionRecordsPairData
+    (
+        c.fieldIOobject("collisionRecordsPairData", IOobject::MUST_READ)
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsPairData);
+
+    labelIOFieldField collisionRecordsWallAccessed
+    (
+        c.fieldIOobject("collisionRecordsWallAccessed", IOobject::MUST_READ)
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsWallAccessed);
+
+    vectorIOFieldField collisionRecordsWallPRel
+    (
+        c.fieldIOobject("collisionRecordsWallPRel", IOobject::MUST_READ)
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsWallPRel);
+
+    wallDataIOFieldField collisionRecordsWallData
+    (
+        c.fieldIOobject("collisionRecordsWallData", IOobject::MUST_READ)
+    );
+    c.checkFieldFieldIOobject(c, collisionRecordsWallData);
+
     label i = 0;
+
     forAllIter(typename Cloud<ParcelType>, c, iter)
     {
         ParcelType& p = iter();
@@ -182,6 +240,17 @@ void Foam::KinematicParcel<ParcelType>::readFields(Cloud<ParcelType>& c)
         p.rho_ = rho[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
+        p.collisionRecords_ = collisionRecordList
+        (
+            collisionRecordsPairAccessed[i],
+            collisionRecordsPairOrigProcOfOther[i],
+            collisionRecordsPairOrigIdOfOther[i],
+            collisionRecordsPairData[i],
+            collisionRecordsWallAccessed[i],
+            collisionRecordsWallPRel[i],
+            collisionRecordsWallData[i]
+        );
+
         i++;
     }
 }
@@ -206,14 +275,56 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
     IOField<vector> f(c.fieldIOobject("f", IOobject::NO_READ), np);
     IOField<vector> angularMomentum
     (
-        c.fieldIOobject("angularMomentum", IOobject::NO_READ), np
+        c.fieldIOobject("angularMomentum", IOobject::NO_READ),
+        np
     );
     IOField<vector> torque(c.fieldIOobject("torque", IOobject::NO_READ), np);
     IOField<scalar> rho(c.fieldIOobject("rho", IOobject::NO_READ), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb", IOobject::NO_READ), np);
     IOField<vector> UTurb(c.fieldIOobject("UTurb", IOobject::NO_READ), np);
 
+    labelIOFieldField collisionRecordsPairAccessed
+    (
+        c.fieldIOobject("collisionRecordsPairAccessed", IOobject::NO_READ),
+        np
+    );
+    labelIOFieldField collisionRecordsPairOrigProcOfOther
+    (
+        c.fieldIOobject
+        (
+            "collisionRecordsPairOrigProcOfOther",
+            IOobject::NO_READ
+        ),
+        np
+    );
+    labelIOFieldField collisionRecordsPairOrigIdOfOther
+    (
+        c.fieldIOobject("collisionRecordsPairOrigIdOfOther", IOobject::NO_READ),
+        np
+    );
+    pairDataIOFieldField collisionRecordsPairData
+    (
+        c.fieldIOobject("collisionRecordsPairData", IOobject::NO_READ),
+        np
+    );
+    labelIOFieldField collisionRecordsWallAccessed
+    (
+        c.fieldIOobject("collisionRecordsWallAccessed", IOobject::NO_READ),
+        np
+    );
+    vectorIOFieldField collisionRecordsWallPRel
+    (
+        c.fieldIOobject("collisionRecordsWallPRel", IOobject::NO_READ),
+        np
+    );
+    wallDataIOFieldField collisionRecordsWallData
+    (
+        c.fieldIOobject("collisionRecordsWallData", IOobject::NO_READ),
+        np
+    );
+
     label i = 0;
+
     forAllConstIter(typename Cloud<ParcelType>, c, iter)
     {
         const KinematicParcel<ParcelType>& p = iter();
@@ -229,6 +340,16 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
         rho[i] = p.rho();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
+        collisionRecordsPairAccessed[i] = p.collisionRecords().pairAccessed();
+        collisionRecordsPairOrigProcOfOther[i] =
+            p.collisionRecords().pairOrigProcOfOther();
+        collisionRecordsPairOrigIdOfOther[i] =
+            p.collisionRecords().pairOrigIdOfOther();
+        collisionRecordsPairData[i] = p.collisionRecords().pairData();
+        collisionRecordsWallAccessed[i] = p.collisionRecords().wallAccessed();
+        collisionRecordsWallPRel[i] = p.collisionRecords().wallPRel();
+        collisionRecordsWallData[i] = p.collisionRecords().wallData();
+
         i++;
     }
 
@@ -243,6 +364,13 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
     rho.write();
     tTurb.write();
     UTurb.write();
+    collisionRecordsPairAccessed.write();
+    collisionRecordsPairOrigProcOfOther.write();
+    collisionRecordsPairOrigIdOfOther.write();
+    collisionRecordsPairData.write();
+    collisionRecordsWallAccessed.write();
+    collisionRecordsWallPRel.write();
+    collisionRecordsWallData.write();
 }
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/defineBasicKinematicParcel.C

@@ -25,6 +25,7 @@ License
 
 #include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
+#include "KinematicParcel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/intermediate/parcels/include/makeParcelCollisionModels.H

@@ -36,6 +36,7 @@ License
 #include "PairSpringSliderDashpot.H"
 
 #include "WallSpringSliderDashpot.H"
+#include "WallLocalSpringSliderDashpot.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,6 +74,13 @@ License
         WallSpringSliderDashpot,                                              \
         KinematicCloud,                                                       \
         ParcelType                                                            \
+    );                                                                        \
+                                                                              \
+    makeWallModelType                                                         \
+    (                                                                         \
+        WallLocalSpringSliderDashpot,                                         \
+        KinematicCloud,                                                       \
+        ParcelType                                                            \
     );
 
 

src/lagrangian/intermediate/parcels/include/makeReactingParcelCollisionModels.H

@@ -37,6 +37,7 @@ License
 #include "PairSpringSliderDashpot.H"
 
 #include "WallSpringSliderDashpot.H"
+#include "WallLocalSpringSliderDashpot.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,9 +100,16 @@ License
         KinematicCloud,                                                       \
         ParcelType,                                                           \
         ThermoType                                                            \
+    );                                                                        \
+                                                                              \
+    makeWallModelThermoType                                                   \
+    (                                                                         \
+        WallLocalSpringSliderDashpot,                                         \
+        KinematicCloud,                                                       \
+        ParcelType,                                                           \
+        ThermoType                                                            \
     );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

@@ -45,7 +45,12 @@ void Foam::PairSpringSliderDashpot<CloudType>::findMinMaxProperties
 
         // Finding minimum diameter to avoid excessive arithmetic
 
-        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+        scalar dEff = p.d();
+
+        if (useEquivalentSize_)
+        {
+            dEff *= cbrt(p.nParticle()*volumeFactor_);
+        }
 
         RMin = min(dEff, RMin);
 
@@ -94,8 +99,14 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
             this->coeffDict().lookup("collisionResolutionSteps")
         )
     ),
-    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
+    volumeFactor_(1.0),
+    useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
 {
+    if (useEquivalentSize_)
+    {
+        volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
+    }
+
     scalar nu = this->owner().constProps().poissonsRatio();
 
     scalar E = this->owner().constProps().youngsModulus();
@@ -158,9 +169,19 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 {
     vector r_AB = (pA.position() - pB.position());
 
-    scalar dAEff = pA.d()*cbrt(pA.nParticle()*volumeFactor_);
+    scalar dAEff = pA.d();
 
-    scalar dBEff = pB.d()*cbrt(pB.nParticle()*volumeFactor_);
+    if (useEquivalentSize_)
+    {
+        dAEff *= cbrt(pA.nParticle()*volumeFactor_);
+    }
+
+    scalar dBEff = pB.d();
+
+    if (useEquivalentSize_)
+    {
+        dBEff *= cbrt(pB.nParticle()*volumeFactor_);
+    }
 
     scalar normalOverlapMag = 0.5*(dAEff + dBEff) - mag(r_AB);
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H

@@ -58,9 +58,6 @@ class PairSpringSliderDashpot
         //  the same Poisson's ratio and Young's modulus
         scalar Gstar_;
 
-        //- Poisson's ratio of both particles
-        scalar sigma_;
-
         //- alpha-coefficient, related to coefficient of restitution
         scalar alpha_;
 
@@ -92,6 +89,10 @@ class PairSpringSliderDashpot
         // factor
         scalar volumeFactor_;
 
+        //- Switch to control use of equivalent size particles.  Used
+        //  because the calculation can be very expensive.
+        bool useEquivalentSize_;
+
 
     // Private Member Functions
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C

@@ -0,0 +1,354 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+     \\/     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 "WallLocalSpringSliderDashpot.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::WallLocalSpringSliderDashpot<CloudType>::findMinMaxProperties
+(
+    scalar& rMin,
+    scalar& rhoMax,
+    scalar& UMagMax
+) const
+{
+    rMin = VGREAT;
+    rhoMax = -VGREAT;
+    UMagMax = -VGREAT;
+
+    forAllConstIter(typename CloudType, this->owner(), iter)
+    {
+        const typename CloudType::parcelType& p = iter();
+
+        // Finding minimum diameter to avoid excessive arithmetic
+
+        scalar dEff = p.d();
+
+        if (useEquivalentSize_)
+        {
+            dEff *= cbrt(p.nParticle()*volumeFactor_);
+        }
+
+        rMin = min(dEff, rMin);
+
+        rhoMax = max(p.rho(), rhoMax);
+
+        UMagMax = max
+        (
+            mag(p.U()) + mag(p.omega())*dEff/2,
+            UMagMax
+        );
+    }
+
+    // Transform the minimum diameter into minimum radius
+    //     rMin = dMin/2
+
+    rMin /= 2.0;
+}
+
+
+template <class CloudType>
+void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
+(
+    typename CloudType::parcelType& p,
+    const point& site,
+    const WallSiteData<vector>& data,
+    scalar pREff
+) const
+{
+    // wall patch index
+    label wPI = patchMap_[data.patchIndex()];
+
+    // data for this patch
+    scalar Estar = Estar_[wPI];
+    scalar Gstar = Gstar_[wPI];
+    scalar alpha = alpha_[wPI];
+    scalar b = b_[wPI];
+    scalar mu = mu_[wPI];
+
+    vector r_PW = p.position() - site;
+
+    vector U_PW = p.U() - data.wallData();
+
+    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
+
+    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
+
+    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
+
+    scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
+
+    vector fN_PW =
+        rHat_PW
+       *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
+
+    p.f() += fN_PW;
+
+    vector USlip_PW =
+        U_PW - (U_PW & rHat_PW)*rHat_PW
+      + (p.omega() ^ (pREff*-rHat_PW));
+
+    scalar deltaT = this->owner().mesh().time().deltaTValue();
+
+    vector& tangentialOverlap_PW =
+        p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
+
+    tangentialOverlap_PW += USlip_PW*deltaT;
+
+    scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
+
+    if (tangentialOverlapMag > VSMALL)
+    {
+        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
+
+        scalar etaT = etaN;
+
+        // Tangential force
+        vector fT_PW;
+
+        if (kT*tangentialOverlapMag > mu*mag(fN_PW))
+        {
+            // Tangential force greater than sliding friction,
+            // particle slips
+
+            fT_PW = -mu*mag(fN_PW)*USlip_PW/mag(USlip_PW);
+
+            tangentialOverlap_PW = vector::zero;
+        }
+        else
+        {
+            fT_PW =
+                -kT*tangentialOverlapMag
+               *tangentialOverlap_PW/tangentialOverlapMag
+              - etaT*USlip_PW;
+        }
+
+        p.f() += fT_PW;
+
+        p.torque() += (pREff*-rHat_PW) ^ fT_PW;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallModel<CloudType>(dict, cloud, typeName),
+    Estar_(),
+    Gstar_(),
+    alpha_(),
+    b_(),
+    mu_(),
+    patchMap_(),
+    maxEstarIndex_(-1),
+    collisionResolutionSteps_
+    (
+        readScalar
+        (
+            this->coeffDict().lookup("collisionResolutionSteps")
+        )
+    ),
+    volumeFactor_(1.0),
+    useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
+{
+    if (useEquivalentSize_)
+    {
+        volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
+    }
+
+    scalar pNu = this->owner().constProps().poissonsRatio();
+
+    scalar pE = this->owner().constProps().youngsModulus();
+
+    const polyMesh& mesh = cloud.mesh();
+
+    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+
+    patchMap_.setSize(bMesh.size(), -1);
+
+    DynamicList<label> wallPatchIndices;
+
+    forAll(bMesh, patchI)
+    {
+        if (isA<wallPolyPatch>(bMesh[patchI]))
+        {
+            wallPatchIndices.append(bMesh[patchI].index());
+        }
+    }
+
+    label nWallPatches = wallPatchIndices.size();
+
+    Estar_.setSize(nWallPatches);
+    Gstar_.setSize(nWallPatches);
+    alpha_.setSize(nWallPatches);
+    b_.setSize(nWallPatches);
+    mu_.setSize(nWallPatches);
+
+    scalar maxEstar = -GREAT;
+
+    forAll(wallPatchIndices, wPI)
+    {
+        const dictionary& patchCoeffDict
+        (
+            this->coeffDict().subDict(bMesh[wallPatchIndices[wPI]].name())
+        );
+
+        patchMap_[wallPatchIndices[wPI]] = wPI;
+
+        scalar nu = dimensionedScalar
+        (
+            patchCoeffDict.lookup("poissonsRatio")
+        ).value();
+
+        scalar E = dimensionedScalar
+        (
+            patchCoeffDict.lookup("youngsModulus")
+        ).value();
+
+        Estar_[wPI] = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
+
+        Gstar_[wPI] = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
+
+        alpha_[wPI] =
+            dimensionedScalar(patchCoeffDict.lookup("alpha")).value();
+
+        b_[wPI] = dimensionedScalar(patchCoeffDict.lookup("b")).value();
+
+        mu_[wPI] = dimensionedScalar(patchCoeffDict.lookup("mu")).value();
+
+        if (Estar_[wPI] > maxEstar)
+        {
+            maxEstarIndex_ = wPI;
+
+            maxEstar = Estar_[wPI];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::WallLocalSpringSliderDashpot<CloudType>::~WallLocalSpringSliderDashpot()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::WallLocalSpringSliderDashpot<CloudType>::pREff
+(
+    const typename CloudType::parcelType& p
+) const
+{
+    if (useEquivalentSize_)
+    {
+        return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
+    }
+    else
+    {
+        return p.d()/2;
+    }
+}
+
+
+template<class CloudType>
+bool Foam::WallLocalSpringSliderDashpot<CloudType>::controlsTimestep() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+Foam::label Foam::WallLocalSpringSliderDashpot<CloudType>::nSubCycles() const
+{
+    if (!(this->owner().size()))
+    {
+        return 1;
+    }
+
+    scalar rMin;
+    scalar rhoMax;
+    scalar UMagMax;
+
+    findMinMaxProperties(rMin, rhoMax, UMagMax);
+
+    // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
+    scalar minCollisionDeltaT =
+        5.429675
+       *rMin
+       *pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + VSMALL), 0.4)
+       /collisionResolutionSteps_;
+
+    return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
+}
+
+
+template<class CloudType>
+void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
+(
+    typename CloudType::parcelType& p,
+    const List<point>& flatSitePoints,
+    const List<WallSiteData<vector> >& flatSiteData,
+    const List<point>& sharpSitePoints,
+    const List<WallSiteData<vector> >& sharpSiteData
+) const
+{
+    scalar pREff = this->pREff(p);
+
+    forAll(flatSitePoints, siteI)
+    {
+        evaluateWall
+        (
+            p,
+            flatSitePoints[siteI],
+            flatSiteData[siteI],
+            pREff
+        );
+    }
+
+    forAll(sharpSitePoints, siteI)
+    {
+        // Treating sharp sites like flat sites
+
+        evaluateWall
+        (
+            p,
+            sharpSitePoints[siteI],
+            sharpSiteData[siteI],
+            pREff
+        );
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+     \\/     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::WallLocalSpringSliderDashpot
+
+Description
+    Forces between particles and walls, interacting with a spring,
+    slider, damper model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallLocalSpringSliderDashpot_H
+#define WallLocalSpringSliderDashpot_H
+
+#include "WallModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                Class WallLocalSpringSliderDashpot Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class WallLocalSpringSliderDashpot
+:
+    public WallModel<CloudType>
+{
+    // Private data
+
+        //- Effective Young's modulus value
+        scalarList Estar_;
+
+        //- Effective shear modulus value
+        scalarList Gstar_;
+
+        //- alpha-coefficient, related to coefficient of restitution
+        scalarList alpha_;
+
+        //- Spring power (b = 1 for linear, b = 3/2 for Hertzian)
+        scalarList b_;
+
+        //- Coefficient of friction in for tangential sliding
+        scalarList mu_;
+
+        //- Mapping the patch index to the model data
+        labelList patchMap_;
+
+        //- Index of the maximum value of Estar_
+        label maxEstarIndex_;
+
+        //- The number of steps over which to resolve the minimum
+        //  harmonic approximation of the collision period
+        scalar collisionResolutionSteps_;
+
+        //- Volume factor for determining the equivalent size of a
+        //  parcel where nParticles is not 1.  The equivalent size of
+        //  the parcel is
+        //      parcelEquivVolume = volumeFactor*nParticles*p.volume()
+        //  so
+        //      parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
+        //  + When volumeFactor = 1, the particles are compressed
+        //    together so that the equivalent volume of the parcel is
+        //    the sum of the constituent particles
+        //  + When volumeFactor = 3*sqrt(2)/pi, the particles are
+        //    close packed, but uncompressed.
+        //  + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
+        //    grouped.
+        // 3*sqrt(2)/pi = 1.350474 is the volume factor for close
+        // packing, i.e pi/(3*sqrt(2)) is the maximum close packing
+        // factor
+        scalar volumeFactor_;
+
+        //- Switch to control use of equivalent size particles.  Used
+        //  because the calculation can be very expensive.
+        bool useEquivalentSize_;
+
+
+    // Private Member Functions
+
+        //- Find the appropriate properties for determining the minimum
+        //- allowable timestep
+        void findMinMaxProperties
+        (
+            scalar& rMin,
+            scalar& rhoMax,
+            scalar& vMagMax
+        ) const;
+
+        //- Calculate the wall interaction for a parcel at a given site
+        void evaluateWall
+        (
+            typename CloudType::parcelType& p,
+            const point& site,
+            const WallSiteData<vector>& data,
+            scalar pREff
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("WallLocalSpringSliderDashpot");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        WallLocalSpringSliderDashpot(const dictionary& dict, CloudType& cloud);
+
+
+    //- Destructor
+    virtual ~WallLocalSpringSliderDashpot();
+
+
+    // Member Functions
+
+        //- Return the volumeFactor
+        inline scalar volumeFactor() const
+        {
+            return volumeFactor_;
+        }
+
+        //- Return the effective radius for a particle for the model
+        virtual scalar pREff(const typename CloudType::parcelType& p) const;
+
+        //- Whether the WallModel has a timestep limit that will
+        //  require subCycling
+        virtual bool controlsTimestep() const;
+
+        //- For WallModels that control the timestep, calculate the
+        //  number of subCycles needed to satisfy the minimum
+        //  allowable timestep
+        virtual label nSubCycles() const;
+
+        //- Calculate the wall interaction for a parcel
+        virtual void evaluateWall
+        (
+            typename CloudType::parcelType& p,
+            const List<point>& flatSitePoints,
+            const List<WallSiteData<vector> >& flatSiteData,
+            const List<point>& sharpSitePoints,
+            const List<WallSiteData<vector> >& sharpSiteData
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallLocalSpringSliderDashpot.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.C

@@ -58,6 +58,14 @@ Foam::WallModel<CloudType>::owner() const
 }
 
 
+template<class CloudType>
+CloudType&
+Foam::WallModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
 template<class CloudType>
 const Foam::dictionary& Foam::WallModel<CloudType>::dict() const
 {

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallModel/WallModel.H

@@ -111,6 +111,9 @@ public:
         //- Return the owner cloud object
         const CloudType& owner() const;
 
+        //- Return non-const access to the owner cloud object
+        CloudType& owner();
+
         //- Return the dictionary
         const dictionary& dict() const;
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C

@@ -45,7 +45,12 @@ void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
 
         // Finding minimum diameter to avoid excessive arithmetic
 
-        scalar dEff = p.d()*cbrt(p.nParticle()*volumeFactor_);
+        scalar dEff = p.d();
+
+        if (useEquivalentSize_)
+        {
+            dEff *= cbrt(p.nParticle()*volumeFactor_);
+        }
 
         rMin = min(dEff, rMin);
 
@@ -64,25 +69,22 @@ void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
     rMin /= 2.0;
 }
 
+
 template <class CloudType>
 void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 (
     typename CloudType::parcelType& p,
     const point& site,
     const WallSiteData<vector>& data,
-    scalar pNu,
-    scalar pE,
     scalar pREff,
-    scalar Estar,
+    scalar kN
-    scalar kN,
-    scalar Gstar
 ) const
 {
     vector r_PW = p.position() - site;
 
     vector U_PW = p.U() - data.wallData();
 
-    scalar normalOverlapMag = pREff - mag(r_PW);
+    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
 
     vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
 
@@ -109,7 +111,7 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 
     if (tangentialOverlapMag > VSMALL)
     {
-        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
+        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar_;
 
         scalar etaT = etaN;
 
@@ -150,8 +152,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
 )
 :
     WallModel<CloudType>(dict, cloud, typeName),
-    E_(dimensionedScalar(this->coeffDict().lookup("youngsModulus")).value()),
+    Estar_(),
-    nu_(dimensionedScalar(this->coeffDict().lookup("poissonsRatio")).value()),
+    Gstar_(),
     alpha_(dimensionedScalar(this->coeffDict().lookup("alpha")).value()),
     b_(dimensionedScalar(this->coeffDict().lookup("b")).value()),
     mu_(dimensionedScalar(this->coeffDict().lookup("mu")).value()),
@@ -162,8 +164,32 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
             this->coeffDict().lookup("collisionResolutionSteps")
         )
     ),
-    volumeFactor_(this->dict().lookupOrDefault("volumeFactor", 1.0))
+    volumeFactor_(1.0),
-{}
+    useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
+{
+    if (useEquivalentSize_)
+    {
+        volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
+    }
+
+    scalar nu = dimensionedScalar
+    (
+        this->coeffDict().lookup("poissonsRatio")
+    ).value();
+
+    scalar E = dimensionedScalar
+    (
+        this->coeffDict().lookup("youngsModulus")
+    ).value();
+
+    scalar pNu = this->owner().constProps().poissonsRatio();
+
+    scalar pE = this->owner().constProps().youngsModulus();
+
+    Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
+
+    Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -181,9 +207,17 @@ Foam::scalar Foam::WallSpringSliderDashpot<CloudType>::pREff
     const typename CloudType::parcelType& p
 ) const
 {
+    if (useEquivalentSize_)
+    {
         return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
+    }
+    else
+    {
+        return p.d()/2;
+    }
 }
 
+
 template<class CloudType>
 bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
 {
@@ -205,17 +239,11 @@ Foam::label Foam::WallSpringSliderDashpot<CloudType>::nSubCycles() const
 
     findMinMaxProperties(rMin, rhoMax, UMagMax);
 
-    scalar pNu = this->owner().constProps().poissonsRatio();
-
-    scalar pE = this->owner().constProps().youngsModulus();
-
-    scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
-
     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
     scalar minCollisionDeltaT =
         5.429675
        *rMin
-       *pow(rhoMax/(Estar*sqrt(UMagMax) + VSMALL), 0.4)
+       *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
        /collisionResolutionSteps_;
 
     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
@@ -232,17 +260,9 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
     const List<WallSiteData<vector> >& sharpSiteData
 ) const
 {
-    scalar pNu = this->owner().constProps().poissonsRatio();
-
-    scalar pE = this->owner().constProps().youngsModulus();
-
     scalar pREff = this->pREff(p);
 
-    scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
+    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar_;
-
-    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
-
-    scalar GStar = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu_ - sqr(nu_))/E_));
 
     forAll(flatSitePoints, siteI)
     {
@@ -251,12 +271,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             p,
             flatSitePoints[siteI],
             flatSiteData[siteI],
-            pNu,
-            pE,
             pREff,
-            Estar,
+            kN
-            kN,
-            GStar
         );
     }
 
@@ -269,12 +285,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
             p,
             sharpSitePoints[siteI],
             sharpSiteData[siteI],
-            pNu,
-            pE,
             pREff,
-            Estar,
+            kN
-            kN,
-            GStar
         );
     }
 }

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H

@@ -50,11 +50,11 @@ class WallSpringSliderDashpot
 {
     // Private data
 
-        //- Young's modulus of the wall
+        //- Effective Young's modulus value
-        scalar E_;
+        scalar Estar_;
 
-        //- Poisson's ratio of the wall
+        //- Effective shear modulus value
-        scalar nu_;
+        scalar Gstar_;
 
         //- alpha-coefficient, related to coefficient of restitution
         scalar alpha_;
@@ -87,6 +87,10 @@ class WallSpringSliderDashpot
         // factor
         scalar volumeFactor_;
 
+        //- Switch to control use of equivalent size particles.  Used
+        //  because the calculation can be very expensive.
+        bool useEquivalentSize_;
+
 
     // Private Member Functions
 
@@ -105,12 +109,8 @@ class WallSpringSliderDashpot
             typename CloudType::parcelType& p,
             const point& site,
             const WallSiteData<vector>& data,
-            scalar pNu,
-            scalar pE,
             scalar pREff,
-            scalar Estar,
+            scalar kN
-            scalar kN,
-            scalar Gstar
         ) const;
 
 

src/parallel/reconstruct/reconstruct/reconstructLagrangian.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "cloud.H"
 #include "polyMesh.H"
 #include "IOobjectList.H"
+#include "IOFieldField.H"
 #include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -67,6 +68,16 @@ tmp<IOField<Type> > reconstructLagrangianField
 );
 
 
+template<class Type>
+tmp<IOFieldField<Field<Type>, Type> > reconstructLagrangianFieldField
+(
+    const word& cloudName,
+    const polyMesh& mesh,
+    const PtrList<fvMesh>& meshes,
+    const word& fieldName
+);
+
+
 template<class Type>
 void reconstructLagrangianFields
 (
@@ -77,6 +88,16 @@ void reconstructLagrangianFields
 );
 
 
+template<class Type>
+void reconstructLagrangianFieldFields
+(
+    const word& cloudName,
+    const polyMesh& mesh,
+    const PtrList<fvMesh>& meshes,
+    const IOobjectList& objects
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/parallel/reconstruct/reconstruct/reconstructLagrangianFields.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "IOField.H"
+#include "IOFieldField.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
@@ -87,6 +88,67 @@ Foam::tmp<Foam::IOField<Type> > Foam::reconstructLagrangianField
 }
 
 
+template<class Type>
+Foam::tmp<Foam::IOFieldField<Foam::Field<Type>, Type> >
+Foam::reconstructLagrangianFieldField
+(
+    const word& cloudName,
+    const polyMesh& mesh,
+    const PtrList<fvMesh>& meshes,
+    const word& fieldName
+)
+{
+    // Construct empty field on mesh
+    tmp<IOFieldField<Field<Type>, Type > > tfield
+    (
+        new IOFieldField<Field<Type>, Type>
+        (
+            IOobject
+            (
+                fieldName,
+                mesh.time().timeName(),
+                cloud::prefix/cloudName,
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            Field<Field<Type> >(0)
+        )
+    );
+    Field<Field<Type> >& field = tfield();
+
+    forAll(meshes, i)
+    {
+        // Check object on local mesh
+        IOobject localIOobject
+        (
+            fieldName,
+            meshes[i].time().timeName(),
+            cloud::prefix/cloudName,
+            meshes[i],
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+
+        if (localIOobject.headerOk())
+        {
+            IOFieldField<Field<Type>, Type> fieldi(localIOobject);
+
+            label offset = field.size();
+            field.setSize(offset + fieldi.size());
+
+            forAll(fieldi, j)
+            {
+                field[offset + j] = fieldi[j];
+            }
+        }
+    }
+
+    return tfield;
+}
+
+
+
 template<class Type>
 void Foam::reconstructLagrangianFields
 (
@@ -122,4 +184,67 @@ void Foam::reconstructLagrangianFields
 }
 
 
+template<class Type>
+void Foam::reconstructLagrangianFieldFields
+(
+    const word& cloudName,
+    const polyMesh& mesh,
+    const PtrList<fvMesh>& meshes,
+    const IOobjectList& objects
+)
+{
+    {
+        const word fieldClassName(IOFieldField<Field<Type>, Type>::typeName);
+
+        IOobjectList fields = objects.lookupClass(fieldClassName);
+
+        if (fields.size())
+        {
+            Info<< "    Reconstructing lagrangian "
+                << fieldClassName << "s\n" << endl;
+
+            forAllConstIter(IOobjectList, fields, fieldIter)
+            {
+                Info<< "        " << fieldIter()->name() << endl;
+                reconstructLagrangianFieldField<Type>
+                (
+                    cloudName,
+                    mesh,
+                    meshes,
+                    fieldIter()->name()
+                )().write();
+            }
+
+            Info<< endl;
+        }
+    }
+
+    {
+        const word fieldClassName(IOField<Field<Type> >::typeName);
+
+        IOobjectList fields = objects.lookupClass(fieldClassName);
+
+        if (fields.size())
+        {
+            Info<< "    Reconstructing lagrangian "
+                << fieldClassName << "s\n" << endl;
+
+            forAllConstIter(IOobjectList, fields, fieldIter)
+            {
+                Info<< "        " << fieldIter()->name() << endl;
+                reconstructLagrangianFieldField<Type>
+                (
+                    cloudName,
+                    mesh,
+                    meshes,
+                    fieldIter()->name()
+                )().write();
+            }
+
+            Info<< endl;
+        }
+    }
+}
+
+
 // ************************************************************************* //