Merge remote-tracking branch 'origin/feature/solidMechanics' into nextRelease

.gitignore

@@ -109,4 +109,12 @@ src/lduSolvers/amg/amgPolicy/samgPolicy.H
 src/lduSolvers/amg/amgPolicy/aamgPolicy.C
 src/lduSolvers/amg/amgPolicy/aamgPolicy.H
 
+# The following files are blacklisted because of a DMCA complaint by ANSYS.
+src/lduSolvers/tools/PriorityArray.C
+src/lduSolvers/tools/PriorityArray.H
+src/lduSolvers/amg/amgPolicy/samgPolicy.C
+src/lduSolvers/amg/amgPolicy/samgPolicy.H
+src/lduSolvers/amg/amgPolicy/aamgPolicy.C
+src/lduSolvers/amg/amgPolicy/aamgPolicy.H
+
 # end-of-file

applications/solvers/solidMechanics/Allwclean

@@ -0,0 +1,20 @@
+#!/bin/sh
+set -x
+
+wclean solidModels
+
+wclean elasticContactSolidFoam
+wclean elasticContactIncrSolidFoam
+wclean elasticContactNonLinULSolidFoam
+wclean elasticGravitySolidFoam
+wclean elasticIncrSolidFoam
+wclean elasticNonLinTLSolidFoam
+wclean elasticNonLinULSolidFoam
+wclean elasticPlasticSolidFoam
+wclean elasticPlasticNonLinULSolidFoam
+wclean elasticSolidFoam
+wclean elasticThermalSolidFoam
+wclean icoFsiElasticNonLinULSolidFoam
+wclean viscoElasticSolidFoam
+
+(cd utilities && ./Allwclean)

applications/solvers/solidMechanics/Allwmake

@@ -0,0 +1,20 @@
+#!/bin/sh
+set -x
+
+wmake libso solidModels
+
+wmake elasticContactSolidFoam
+wmake elasticContactIncrSolidFoam
+wmake elasticContactNonLinULSolidFoam
+wmake elasticGravitySolidFoam
+wmake elasticIncrSolidFoam
+wmake elasticNonLinTLSolidFoam
+wmake elasticNonLinULSolidFoam
+wmake elasticPlasticNonLinULSolidFoam
+wmake elasticPlasticSolidFoam
+wmake elasticSolidFoam
+wmake elasticThermalSolidFoam
+wmake icoFsiElasticNonLinULSolidFoam
+wmake viscoElasticSolidFoam
+
+(cd utilities && ./Allwmake)

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude \
+    -I$(FOAM_SRC)/VectorN/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian \
+    -lVectorN

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/Make/options_myLib

@@ -0,0 +1,15 @@
+EXE_INC = \
+//    -I../../../myLibraries/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I../materialModels/lnInclude \
+    -I./solidInterface
+//    -I./patchToPatchInterpolation
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+//    -lfiniteVolume_philipc \
+//    -lfiniteVolume \
+    -lmaterialModels_philipc \
+    -llduSolvers \
+    -lmeshTools

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/Make/options_old

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/dynamicMesh/lnInclude \
+    -I../materialModels/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+    -lfiniteVolume \
+    -lmaterialModels_philipc

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/calculateDEpsilonDSigma.H

@@ -0,0 +1,3 @@
+DEpsilon = symm(gradDU);
+
+DSigma = 2*mu*DEpsilon + lambda*(I*tr(DEpsilon));

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/calculateDivDSigmaExp.H

@@ -0,0 +1,47 @@
+if(divDSigmaExpMethod == "standard")
+  {
+    divDSigmaExp = fvc::div
+      (
+       mu*gradDU.T() + lambda*(I*tr(gradDU)) - (mu + lambda)*gradDU,
+       "div(sigma)"
+       );
+  }
+ else if(divDSigmaExpMethod == "surface")
+   { 
+     divDSigmaExp = fvc::div
+       (
+	muf*(mesh.Sf() & fvc::interpolate(gradDU.T()))
+	+ lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradDU)))
+	- (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradDU))
+	);
+   }
+ else if(divDSigmaExpMethod == "decompose")
+   {
+     surfaceTensorField shearGradDU =
+       ((I - n*n)&fvc::interpolate(gradDU));
+     
+     divDSigmaExp = fvc::div
+       (
+	mesh.magSf()
+	*(
+	  - (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
+	  + lambdaf*tr(shearGradDU&(I - n*n))*n
+	  + muf*(shearGradDU&n)
+	  )
+	);
+   }
+ else if(divDSigmaExpMethod == "laplacian")
+   {
+     divDSigmaExp =
+       - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+       + fvc::div
+       (
+	mu*gradDU.T()
+	+ lambda*(I*tr(gradDU)),
+	"div(sigma)"
+	);
+   }
+ else
+   {
+     FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
+   }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magDU = mag(DU.internalField());
+
+    forAll(magDU, cellI)
+    {
+        if (magDU[cellI] < SMALL)
+        {
+            magDU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                DU.internalField()
+              - DU.prevIter().internalField()
+            )
+           /magDU
+        );
+}

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/correctGlobalFaceZoneMesh.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+correctGlobalFaceZoneMesh.H
+
+When there is a globalFaceZone and the mesh is moved by interpolating U to the
+vertices with volPointInterpolation, then there are two problems:
+    -some points on the patch with the faceZone are moved incorrectly, I think
+     it is because the faceZone has no U and causes an incorrect interpolation,
+    -the faceZones points not on the proc cells are not moved at all because
+     they have no U.
+
+So the points on the patch with the faceZone need to be fixed and also all the
+faceZone points need to be moved and synchronised so each proc has the same
+full faceZone mesh.
+
+The mapping of procs faceZone order of points to the master procs faceZone point
+order is kept in procToGlobalFZmap, which is calculated at the start of the run
+in the createGlobalToLocalFaceZonePointMap.H header file.
+
+Note: DU is used for updated Lagrangian solver instead of U
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
+    contactPatchPairList& contacts = contact;
+    
+    forAll(contacts, contactI)
+      {
+	label masterID = contacts[contactI].masterPatch().index();
+	label slaveID = contacts[contactI].slavePatch().index();
+	
+	primitivePatchInterpolation masterInterpolator
+	  (mesh.boundaryMesh()[masterID]);
+	primitivePatchInterpolation slaveInterpolator
+	  (mesh.boundaryMesh()[slaveID]);
+
+	//- U must be interpolated to the vertices, this ignores the faceZone
+	//- points with no U (unlike volPointInterpolation)
+	vectorField correctMasterPointU =
+	  masterInterpolator.faceToPointInterpolate<vector>
+	  (
+	   U.boundaryField()[masterID]
+	   );
+	vectorField correctSlavePointU =
+	slaveInterpolator.faceToPointInterpolate<vector>
+	  (
+	   U.boundaryField()[slaveID]
+	   );
+	
+	vectorField oldMasterPoints =
+	  mesh.boundaryMesh()[masterID].localPoints();
+	vectorField oldSlavePoints =
+	  mesh.boundaryMesh()[slaveID].localPoints();
+	
+	labelList masterPointLabels =
+	  mesh.boundaryMesh()[masterID].meshPoints();
+	labelList slavePointLabels =
+	  mesh.boundaryMesh()[slaveID].meshPoints();
+	
+	//- correct the patch newPoints
+	forAll(masterPointLabels, pointI)
+	  {
+	    label pointGlobalLabel = masterPointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldMasterPoints[pointI]
+	      +
+	      correctMasterPointU[pointI];
+	  }
+	forAll(slavePointLabels, pointI)
+	  {
+	    label pointGlobalLabel = slavePointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldSlavePoints[pointI]
+	      +
+	      correctSlavePointU[pointI];
+	  }
+      }
+
+    
+
+
+    //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
+
+    forAll(mesh.faceZones(), faceZoneI)
+      {
+	//- find the patch corresponding to this faceZone
+	//- assuming that the FZ is called <patch_name>FaceZone
+	string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	//- remove the string FaceZone from the end of the face zone name to get the patch name
+	string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	if(patchID == -1)
+	  {
+	    FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+		       << faceZoneName << exit(FatalError);
+	  }
+	    
+	vectorField globalFZpoints =
+	  mesh.faceZones()[faceZoneI]().localPoints();
+	
+	//- new points for the face zone
+	vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	//- inter-proc points are shared by multiple procs
+	//- pointNumProc is the number of procs which a point lies on
+	scalarField pointNumProcs(globalFZpoints.size(), 0.0);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
+	    if(pointOnLocalProcPatch[faceZoneI][localPoint])
+	      {
+		label procPoint =
+		  mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
+		globalFZnewPoints[globalPointI] =
+		  newPoints[procPoint];
+		pointNumProcs[globalPointI] = 1;
+	      }		
+	  }
+	
+	reduce(globalFZnewPoints, sumOp<vectorField>());
+	reduce(pointNumProcs, sumOp<scalarField>());
+	
+	//- now average the newPoints between all procs
+	if(min(pointNumProcs) < 1)
+	  {
+	    FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
+	  }
+	globalFZnewPoints /= pointNumProcs;
+	
+	//- the globalFZnewPoints now contains the correct FZ new points in
+	//- a global order, now convert them back into the local proc order
+	
+	vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    procFZnewPoints[localPoint] =
+	      globalFZnewPoints[globalPointI];
+	  }
+	
+	//- now fix the newPoints points on the globalFaceZones
+	labelList procFZmeshPoints =
+	  mesh.faceZones()[faceZoneI]().meshPoints();
+	
+	forAll(procFZmeshPoints, pointI)
+	  {
+	    label procPoint = procFZmeshPoints[pointI];
+	    newPoints[procPoint] =
+	      procFZnewPoints[pointI];
+	  }
+      }
+  }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/createFields.H

@@ -0,0 +1,115 @@
+    Info<< "Reading field DU\n" << endl;
+    volVectorField DU
+    (
+        IOobject
+        (
+            "DU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volTensorField gradDU = fvc::grad(DU);
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+	dimensionedVector("zero", dimLength, vector::zero)
+    );
+
+    volSymmTensorField DEpsilon
+    (
+        IOobject
+        (
+            "DEpsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField DSigma
+    (
+        IOobject
+        (
+            "DSigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+      volVectorField divDSigmaExp
+      (
+        IOobject
+        (
+	 "divDSigmaExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+
+   // read rheology properties
+    rheologyModel rheology(sigma);
+
+    volScalarField rho = rheology.rho();
+
+    volScalarField mu = rheology.mu();
+    volScalarField lambda = rheology.lambda();
+    surfaceScalarField muf = fvc::interpolate(mu, "mu");
+    surfaceScalarField lambdaf = fvc::interpolate(lambda, "lambda");
+
+    surfaceVectorField n = mesh.Sf()/mesh.magSf();
+
+   //- create contact problem
+   contactProblem contact(DU);

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/createGlobalToLocalFaceZonePointMap.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+createGlobalToLocalFaceZonePointMap.H
+
+I could not figure out the order of the points on a globalFaceZone, they is a
+different order for each processor where the processors own points come first.
+
+So I decide that the master proc order is the global order and I find the map
+from this order to each proc order here by just commparing actual point
+coordinates.
+
+This map is then used to correct and synchronise the globalFaceZone points on
+all the procs when the mesh is moved, in header file correctGlobalFaceZoneMesh.H.
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- procToGlobalFZmap is a map from the current proc faceZone point order to the
+//- master proc point order
+//- these are read if present to allow restarting of contact cases
+IOList<labelList> procToGlobalFZmap
+(
+ IOobject
+ (
+  "procToGlobalFZmap",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+IOList<labelList> pointOnLocalProcPatch
+(
+ IOobject
+ (
+  "pointOnLocalProcPatch",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+//- if they have been read then don't recalculate it
+bool globalFaceZoneMappingSet = false;
+if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
+  {
+    Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
+	 << endl; 
+    globalFaceZoneMappingSet = true;
+  }
+ else
+   {
+     Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl
+	  << "this message should only appear starting a new analysis" << endl;
+   }
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    if(!globalFaceZoneMappingSet)
+      {
+	forAll(mesh.faceZones(), faceZoneI)
+	  {
+	    vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
+	
+	    //- set all slave points to zero because only the master order is used
+	    if(!Pstream::master())
+	      globalFZpoints *= 0.0;
+	    
+	    //- pass points to all procs
+	    reduce(globalFZpoints, sumOp<vectorField>());
+	    
+	    
+	    //- now every proc has the master's list of FZ points
+	    //- every proc must now find the mapping from their local FZpoints to
+	    //- the globalFZpoints
+	    
+	    vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    forAll(globalFZpoints, globalPointI)
+	      {
+		forAll(procFZpoints, procPointI)
+		  {
+		    if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
+		      {
+			procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
+			break;
+		      }
+		  }
+	      }
+	    //- procToGlobalFZmap now contains the local FZpoint label for each
+	    //- global FZ point label - for each faceZone
+	    
+	    //- check what points are on the current proc patch
+	    pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
+	    
+	    //- find corresponding patch
+	    string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	    //- remove the string FaceZone from the end of the face zone name to get the patch name
+	    string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	    label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	    if(patchID == -1)
+	      {
+		FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+			   << faceZoneName << exit(FatalError);
+	      }
+	    
+	    forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
+	      {
+		forAll(mesh.boundaryMesh()[patchID].localPoints(), pi)
+		  {
+		    if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi])
+		      {
+			pointOnLocalProcPatch[faceZoneI][fzpi] = 1;
+			break;
+		      }
+		  }
+	      }
+	  }
+      } //- end if(!globalFaceZoneMappingSet)
+  }
+
+//- write to disk to allow restart of cases
+//- because it is not possible to calculate the 
+//- mapping after the meshes have moved
+if(!globalFaceZoneMappingSet && Pstream::parRun())
+  {
+    procToGlobalFZmap.write();
+    pointOnLocalProcPatch.write();
+  }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/createSolidInterface.H

@@ -0,0 +1,25 @@
+Switch solidInterfaceCorr(false);
+
+solidInterface* solidInterfacePtr(NULL);
+
+{
+    const dictionary& stressControl =
+      mesh.solutionDict().subDict("stressedFoam");
+
+    solidInterfaceCorr = Switch(stressControl.lookup("solidInterface"));
+
+    if(solidInterfaceCorr)
+      {
+	Info << "Creating solid interface correction" << endl;
+	solidInterfacePtr = new solidInterface(mesh, rheology);
+	solidInterfacePtr->modifyProperties(muf, lambdaf);
+	gradDU = solidInterfacePtr->grad(DU);
+	
+	//- solidInterface needs muf and lambdaf to be used for divSigmaExp
+	if(divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose")
+	  {
+	    FatalError << "divDSigmaExp must be decompose or surface when solidInterface is on"
+		       << exit(FatalError);
+	  }
+      }
+}

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/elasticContactIncrSolidFoam.C

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    elasticContactIncrSolidFoam
+
+Description
+    Transient/steady-state segregated finite-volume solver for small strain
+    elastic solid bodies in contact, using an incremental total Lagrangian
+    approach.
+    Works in parallel.
+
+    Solves for the displacement increment vector field DU, also generating the
+    stress tensor field sigma.
+
+    It is only for frictionless contact, friction not implemented yet.
+
+Author
+    Philip Cardiff
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "rheologyModel.H"
+#include "contactProblem.H"
+#include "solidInterface.H"
+
+#include "volPointInterpolation.H"
+#include "pointPatchInterpolation.H"
+#include "primitivePatchInterpolation.H"
+#include "fixedValuePointPatchFields.H"
+#include "pointFields.H"
+#include "pointMesh.H"
+#include "pointBoundaryMesh.H"
+#include "primitivePatchInterpolation.H"
+#include "twoDPointCorrector.H"
+//#include "leastSquaresVolPointInterpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+
+# include "createTime.H"
+
+# include "createMesh.H"
+
+# include "createFields.H"
+
+# include "readDivDSigmaExpMethod.H"
+
+# include "createGlobalToLocalFaceZonePointMap.H"
+
+# include "createSolidInterface.H"
+
+  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+  
+  Info<< "\nStarting time loop\n" << endl;
+  
+  for (runTime++; !runTime.end(); runTime++)
+    {
+      Info<< "Time: " << runTime.timeName() << endl;
+      
+#     include "readContactControls.H"
+
+#     include "readStressedFoamControls.H"
+  
+      //-- for moving the mesh and then back again
+      vectorField oldMeshPoints = mesh.allPoints();
+      
+      int iCorr = 0;
+      lduMatrix::solverPerformance solverPerf;
+      word solverName;
+      lduMatrix::debug = 0;
+      scalar residual = GREAT;
+      scalar initialResidual = 0;
+      scalar relativeResidual = GREAT;
+
+      //- reset DU to zero at the start of the time-step if
+      //- a predictor is not required
+      if(!predictor)
+	DU = dimensionedVector("zero", dimLength, vector::zero);
+
+      do //- start of momentum loop
+	{
+	  DU.storePrevIter();
+	  
+	  //- correct the contact boundaries 
+         if(iCorr % uEqnContactCorrFreq == 0)
+	    {
+	      Info << "\t\tCorrecting contact in the momentum loop "
+		   << "iteration: " << iCorr
+		   << ", residual: " << residual
+		   << endl;
+	      //#                 include "moveMeshLeastSquares.H"
+#             include "moveSolidMesh.H"
+	      contact.correct();
+	      mesh.movePoints(oldMeshPoints);
+	    }
+  
+#         include "calculateDivDSigmaExp.H"
+
+	  fvVectorMatrix DUEqn
+	    (
+	     fvm::d2dt2(rho, DU)
+	     ==
+	     fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
+	     + divDSigmaExp
+	     );
+	  
+          if(solidInterfaceCorr)
+            {
+              solidInterfacePtr->correct(DUEqn);
+            }
+
+	  solverPerf = DUEqn.solve();
+	  
+	  DU.relax();
+	  
+	  solverName = solverPerf.solverName();
+	  
+          if(solidInterfaceCorr)
+            {
+              gradDU = solidInterfacePtr->grad(DU);
+            }
+          else
+            {
+              gradDU = fvc::grad(DU);
+            }
+	  
+	  U = U.oldTime() + DU;
+
+	  residual = solverPerf.initialResidual();
+	  
+	  //****************************************************//
+	  // The contact residual is the initial residual for the
+	  // first iteration of the momentum equation
+	  //****************************************************//
+	  if(iCorr == 0)
+	    {
+	      initialResidual = solverPerf.initialResidual();
+	    }
+	  
+#         include "calculateRelativeResidual.H"	      
+
+	  Info << "\tTime " << runTime.value()
+	       << ", Corrector " << iCorr
+	       << ", Solving for " << DU.name()
+	       << " using " << solverPerf.solverName()
+	       << ", residual = " << solverPerf.initialResidual()
+	       << ", relative residual = " << relativeResidual << endl;
+	} //- end of momentum loop
+      while
+	(
+	 relativeResidual > convergenceTolerance
+	 //residual > convergenceTolerance
+	 &&
+	 ++iCorr < nCorr
+	 );
+      
+      // Print out info per contact iteration
+      Info << "\t\tSolving for " << DU.name()
+	   << " using " << solverName
+	   << ", Initial residual = " << initialResidual
+	   << ", Final residual = " << solverPerf.initialResidual()
+	   << ", No outer iterations " << iCorr << endl;
+      
+      lduMatrix::debug = 1;
+      
+#     include "calculateDEpsilonDSigma.H"
+            
+      epsilon += DEpsilon;
+      
+      sigma += DSigma;
+
+#     include "writeFields.H"
+  
+      //#     include "writeBoundaryNetForces.H"
+      
+      //#     include "moveMeshLeastSquares.H"
+      //#     include "moveSolidMesh.H"
+      //#     include "printContactResults.H"
+      //mesh.movePoints(oldMeshPoints);
+      
+      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+	  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+	  << endl << endl;
+    }
+  
+  Info<< "End\n" << endl;
+  
+  return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/moveMeshLeastSquares.H

@@ -0,0 +1,56 @@
+//--------------------------------------------------//
+//- move mesh
+//--------------------------------------------------//
+if(min(J.internalField()) > 0)
+  {
+    Info << "Moving mesh using least squares interpolation" << endl;
+
+    leastSquaresVolPointInterpolation pointInterpolation(mesh);
+  
+    // Create point mesh
+    pointMesh pMesh(mesh);
+    
+    wordList types
+      (
+       pMesh.boundary().size(),
+       calculatedFvPatchVectorField::typeName
+       );
+    
+    pointVectorField pointDU
+      (
+       IOobject
+       (
+        "pointDU",
+        runTime.timeName(),
+        mesh
+        ),
+       pMesh,
+       dimensionedVector("zero", dimLength, vector::zero),
+       types
+       );
+    
+    pointInterpolation.interpolate(DU, pointDU);
+
+    const vectorField& pointDUI =
+      pointDU.internalField();
+
+    //- Move mesh
+    vectorField newPoints = mesh.allPoints();
+
+    forAll (pointDUI, pointI)
+      {
+        newPoints[pointI] += pointDUI[pointI];
+      }
+
+    twoDPointCorrector twoDCorrector(mesh);
+    twoDCorrector.correctPoints(newPoints);
+    mesh.movePoints(newPoints);
+    mesh.V00();
+    mesh.moving(false);
+  }
+  else
+    {
+      FatalErrorIn(args.executable())
+	<< "Negative Jacobian"
+	<< exit(FatalError);
+    }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/moveSolidMesh.H

@@ -0,0 +1,28 @@
+{  
+  //- move mesh for the contact correction
+
+  // Create point interpolation
+  volPointInterpolation pointInterpolation(mesh);
+
+  // Calculate mesh points displacement
+  pointVectorField pointU = pointInterpolation.interpolate(U);
+
+  const vectorField& pointUI = pointU.internalField();
+  
+  // Move mesh
+  vectorField newPoints = mesh.allPoints();
+  
+  forAll (pointUI, pointI)
+    {
+      newPoints[pointI] += pointUI[pointI];
+    }
+
+# include "correctGlobalFaceZoneMesh.H"
+ 
+  twoDPointCorrector twoDCorrector(mesh);
+  twoDCorrector.correctPoints(newPoints);  
+  
+  mesh.movePoints(newPoints);
+  mesh.V00();
+  mesh.moving(false);
+}

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/printContactResults.H

@@ -0,0 +1,55 @@
+if (runTime.outputTime())
+  {
+    // FAILS IN PARALLEL - FIX
+    //    Info << "Print contact area" << endl;
+    //volScalarField ca = contact.contactArea();
+    //ca.write();
+    
+    //-------------------------------------------------------------//
+    // I couldn't get tmp to return the pointScalarField correctly //
+    // so I had to make the pointScalarField here and pass it to   //
+    // contactGapPoints and pointContactForce to populate          //
+    //-------------------------------------------------------------//
+    //This is the point distance for each contact vertex
+    pointScalarField cGapPoints 
+      (
+       IOobject
+       (
+	"pointContactGap",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedScalar("scalar", dimLength, 0.0),
+       "calculated"
+       );
+
+    contact.contactGapPoints(cGapPoints);
+    cGapPoints.write();
+	  
+    
+    //- This is the point distance for each contact vertex
+    pointVectorField cPointForce 
+      (
+       IOobject
+       (
+	"pointContactForce",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedVector("vector", dimForce, vector::zero),
+       "calculated"
+       );
+    contact.contactPointForce(cPointForce);
+    cPointForce.write();
+    
+    //- this is the actual (sigma&n)&n) on the contact patches
+    //- SHOULD THIS BE A REF TO A TMP...?
+    volScalarField cPressure = contact.contactPressure();
+    cPressure.write();
+  }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/readContactControls.H

@@ -0,0 +1,16 @@
+//****************************************************//
+// Read the contact tol and max contact correctors
+//****************************************************//
+const IOdictionary& contactControl =
+  IOobject
+  (
+   "contactProperties",
+   U.time().constant(),
+   U.db(),
+   IOobject::MUST_READ,
+   IOobject::NO_WRITE
+   );
+
+
+//- The frequnecy the conatct is corrected inside the momentum loop
+int uEqnContactCorrFreq(readInt(contactControl.lookup("innerContactCorrectFreq")));

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/readDivDSigmaExpMethod.H

@@ -0,0 +1,9 @@
+//- how explicit component of sigma is to be calculated
+word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp"));
+Info << divDSigmaExpMethod << " method chosen for calculation of sigmaExp" << endl;
+if(divDSigmaExpMethod != "standard" && divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose" && divDSigmaExpMethod != "laplacian")
+  {
+    FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
+	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	       << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/readStressedFoamControls.H

@@ -0,0 +1,7 @@
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("stressedFoam");
+
+    int nCorr(readInt(stressControl.lookup("nCorrectors")));
+    scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));
+
+Switch predictor(stressControl.lookup("predictor"));

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/writeBoundaryNetForces.H

@@ -0,0 +1,14 @@
+// * * * * * * * * * * * * * * * * NET FORCES * * * * * * * * * * * * * * *  //
+
+vectorField netForces(mesh.boundary().size(), vector::zero);
+
+Info << nl;
+forAll(netForces, patchI)
+{
+  netForces[patchI] = gSum(mesh.Sf().boundaryField()[patchI] & sigma.boundaryField()[patchI]);
+  
+  Info << "patch\t" << mesh.boundary()[patchI].name() << "\t\tnet force is\t"
+       << netForces[patchI] << " N" << endl;
+}   
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/solidMechanics/elasticContactIncrSolidFoam/writeFields.H

@@ -0,0 +1,77 @@
+if (runTime.outputTime())
+  {
+    volScalarField epsilonEq
+      (
+       IOobject
+       (
+	"epsilonEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+       );
+
+    Info<< "Max epsilonEq = " << max(epsilonEq).value()
+	<< endl;
+
+    volScalarField sigmaEq
+      (
+       IOobject
+       (
+	"sigmaEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((3.0/2.0)*magSqr(dev(sigma)))
+       );
+    
+    Info<< "Max sigmaEq = " << max(sigmaEq).value()
+	<< endl;
+
+    volScalarField pressure
+      (
+       IOobject
+       (
+        "pressure",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+        ),
+       tr(sigma)/3.0
+       );
+    //- boundary surface pressure
+    forAll(pressure.boundaryField(), patchi)
+      {
+        const vectorField& nb = n.boundaryField()[patchi];
+        pressure.boundaryField()[patchi] =
+          -(nb & ( nb & sigma.boundaryField()[patchi] ));
+      }
+
+    //- contact slave penetration
+    #     include "moveSolidMesh.H"
+    pointMesh pMesh(mesh);
+    pointScalarField cGapPoints
+      (
+       IOobject
+       (
+        "pointContactGap",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+        ),
+       pMesh,
+       dimensionedScalar("scalar", dimLength, 0.0),
+       "calculated"
+       );
+    contact.contactGapPoints(cGapPoints);
+    cGapPoints.write();
+    mesh.movePoints(oldMeshPoints);  
+
+    runTime.write();
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude \
+    -I$(FOAM_SRC)/VectorN/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian \
+    -lVectorN

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/Make/options_myLib

@@ -0,0 +1,15 @@
+EXE_INC = \
+//    -I../../../myLibraries/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I../materialModels/lnInclude \
+    -I./solidInterface
+//    -I./patchToPatchInterpolation
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+//    -lfiniteVolume_philipc \
+//    -lfiniteVolume \
+    -lmaterialModels_philipc \
+    -llduSolvers \
+    -lmeshTools

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/Make/options_old

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/dynamicMesh/lnInclude \
+    -I../materialModels/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+    -lfiniteVolume \
+    -lmaterialModels_philipc

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculateDEpsilonDSigma.H

@@ -0,0 +1,3 @@
+DEpsilon = symm(gradDU) + 0.5*symm(gradDU & gradDU.T());
+
+DSigma = 2*mu*DEpsilon + lambda*(I*tr(DEpsilon));

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculateDivDSigmaExp.H

@@ -0,0 +1,47 @@
+if(divDSigmaExpMethod == "standard")
+  {
+    divDSigmaExp = fvc::div
+      (
+       mu*gradDU.T() + lambda*(I*tr(gradDU)) - (mu + lambda)*gradDU,
+       "div(sigma)"
+       );
+  }
+ else if(divDSigmaExpMethod == "surface")
+   { 
+     divDSigmaExp = fvc::div
+       (
+	muf*(mesh.Sf() & fvc::interpolate(gradDU.T()))
+	+ lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradDU)))
+	- (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradDU))
+	);
+   }
+ else if(divDSigmaExpMethod == "decompose")
+   {
+     surfaceTensorField shearGradDU =
+       ((I - n*n)&fvc::interpolate(gradDU));
+     
+     divDSigmaExp = fvc::div
+       (
+	mesh.magSf()
+	*(
+	  - (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
+	  + lambdaf*tr(shearGradDU&(I - n*n))*n
+	  + muf*(shearGradDU&n)
+	  )
+	);
+   }
+ else if(divDSigmaExpMethod == "laplacian")
+   {
+     divDSigmaExp =
+       - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+       + fvc::div
+       (
+	mu*gradDU.T()
+	+ lambda*(I*tr(gradDU)),
+	"div(sigma)"
+	);
+   }
+ else
+   {
+     FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
+   }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculateDivDSigmaExpLargeStrain.H

@@ -0,0 +1,34 @@
+//----------------------------------------------------//
+//- sigma explicit large strain explicit terms
+//----------------------------------------------------//
+if(divDSigmaLargeStrainExpMethod == "standard")
+  {
+    divDSigmaLargeStrainExp =
+      fvc::div
+      (
+       mu*(gradDU & gradDU.T())	
+       + 0.5*lambda*(gradDU && gradDU)*I //- equivalent to 0.5*lambda*(I*tr(gradDU & gradDU.T()))
+       + ((sigma + DSigma) & DF.T()),
+       "div(sigma)"
+       );
+  }
+ else if(divDSigmaLargeStrainExpMethod == "surface")
+   {
+     divDSigmaLargeStrainExp =
+       fvc::div
+       (
+	muf * (mesh.Sf() & fvc::interpolate(gradDU & gradDU.T()))
+	+ 0.5*lambdaf * (mesh.Sf() & (fvc::interpolate(gradDU && gradDU)*I))
+	+ (mesh.Sf() & fvc::interpolate( sigma & DF.T() ))
+	+ (mesh.Sf() & fvc::interpolate(DSigma & DF.T() ))
+	);
+   }
+ else
+   {
+     FatalError
+       << "divDSigmaLargeStrainExp not found!"
+       << exit(FatalError);
+   }
+
+//- relax large strain component
+divDSigmaLargeStrainExp.relax();

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculateExtrapolationVectors.H

@@ -0,0 +1,73 @@
+//void Foam::edgeCorrectedVolPointInterpolation::makeExtrapolationWeights() const
+//{
+// Interpolation:
+// For each point to be corrected calculate the extrapolated value
+// for each face around it and combine them using the inverse
+// distance weighting factors
+
+//const labelList& ptc = boundaryPoints();
+
+// Calculate the correction vectors
+
+//extrapolationVectorsPtr_ =
+//  new FieldField<Field, vector>(ptc.size());
+//FieldField<Field, vector>& extraVecs = *extrapolationVectorsPtr_;
+FieldField<Field, vector> extraVecs(ptc.size());
+
+{
+  const labelListList& pfaces = mesh.pointFaces();
+
+  const volVectorField& centres = mesh.C();
+
+  const fvBoundaryMesh& bm = mesh.boundary();
+
+  forAll (ptc, pointI)
+    {
+      const label curPoint = ptc[pointI];
+      
+      const labelList& curFaces = pfaces[curPoint];
+      
+      // extraVecs.hook(new vectorField(curFaces.size())); //- no hook function
+       extraVecs.set
+       (
+	pointI,
+	new vectorField(curFaces.size())
+	);
+
+      vectorField& curExtraVectors = extraVecs[pointI];
+      
+      label nFacesAroundPoint = 0;
+
+      const vector& pointLoc = mesh.points()[curPoint];
+      
+      // Go through all the faces
+      forAll (curFaces, faceI)
+	{
+	  if (!mesh.isInternalFace(curFaces[faceI]))
+	    {
+	      // This is a boundary face.  If not in the empty patch
+	      // or coupled calculate the extrapolation vector
+	      label patchID =
+		mesh.boundaryMesh().whichPatch(curFaces[faceI]);
+	      
+	      if
+		(
+		 !isA<emptyFvPatch>(bm[patchID])
+		 && !bm[patchID].coupled()
+		 )
+		{
+		  // Found a face for extrapolation
+		  curExtraVectors[nFacesAroundPoint] =
+		    pointLoc
+		    - centres.boundaryField()[patchID]
+		    [bm[patchID].patch().whichFace(curFaces[faceI])];
+		  
+		  nFacesAroundPoint++;
+		}
+	    }
+	}
+      
+      curExtraVectors.setSize(nFacesAroundPoint);
+    }
+
+}

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculatePointBoundaryWeights.H

@@ -0,0 +1,121 @@
+//void volPointInterpolation::makeBoundaryWeights() const
+//{
+// const labelList& ptc = boundaryPoints();
+
+// Calculate the correction vectors and weighting factors
+//pointBoundaryWeightsPtr_ = new FieldField<Field, scalar>(ptc.size());
+//FieldField<Field, scalar>& w = *pointBoundaryWeightsPtr_;
+FieldField<Field, scalar> w(ptc.size());
+
+{
+  const labelListList& pf = mesh.pointFaces();
+
+  const volVectorField& centres = mesh.C();
+
+  const fvBoundaryMesh& bm = mesh.boundary();
+  
+  pointScalarField volPointSumWeights
+    (
+     IOobject
+     (
+      "volPointSumWeights",
+      mesh.polyMesh::instance(),
+      mesh
+      ),
+     pMesh,
+     dimensionedScalar("zero", dimless, 0)
+     );
+  
+  forAll (ptc, pointI)
+    {
+      const label curPoint = ptc[pointI];
+      
+      const labelList& curFaces = pf[curPoint];
+      
+      //w.hook(new scalarField(curFaces.size())); //philipc no hook function
+      w.set
+       (
+	pointI,
+	new scalarField(curFaces.size())
+	);
+
+      scalarField& curWeights = w[pointI];
+      
+      label nFacesAroundPoint = 0;
+      
+      const vector& pointLoc = mesh.points()[curPoint];
+      
+      // Go through all the faces
+      forAll (curFaces, faceI)
+	{
+	  if (!mesh.isInternalFace(curFaces[faceI]))
+	    {
+	      // This is a boundary face.  If not in the empty patch
+	      // or coupled calculate the extrapolation vector
+	      label patchID =
+		mesh.boundaryMesh().whichPatch(curFaces[faceI]);
+	      
+	      if
+		(
+		 !isA<emptyFvPatch>(bm[patchID])
+		 && !(
+		      bm[patchID].coupled()
+		      //&& Pstream::parRun()
+		      //&& !mesh.parallelData().cyclicParallel()
+		      )
+		 )
+		{
+		  curWeights[nFacesAroundPoint] =
+		    1.0/mag
+		    (
+		     pointLoc
+		     - centres.boundaryField()[patchID]
+		     [
+		      bm[patchID].patch().whichFace(curFaces[faceI])
+		      ]
+		     );
+		  
+		  nFacesAroundPoint++;
+		}
+	    }
+	}
+      
+      // Reset the sizes of the local weights
+      curWeights.setSize(nFacesAroundPoint);
+      
+      // Collect the sum of weights for parallel correction
+      volPointSumWeights[curPoint] += sum(curWeights);
+    }
+  
+  // Do parallel correction of weights
+  
+  // Update coupled boundaries
+  // Work-around for cyclic parallels.  
+  /*if (Pstream::parRun() && !mesh.parallelData().cyclicParallel())
+    {
+    forAll (volPointSumWeights.boundaryField(), patchI)
+    {
+    if (volPointSumWeights.boundaryField()[patchI].coupled())
+    {
+    volPointSumWeights.boundaryField()[patchI].initAddField();
+    }
+    }
+    
+    forAll (volPointSumWeights.boundaryField(), patchI)
+    {
+    if (volPointSumWeights.boundaryField()[patchI].coupled())
+    {
+    volPointSumWeights.boundaryField()[patchI].addField
+    (
+    volPointSumWeights.internalField()
+    );
+    }
+    }
+    }*/
+  
+  // Re-scale the weights for the current point
+  forAll (ptc, pointI)
+    {
+      w[pointI] /= volPointSumWeights[ptc[pointI]];
+    }
+}

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magDU = mag(DU.internalField());
+
+    forAll(magDU, cellI)
+    {
+        if (magDU[cellI] < SMALL)
+        {
+            magDU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                DU.internalField()
+              - DU.prevIter().internalField()
+            )
+           /magDU
+        );
+}

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/correctGlobalFaceZoneMesh.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+correctGlobalFaceZoneMesh.H
+
+When there is a globalFaceZone and the mesh is moved by interpolating U to the
+vertices with volPointInterpolation, then there are two problems:
+    -some points on the patch with the faceZone are moved incorrectly, I think
+     it is because the faceZone has no U and causes an incorrect interpolation,
+    -the faceZones points not on the proc cells are not moved at all because
+     they have no U.
+
+So the points on the patch with the faceZone need to be fixed and also all the
+faceZone points need to be moved and synchronised so each proc has the same
+full faceZone mesh.
+
+The mapping of procs faceZone order of points to the master procs faceZone point
+order is kept in procToGlobalFZmap, which is calculated at the start of the run
+in the createGlobalToLocalFaceZonePointMap.H header file.
+
+Note: DU is used for updated Lagrangian solver instead of U
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
+    contactPatchPairList& contacts = contact;
+    
+    forAll(contacts, contactI)
+      {
+	label masterID = contacts[contactI].masterPatch().index();
+	label slaveID = contacts[contactI].slavePatch().index();
+	
+	primitivePatchInterpolation masterInterpolator
+	  (mesh.boundaryMesh()[masterID]);
+	primitivePatchInterpolation slaveInterpolator
+	  (mesh.boundaryMesh()[slaveID]);
+
+	//- DU must be interpolated to the vertices, this ignores the faceZone
+	//- points with no DU (unlike volPointInterpolation)
+	vectorField correctMasterPointDU =
+	  masterInterpolator.faceToPointInterpolate<vector>
+	  (
+	   DU.boundaryField()[masterID]
+	   );
+	vectorField correctSlavePointDU =
+	slaveInterpolator.faceToPointInterpolate<vector>
+	  (
+	   DU.boundaryField()[slaveID]
+	   );
+	
+	vectorField oldMasterPoints =
+	  mesh.boundaryMesh()[masterID].localPoints();
+	vectorField oldSlavePoints =
+	  mesh.boundaryMesh()[slaveID].localPoints();
+	
+	labelList masterPointLabels =
+	  mesh.boundaryMesh()[masterID].meshPoints();
+	labelList slavePointLabels =
+	  mesh.boundaryMesh()[slaveID].meshPoints();
+	
+	//- correct the patch newPoints
+	forAll(masterPointLabels, pointI)
+	  {
+	    label pointGlobalLabel = masterPointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldMasterPoints[pointI]
+	      +
+	      correctMasterPointDU[pointI];
+	  }
+	forAll(slavePointLabels, pointI)
+	  {
+	    label pointGlobalLabel = slavePointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldSlavePoints[pointI]
+	      +
+	      correctSlavePointDU[pointI];
+	  }
+      }
+
+    
+
+
+    //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
+
+    forAll(mesh.faceZones(), faceZoneI)
+      {
+	//- find the patch corresponding to this faceZone
+	//- assuming that the FZ is called <patch_name>FaceZone
+	string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	//- remove the string FaceZone from the end of the face zone name to get the patch name
+	string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	if(patchID == -1)
+	  {
+	    FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+		       << faceZoneName << exit(FatalError);
+	  }
+	    
+	vectorField globalFZpoints =
+	  mesh.faceZones()[faceZoneI]().localPoints();
+	
+	//- new points for the face zone
+	vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	//- inter-proc points are shared by multiple procs
+	//- pointNumProc is the number of procs which a point lies on
+	scalarField pointNumProcs(globalFZpoints.size(), 0.0);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
+	    if(pointOnLocalProcPatch[faceZoneI][localPoint])
+	      {
+		label procPoint =
+		  mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
+		globalFZnewPoints[globalPointI] =
+		  newPoints[procPoint];
+		pointNumProcs[globalPointI] = 1;
+	      }		
+	  }
+	
+	reduce(globalFZnewPoints, sumOp<vectorField>());
+	reduce(pointNumProcs, sumOp<scalarField>());
+	
+	//- now average the newPoints between all procs
+	if(min(pointNumProcs) < 1)
+	  {
+	    FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
+	  }
+	globalFZnewPoints /= pointNumProcs;
+	
+	//- the globalFZnewPoints now contains the correct FZ new points in
+	//- a global order, now convert them back into the local proc order
+	
+	vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    procFZnewPoints[localPoint] =
+	      globalFZnewPoints[globalPointI];
+	  }
+	
+	//- now fix the newPoints points on the globalFaceZones
+	labelList procFZmeshPoints =
+	  mesh.faceZones()[faceZoneI]().meshPoints();
+	
+	forAll(procFZmeshPoints, pointI)
+	  {
+	    label procPoint = procFZmeshPoints[pointI];
+	    newPoints[procPoint] =
+	      procFZnewPoints[pointI];
+	  }
+      }
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/createFields.H

@@ -0,0 +1,132 @@
+    Info<< "Reading field DU\n" << endl;
+    volVectorField DU
+    (
+        IOobject
+        (
+            "DU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volTensorField gradDU = fvc::grad(DU);
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+	dimensionedVector("zero", dimLength, vector::zero)
+    );
+
+    volSymmTensorField DEpsilon
+    (
+        IOobject
+        (
+            "DEpsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField DSigma
+    (
+        IOobject
+        (
+            "DSigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+      volVectorField divDSigmaExp
+      (
+        IOobject
+        (
+	 "divDSigmaExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+      volVectorField divDSigmaLargeStrainExp
+      (
+        IOobject
+        (
+	 "divDSigmaLargeStrainExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+   // read rheology properties
+    rheologyModel rheology(sigma);
+
+    volScalarField rho = rheology.rho();
+
+    volScalarField mu = rheology.mu();
+    volScalarField lambda = rheology.lambda();
+    surfaceScalarField muf = fvc::interpolate(mu, "mu");
+    surfaceScalarField lambdaf = fvc::interpolate(lambda, "lambda");
+
+    surfaceVectorField n = mesh.Sf()/mesh.magSf();
+
+    volTensorField DF = gradDU.T();
+    volTensorField F = I + DF;
+    volScalarField J = det(F);
+
+   //- create contact problem
+   contactProblem contact(DU);

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/createGlobalToLocalFaceZonePointMap.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+createGlobalToLocalFaceZonePointMap.H
+
+I could not figure out the order of the points on a globalFaceZone, they is a
+different order for each processor where the processors own points come first.
+
+So I decide that the master proc order is the global order and I find the map
+from this order to each proc order here by just commparing actual point
+coordinates.
+
+This map is then used to correct and synchronise the globalFaceZone points on
+all the procs when the mesh is moved, in header file correctGlobalFaceZoneMesh.H.
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- procToGlobalFZmap is a map from the current proc faceZone point order to the
+//- master proc point order
+//- these are read if present to allow restarting of contact cases
+IOList<labelList> procToGlobalFZmap
+(
+ IOobject
+ (
+  "procToGlobalFZmap",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+IOList<labelList> pointOnLocalProcPatch
+(
+ IOobject
+ (
+  "pointOnLocalProcPatch",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+//- if they have been read then don't recalculate it
+bool globalFaceZoneMappingSet = false;
+if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
+  {
+    Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
+	 << endl; 
+    globalFaceZoneMappingSet = true;
+  }
+ else
+   {
+     Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl
+	  << "this message should only appear starting a new analysis" << endl;
+   }
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    if(!globalFaceZoneMappingSet)
+      {
+	forAll(mesh.faceZones(), faceZoneI)
+	  {
+	    vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
+	
+	    //- set all slave points to zero because only the master order is used
+	    if(!Pstream::master())
+	      globalFZpoints *= 0.0;
+	    
+	    //- pass points to all procs
+	    reduce(globalFZpoints, sumOp<vectorField>());
+	    
+	    
+	    //- now every proc has the master's list of FZ points
+	    //- every proc must now find the mapping from their local FZpoints to
+	    //- the globalFZpoints
+	    
+	    vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    forAll(globalFZpoints, globalPointI)
+	      {
+		forAll(procFZpoints, procPointI)
+		  {
+		    if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
+		      {
+			procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
+			break;
+		      }
+		  }
+	      }
+	    //- procToGlobalFZmap now contains the local FZpoint label for each
+	    //- global FZ point label - for each faceZone
+	    
+	    //- check what points are on the current proc patch
+	    pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
+	    
+	    //- find corresponding patch
+	    string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	    //- remove the string FaceZone from the end of the face zone name to get the patch name
+	    string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	    label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	    if(patchID == -1)
+	      {
+		FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+			   << faceZoneName << exit(FatalError);
+	      }
+	    
+	    forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
+	      {
+		forAll(mesh.boundaryMesh()[patchID].localPoints(), pi)
+		  {
+		    if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi])
+		      {
+			pointOnLocalProcPatch[faceZoneI][fzpi] = 1;
+			break;
+		      }
+		  }
+	      }
+	  }
+      } //- end if(!globalFaceZoneMappingSet)
+  }
+
+//- write to disk to allow restart of cases
+//- because it is not possible to calculate the 
+//- mapping after the meshes have moved
+if(!globalFaceZoneMappingSet)
+  {
+    procToGlobalFZmap.write();
+    pointOnLocalProcPatch.write();
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/elasticContactNonLinULSolidFoam.C

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    elasticContactNonLinSolidFoam
+
+Description
+    Transient/steady-state segregated finite-volume solver for large strain
+    elastic solid bodies in contact, using an incremental updated Lagrangian
+    approach.
+
+    Works in parallel but mesh.movePoints sometimes fails for some unknown
+    reason depending on the decomposition.
+
+    Solves for the displacement increment vector field DU, also generating the
+    stress tensor field sigma.
+
+    It is only for frictionless contact, friction not implemented yet.
+
+Author
+    Philip Cardiff
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "rheologyModel.H"
+#include "contactProblem.H"
+
+#include "volPointInterpolation.H"
+#include "leastSquaresVolPointInterpolation.H"
+#include "pointPatchInterpolation.H"
+#include "primitivePatchInterpolation.H"
+#include "fixedValuePointPatchFields.H"
+#include "pointFields.H"
+#include "pointMesh.H"
+#include "pointBoundaryMesh.H"
+#include "primitivePatchInterpolation.H"
+#include "twoDPointCorrector.H"
+#include "plane.H"
+#include "meshSearch.H"
+
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+
+# include "createTime.H"
+
+# include "createMesh.H"
+
+# include "createFields.H"
+
+# include "readDivDSigmaExpMethod.H"
+
+# include "readDivDSigmaLargeStrainMethod.H"
+
+# include "readMoveMeshMethod.H"
+
+# include "createGlobalToLocalFaceZonePointMap.H"
+
+  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+  
+  Info<< "\nStarting time loop\n" << endl;
+  
+  for (runTime++; !runTime.end(); runTime++)
+    {
+      Info<< "Time: " << runTime.timeName() << endl;
+      
+#     include "readContactControls.H"
+
+#     include "readStressedFoamControls.H"
+  
+      //-- for moving the mesh and then back again
+      vectorField oldMeshPoints = mesh.allPoints();
+      
+      int iCorr = 0;
+      lduMatrix::solverPerformance solverPerf;
+      word solverName;
+      lduMatrix::debug = 0;
+      scalar residual = GREAT;
+      scalar initialResidual = 0;
+      scalar relativeResidual = GREAT;
+            
+      do //- start of momentum loop
+	{
+	  DU.storePrevIter();
+	  
+	  divDSigmaLargeStrainExp.storePrevIter();
+
+	  //- correct the contact boundaries 
+	  if(iCorr % uEqnContactCorrFreq == 0)
+	    {
+	      Info << "\t\tCorrecting contact in the momentum loop "
+		   << "iteration: " << iCorr
+		   << ", residual: " << residual
+		   << endl;
+	      //#                 include "moveMeshLeastSquares.H"
+#             include "moveSolidMeshForContact.H"
+	      contact.correct();
+	      mesh.movePoints(oldMeshPoints);
+	    }
+	  
+#         include "calculateDivDSigmaExp.H"
+	  
+#         include "calculateDivDSigmaExpLargeStrain.H"
+
+	  fvVectorMatrix DUEqn
+	    (
+	     fvm::d2dt2(rho, DU)
+	     ==
+	     fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
+	     + divDSigmaExp
+	     + divDSigmaLargeStrainExp
+
+	     );
+	  
+	  solverPerf = DUEqn.solve();
+	  
+	  DU.relax();
+	  
+	  solverName = solverPerf.solverName();
+	  
+	  gradDU = fvc::grad(DU);
+	  
+	  DF = gradDU.T();
+
+#         include "calculateDEpsilonDSigma.H"
+
+	  residual = solverPerf.initialResidual();
+	  
+	  if(iCorr == 0)
+	    {
+	      initialResidual = solverPerf.initialResidual();
+	    }
+	  
+#         include "calculateRelativeResidual.H"	      
+	  
+	  Info << "\tTime " << runTime.value()
+	       << ", Corrector " << iCorr
+	       << ", Solving for " << DU.name()
+	       << " using " << solverPerf.solverName()
+	       << ", residual = " << solverPerf.initialResidual()
+	       << ", relative residual = " << relativeResidual << endl;
+	} //- end of momentum loop
+      while
+	(
+	 relativeResidual > convergenceTolerance
+	 //residual > convergenceTolerance
+	 &&
+	 ++iCorr < nCorr
+	 );
+      
+      // Print out info per contact iteration
+      Info << "\t\tSolving for " << DU.name()
+	   << " using " << solverName
+	   << ", Initial residual = " << initialResidual 
+	   << ", Final residual = " << solverPerf.initialResidual()
+	   << ", No outer iterations " << iCorr << endl;
+      
+      lduMatrix::debug = 1;
+      
+#     include "rotateFields.H"
+
+#     include "moveMesh.H"
+
+#     include "writeFields.H"
+        
+      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+	  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+	  << endl << endl;
+    }
+  
+  Info<< "End\n" << endl;
+  
+  return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/findBoundaryPoints.H

@@ -0,0 +1,39 @@
+//void volPointInterpolation::makeBoundaryAddressing() const
+//{
+
+// Go through all patches and mark up the external edge points
+labelHashSet pointsCorrectionMap(max(mesh.nPoints()/10, 100));
+
+const fvBoundaryMesh& bm = mesh.boundary();
+
+forAll (bm, patchI)
+{
+  // If the patch is empty, skip it
+  // If the patch is coupled, and there are no cyclic parallels, skip it
+  if
+    (
+     !isA<emptyFvPatch>(bm[patchI])
+     && !(
+	  bm[patchI].coupled()
+	  //&& Pstream::parRun()
+	  //&& !mesh.parallelData().cyclicParallel()
+	  )
+     )
+    {
+      const labelList& bp = bm[patchI].patch().boundaryPoints();
+      
+      const labelList& meshPoints = bm[patchI].patch().meshPoints();
+
+      forAll (bp, pointI)
+	{
+	  pointsCorrectionMap.insert(meshPoints[bp[pointI]]);
+	}
+    }
+}
+
+// Grab the points to correct
+//boundaryPointsPtr_ = new labelList(pointsCorrectionMap.toc());
+//labelList& ptc = *boundaryPointsPtr_;
+labelList ptc(pointsCorrectionMap.toc());
+
+sort(ptc);

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/moveMesh.H

@@ -0,0 +1,15 @@
+if(moveMeshMethod == "inverseDistance")
+  {
+#   include "moveMeshInverseDistance.H"
+  }
+ else if(moveMeshMethod == "leastSquares")
+   {
+#   include "moveMeshLeastSquares.H"
+   }
+ else
+   {
+     FatalError << "move mesh method " << moveMeshMethod << " not recognised" << nl
+		<< "available methods are:" << nl
+		<< "inverseDistance" << nl
+		<< "leastSquares" << exit(FatalError);
+   }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/moveMeshInverseDistance.H

@@ -0,0 +1,66 @@
+//--------------------------------------------------//
+//- move mesh
+//--------------------------------------------------//
+if(min(J.internalField()) > 0)
+  {
+    Info << "Move solid mesh using inverse distance interpolation" << endl;
+    
+    // Create point mesh
+    pointMesh pMesh(mesh);
+
+    // Create point interpolation
+    volPointInterpolation pointInterpolation(mesh);
+    
+    wordList types
+      (
+       pMesh.boundary().size(),
+       //fixedValueFvPatchVectorField::typeName
+       calculatedFvPatchVectorField::typeName
+       );
+
+    pointVectorField pointDU
+      (
+       IOobject
+       (
+	"pointDU",
+	runTime.timeName(),
+	mesh
+	),
+       pMesh,
+       dimensionedVector("zero", dimLength, vector::zero),
+       types
+       );
+
+    // Calculate mesh points displacement
+    pointInterpolation.interpolate(DU, pointDU);
+
+    //- correct edge interpolation
+    //- this is the stuff from edgeCorrectedVolPointInterpolation but
+    //- that class no longer works
+#   include "performEdgeCorrectedVolPointInterpolation.H"
+
+    //pointDU.write();
+
+    const vectorField& pointDUI =
+      pointDU.internalField();
+
+    // Move mesh
+    vectorField newPoints = mesh.allPoints();
+    
+     forAll (pointDUI, pointI)
+       {
+ 	newPoints[pointI] += pointDUI[pointI];
+      }
+  
+    twoDPointCorrector twoDCorrector(mesh);
+    twoDCorrector.correctPoints(newPoints);
+    mesh.movePoints(newPoints);
+    mesh.V00();
+    mesh.moving(false);
+  }
+  else
+    {
+      FatalErrorIn(args.executable())
+	<< "Negative Jacobian"
+	<< exit(FatalError);
+    }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/moveMeshLeastSquares.H

@@ -0,0 +1,56 @@
+//--------------------------------------------------//
+//- move mesh
+//--------------------------------------------------//
+if(min(J.internalField()) > 0)
+  {
+    Info << "Moving mesh using least squares interpolation" << endl;
+
+    leastSquaresVolPointInterpolation pointInterpolation(mesh);
+  
+    // Create point mesh
+    pointMesh pMesh(mesh);
+    
+    wordList types
+      (
+       pMesh.boundary().size(),
+       calculatedFvPatchVectorField::typeName
+       );
+    
+    pointVectorField pointDU
+      (
+       IOobject
+       (
+        "pointDU",
+        runTime.timeName(),
+        mesh
+        ),
+       pMesh,
+       dimensionedVector("zero", dimLength, vector::zero),
+       types
+       );
+    
+    pointInterpolation.interpolate(DU, pointDU);
+
+    const vectorField& pointDUI =
+      pointDU.internalField();
+
+    //- Move mesh
+    vectorField newPoints = mesh.allPoints();
+
+    forAll (pointDUI, pointI)
+      {
+        newPoints[pointI] += pointDUI[pointI];
+      }
+
+    twoDPointCorrector twoDCorrector(mesh);
+    twoDCorrector.correctPoints(newPoints);
+    mesh.movePoints(newPoints);
+    mesh.V00();
+    mesh.moving(false);
+  }
+  else
+    {
+      FatalErrorIn(args.executable())
+	<< "Negative Jacobian"
+	<< exit(FatalError);
+    }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/moveSolidMeshForContact.H

@@ -0,0 +1,27 @@
+{  
+  // Create point interpolation
+  volPointInterpolation pointInterpolation(mesh);
+
+  // Calculate mesh points displacement
+  pointVectorField pointDU = pointInterpolation.interpolate(DU);
+
+  vectorField pointDUI = pointDU.internalField();
+  
+  // Move mesh
+  vectorField newPoints = mesh.allPoints();
+  
+  forAll (pointDUI, pointI)
+    {
+      newPoints[pointI] += pointDUI[pointI];
+    }
+
+# include "correctGlobalFaceZoneMesh.H"
+ 
+  twoDPointCorrector twoDCorrector(mesh);
+  twoDCorrector.correctPoints(newPoints);  
+  
+  mesh.movePoints(newPoints);
+  //  pMesh.movePoints(newPoints);
+  mesh.V00();
+  mesh.moving(false);
+}

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/performEdgeCorrectedVolPointInterpolation.H

@@ -0,0 +1,133 @@
+// Do "normal" interpolation
+//volPointInterpolation::interpolate(vf, pf);
+
+//- interpolation is done just before this file
+pointVectorField& pf = pointDU;
+
+
+// Do the correction
+//GeometricField<Type, pointPatchField, pointMesh>  pfCorr
+/*pointVectorField pfCorr
+(
+ IOobject
+ (
+  //  "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
+  "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
+  //vf.instance(),
+  DU,
+  pMesh,
+  IOobject::NO_READ,
+  IOobject::NO_WRITE
+  ),
+ pMesh,
+ //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
+ dimensionedVector("zero", pf.dimensions(), vector::zero),
+ pf.boundaryField().types()
+ );*/
+
+pointVectorField pfCorr 
+(
+ IOobject
+ (
+  "pointDUcorr",
+  runTime.timeName(),
+  mesh,
+  IOobject::NO_READ,
+  IOobject::NO_WRITE
+  ),
+ pMesh,
+ dimensionedVector("vector", dimLength, vector::zero),
+ "calculated"
+ );
+
+//const labelList& ptc = boundaryPoints();
+#include "findBoundaryPoints.H"
+
+//const FieldField<Field, vector>& extraVecs = extrapolationVectors(); ///*********
+#include "calculateExtrapolationVectors.H"
+
+//const FieldField<Field, scalar>& w = pointBoundaryWeights(); ///*********
+#include "calculatePointBoundaryWeights.H"
+
+//Info << "w: " << w << endl;
+
+const labelListList& PointFaces = mesh.pointFaces();
+
+//const fvBoundaryMesh& bm = mesh.boundary(); // declared in findBoundaryPoints.H
+
+forAll (ptc, pointI)
+{
+  const label curPoint = ptc[pointI];
+  
+  const labelList& curFaces = PointFaces[curPoint];
+  
+  label fI = 0;
+  
+  // Go through all the faces
+  forAll (curFaces, faceI)
+    {
+      if (!mesh.isInternalFace(curFaces[faceI]))
+	{
+	  // This is a boundary face.  If not in the empty patch
+	  // or coupled calculate the extrapolation vector
+	  label patchID =
+	    mesh.boundaryMesh().whichPatch(curFaces[faceI]);
+	  
+	  if
+	    (
+	     !isA<emptyFvPatch>(mesh.boundary()[patchID])
+	     && !mesh.boundary()[patchID].coupled()
+	     )
+	    {
+	      label faceInPatchID =
+		bm[patchID].patch().whichFace(curFaces[faceI]);
+	      
+	      pfCorr[curPoint] +=
+		w[pointI][fI]*
+		(
+		 extraVecs[pointI][fI]
+		 & gradDU.boundaryField()[patchID][faceInPatchID]
+		 );
+	      
+	      fI++;
+	    }
+	}
+    }
+}
+
+// Update coupled boundaries
+/*forAll (pfCorr.boundaryField(), patchI)
+{
+  if (pfCorr.boundaryField()[patchI].coupled())
+    {
+      pfCorr.boundaryField()[patchI].initAddField();
+    }
+    }*/
+
+ /*forAll (pfCorr.boundaryField(), patchI)
+{
+  if (pfCorr.boundaryField()[patchI].coupled())
+    {
+      pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
+    }
+    }*/
+
+
+  //Info << "pfCorr: " << pfCorr << endl;
+ pfCorr.correctBoundaryConditions();
+
+//pfCorr.write();
+
+//Info << "pf: " << pf << endl;
+pf += pfCorr;
+
+//- this was already commented
+// Replace extrapolated values in pf
+//     forAll (ptc, pointI)
+//     {
+//         pf[ptc[pointI]] = pfCorr[ptc[pointI]];
+//     }
+
+pf.correctBoundaryConditions();
+
+//pf.write();

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/printContactResults.H

@@ -0,0 +1,55 @@
+if (runTime.outputTime())
+  {
+    // FAILS IN PARALLEL - FIX
+    //    Info << "Print contact area" << endl;
+    //volScalarField ca = contact.contactArea();
+    //ca.write();
+    
+    //-------------------------------------------------------------//
+    // I couldn't get tmp to return the pointScalarField correctly //
+    // so I had to make the pointScalarField here and pass it to   //
+    // contactGapPoints and pointContactForce to populate          //
+    //-------------------------------------------------------------//
+    //This is the point distance for each contact vertex
+    pointScalarField cGapPoints 
+      (
+       IOobject
+       (
+	"pointContactGap",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedScalar("scalar", dimLength, 0.0),
+       "calculated"
+       );
+
+    contact.contactGapPoints(cGapPoints);
+    cGapPoints.write();
+	  
+    
+    //- This is the point distance for each contact vertex
+    pointVectorField cPointForce 
+      (
+       IOobject
+       (
+	"pointContactForce",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedVector("vector", dimForce, vector::zero),
+       "calculated"
+       );
+    contact.contactPointForce(cPointForce);
+    cPointForce.write();
+    
+    //- this is the actual (sigma&n)&n) on the contact patches
+    //- SHOULD THIS BE A REF TO A TMP...?
+    volScalarField cPressure = contact.contactPressure();
+    cPressure.write();
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/readContactControls.H

@@ -0,0 +1,16 @@
+//****************************************************//
+// Read the contact tol and max contact correctors
+//****************************************************//
+const IOdictionary& contactControl =
+  IOobject
+  (
+   "contactProperties",
+   U.time().constant(),
+   U.db(),
+   IOobject::MUST_READ,
+   IOobject::NO_WRITE
+   );
+
+
+//- The frequnecy the conatct is corrected inside the momentum loop
+int uEqnContactCorrFreq(readInt(contactControl.lookup("innerContactCorrectFreq")));

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/readDivDSigmaExpMethod.H

@@ -0,0 +1,9 @@
+//- how explicit component of sigma is to be calculated
+word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp"));
+Info << divDSigmaExpMethod << " method chosen for calculation of DSigmaExp" << endl;
+if(divDSigmaExpMethod != "standard" && divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose" && divDSigmaExpMethod != "laplacian")
+  {
+    FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
+	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	       << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/readDivDSigmaLargeStrainMethod.H

@@ -0,0 +1,10 @@
+//-  the method used to calculate the explicit component of sigma
+word divDSigmaLargeStrainExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaLargeStrainExp"));
+Info << "divDSigmaLargeStrainExp calculation method: " << divDSigmaLargeStrainExpMethod << endl;
+if(divDSigmaLargeStrainExpMethod != "standard"
+   && divDSigmaLargeStrainExpMethod != "surface")
+  {
+    FatalError << "divDSigmaLargeStrainExp method " << divDSigmaLargeStrainExpMethod << " not found!" << nl
+               << "valid methods are:\nstandard\nsurface"
+               << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/readMoveMeshMethod.H

@@ -0,0 +1,10 @@
+//-  the method used to interpolate in cell centre displacements to the vertices to move the mesh
+word moveMeshMethod(mesh.solutionDict().subDict("stressedFoam").lookup("moveMeshMethod"));
+Info << "Move Mesh method: " << moveMeshMethod << endl;
+if(moveMeshMethod != "inverseDistance"
+   && moveMeshMethod != "leastSquares")
+  {
+    FatalError << "moveMeshMethod " << moveMeshMethod << " not found!" << nl
+               << "valid methods are:\ninvreseDistance\nleastSquares"
+               << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/readStressedFoamControls.H

@@ -0,0 +1,5 @@
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("stressedFoam");
+
+    int nCorr(readInt(stressControl.lookup("nCorrectors")));
+    scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/rotateFields.H

@@ -0,0 +1,24 @@
+//--------------------------------------------------//
+//- rotate fields
+//--------------------------------------------------//
+{
+  Info << "Rotating fields" << endl;
+
+  F = I + DF;
+  
+  U += DU;
+  
+  epsilon += DEpsilon;
+
+  sigma += DSigma;
+
+  volTensorField Finv = inv(F);
+
+  J = det(F);
+  
+  rho = rho/J;
+ 
+  epsilon = symm(Finv.T() & epsilon & Finv);
+
+  sigma  = 1/J * symm(F & sigma & F.T());
+}

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/writeBoundaryNetForces.H

@@ -0,0 +1,14 @@
+// * * * * * * * * * * * * * * * * NET FORCES * * * * * * * * * * * * * * *  //
+
+vectorField netForces(mesh.boundary().size(), vector::zero);
+
+Info << nl;
+forAll(netForces, patchI)
+{
+  netForces[patchI] = gSum(mesh.Sf().boundaryField()[patchI] & sigma.boundaryField()[patchI]);
+  
+  Info << "patch " << mesh.boundary()[patchI].name() << " net force is "
+       << netForces[patchI] << " N" << endl;
+}   
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/solidMechanics/elasticContactNonLinULSolidFoam/writeFields.H

@@ -0,0 +1,36 @@
+if (runTime.outputTime())
+  {
+    volScalarField epsilonEq
+      (
+       IOobject
+       (
+	"epsilonEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+       );
+
+    Info<< "Max epsilonEq = " << max(epsilonEq).value()
+	<< endl;
+
+    volScalarField sigmaEq
+      (
+       IOobject
+       (
+	"sigmaEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((3.0/2.0)*magSqr(dev(sigma)))
+       );
+    
+    Info<< "Max sigmaEq = " << max(sigmaEq).value()
+	<< endl;
+    
+    runTime.write();
+  }

applications/solvers/solidMechanics/elasticContactSolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticContactSolidFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude \
+    -I$(FOAM_SRC)/VectorN/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian \
+    -lVectorN

applications/solvers/solidMechanics/elasticContactSolidFoam/Make/options_myLib

@@ -0,0 +1,15 @@
+EXE_INC = \
+//    -I../../../myLibraries/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I../materialModels/lnInclude \
+    -I./solidInterface
+//    -I./patchToPatchInterpolation
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+//    -lfiniteVolume_philipc \
+//    -lfiniteVolume \
+    -lmaterialModels_philipc \
+    -llduSolvers \
+    -lmeshTools

applications/solvers/solidMechanics/elasticContactSolidFoam/Make/options_old

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/dynamicMesh/lnInclude \
+    -I../materialModels/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+    -lfiniteVolume \
+    -lmaterialModels_philipc

applications/solvers/solidMechanics/elasticContactSolidFoam/calculateDivSigmaExp.H

@@ -0,0 +1,47 @@
+if(divSigmaExpMethod == "standard")
+  {
+    divSigmaExp = fvc::div
+      (
+       mu*gradU.T() + lambda*(I*tr(gradU)) - (mu + lambda)*gradU,
+       "div(sigma)"
+       );
+  }
+ else if(divSigmaExpMethod == "surface")
+   { 
+     divSigmaExp = fvc::div
+       (
+	muf*(mesh.Sf() & fvc::interpolate(gradU.T()))
+	+ lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradU)))
+	- (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradU))
+	);
+   }
+ else if(divSigmaExpMethod == "decompose")
+   {
+     surfaceTensorField shearGradU =
+       ((I - n*n)&fvc::interpolate(gradU));
+     
+     divSigmaExp = fvc::div
+       (
+	mesh.magSf()
+	*(
+	  - (muf + lambdaf)*(snGradU&(I - n*n))
+	  + lambdaf*tr(shearGradU&(I - n*n))*n
+	  + muf*(shearGradU&n)
+	  )
+	);
+   }
+ else if(divSigmaExpMethod == "expLaplacian")
+   {
+     divSigmaExp =
+       - fvc::laplacian(mu + lambda, U, "laplacian(U,U)")
+       + fvc::div
+       (
+	mu*gradU.T()
+	+ lambda*(I*tr(gradU)),
+	"div(sigma)"
+	);
+   }
+ else
+   {
+     FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
+   }

applications/solvers/solidMechanics/elasticContactSolidFoam/calculateEpsilonSigma.H

@@ -0,0 +1,3 @@
+epsilon = symm(gradU);
+
+sigma = 2*mu*epsilon + lambda*(I*tr(epsilon));

applications/solvers/solidMechanics/elasticContactSolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magU = mag(U.internalField() - U.oldTime().internalField());
+
+    forAll(magU, cellI)
+    {
+        if (magU[cellI] < SMALL)
+        {
+            magU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                U.internalField()
+              - U.prevIter().internalField()
+            )
+           /magU
+        );
+}

applications/solvers/solidMechanics/elasticContactSolidFoam/correctGlobalFaceZoneMesh.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+correctGlobalFaceZoneMesh.H
+
+When there is a globalFaceZone and the mesh is moved by interpolating U to the
+vertices with volPointInterpolation, then there are two problems:
+    -some points on the patch with the faceZone are moved incorrectly, I think
+     it is because the faceZone has no U and causes an incorrect interpolation,
+    -the faceZones points not on the proc cells are not moved at all because
+     they have no U.
+
+So the points on the patch with the faceZone need to be fixed and also all the
+faceZone points need to be moved and synchronised so each proc has the same
+full faceZone mesh.
+
+The mapping of procs faceZone order of points to the master procs faceZone point
+order is kept in procToGlobalFZmap, which is calculated at the start of the run
+in the createGlobalToLocalFaceZonePointMap.H header file.
+
+Note: DU is used for updated Lagrangian solver instead of U
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
+    contactPatchPairList& contacts = contact;
+    
+    forAll(contacts, contactI)
+      {
+	label masterID = contacts[contactI].masterPatch().index();
+	label slaveID = contacts[contactI].slavePatch().index();
+	
+	primitivePatchInterpolation masterInterpolator
+	  (mesh.boundaryMesh()[masterID]);
+	primitivePatchInterpolation slaveInterpolator
+	  (mesh.boundaryMesh()[slaveID]);
+
+	//- U must be interpolated to the vertices, this ignores the faceZone
+	//- points with no U (unlike volPointInterpolation)
+	vectorField correctMasterPointU =
+	  masterInterpolator.faceToPointInterpolate<vector>
+	  (
+	   U.boundaryField()[masterID]
+	   );
+	vectorField correctSlavePointU =
+	slaveInterpolator.faceToPointInterpolate<vector>
+	  (
+	   U.boundaryField()[slaveID]
+	   );
+	
+	vectorField oldMasterPoints =
+	  mesh.boundaryMesh()[masterID].localPoints();
+	vectorField oldSlavePoints =
+	  mesh.boundaryMesh()[slaveID].localPoints();
+	
+	labelList masterPointLabels =
+	  mesh.boundaryMesh()[masterID].meshPoints();
+	labelList slavePointLabels =
+	  mesh.boundaryMesh()[slaveID].meshPoints();
+	
+	//- correct the patch newPoints
+	forAll(masterPointLabels, pointI)
+	  {
+	    label pointGlobalLabel = masterPointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldMasterPoints[pointI]
+	      +
+	      correctMasterPointU[pointI];
+	  }
+	forAll(slavePointLabels, pointI)
+	  {
+	    label pointGlobalLabel = slavePointLabels[pointI];
+	    newPoints[pointGlobalLabel] =
+	      oldSlavePoints[pointI]
+	      +
+	      correctSlavePointU[pointI];
+	  }
+      }
+
+    
+
+
+    //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
+
+    forAll(mesh.faceZones(), faceZoneI)
+      {
+	//- find the patch corresponding to this faceZone
+	//- assuming that the FZ is called <patch_name>FaceZone
+	string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	//- remove the string FaceZone from the end of the face zone name to get the patch name
+	string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	if(patchID == -1)
+	  {
+	    FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+		       << faceZoneName << exit(FatalError);
+	  }
+	    
+	vectorField globalFZpoints =
+	  mesh.faceZones()[faceZoneI]().localPoints();
+	
+	//- new points for the face zone
+	vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	//- inter-proc points are shared by multiple procs
+	//- pointNumProc is the number of procs which a point lies on
+	scalarField pointNumProcs(globalFZpoints.size(), 0.0);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
+	    if(pointOnLocalProcPatch[faceZoneI][localPoint])
+	      {
+		label procPoint =
+		  mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
+		globalFZnewPoints[globalPointI] =
+		  newPoints[procPoint];
+		pointNumProcs[globalPointI] = 1;
+	      }		
+	  }
+	
+	reduce(globalFZnewPoints, sumOp<vectorField>());
+	reduce(pointNumProcs, sumOp<scalarField>());
+	
+	//- now average the newPoints between all procs
+	if(min(pointNumProcs) < 1)
+	  {
+	    FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
+	  }
+	globalFZnewPoints /= pointNumProcs;
+	
+	//- the globalFZnewPoints now contains the correct FZ new points in
+	//- a global order, now convert them back into the local proc order
+	
+	vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
+	
+	forAll(globalFZnewPoints, globalPointI)
+	  {
+	    label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
+	    
+	    procFZnewPoints[localPoint] =
+	      globalFZnewPoints[globalPointI];
+	  }
+	
+	//- now fix the newPoints points on the globalFaceZones
+	labelList procFZmeshPoints =
+	  mesh.faceZones()[faceZoneI]().meshPoints();
+	
+	forAll(procFZmeshPoints, pointI)
+	  {
+	    label procPoint = procFZmeshPoints[pointI];
+	    newPoints[procPoint] =
+	      procFZnewPoints[pointI];
+	  }
+      }
+  }

applications/solvers/solidMechanics/elasticContactSolidFoam/createFields.H

@@ -0,0 +1,91 @@
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volTensorField gradU = fvc::grad(U);
+    surfaceVectorField snGradU = fvc::snGrad(U);
+
+    volVectorField V
+    (
+        IOobject
+        (
+            "V",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fvc::ddt(U)
+    );
+
+   volTensorField gradV = fvc::ddt(gradU);
+   surfaceVectorField snGradV = (snGradU - snGradU.oldTime())/runTime.deltaT();
+
+    volSymmTensorField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+      volVectorField divSigmaExp
+      (
+        IOobject
+        (
+	 "divSigmaExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+
+   // read rheology properties
+    rheologyModel rheology(sigma);
+
+    volScalarField rho = rheology.rho();
+
+    volScalarField mu = rheology.mu();
+    volScalarField lambda = rheology.lambda();
+    surfaceScalarField muf = fvc::interpolate(mu, "mu");
+    surfaceScalarField lambdaf = fvc::interpolate(lambda, "lambda");
+
+    surfaceVectorField n = mesh.Sf()/mesh.magSf();
+
+
+   //- create contact problem
+   contactProblem contact(U);

applications/solvers/solidMechanics/elasticContactSolidFoam/createGlobalToLocalFaceZonePointMap.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+createGlobalToLocalFaceZonePointMap.H
+
+I could not figure out the order of the points on a globalFaceZone, they is a
+different order for each processor where the processors own points come first.
+
+So I decide that the master proc order is the global order and I find the map
+from this order to each proc order here by just commparing actual point
+coordinates.
+
+This map is then used to correct and synchronise the globalFaceZone points on
+all the procs when the mesh is moved, in header file correctGlobalFaceZoneMesh.H.
+
+philipc
+\*---------------------------------------------------------------------------*/
+
+//- procToGlobalFZmap is a map from the current proc faceZone point order to the
+//- master proc point order
+//- these are read if present to allow restarting of contact cases
+IOList<labelList> procToGlobalFZmap
+(
+ IOobject
+ (
+  "procToGlobalFZmap",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+IOList<labelList> pointOnLocalProcPatch
+(
+ IOobject
+ (
+  "pointOnLocalProcPatch",
+  runTime.timeName(),
+  mesh,
+  IOobject::READ_IF_PRESENT,
+  IOobject::AUTO_WRITE
+  ),
+ mesh.faceZones().size()
+ );
+
+//- if they have been read then don't recalculate it
+bool globalFaceZoneMappingSet = false;
+if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
+  {
+    Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
+	 << endl; 
+    globalFaceZoneMappingSet = true;
+  }
+ else
+   {
+     Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl
+	  << "this message should only appear starting a new analysis" << endl;
+   }
+
+//- this is only needed in a parallel runs
+if(Pstream::parRun())
+  {
+    if(!globalFaceZoneMappingSet)
+      {
+	forAll(mesh.faceZones(), faceZoneI)
+	  {
+	    vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
+	
+	    //- set all slave points to zero because only the master order is used
+	    if(!Pstream::master())
+	      globalFZpoints *= 0.0;
+	    
+	    //- pass points to all procs
+	    reduce(globalFZpoints, sumOp<vectorField>());
+	    
+	    
+	    //- now every proc has the master's list of FZ points
+	    //- every proc must now find the mapping from their local FZpoints to
+	    //- the globalFZpoints
+	    
+	    vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
+	    
+	    forAll(globalFZpoints, globalPointI)
+	      {
+		forAll(procFZpoints, procPointI)
+		  {
+		    if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
+		      {
+			procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
+			break;
+		      }
+		  }
+	      }
+	    //- procToGlobalFZmap now contains the local FZpoint label for each
+	    //- global FZ point label - for each faceZone
+	    
+	    //- check what points are on the current proc patch
+	    pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
+	    
+	    //- find corresponding patch
+	    string faceZoneName = mesh.faceZones().names()[faceZoneI];
+	    //- remove the string FaceZone from the end of the face zone name to get the patch name
+	    string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
+	    label patchID = mesh.boundaryMesh().findPatchID(patchName);
+	    if(patchID == -1)
+	      {
+		FatalError << "Patch " << patchName << " not found corresponding for faceZone"
+			   << faceZoneName << exit(FatalError);
+	      }
+	    
+	    forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
+	      {
+		forAll(mesh.boundaryMesh()[patchID].localPoints(), pi)
+		  {
+		    if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi])
+		      {
+			pointOnLocalProcPatch[faceZoneI][fzpi] = 1;
+			break;
+		      }
+		  }
+	      }
+	  }
+      } //- end if(!globalFaceZoneMappingSet)
+  }
+
+//- write to disk to allow restart of cases
+//- because it is not possible to calculate the 
+//- mapping after the meshes have moved
+if(!globalFaceZoneMappingSet)
+  {
+    procToGlobalFZmap.write();
+    pointOnLocalProcPatch.write();
+  }

applications/solvers/solidMechanics/elasticContactSolidFoam/elasticContactSolidFoam.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    elasticContactSolidFoam
+
+Description
+    Transient/steady-state segregated finite-volume solver for small strain
+    elastic solid bodies in contact, using an total strain total Lagrangian
+    approach.
+
+    Works in parallel but mesh.movePoints sometimes fails for some unknown
+    reason depending on the decomposition.
+
+    Solves for the displacement increment vector field DU, also generating the
+    stress tensor field sigma.
+
+    It is only for frictionless contact, friction not implemented yet.
+
+Author
+    Philip Cardiff
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "rheologyModel.H"
+#include "contactProblem.H"
+
+#include "volPointInterpolation.H"
+#include "pointPatchInterpolation.H"
+#include "primitivePatchInterpolation.H"
+#include "fixedValuePointPatchFields.H"
+#include "pointFields.H"
+#include "pointMesh.H"
+#include "pointBoundaryMesh.H"
+#include "primitivePatchInterpolation.H"
+#include "twoDPointCorrector.H"
+#include "plane.H"
+#include "meshSearch.H"
+//#include "leastSquaresVolPointInterpolation.H"
+
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+
+# include "createTime.H"
+
+# include "createMesh.H"
+
+# include "createFields.H"
+
+# include "readDivSigmaExpMethod.H"
+
+# include "createGlobalToLocalFaceZonePointMap.H"
+
+  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+  
+  Info<< "\nStarting time loop\n" << endl;
+  
+  for (runTime++; !runTime.end(); runTime++)
+    {
+      Info<< "Time: " << runTime.timeName() << endl;
+      
+#     include "readContactControls.H"
+
+#     include "readStressedFoamControls.H"
+  
+      //-- for moving the mesh and then back again
+      vectorField oldMeshPoints = mesh.allPoints();
+      
+      int iCorr = 0;
+      lduMatrix::solverPerformance solverPerf;
+      word solverName;
+      lduMatrix::debug = 0;
+      scalar residual = GREAT;
+      scalar initialResidual = 0;
+      scalar relativeResidual = GREAT;
+
+      //- Predictor step
+      if (predictor)
+        {
+	  Info << "\nPredicting U, gradU and snGradU based on V, gradV and snGradV\n" << endl;
+	  U += V*runTime.deltaT();
+	  gradU += gradV*runTime.deltaT();
+	  snGradU += snGradV*runTime.deltaT();
+        }
+
+      do //- start of momentum loop
+	{
+	  U.storePrevIter();
+	  
+	  //- correct the contact boundaries 
+          if(iCorr % uEqnContactCorrFreq == 0)
+	    {
+	      Info << "\t\tCorrecting contact in the momentum loop "
+		   << "iteration: " << iCorr
+		   << ", residual: " << residual
+		   << endl;
+	      //#                 include "moveMeshLeastSquares.H"
+#             include "moveSolidMesh.H"
+	      contact.correct();
+	      mesh.movePoints(oldMeshPoints);
+	    }
+  
+#         include "calculateDivSigmaExp.H"
+
+	  fvVectorMatrix UEqn
+	    (
+	     fvm::d2dt2(rho, U)
+	     ==
+	     fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)")
+	     + divSigmaExp
+	     );
+	  
+	  solverPerf = UEqn.solve();
+	  
+	  U.relax();
+	  
+	  solverName = solverPerf.solverName();
+	  
+	  gradU = fvc::grad(U);
+	  snGradU = fvc::snGrad(U);
+
+	  residual = solverPerf.initialResidual();
+	  
+	  if(iCorr == 0)
+	    {
+	      initialResidual = solverPerf.initialResidual();
+	    }
+	  
+#         include "calculateRelativeResidual.H"	      
+
+	  Info << "\tTime " << runTime.value()
+	       << ", Corrector " << iCorr
+	       << ", Solving for " << U.name()
+	       << " using " << solverPerf.solverName()
+	       << ", residual = " << solverPerf.initialResidual()
+	       << ", relative residual = " << relativeResidual << endl;
+	} //- end of momentum loop
+      while
+	(
+	 //relativeResidual > convergenceTolerance
+	 residual > convergenceTolerance
+	 &&
+	 ++iCorr < nCorr
+	 );
+      
+      // Print out info per contact iteration
+      Info << "\t\tSolving for " << U.name()
+	   << " using " << solverName
+	   << ", Initial residual = " << initialResidual
+	   << ", Final residual = " << solverPerf.initialResidual()
+	   << ", No outer iterations " << iCorr << endl;
+      
+      lduMatrix::debug = 1;
+      
+      V = fvc::ddt(U);
+      gradV = fvc::ddt(gradU);
+      snGradV = (snGradU - snGradU.oldTime())/runTime.deltaT();
+
+#     include "calculateEpsilonSigma.H"
+      
+#     include "writeFields.H"
+        
+      //#     include "moveMeshLeastSquares.H"
+      //#     include "moveSolidMesh.H"
+      //#     include "printContactResults.H"
+      //mesh.movePoints(oldMeshPoints);
+      
+      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+	  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+	  << endl << endl;
+    }
+  
+  Info<< "End\n" << endl;
+  
+  return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/solidMechanics/elasticContactSolidFoam/moveMeshLeastSquares.H

@@ -0,0 +1,56 @@
+//--------------------------------------------------//
+//- move mesh
+//--------------------------------------------------//
+if(min(J.internalField()) > 0)
+  {
+    Info << "Moving mesh using least squares interpolation" << endl;
+
+    leastSquaresVolPointInterpolation pointInterpolation(mesh);
+  
+    // Create point mesh
+    pointMesh pMesh(mesh);
+    
+    wordList types
+      (
+       pMesh.boundary().size(),
+       calculatedFvPatchVectorField::typeName
+       );
+    
+    pointVectorField pointDU
+      (
+       IOobject
+       (
+        "pointDU",
+        runTime.timeName(),
+        mesh
+        ),
+       pMesh,
+       dimensionedVector("zero", dimLength, vector::zero),
+       types
+       );
+    
+    pointInterpolation.interpolate(DU, pointDU);
+
+    const vectorField& pointDUI =
+      pointDU.internalField();
+
+    //- Move mesh
+    vectorField newPoints = mesh.allPoints();
+
+    forAll (pointDUI, pointI)
+      {
+        newPoints[pointI] += pointDUI[pointI];
+      }
+
+    twoDPointCorrector twoDCorrector(mesh);
+    twoDCorrector.correctPoints(newPoints);
+    mesh.movePoints(newPoints);
+    mesh.V00();
+    mesh.moving(false);
+  }
+  else
+    {
+      FatalErrorIn(args.executable())
+	<< "Negative Jacobian"
+	<< exit(FatalError);
+    }

applications/solvers/solidMechanics/elasticContactSolidFoam/moveSolidMesh.H

@@ -0,0 +1,30 @@
+{  
+  // Create point interpolation
+  volPointInterpolation pointInterpolation(mesh);
+
+
+  // Calculate mesh points displacement
+  pointVectorField pointU = pointInterpolation.interpolate(U);
+
+  //  const vectorField& pointUI = pointU.internalField();
+  vectorField pointUI = pointU.internalField();
+  
+
+  // Move mesh
+  vectorField newPoints = mesh.allPoints();
+  
+  forAll (pointUI, pointI)
+    {
+      newPoints[pointI] += pointUI[pointI];
+    }
+
+# include "correctGlobalFaceZoneMesh.H"
+ 
+  twoDPointCorrector twoDCorrector(mesh);
+  twoDCorrector.correctPoints(newPoints);  
+  
+  mesh.movePoints(newPoints);
+  //  pMesh.movePoints(newPoints);
+  mesh.V00();
+  mesh.moving(false);
+}

applications/solvers/solidMechanics/elasticContactSolidFoam/printContactResults.H

@@ -0,0 +1,55 @@
+if (runTime.outputTime())
+  {
+    // FAILS IN PARALLEL - FIX
+    //    Info << "Print contact area" << endl;
+    //volScalarField ca = contact.contactArea();
+    //ca.write();
+    
+    //-------------------------------------------------------------//
+    // I couldn't get tmp to return the pointScalarField correctly //
+    // so I had to make the pointScalarField here and pass it to   //
+    // contactGapPoints and pointContactForce to populate          //
+    //-------------------------------------------------------------//
+    //This is the point distance for each contact vertex
+    pointScalarField cGapPoints 
+      (
+       IOobject
+       (
+	"pointContactGap",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedScalar("scalar", dimLength, 0.0),
+       "calculated"
+       );
+
+    contact.contactGapPoints(cGapPoints);
+    cGapPoints.write();
+	  
+    
+    //- This is the point distance for each contact vertex
+    pointVectorField cPointForce 
+      (
+       IOobject
+       (
+	"pointContactForce",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       pMesh,
+       dimensionedVector("vector", dimForce, vector::zero),
+       "calculated"
+       );
+    contact.contactPointForce(cPointForce);
+    cPointForce.write();
+    
+    //- this is the actual (sigma&n)&n) on the contact patches
+    //- SHOULD THIS BE A REF TO A TMP...?
+    volScalarField cPressure = contact.contactPressure();
+    cPressure.write();
+  }

applications/solvers/solidMechanics/elasticContactSolidFoam/readContactControls.H

@@ -0,0 +1,16 @@
+//****************************************************//
+// Read the contact tol and max contact correctors
+//****************************************************//
+const IOdictionary& contactControl =
+  IOobject
+  (
+   "contactProperties",
+   U.time().constant(),
+   U.db(),
+   IOobject::MUST_READ,
+   IOobject::NO_WRITE
+   );
+
+
+//- The frequnecy the conatct is corrected inside the momentum loop
+int uEqnContactCorrFreq(readInt(contactControl.lookup("innerContactCorrectFreq")));

applications/solvers/solidMechanics/elasticContactSolidFoam/readDivSigmaExpMethod.H

@@ -0,0 +1,9 @@
+//- how explicit component of sigma is to be calculated
+word divSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divSigmaExp"));
+Info << divSigmaExpMethod << " method chosen for calculation of sigmaExp" << endl;
+if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian")
+  {
+    FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl
+	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	       << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticContactSolidFoam/readStressedFoamControls.H

@@ -0,0 +1,7 @@
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("stressedFoam");
+
+    int nCorr(readInt(stressControl.lookup("nCorrectors")));
+    scalar convergenceTolerance(readScalar(stressControl.lookup("U")));
+
+    Switch predictor(stressControl.lookup("predictor"));

applications/solvers/solidMechanics/elasticContactSolidFoam/writeBoundaryNetForces.H

@@ -0,0 +1,14 @@
+// * * * * * * * * * * * * * * * * NET FORCES * * * * * * * * * * * * * * *  //
+
+vectorField netForces(mesh.boundary().size(), vector::zero);
+
+Info << nl;
+forAll(netForces, patchI)
+{
+  netForces[patchI] = gSum(mesh.Sf().boundaryField()[patchI] & sigma.boundaryField()[patchI]);
+  
+  Info << "patch\t" << mesh.boundary()[patchI].name() << "\t\tnet force is\t"
+       << netForces[patchI] << " N" << endl;
+}   
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/solidMechanics/elasticContactSolidFoam/writeFields.H

@@ -0,0 +1,36 @@
+if (runTime.outputTime())
+  {
+    volScalarField epsilonEq
+      (
+       IOobject
+       (
+	"epsilonEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+       );
+
+    Info<< "Max epsilonEq = " << max(epsilonEq).value()
+	<< endl;
+
+    volScalarField sigmaEq
+      (
+       IOobject
+       (
+	"sigmaEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((3.0/2.0)*magSqr(dev(sigma)))
+       );
+    
+    Info<< "Max sigmaEq = " << max(sigmaEq).value()
+	<< endl;
+    
+    runTime.write();
+  }

applications/solvers/solidMechanics/elasticGravitySolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticGravitySolidFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian

applications/solvers/solidMechanics/elasticGravitySolidFoam/calculateDivSigmaExp.H

@@ -0,0 +1,47 @@
+if(divSigmaExpMethod == "standard")
+  {
+    divSigmaExp = fvc::div
+      (
+       mu*gradU.T() + lambda*(I*tr(gradU)) - (mu + lambda)*gradU,
+       "div(sigma)"
+       );
+  }
+ else if(divSigmaExpMethod == "surface")
+   { 
+     divSigmaExp = fvc::div
+       (
+	muf*(mesh.Sf() & fvc::interpolate(gradU.T()))
+	+ lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradU)))
+	- (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradU))
+	);
+   }
+ else if(divSigmaExpMethod == "decompose")
+   {
+     surfaceTensorField shearGradU =
+       ((I - n*n)&fvc::interpolate(gradU));
+     
+     divSigmaExp = fvc::div
+       (
+	mesh.magSf()
+	*(
+	  - (muf + lambdaf)*(fvc::snGrad(U)&(I - n*n))
+	  + lambdaf*tr(shearGradU&(I - n*n))*n
+	  + muf*(shearGradU&n)
+	  )
+	);
+   }
+ else if(divSigmaExpMethod == "expLaplacian")
+   {
+     divSigmaExp =
+       - fvc::laplacian(mu + lambda, U, "laplacian(DU,U)")
+       + fvc::div
+       (
+	mu*gradU.T()
+	+ lambda*(I*tr(gradU)),
+	"div(sigma)"
+	);
+   }
+ else
+   {
+     FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
+   }

applications/solvers/solidMechanics/elasticGravitySolidFoam/calculateEpsilonSigma.H

@@ -0,0 +1,3 @@
+epsilon = symm(gradU);
+
+sigma = 2*mu*epsilon + lambda*(I*tr(epsilon));

applications/solvers/solidMechanics/elasticGravitySolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magU = mag(U.internalField() - U.oldTime().internalField());
+
+    forAll(magU, cellI)
+    {
+        if (magU[cellI] < SMALL)
+        {
+            magU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                U.internalField()
+              - U.prevIter().internalField()
+            )
+           /magU
+        );
+}

applications/solvers/solidMechanics/elasticGravitySolidFoam/createFields.H

@@ -0,0 +1,68 @@
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volTensorField gradU = fvc::grad(U);
+
+    volSymmTensorField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+      volVectorField divSigmaExp
+      (
+        IOobject
+        (
+	 "divSigmaExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+    rheologyModel rheology(sigma);
+
+    volScalarField rho = rheology.rho();
+
+    volScalarField mu = rheology.mu();
+    volScalarField lambda = rheology.lambda();
+    surfaceScalarField muf = fvc::interpolate(mu, "mu");
+    surfaceScalarField lambdaf = fvc::interpolate(lambda, "lambda");
+
+    surfaceVectorField n = mesh.Sf()/mesh.magSf();

applications/solvers/solidMechanics/elasticGravitySolidFoam/createSolidInterface.H

@@ -0,0 +1,24 @@
+Switch solidInterfaceCorr(false);
+
+solidInterface* solidInterfacePtr(NULL);
+
+{
+    const dictionary& stressControl =
+      mesh.solutionDict().subDict("stressedFoam");
+
+    solidInterfaceCorr = Switch(stressControl.lookup("solidInterface"));
+
+    if(solidInterfaceCorr)
+      {
+	Info << "Creating solid interface correction" << endl;
+	solidInterfacePtr = new solidInterface(mesh, rheology);
+	solidInterfacePtr->modifyProperties(muf, lambdaf);
+
+	//- solidInterface needs muf and lambdaf to be used for divSigmaExp
+	if(divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose")
+	  {
+	    FatalError << "divSigmaExp must be decompose or surface when solidInterface is on"
+		       << exit(FatalError);
+	  }
+      }
+}

applications/solvers/solidMechanics/elasticGravitySolidFoam/elasticGravitySolidFoam.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    elasticGravitySolidFoam
+
+Description
+    Transient/steady-state segregated finite-volume solver for small strain
+    elastic solid bodies.
+
+    Displacement field U is solved for using a total Lagrangian approach,
+    also generating the strain tensor field epsilon and stress tensor
+    field sigma.
+
+    With optional multi-material solid interface correction ensuring
+    correct tractions on multi-material interfaces.
+
+    Graivty added as a body force.
+
+Author
+    Philip Cardiff
+    multi-material by Tukovic et al. 2012
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "rheologyModel.H"
+#include "solidInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+
+# include "createTime.H"
+
+# include "createMesh.H"
+
+# include "createFields.H"
+
+# include "readDivSigmaExpMethod.H"
+
+# include "createSolidInterface.H"
+
+# include "readGravity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+  Info<< "\nCalculating displacement field\n" << endl;
+
+  while(runTime.loop())
+    {
+      Info<< "Time: " << runTime.timeName() << nl << endl;
+      
+#     include "readStressedFoamControls.H"
+
+      int iCorr = 0;
+      scalar initialResidual = 0;
+      lduMatrix::solverPerformance solverPerf;
+      scalar relativeResidual = GREAT;
+
+      lduMatrix::debug=0;
+
+      do
+        {
+	  U.storePrevIter();
+	  
+#         include "calculateDivSigmaExp.H"
+
+	  //- linear momentum equation
+	  fvVectorMatrix UEqn
+            (
+	     fvm::d2dt2(rho, U)
+	     ==
+	     fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+	     + divSigmaExp
+	     + rho*gravity
+	     );
+
+	  if(solidInterfaceCorr)
+	    {
+	      solidInterfacePtr->correct(UEqn);
+	    }
+
+	  solverPerf = UEqn.solve();
+
+	  if(iCorr == 0)
+	    {
+	      initialResidual = solverPerf.initialResidual();
+	    }
+	  
+	  U.relax();
+
+	  if(solidInterfaceCorr)
+	    {
+	      gradU = solidInterfacePtr->grad(U);
+	    }
+	  else
+	    {
+	      gradU = fvc::grad(U);
+	    }
+
+#         include "calculateRelativeResidual.H"
+	  
+	  Info << "\tTime " << runTime.value()
+	       << ", Corrector " << iCorr
+	       << ", Solving for " << U.name()
+	       << " using " << solverPerf.solverName()
+	       << ", residual = " << solverPerf.initialResidual()
+	       << ", relative residual = " << relativeResidual << endl;
+        }
+	while
+	  (
+	   //solverPerf.initialResidual() > convergenceTolerance
+	   relativeResidual > convergenceTolerance
+	   &&
+	   ++iCorr < nCorr
+	   );
+	
+      Info << nl << "Time " << runTime.value() << ", Solving for " << U.name() 
+	   << ", Initial residual = " << initialResidual 
+	   << ", Final residual = " << solverPerf.initialResidual()
+	   << ", Relative residual = " << relativeResidual
+	   << ", No outer iterations " << iCorr
+	   << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+	   << "  ClockTime = " << runTime.elapsedClockTime() << " s" 
+	   << endl;
+      
+      lduMatrix::debug=0;
+
+#     include "calculateEpsilonSigma.H"
+
+#     include "writeFields.H"
+
+      Info<< "ExecutionTime = "
+	  << runTime.elapsedCpuTime()
+	  << " s\n\n" << endl;
+    }
+  
+  Info<< "End\n" << endl;
+  
+  return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/solidMechanics/elasticGravitySolidFoam/readDivSigmaExpMethod.H

@@ -0,0 +1,9 @@
+//- how explicit component of sigma is to be calculated
+word divSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divSigmaExp"));
+Info << "Selecting divSigmaExp calculation method " << divSigmaExpMethod <<  endl;
+if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian")
+  {
+    FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl
+	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	       << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticGravitySolidFoam/readGravity.H

@@ -0,0 +1 @@
+dimensionedVector gravity(mesh.solutionDict().subDict("stressedFoam").lookup("gravity"));

applications/solvers/solidMechanics/elasticGravitySolidFoam/readStressedFoamControls.H

@@ -0,0 +1,5 @@
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("stressedFoam");
+
+    int nCorr(readInt(stressControl.lookup("nCorrectors")));
+    scalar convergenceTolerance(readScalar(stressControl.lookup("U")));

applications/solvers/solidMechanics/elasticGravitySolidFoam/writeFields.H

@@ -0,0 +1,36 @@
+if (runTime.outputTime())
+  {
+    volScalarField epsilonEq
+      (
+       IOobject
+       (
+	"epsilonEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+       );
+
+    Info<< "Max epsilonEq = " << max(epsilonEq).value()
+	<< endl;
+
+    volScalarField sigmaEq
+      (
+       IOobject
+       (
+	"sigmaEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((3.0/2.0)*magSqr(dev(sigma)))
+       );
+    
+    Info<< "Max sigmaEq = " << max(sigmaEq).value()
+	<< endl;
+    
+    runTime.write();
+  }

applications/solvers/solidMechanics/elasticIncrSolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticIncrSolidFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude \
+    -I$(FOAM_SRC)/VectorN/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian \
+    -lVectorN

applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateDEpsilonDSigma.H

@@ -0,0 +1,3 @@
+DEpsilon = symm(gradDU);
+
+DSigma = 2*mu*DEpsilon + lambda*(I*tr(DEpsilon));

applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateDivDSigmaExp.H

@@ -0,0 +1,47 @@
+if(divDSigmaExpMethod == "standard")
+  {
+    divDSigmaExp = fvc::div
+      (
+       mu*gradDU.T() + lambda*(I*tr(gradDU)) - (mu + lambda)*gradDU,
+       "div(sigma)"
+       );
+  }
+ else if(divDSigmaExpMethod == "surface")
+   { 
+     divDSigmaExp = fvc::div
+       (
+	muf*(mesh.Sf() & fvc::interpolate(gradDU.T()))
+	+ lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradDU)))
+	- (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradDU))
+	);
+   }
+ else if(divDSigmaExpMethod == "decompose")
+   {
+     surfaceTensorField shearGradDU =
+       ((I - n*n)&fvc::interpolate(gradDU));
+     
+     divDSigmaExp = fvc::div
+       (
+	mesh.magSf()
+	*(
+	  - (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
+	  + lambdaf*tr(shearGradDU&(I - n*n))*n
+	  + muf*(shearGradDU&n)
+	  )
+	);
+   }
+ else if(divDSigmaExpMethod == "laplacian")
+   {
+     divDSigmaExp =
+       - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+       + fvc::div
+       (
+	mu*gradDU.T()
+	+ lambda*(I*tr(gradDU)),
+	"div(sigma)"
+	);
+   }
+ else
+   {
+     FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
+   }

applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateNetForces.H

@@ -0,0 +1,12 @@
+//- net forces on boundary patches
+
+vectorField netForces(mesh.boundary().size(), vector::zero);
+
+Info << nl;
+forAll(netForces, patchI)
+{
+  netForces[patchI] = gSum(mesh.Sf().boundaryField()[patchI] & sigma.boundaryField()[patchI]);
+  
+  Info << "patch\t" << mesh.boundary()[patchI].name() << "\t\tnet force is\t"
+       << netForces[patchI] << " N" << endl;
+}   

applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magDU = mag(DU.internalField());
+
+    forAll(magDU, cellI)
+    {
+        if (magDU[cellI] < SMALL)
+        {
+            magDU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                DU.internalField()
+              - DU.prevIter().internalField()
+            )
+           /magDU
+        );
+}

applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H

@@ -0,0 +1,109 @@
+    Info<< "Reading field DU\n" << endl;
+    volVectorField DU
+    (
+        IOobject
+        (
+            "DU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    volTensorField gradDU = fvc::grad(DU);
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimLength, vector::zero)	
+    );
+
+    volSymmTensorField DEpsilon
+    (
+        IOobject
+        (
+            "DEpsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField epsilon
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    );
+
+    volSymmTensorField DSigma
+    (
+        IOobject
+        (
+            "DSigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+      volVectorField divDSigmaExp
+      (
+        IOobject
+        (
+	 "divDSigmaExp",
+	 runTime.timeName(),
+	 mesh,
+	 IOobject::NO_READ,
+	 IOobject::NO_WRITE
+	 ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+       );
+
+    rheologyModel rheology(sigma);
+
+    volScalarField rho = rheology.rho();
+
+    volScalarField mu = rheology.mu();
+    volScalarField lambda = rheology.lambda();
+    surfaceScalarField muf = fvc::interpolate(mu, "mu");
+    surfaceScalarField lambdaf = fvc::interpolate(lambda, "lambda");
+
+    surfaceVectorField n = mesh.Sf()/mesh.magSf();

applications/solvers/solidMechanics/elasticIncrSolidFoam/createSolidInterface.H

@@ -0,0 +1,24 @@
+Switch solidInterfaceCorr(false);
+
+solidInterface* solidInterfacePtr(NULL);
+
+{
+    const dictionary& stressControl =
+      mesh.solutionDict().subDict("stressedFoam");
+
+    solidInterfaceCorr = Switch(stressControl.lookup("solidInterface"));
+
+    if(solidInterfaceCorr)
+      {
+	Info << "Creating solid interface correction" << endl;
+	solidInterfacePtr = new solidInterface(mesh, rheology);
+	solidInterfacePtr->modifyProperties(muf, lambdaf);
+
+	//- solidInterface needs muf and lambdaf to be used for divDSigmaExp
+	if(divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose")
+	  {
+	    FatalError << "divDSigmaExp must be decompose or surface when solidInterface is on"
+		       << exit(FatalError);
+	  }
+      }
+}

applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    elasticIncrSolidFoam
+
+Description
+    Transient/steady-state segregated finite-volume solver for small strain
+    elastic solid bodies.
+
+    Displacement Increment field DU is solved for using a incremental total
+    Lagrangian approach,
+    also generating the strain tensor field epsilon and stress tensor
+    field sigma.
+
+    With optional multi-material solid interface correction ensuring
+    correct tractions on multi-material interfaces
+
+Author
+    Philip Cardiff
+    multi-material by Tukovic et al. 2012
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "rheologyModel.H"
+#include "solidInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+
+# include "createTime.H"
+
+# include "createMesh.H"
+
+# include "createFields.H"
+
+# include "readDivDSigmaExpMethod.H"
+
+# include "createSolidInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+  Info<< "\nCalculating displacement field\n" << endl;
+
+  while(runTime.loop())
+    {
+      Info<< "Time: " << runTime.timeName() << nl << endl;
+      
+#     include "readStressedFoamControls.H"
+
+      int iCorr = 0;
+      scalar initialResidual = 0;
+      scalar relativeResidual = GREAT;
+      lduMatrix::solverPerformance solverPerf;
+      lduMatrix::debug = 0;
+
+      do
+        {
+	  DU.storePrevIter();
+
+#         include "calculateDivDSigmaExp.H"
+
+	  //- linear momentum equation
+	  fvVectorMatrix DUEqn
+            (
+	     fvm::d2dt2(rho, DU)
+	     ==
+	     fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+	     + divDSigmaExp
+	     );
+
+	  if(solidInterfaceCorr)
+	    {
+	      solidInterfacePtr->correct(DUEqn);
+	    }
+
+	  solverPerf = DUEqn.solve();
+
+	  if(iCorr == 0)
+	    {
+	      initialResidual = solverPerf.initialResidual();
+	    }
+	  
+	  DU.relax();
+
+	  if(solidInterfaceCorr)
+	    {
+	      gradDU = solidInterfacePtr->grad(DU);
+	    }
+	  else
+	    {
+	      gradDU = fvc::grad(DU);
+	    }
+
+#         include "calculateRelativeResidual.H"
+	  
+	  Info << "\tTime " << runTime.value()
+	       << ", Corrector " << iCorr
+	       << ", Solving for " << DU.name()
+	       << " using " << solverPerf.solverName()
+	       << ", residual = " << solverPerf.initialResidual()
+	       << ", relative residual = " << relativeResidual << endl;
+        }
+	while
+	  (
+	   //solverPerf.initialResidual() > convergenceTolerance
+	   relativeResidual > convergenceTolerance
+	   &&
+	   ++iCorr < nCorr
+	   );
+	
+      Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name() 
+	   << ", Initial residual = " << initialResidual 
+	   << ", Final residual = " << solverPerf.initialResidual()
+	   << ", Relative residual = " << relativeResidual
+	   << ", No outer iterations " << iCorr
+	   << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+	   << "  ClockTime = " << runTime.elapsedClockTime() << " s" 
+	   << endl;
+      
+      lduMatrix::debug=0;
+
+#     include "calculateDEpsilonDSigma.H"
+
+      U += DU;
+      sigma += DSigma;
+      epsilon += DEpsilon;
+
+#     include "writeFields.H"
+
+#     include "calculateNetForces.H"
+
+      Info<< "ExecutionTime = "
+	  << runTime.elapsedCpuTime()
+	  << " s\n\n" << endl;
+    }
+  
+  Info<< "End\n" << endl;
+  
+  return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H

@@ -0,0 +1,9 @@
+//- how explicit component of sigma is to be calculated
+word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp"));
+Info << "Selecting divDSigmaExp calculation method " << divDSigmaExpMethod <<  endl;
+if(divDSigmaExpMethod != "standard" && divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose" && divDSigmaExpMethod != "laplacian")
+  {
+    FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
+	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	       << exit(FatalError);
+  }

applications/solvers/solidMechanics/elasticIncrSolidFoam/readStressedFoamControls.H

@@ -0,0 +1,5 @@
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("stressedFoam");
+
+    int nCorr(readInt(stressControl.lookup("nCorrectors")));
+    scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));

applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H

@@ -0,0 +1,36 @@
+if (runTime.outputTime())
+  {
+    volScalarField epsilonEq
+      (
+       IOobject
+       (
+	"epsilonEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+       );
+
+    Info<< "Max epsilonEq = " << max(epsilonEq).value()
+	<< endl;
+
+    volScalarField sigmaEq
+      (
+       IOobject
+       (
+	"sigmaEq",
+	runTime.timeName(),
+	mesh,
+	IOobject::NO_READ,
+	IOobject::AUTO_WRITE
+	),
+       sqrt((3.0/2.0)*magSqr(dev(sigma)))
+       );
+    
+    Info<< "Max sigmaEq = " << max(sigmaEq).value()
+	<< endl;
+    
+    runTime.write();
+  }

applications/solvers/solidMechanics/elasticIncrSolidFoam/writeForceDisplacement.H

@@ -0,0 +1,48 @@
+//- write force displacement to file
+
+label leftPatchID = mesh.boundaryMesh().findPatchID("leftClamp");
+if(leftPatchID == -1)
+  {
+    FatalError << "Cannot find patch left for calculating force" << endl;
+  }
+
+//- calculate force in x direction on leftClamp patch
+scalar leftForce = gSum(
+			vector(1, 0, 0) &
+			(mesh.boundary()[leftPatchID].Sf() & sigma.boundaryField()[leftPatchID])
+			);
+
+//- patchIntegrate utility integrates it this way but this is worng because the sigma tensor should
+//- be dotted with the surface normal to give the actual traction/force
+//- you cannot just take the component of the sigma tensor
+//scalar leftForcePatchIntegrateMethod = gSum(
+//					    mesh.magSf().boundaryField()[leftPatchID]
+//					    *sigma.boundaryField()[leftPatchID].component(symmTensor::XY)
+//					    );
+
+vector gaugeU1 = vector::zero;
+vector gaugeU2 = vector::zero;
+if(gaugeFaceID1 != -1)
+  {
+    gaugeU1 = U.boundaryField()[gaugeFacePatchID1][gaugeFaceID1];
+  }
+if(gaugeFaceID2 != -1)
+  {
+    gaugeU2 = U.boundaryField()[gaugeFacePatchID2][gaugeFaceID2];
+  }
+
+//- reduce across procs
+reduce(gaugeU1, sumOp<vector>());
+reduce(gaugeU2, sumOp<vector>());
+
+Pout << "gaugeU1 is " << gaugeU1 << nl
+     << "gaugeU2 is " << gaugeU2 << endl;
+
+scalar gaugeDisp = mag(gaugeU1 - gaugeU2);
+
+//- write to file
+if(Pstream::master())
+  {
+    OFstream& forceDispFile = *filePtr;
+    forceDispFile << 1000*gaugeDisp << "\t" << -1*leftForce << endl;
+  }

applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/Make/files

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

applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I../solidModels/lnInclude \
+    -I$(FOAM_SRC)/VectorN/lnInclude
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) -lsolidModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -lmeshTools \
+    -llagrangian \
+    -lVectorN
+

applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/calculateEpsilonSigma.H

@@ -0,0 +1,5 @@
+//- Green finite strain tensor
+epsilon = symm(gradU) + 0.5*symm(gradU & gradU.T());
+
+//- second Piola-Kirchhoff stress tensor
+sigma = 2*mu*epsilon + lambda*(I*tr(epsilon));

applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/calculateRelativeResidual.H

@@ -0,0 +1,22 @@
+{
+    scalarField magU = mag(U.internalField() - U.oldTime().internalField());
+
+    forAll(magU, cellI)
+    {
+        if (magU[cellI] < SMALL)
+        {
+            magU[cellI] = SMALL;
+        }
+    }
+
+    relativeResidual =
+        gMax
+        (
+            mag
+            (
+                U.internalField()
+              - U.prevIter().internalField()
+            )
+           /magU
+        );
+}