--HG--
branch : bgschaid/minorAdditionsBranch
This commit is contained in:
Bernhard F.W. Gschaider 2013-07-19 00:36:27 +02:00
commit 72c243a96f
1416 changed files with 18017 additions and 18012 deletions

View file

@ -49,7 +49,7 @@ IF (NOT $ENV{CDASH_SUBMIT_LOCAL_HOST_ID} STREQUAL "")
# $CDASH_SUBMIT_LOCAL_HOST_ID # $CDASH_SUBMIT_LOCAL_HOST_ID
SET( SET(
SITENAME $ENV{CDASH_SUBMIT_LOCAL_HOST_ID} SITENAME $ENV{CDASH_SUBMIT_LOCAL_HOST_ID}
CACHE STRING "Name of the local site" CACHE STRING "Name of the local site"
) )
ELSE (NOT $ENV{CDASH_SUBMIT_LOCAL_HOST_ID} STREQUAL "") ELSE (NOT $ENV{CDASH_SUBMIT_LOCAL_HOST_ID} STREQUAL "")
# Grab the hostname FQN; will be used for the sitename # Grab the hostname FQN; will be used for the sitename
@ -128,7 +128,7 @@ if(GIT_FOUND)
if (GIT_BRANCH_NAME STREQUAL "") if (GIT_BRANCH_NAME STREQUAL "")
message("No git-branch. Mercurial?") message("No git-branch. Mercurial?")
EXEC_PROGRAM(hg EXEC_PROGRAM(hg
ARGS branch ARGS branch
OUTPUT_VARIABLE GIT_BRANCH_NAME OUTPUT_VARIABLE GIT_BRANCH_NAME
) )
message("Git branch (mercurial): ${GIT_BRANCH_NAME}") message("Git branch (mercurial): ${GIT_BRANCH_NAME}")

14
COPYING
View file

@ -1,12 +1,12 @@
GNU GENERAL PUBLIC LICENSE GNU GENERAL PUBLIC LICENSE
Version 2, June 1991 Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc. Copyright (C) 1989, 1991 Free Software Foundation, Inc.
51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed. of this license document, but changing it is not allowed.
Preamble Preamble
The licenses for most software are designed to take away your The licenses for most software are designed to take away your
freedom to share and change it. By contrast, the GNU General Public freedom to share and change it. By contrast, the GNU General Public
@ -56,7 +56,7 @@ patent must be licensed for everyone's free use or not licensed at all.
The precise terms and conditions for copying, distribution and The precise terms and conditions for copying, distribution and
modification follow. modification follow.
GNU GENERAL PUBLIC LICENSE GNU GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License applies to any program or other work which contains 0. This License applies to any program or other work which contains
@ -255,7 +255,7 @@ make exceptions for this. Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally. of promoting the sharing and reuse of software generally.
NO WARRANTY NO WARRANTY
11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
@ -277,9 +277,9 @@ YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES. POSSIBILITY OF SUCH DAMAGES.
END OF TERMS AND CONDITIONS END OF TERMS AND CONDITIONS
How to Apply These Terms to Your New Programs How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it possible use to the public, the best way to achieve this is to make it

8
README
View file

@ -59,7 +59,7 @@
Then update the environment variables by sourcing the $HOME/.bashrc file by Then update the environment variables by sourcing the $HOME/.bashrc file by
typing in the terminal: typing in the terminal:
+ . $HOME/.bashrc + . $HOME/.bashrc
2) OR, if running tcsh or csh, source the etc/cshrc file by adding the 2) OR, if running tcsh or csh, source the etc/cshrc file by adding the
following line to the end of your $HOME/.cshrc file: following line to the end of your $HOME/.cshrc file:
@ -69,7 +69,7 @@
Then update the environment variables by sourcing the $HOME/.cshrc file by Then update the environment variables by sourcing the $HOME/.cshrc file by
typing in the terminal: typing in the terminal:
+ source $HOME/.cshrc + source $HOME/.cshrc
*** Installation in alternative locations *** Installation in alternative locations
OpenFOAM may also be installed in alternative locations. However, the OpenFOAM may also be installed in alternative locations. However, the
@ -79,13 +79,13 @@
The environment variable 'FOAM_INST_DIR' can be used to find and source the The environment variable 'FOAM_INST_DIR' can be used to find and source the
appropriate resource file. Here is a bash/ksh/sh example: appropriate resource file. Here is a bash/ksh/sh example:
+ export FOAM_INST_DIR=/data/app/OpenFOAM + export FOAM_INST_DIR=/data/app/OpenFOAM
+ foamDotFile=$FOAM_INST_DIR/OpenFOAM-<VERSION>/etc/bashrc + foamDotFile=$FOAM_INST_DIR/OpenFOAM-<VERSION>/etc/bashrc
+ [ -f $foamDotFile ] && . $foamDotFile + [ -f $foamDotFile ] && . $foamDotFile
and a csh/tcsh example: and a csh/tcsh example:
+ setenv FOAM_INST_DIR /data/app/OpenFOAM + setenv FOAM_INST_DIR /data/app/OpenFOAM
+ set foamDotFile=$FOAM_INST_DIR/OpenFOAM-<VERSION>/etc/cshrc + set foamDotFile=$FOAM_INST_DIR/OpenFOAM-<VERSION>/etc/cshrc
+ if ( -f $foamDotFile ) source $foamDotFile + if ( -f $foamDotFile ) source $foamDotFile

View file

@ -3,9 +3,9 @@
// Creates the porosity field for MULES // Creates the porosity field for MULES
volScalarField porosity volScalarField porosity
( (
IOobject IOobject
( (
"porosity", "porosity",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
@ -26,7 +26,7 @@
forAll( cells, cellI ) forAll( cells, cellI )
{ {
porosity[cells[cellI]] = zonePorosity; porosity[cells[cellI]] = zonePorosity;
} }
} }

View file

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

View file

@ -23,138 +23,131 @@ philipc
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
//***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****// //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
contactPatchPairList& contacts = contact; contactPatchPairList& contacts = contact;
forAll(contacts, contactI) forAll(contacts, contactI)
{ {
label masterID = contacts[contactI].masterPatch().index(); label masterID = contacts[contactI].masterPatch().index();
label slaveID = contacts[contactI].slavePatch().index(); label slaveID = contacts[contactI].slavePatch().index();
primitivePatchInterpolation masterInterpolator primitivePatchInterpolation masterInterpolator
(mesh.boundaryMesh()[masterID]); (
primitivePatchInterpolation slaveInterpolator mesh.boundaryMesh()[masterID]
(mesh.boundaryMesh()[slaveID]); );
primitivePatchInterpolation slaveInterpolator
(
mesh.boundaryMesh()[slaveID]
);
//- U must be interpolated to the vertices, this ignores the faceZone //- U must be interpolated to the vertices, this ignores the faceZone
//- points with no U (unlike volPointInterpolation) //- points with no U (unlike volPointInterpolation)
vectorField correctMasterPointU = vectorField correctMasterPointU =
masterInterpolator.faceToPointInterpolate<vector> masterInterpolator.faceToPointInterpolate<vector>
( (
U.boundaryField()[masterID] U.boundaryField()[masterID]
); );
vectorField correctSlavePointU = vectorField correctSlavePointU =
slaveInterpolator.faceToPointInterpolate<vector> slaveInterpolator.faceToPointInterpolate<vector>
( (
U.boundaryField()[slaveID] U.boundaryField()[slaveID]
); );
vectorField oldMasterPoints = vectorField oldMasterPoints =
mesh.boundaryMesh()[masterID].localPoints(); mesh.boundaryMesh()[masterID].localPoints();
vectorField oldSlavePoints = vectorField oldSlavePoints =
mesh.boundaryMesh()[slaveID].localPoints(); 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];
}
}
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 *****// //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
//- find the patch corresponding to this faceZone //- find the patch corresponding to this faceZone
//- assuming that the FZ is called <patch_name>FaceZone //- assuming that the FZ is called <patch_name>FaceZone
string faceZoneName = mesh.faceZones().names()[faceZoneI]; string faceZoneName = mesh.faceZones().names()[faceZoneI];
//- remove the string FaceZone from the end of the face zone name to get the patch name //- 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)); string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
label patchID = mesh.boundaryMesh().findPatchID(patchName); label patchID = mesh.boundaryMesh().findPatchID(patchName);
if(patchID == -1) if(patchID == -1)
{ {
FatalError << "Patch " << patchName << " not found corresponding for faceZone" FatalError << "Patch " << patchName << " not found corresponding for faceZone"
<< faceZoneName << exit(FatalError); << faceZoneName << exit(FatalError);
} }
vectorField globalFZpoints = vectorField globalFZpoints =
mesh.faceZones()[faceZoneI]().localPoints(); mesh.faceZones()[faceZoneI]().localPoints();
//- new points for the face zone //- new points for the face zone
vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero); vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
//- inter-proc points are shared by multiple procs //- inter-proc points are shared by multiple procs
//- pointNumProc is the number of procs which a point lies on //- pointNumProc is the number of procs which a point lies on
scalarField pointNumProcs(globalFZpoints.size(), 0.0); scalarField pointNumProcs(globalFZpoints.size(), 0.0);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
//if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size()) //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
if(pointOnLocalProcPatch[faceZoneI][localPoint]) if(pointOnLocalProcPatch[faceZoneI][localPoint])
{ {
label procPoint = label procPoint =
mesh.faceZones()[faceZoneI]().meshPoints()[localPoint]; mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
globalFZnewPoints[globalPointI] = globalFZnewPoints[globalPointI] = newPoints[procPoint];
newPoints[procPoint]; pointNumProcs[globalPointI] = 1;
pointNumProcs[globalPointI] = 1; }
} }
}
reduce(globalFZnewPoints, sumOp<vectorField>()); reduce(globalFZnewPoints, sumOp<vectorField>());
reduce(pointNumProcs, sumOp<scalarField>()); reduce(pointNumProcs, sumOp<scalarField>());
//- now average the newPoints between all procs //- now average the newPoints between all procs
if(min(pointNumProcs) < 1) if(min(pointNumProcs) < 1)
{ {
FatalError << "pointNumProc has not been set for all points" << exit(FatalError); FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
} }
globalFZnewPoints /= pointNumProcs; globalFZnewPoints /= pointNumProcs;
//- the globalFZnewPoints now contains the correct FZ new points in //- the globalFZnewPoints now contains the correct FZ new points in
//- a global order, now convert them back into the local proc order //- a global order, now convert them back into the local proc order
vectorField procFZnewPoints(globalFZpoints.size(), vector::zero); vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
procFZnewPoints[localPoint] = globalFZnewPoints[globalPointI];
}
procFZnewPoints[localPoint] = //- now fix the newPoints points on the globalFaceZones
globalFZnewPoints[globalPointI]; labelList procFZmeshPoints = mesh.faceZones()[faceZoneI]().meshPoints();
}
//- now fix the newPoints points on the globalFaceZones forAll(procFZmeshPoints, pointI)
labelList procFZmeshPoints = {
mesh.faceZones()[faceZoneI]().meshPoints(); label procPoint = procFZmeshPoints[pointI];
newPoints[procPoint] = procFZnewPoints[pointI];
forAll(procFZmeshPoints, pointI) }
{ }
label procPoint = procFZmeshPoints[pointI]; }
newPoints[procPoint] =
procFZnewPoints[pointI];
}
}
}

View file

@ -25,7 +25,7 @@
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimLength, vector::zero) dimensionedVector("zero", dimLength, vector::zero)
); );
volSymmTensorField DEpsilon volSymmTensorField DEpsilon
@ -84,22 +84,22 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
// read rheology properties // read rheology properties
rheologyModel rheology(sigma); rheologyModel rheology(sigma);
volScalarField rho = rheology.rho(); volScalarField rho = rheology.rho();
@ -111,5 +111,5 @@
surfaceVectorField n = mesh.Sf()/mesh.magSf(); surfaceVectorField n = mesh.Sf()/mesh.magSf();
//- create contact problem //- create contact problem
contactProblem contact(DU); contactProblem contact(DU);

View file

@ -19,117 +19,120 @@ philipc
//- these are read if present to allow restarting of contact cases //- these are read if present to allow restarting of contact cases
IOList<labelList> procToGlobalFZmap IOList<labelList> procToGlobalFZmap
( (
IOobject IOobject
( (
"procToGlobalFZmap", "procToGlobalFZmap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
IOList<labelList> pointOnLocalProcPatch IOList<labelList> pointOnLocalProcPatch
( (
IOobject IOobject
( (
"pointOnLocalProcPatch", "pointOnLocalProcPatch",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
//- if they have been read then don't recalculate it //- if they have been read then don't recalculate it
bool globalFaceZoneMappingSet = false; bool globalFaceZoneMappingSet = false;
if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0) if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
{ {
Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases" Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
<< endl; << endl;
globalFaceZoneMappingSet = true; globalFaceZoneMappingSet = true;
} }
else else
{ {
Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl 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 message should only appear starting a new analysis" << endl;
} }
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
if(!globalFaceZoneMappingSet) if(!globalFaceZoneMappingSet)
{ {
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0); procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
//- set all slave points to zero because only the master order is used //- set all slave points to zero because only the master order is used
if(!Pstream::master()) if(!Pstream::master())
globalFZpoints *= 0.0; {
globalFZpoints *= 0.0;
}
//- pass points to all procs //- pass points to all procs
reduce(globalFZpoints, sumOp<vectorField>()); reduce(globalFZpoints, sumOp<vectorField>());
//- now every proc has the master's list of FZ points //- now every proc has the master's list of FZ points
//- every proc must now find the mapping from their local FZpoints to //- every proc must now find the mapping from their local FZpoints to
//- the globalFZpoints //- the globalFZpoints
vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
forAll(globalFZpoints, globalPointI) forAll(globalFZpoints, globalPointI)
{ {
forAll(procFZpoints, procPointI) forAll(procFZpoints, procPointI)
{ {
if(procFZpoints[procPointI] == globalFZpoints[globalPointI]) if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
{ {
procToGlobalFZmap[faceZoneI][globalPointI] = procPointI; procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
break; 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 //- procToGlobalFZmap now contains the local FZpoint label for each
pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0); //- global FZ point label - for each faceZone
//- find corresponding patch //- check what points are on the current proc patch
string faceZoneName = mesh.faceZones().names()[faceZoneI]; pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
//- 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) //- find corresponding patch
{ string faceZoneName = mesh.faceZones().names()[faceZoneI];
forAll(mesh.boundaryMesh()[patchID].localPoints(), pi) //- 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));
if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi]) label patchID = mesh.boundaryMesh().findPatchID(patchName);
{ if(patchID == -1)
pointOnLocalProcPatch[faceZoneI][fzpi] = 1; {
break; FatalError << "Patch " << patchName << " not found corresponding for faceZone"
} << faceZoneName << exit(FatalError);
} }
}
} forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
} //- end if(!globalFaceZoneMappingSet) {
} 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 //- write to disk to allow restart of cases
//- because it is not possible to calculate the //- because it is not possible to calculate the
//- mapping after the meshes have moved //- mapping after the meshes have moved
if(!globalFaceZoneMappingSet && Pstream::parRun()) if(!globalFaceZoneMappingSet && Pstream::parRun())
{ {
procToGlobalFZmap.write(); procToGlobalFZmap.write();
pointOnLocalProcPatch.write(); pointOnLocalProcPatch.write();
} }

View file

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

View file

@ -61,21 +61,21 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "createGlobalToLocalFaceZonePointMap.H" # include "createGlobalToLocalFaceZonePointMap.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl; Info<< "\nStarting time loop\n" << endl;
@ -101,91 +101,93 @@ int main(int argc, char *argv[])
//- reset DU to zero at the start of the time-step if //- reset DU to zero at the start of the time-step if
//- a predictor is not required //- a predictor is not required
if(!predictor) if(!predictor)
DU = dimensionedVector("zero", dimLength, vector::zero); {
DU = dimensionedVector("zero", dimLength, vector::zero);
}
do //- start of momentum loop do //- start of momentum loop
{ {
DU.storePrevIter(); DU.storePrevIter();
//- correct the contact boundaries //- correct the contact boundaries
if(iCorr % uEqnContactCorrFreq == 0) if(iCorr % uEqnContactCorrFreq == 0)
{ {
Info << "\t\tCorrecting contact in the momentum loop " Info << "\t\tCorrecting contact in the momentum loop "
<< "iteration: " << iCorr << "iteration: " << iCorr
<< ", residual: " << residual << ", residual: " << residual
<< endl; << endl;
//# include "moveMeshLeastSquares.H" //# include "moveMeshLeastSquares.H"
# include "moveSolidMesh.H" # include "moveSolidMesh.H"
contact.correct(); contact.correct();
mesh.movePoints(oldMeshPoints); mesh.movePoints(oldMeshPoints);
} }
# include "calculateDivDSigmaExp.H" # include "calculateDivDSigmaExp.H"
fvVectorMatrix DUEqn fvVectorMatrix DUEqn
( (
fvm::d2dt2(rho, DU) fvm::d2dt2(rho, DU)
== ==
fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)") fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
+ divDSigmaExp + divDSigmaExp
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->correct(DUEqn); solidInterfacePtr->correct(DUEqn);
} }
solverPerf = DUEqn.solve(); solverPerf = DUEqn.solve();
DU.relax(); DU.relax();
solverName = solverPerf.solverName(); solverName = solverPerf.solverName();
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
gradDU = solidInterfacePtr->grad(DU); gradDU = solidInterfacePtr->grad(DU);
} }
else else
{ {
gradDU = fvc::grad(DU); gradDU = fvc::grad(DU);
} }
U = U.oldTime() + DU; U = U.oldTime() + DU;
residual = solverPerf.initialResidual(); residual = solverPerf.initialResidual();
//****************************************************// //****************************************************//
// The contact residual is the initial residual for the // The contact residual is the initial residual for the
// first iteration of the momentum equation // first iteration of the momentum equation
//****************************************************// //****************************************************//
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << DU.name() << ", Solving for " << DU.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} //- end of momentum loop } //- end of momentum loop
while while
( (
relativeResidual > convergenceTolerance relativeResidual > convergenceTolerance
//residual > convergenceTolerance //residual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
// Print out info per contact iteration // Print out info per contact iteration
Info << "\t\tSolving for " << DU.name() Info << "\t\tSolving for " << DU.name()
<< " using " << solverName << " using " << solverName
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl; << ", No outer iterations " << iCorr << endl;
lduMatrix::debug = 1; lduMatrix::debug = 1;
@ -205,8 +207,8 @@ int main(int argc, char *argv[])
//mesh.movePoints(oldMeshPoints); //mesh.movePoints(oldMeshPoints);
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl << endl; << endl << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;

View file

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

View file

@ -1,5 +1,5 @@
if (runTime.outputTime()) if (runTime.outputTime())
{ {
// FAILS IN PARALLEL - FIX // FAILS IN PARALLEL - FIX
// Info << "Print contact area" << endl; // Info << "Print contact area" << endl;
//volScalarField ca = contact.contactArea(); //volScalarField ca = contact.contactArea();
@ -12,39 +12,39 @@ if (runTime.outputTime())
//-------------------------------------------------------------// //-------------------------------------------------------------//
//This is the point distance for each contact vertex //This is the point distance for each contact vertex
pointScalarField cGapPoints pointScalarField cGapPoints
( (
IOobject IOobject
( (
"pointContactGap", "pointContactGap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedScalar("scalar", dimLength, 0.0), dimensionedScalar("scalar", dimLength, 0.0),
"calculated" "calculated"
); );
contact.contactGapPoints(cGapPoints); contact.contactGapPoints(cGapPoints);
cGapPoints.write(); cGapPoints.write();
//- This is the point distance for each contact vertex //- This is the point distance for each contact vertex
pointVectorField cPointForce pointVectorField cPointForce
( (
IOobject IOobject
( (
"pointContactForce", "pointContactForce",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimForce, vector::zero), dimensionedVector("vector", dimForce, vector::zero),
"calculated" "calculated"
); );
contact.contactPointForce(cPointForce); contact.contactPointForce(cPointForce);
cPointForce.write(); cPointForce.write();
@ -52,4 +52,4 @@ if (runTime.outputTime())
//- SHOULD THIS BE A REF TO A TMP...? //- SHOULD THIS BE A REF TO A TMP...?
volScalarField cPressure = contact.contactPressure(); volScalarField cPressure = contact.contactPressure();
cPressure.write(); cPressure.write();
} }

View file

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

View file

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

View file

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

View file

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

View file

@ -15,59 +15,58 @@
FieldField<Field, vector> extraVecs(ptc.size()); FieldField<Field, vector> extraVecs(ptc.size());
{ {
const labelListList& pfaces = mesh.pointFaces(); const labelListList& pfaces = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = pfaces[curPoint]; const labelList& curFaces = pfaces[curPoint];
// extraVecs.hook(new vectorField(curFaces.size())); //- no hook function // extraVecs.hook(new vectorField(curFaces.size())); //- no hook function
extraVecs.set extraVecs.set
( (
pointI, pointI,
new vectorField(curFaces.size()) new vectorField(curFaces.size())
); );
vectorField& curExtraVectors = extraVecs[pointI]; vectorField& curExtraVectors = extraVecs[pointI];
label nFacesAroundPoint = 0; label nFacesAroundPoint = 0;
const vector& pointLoc = mesh.points()[curPoint]; const vector& pointLoc = mesh.points()[curPoint];
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(bm[patchID]) !isA<emptyFvPatch>(bm[patchID])
&& !bm[patchID].coupled() && !bm[patchID].coupled()
) )
{ {
// Found a face for extrapolation // Found a face for extrapolation
curExtraVectors[nFacesAroundPoint] = curExtraVectors[nFacesAroundPoint] =
pointLoc pointLoc
- centres.boundaryField()[patchID] - centres.boundaryField()[patchID]
[bm[patchID].patch().whichFace(curFaces[faceI])]; [bm[patchID].patch().whichFace(curFaces[faceI])];
nFacesAroundPoint++; nFacesAroundPoint++;
} }
} }
} }
curExtraVectors.setSize(nFacesAroundPoint); curExtraVectors.setSize(nFacesAroundPoint);
} }
} }

View file

@ -8,114 +8,116 @@
FieldField<Field, scalar> w(ptc.size()); FieldField<Field, scalar> w(ptc.size());
{ {
const labelListList& pf = mesh.pointFaces(); const labelListList& pf = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
pointScalarField volPointSumWeights pointScalarField volPointSumWeights
( (
IOobject IOobject
( (
"volPointSumWeights", "volPointSumWeights",
mesh.polyMesh::instance(), mesh.polyMesh::instance(),
mesh mesh
), ),
pMesh, pMesh,
dimensionedScalar("zero", dimless, 0) 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)
forAll (ptc, pointI)
{ {
w[pointI] /= volPointSumWeights[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]];
} }
} }

View file

@ -23,138 +23,131 @@ philipc
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
//***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****// //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
contactPatchPairList& contacts = contact; contactPatchPairList& contacts = contact;
forAll(contacts, contactI) forAll(contacts, contactI)
{ {
label masterID = contacts[contactI].masterPatch().index(); label masterID = contacts[contactI].masterPatch().index();
label slaveID = contacts[contactI].slavePatch().index(); label slaveID = contacts[contactI].slavePatch().index();
primitivePatchInterpolation masterInterpolator primitivePatchInterpolation masterInterpolator
(mesh.boundaryMesh()[masterID]); (
primitivePatchInterpolation slaveInterpolator mesh.boundaryMesh()[masterID]
(mesh.boundaryMesh()[slaveID]); );
primitivePatchInterpolation slaveInterpolator
(
mesh.boundaryMesh()[slaveID]
);
//- DU must be interpolated to the vertices, this ignores the faceZone //- DU must be interpolated to the vertices, this ignores the faceZone
//- points with no DU (unlike volPointInterpolation) //- points with no DU (unlike volPointInterpolation)
vectorField correctMasterPointDU = vectorField correctMasterPointDU =
masterInterpolator.faceToPointInterpolate<vector> masterInterpolator.faceToPointInterpolate<vector>
( (
DU.boundaryField()[masterID] DU.boundaryField()[masterID]
); );
vectorField correctSlavePointDU = vectorField correctSlavePointDU =
slaveInterpolator.faceToPointInterpolate<vector> slaveInterpolator.faceToPointInterpolate<vector>
( (
DU.boundaryField()[slaveID] DU.boundaryField()[slaveID]
); );
vectorField oldMasterPoints = vectorField oldMasterPoints =
mesh.boundaryMesh()[masterID].localPoints(); mesh.boundaryMesh()[masterID].localPoints();
vectorField oldSlavePoints = vectorField oldSlavePoints =
mesh.boundaryMesh()[slaveID].localPoints(); 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];
}
}
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 *****// //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
//- find the patch corresponding to this faceZone //- find the patch corresponding to this faceZone
//- assuming that the FZ is called <patch_name>FaceZone //- assuming that the FZ is called <patch_name>FaceZone
string faceZoneName = mesh.faceZones().names()[faceZoneI]; string faceZoneName = mesh.faceZones().names()[faceZoneI];
//- remove the string FaceZone from the end of the face zone name to get the patch name //- 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)); string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
label patchID = mesh.boundaryMesh().findPatchID(patchName); label patchID = mesh.boundaryMesh().findPatchID(patchName);
if(patchID == -1) if(patchID == -1)
{ {
FatalError << "Patch " << patchName << " not found corresponding for faceZone" FatalError << "Patch " << patchName << " not found corresponding for faceZone"
<< faceZoneName << exit(FatalError); << faceZoneName << exit(FatalError);
} }
vectorField globalFZpoints = vectorField globalFZpoints =
mesh.faceZones()[faceZoneI]().localPoints(); mesh.faceZones()[faceZoneI]().localPoints();
//- new points for the face zone //- new points for the face zone
vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero); vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
//- inter-proc points are shared by multiple procs //- inter-proc points are shared by multiple procs
//- pointNumProc is the number of procs which a point lies on //- pointNumProc is the number of procs which a point lies on
scalarField pointNumProcs(globalFZpoints.size(), 0.0); scalarField pointNumProcs(globalFZpoints.size(), 0.0);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
//if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size()) //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
if(pointOnLocalProcPatch[faceZoneI][localPoint]) if(pointOnLocalProcPatch[faceZoneI][localPoint])
{ {
label procPoint = label procPoint =
mesh.faceZones()[faceZoneI]().meshPoints()[localPoint]; mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
globalFZnewPoints[globalPointI] = globalFZnewPoints[globalPointI] = newPoints[procPoint];
newPoints[procPoint]; pointNumProcs[globalPointI] = 1;
pointNumProcs[globalPointI] = 1; }
} }
}
reduce(globalFZnewPoints, sumOp<vectorField>()); reduce(globalFZnewPoints, sumOp<vectorField>());
reduce(pointNumProcs, sumOp<scalarField>()); reduce(pointNumProcs, sumOp<scalarField>());
//- now average the newPoints between all procs //- now average the newPoints between all procs
if(min(pointNumProcs) < 1) if(min(pointNumProcs) < 1)
{ {
FatalError << "pointNumProc has not been set for all points" << exit(FatalError); FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
} }
globalFZnewPoints /= pointNumProcs; globalFZnewPoints /= pointNumProcs;
//- the globalFZnewPoints now contains the correct FZ new points in //- the globalFZnewPoints now contains the correct FZ new points in
//- a global order, now convert them back into the local proc order //- a global order, now convert them back into the local proc order
vectorField procFZnewPoints(globalFZpoints.size(), vector::zero); vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
procFZnewPoints[localPoint] = globalFZnewPoints[globalPointI];
}
procFZnewPoints[localPoint] = //- now fix the newPoints points on the globalFaceZones
globalFZnewPoints[globalPointI]; labelList procFZmeshPoints = mesh.faceZones()[faceZoneI]().meshPoints();
}
//- now fix the newPoints points on the globalFaceZones forAll(procFZmeshPoints, pointI)
labelList procFZmeshPoints = {
mesh.faceZones()[faceZoneI]().meshPoints(); label procPoint = procFZmeshPoints[pointI];
newPoints[procPoint] = procFZnewPoints[pointI];
forAll(procFZmeshPoints, pointI) }
{ }
label procPoint = procFZmeshPoints[pointI]; }
newPoints[procPoint] =
procFZnewPoints[pointI];
}
}
}

View file

@ -25,7 +25,7 @@
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimLength, vector::zero) dimensionedVector("zero", dimLength, vector::zero)
); );
volSymmTensorField DEpsilon volSymmTensorField DEpsilon
@ -84,35 +84,35 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
volVectorField divDSigmaLargeStrainExp volVectorField divDSigmaLargeStrainExp
( (
IOobject IOobject
( (
"divDSigmaLargeStrainExp", "divDSigmaLargeStrainExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
// read rheology properties // read rheology properties
rheologyModel rheology(sigma); rheologyModel rheology(sigma);
volScalarField rho = rheology.rho(); volScalarField rho = rheology.rho();

View file

@ -19,117 +19,117 @@ philipc
//- these are read if present to allow restarting of contact cases //- these are read if present to allow restarting of contact cases
IOList<labelList> procToGlobalFZmap IOList<labelList> procToGlobalFZmap
( (
IOobject IOobject
( (
"procToGlobalFZmap", "procToGlobalFZmap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
IOList<labelList> pointOnLocalProcPatch IOList<labelList> pointOnLocalProcPatch
( (
IOobject IOobject
( (
"pointOnLocalProcPatch", "pointOnLocalProcPatch",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
//- if they have been read then don't recalculate it //- if they have been read then don't recalculate it
bool globalFaceZoneMappingSet = false; bool globalFaceZoneMappingSet = false;
if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0) if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
{ {
Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases" Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
<< endl; << endl;
globalFaceZoneMappingSet = true; globalFaceZoneMappingSet = true;
} }
else else
{ {
Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl 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 message should only appear starting a new analysis" << endl;
} }
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
if(!globalFaceZoneMappingSet) if(!globalFaceZoneMappingSet)
{ {
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0); procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
//- set all slave points to zero because only the master order is used //- set all slave points to zero because only the master order is used
if(!Pstream::master()) if(!Pstream::master())
globalFZpoints *= 0.0; {
globalFZpoints *= 0.0;
}
//- pass points to all procs //- pass points to all procs
reduce(globalFZpoints, sumOp<vectorField>()); reduce(globalFZpoints, sumOp<vectorField>());
//- now every proc has the master's list of FZ points //- now every proc has the master's list of FZ points
//- every proc must now find the mapping from their local FZpoints to //- every proc must now find the mapping from their local FZpoints to
//- the globalFZpoints //- the globalFZpoints
vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
forAll(globalFZpoints, globalPointI) forAll(globalFZpoints, globalPointI)
{ {
forAll(procFZpoints, procPointI) forAll(procFZpoints, procPointI)
{ {
if(procFZpoints[procPointI] == globalFZpoints[globalPointI]) if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
{ {
procToGlobalFZmap[faceZoneI][globalPointI] = procPointI; procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
break; 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 //- check what points are on the current proc patch
pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0); pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
//- find corresponding patch //- find corresponding patch
string faceZoneName = mesh.faceZones().names()[faceZoneI]; string faceZoneName = mesh.faceZones().names()[faceZoneI];
//- remove the string FaceZone from the end of the face zone name to get the patch name //- 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)); string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
label patchID = mesh.boundaryMesh().findPatchID(patchName); label patchID = mesh.boundaryMesh().findPatchID(patchName);
if(patchID == -1) if(patchID == -1)
{ {
FatalError << "Patch " << patchName << " not found corresponding for faceZone" FatalError << "Patch " << patchName << " not found corresponding for faceZone"
<< faceZoneName << exit(FatalError); << faceZoneName << exit(FatalError);
} }
forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi) forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
{ {
forAll(mesh.boundaryMesh()[patchID].localPoints(), pi) forAll(mesh.boundaryMesh()[patchID].localPoints(), pi)
{ {
if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi]) if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi])
{ {
pointOnLocalProcPatch[faceZoneI][fzpi] = 1; pointOnLocalProcPatch[faceZoneI][fzpi] = 1;
break; break;
} }
} }
} }
} }
} //- end if(!globalFaceZoneMappingSet) } //- end if(!globalFaceZoneMappingSet)
} }
//- write to disk to allow restart of cases //- write to disk to allow restart of cases
//- because it is not possible to calculate the //- because it is not possible to calculate the
//- mapping after the meshes have moved //- mapping after the meshes have moved
if(!globalFaceZoneMappingSet) if(!globalFaceZoneMappingSet)
{ {
procToGlobalFZmap.write(); procToGlobalFZmap.write();
pointOnLocalProcPatch.write(); pointOnLocalProcPatch.write();
} }

View file

@ -66,137 +66,136 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "readDivDSigmaLargeStrainMethod.H" # include "readDivDSigmaLargeStrainMethod.H"
# include "readMoveMeshMethod.H" # include "readMoveMeshMethod.H"
# include "createGlobalToLocalFaceZonePointMap.H" # include "createGlobalToLocalFaceZonePointMap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl; Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++) for (runTime++; !runTime.end(); runTime++)
{ {
Info<< "Time: " << runTime.timeName() << endl; Info<< "Time: " << runTime.timeName() << endl;
# include "readContactControls.H" # include "readContactControls.H"
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
//-- for moving the mesh and then back again //-- for moving the mesh and then back again
vectorField oldMeshPoints = mesh.allPoints(); vectorField oldMeshPoints = mesh.allPoints();
int iCorr = 0; int iCorr = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
word solverName; word solverName;
lduMatrix::debug = 0; lduMatrix::debug = 0;
scalar residual = GREAT; scalar residual = GREAT;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
do //- start of momentum loop do //- start of momentum loop
{ {
DU.storePrevIter(); DU.storePrevIter();
divDSigmaLargeStrainExp.storePrevIter(); divDSigmaLargeStrainExp.storePrevIter();
//- correct the contact boundaries //- correct the contact boundaries
if(iCorr % uEqnContactCorrFreq == 0) if(iCorr % uEqnContactCorrFreq == 0)
{ {
Info << "\t\tCorrecting contact in the momentum loop " Info << "\t\tCorrecting contact in the momentum loop "
<< "iteration: " << iCorr << "iteration: " << iCorr
<< ", residual: " << residual << ", residual: " << residual
<< endl; << endl;
//# include "moveMeshLeastSquares.H" //# include "moveMeshLeastSquares.H"
# include "moveSolidMeshForContact.H" # include "moveSolidMeshForContact.H"
contact.correct(); contact.correct();
mesh.movePoints(oldMeshPoints); mesh.movePoints(oldMeshPoints);
} }
# include "calculateDivDSigmaExp.H" # include "calculateDivDSigmaExp.H"
# include "calculateDivDSigmaExpLargeStrain.H" # include "calculateDivDSigmaExpLargeStrain.H"
fvVectorMatrix DUEqn fvVectorMatrix DUEqn
( (
fvm::d2dt2(rho, DU) fvm::d2dt2(rho, DU)
== ==
fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)") fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
+ divDSigmaExp + divDSigmaExp
+ divDSigmaLargeStrainExp + divDSigmaLargeStrainExp
);
); solverPerf = DUEqn.solve();
solverPerf = DUEqn.solve(); DU.relax();
DU.relax(); solverName = solverPerf.solverName();
solverName = solverPerf.solverName(); gradDU = fvc::grad(DU);
gradDU = fvc::grad(DU); DF = gradDU.T();
DF = gradDU.T(); # include "calculateDEpsilonDSigma.H"
# include "calculateDEpsilonDSigma.H" residual = solverPerf.initialResidual();
residual = solverPerf.initialResidual(); if(iCorr == 0)
{
initialResidual = solverPerf.initialResidual();
}
if(iCorr == 0) # include "calculateRelativeResidual.H"
{
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
);
Info << "\tTime " << runTime.value() // Print out info per contact iteration
<< ", Corrector " << iCorr Info << "\t\tSolving for " << DU.name()
<< ", Solving for " << DU.name() << " using " << solverName
<< " using " << solverPerf.solverName() << ", Initial residual = " << initialResidual
<< ", residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", No outer iterations " << iCorr << endl;
} //- end of momentum loop
while
(
relativeResidual > convergenceTolerance
//residual > convergenceTolerance
&&
++iCorr < nCorr
);
// Print out info per contact iteration lduMatrix::debug = 1;
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 "rotateFields.H" # include "moveMesh.H"
# include "moveMesh.H" # include "writeFields.H"
# include "writeFields.H" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << endl << endl;
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

@ -8,26 +8,26 @@ const fvBoundaryMesh& bm = mesh.boundary();
forAll (bm, patchI) forAll (bm, patchI)
{ {
// If the patch is empty, skip it // If the patch is empty, skip it
// If the patch is coupled, and there are no cyclic parallels, skip it // If the patch is coupled, and there are no cyclic parallels, skip it
if if
( (
!isA<emptyFvPatch>(bm[patchI]) !isA<emptyFvPatch>(bm[patchI])
&& !( && !(
bm[patchI].coupled() bm[patchI].coupled()
//&& Pstream::parRun() //&& Pstream::parRun()
//&& !mesh.parallelData().cyclicParallel() //&& !mesh.parallelData().cyclicParallel()
) )
) )
{ {
const labelList& bp = bm[patchI].patch().boundaryPoints(); const labelList& bp = bm[patchI].patch().boundaryPoints();
const labelList& meshPoints = bm[patchI].patch().meshPoints(); const labelList& meshPoints = bm[patchI].patch().meshPoints();
forAll (bp, pointI) forAll (bp, pointI)
{ {
pointsCorrectionMap.insert(meshPoints[bp[pointI]]); pointsCorrectionMap.insert(meshPoints[bp[pointI]]);
} }
} }
} }

View file

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

View file

@ -2,7 +2,7 @@
//- move mesh //- move mesh
//--------------------------------------------------// //--------------------------------------------------//
if(min(J.internalField()) > 0) if(min(J.internalField()) > 0)
{ {
Info << "Move solid mesh using inverse distance interpolation" << endl; Info << "Move solid mesh using inverse distance interpolation" << endl;
// Create point mesh // Create point mesh
@ -12,24 +12,24 @@ if(min(J.internalField()) > 0)
volPointInterpolation pointInterpolation(mesh); volPointInterpolation pointInterpolation(mesh);
wordList types wordList types
( (
pMesh.boundary().size(), pMesh.boundary().size(),
//fixedValueFvPatchVectorField::typeName //fixedValueFvPatchVectorField::typeName
calculatedFvPatchVectorField::typeName calculatedFvPatchVectorField::typeName
); );
pointVectorField pointDU pointVectorField pointDU
( (
IOobject IOobject
( (
"pointDU", "pointDU",
runTime.timeName(), runTime.timeName(),
mesh mesh
), ),
pMesh, pMesh,
dimensionedVector("zero", dimLength, vector::zero), dimensionedVector("zero", dimLength, vector::zero),
types types
); );
// Calculate mesh points displacement // Calculate mesh points displacement
pointInterpolation.interpolate(DU, pointDU); pointInterpolation.interpolate(DU, pointDU);
@ -41,26 +41,25 @@ if(min(J.internalField()) > 0)
//pointDU.write(); //pointDU.write();
const vectorField& pointDUI = const vectorField& pointDUI = pointDU.internalField();
pointDU.internalField();
// Move mesh // Move mesh
vectorField newPoints = mesh.allPoints(); vectorField newPoints = mesh.allPoints();
forAll (pointDUI, pointI) forAll (pointDUI, pointI)
{ {
newPoints[pointI] += pointDUI[pointI]; newPoints[pointI] += pointDUI[pointI];
} }
twoDPointCorrector twoDCorrector(mesh); twoDPointCorrector twoDCorrector(mesh);
twoDCorrector.correctPoints(newPoints); twoDCorrector.correctPoints(newPoints);
mesh.movePoints(newPoints); mesh.movePoints(newPoints);
mesh.V00(); mesh.V00();
mesh.moving(false); mesh.moving(false);
} }
else else
{ {
FatalErrorIn(args.executable()) FatalErrorIn(args.executable())
<< "Negative Jacobian" << "Negative Jacobian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

@ -7,38 +7,40 @@ pointVectorField& pf = pointDU;
// Do the correction // Do the correction
//GeometricField<Type, pointPatchField, pointMesh> pfCorr //GeometricField<Type, pointPatchField, pointMesh> pfCorr
/*pointVectorField pfCorr /*
pointVectorField pfCorr
( (
IOobject IOobject
( (
// "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr", // "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
"edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr", "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
//vf.instance(), //vf.instance(),
DU, DU,
pMesh, pMesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
//dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero), //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
dimensionedVector("zero", pf.dimensions(), vector::zero), dimensionedVector("zero", pf.dimensions(), vector::zero),
pf.boundaryField().types() pf.boundaryField().types()
);*/ );
*/
pointVectorField pfCorr pointVectorField pfCorr
( (
IOobject IOobject
( (
"pointDUcorr", "pointDUcorr",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimLength, vector::zero), dimensionedVector("vector", dimLength, vector::zero),
"calculated" "calculated"
); );
//const labelList& ptc = boundaryPoints(); //const labelList& ptc = boundaryPoints();
#include "findBoundaryPoints.H" #include "findBoundaryPoints.H"
@ -57,64 +59,68 @@ const labelListList& PointFaces = mesh.pointFaces();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = PointFaces[curPoint]; const labelList& curFaces = PointFaces[curPoint];
label fI = 0; label fI = 0;
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(mesh.boundary()[patchID]) !isA<emptyFvPatch>(mesh.boundary()[patchID])
&& !mesh.boundary()[patchID].coupled() && !mesh.boundary()[patchID].coupled()
) )
{ {
label faceInPatchID = label faceInPatchID =
bm[patchID].patch().whichFace(curFaces[faceI]); bm[patchID].patch().whichFace(curFaces[faceI]);
pfCorr[curPoint] += pfCorr[curPoint] +=
w[pointI][fI]* w[pointI][fI]*
( (
extraVecs[pointI][fI] extraVecs[pointI][fI]
& gradDU.boundaryField()[patchID][faceInPatchID] & gradDU.boundaryField()[patchID][faceInPatchID]
); );
fI++; fI++;
} }
} }
} }
} }
// Update coupled boundaries // Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].initAddField(); pfCorr.boundaryField()[patchI].initAddField();
} }
}*/ }
*/
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].addField(pfCorr.internalField()); pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
} }
}*/ }
*/
//Info << "pfCorr: " << pfCorr << endl; //Info << "pfCorr: " << pfCorr << endl;
pfCorr.correctBoundaryConditions(); pfCorr.correctBoundaryConditions();
//pfCorr.write(); //pfCorr.write();

View file

@ -1,5 +1,5 @@
if (runTime.outputTime()) if (runTime.outputTime())
{ {
// FAILS IN PARALLEL - FIX // FAILS IN PARALLEL - FIX
// Info << "Print contact area" << endl; // Info << "Print contact area" << endl;
//volScalarField ca = contact.contactArea(); //volScalarField ca = contact.contactArea();
@ -12,39 +12,39 @@ if (runTime.outputTime())
//-------------------------------------------------------------// //-------------------------------------------------------------//
//This is the point distance for each contact vertex //This is the point distance for each contact vertex
pointScalarField cGapPoints pointScalarField cGapPoints
( (
IOobject IOobject
( (
"pointContactGap", "pointContactGap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedScalar("scalar", dimLength, 0.0), dimensionedScalar("scalar", dimLength, 0.0),
"calculated" "calculated"
); );
contact.contactGapPoints(cGapPoints); contact.contactGapPoints(cGapPoints);
cGapPoints.write(); cGapPoints.write();
//- This is the point distance for each contact vertex //- This is the point distance for each contact vertex
pointVectorField cPointForce pointVectorField cPointForce
( (
IOobject IOobject
( (
"pointContactForce", "pointContactForce",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimForce, vector::zero), dimensionedVector("vector", dimForce, vector::zero),
"calculated" "calculated"
); );
contact.contactPointForce(cPointForce); contact.contactPointForce(cPointForce);
cPointForce.write(); cPointForce.write();
@ -52,4 +52,4 @@ if (runTime.outputTime())
//- SHOULD THIS BE A REF TO A TMP...? //- SHOULD THIS BE A REF TO A TMP...?
volScalarField cPressure = contact.contactPressure(); volScalarField cPressure = contact.contactPressure();
cPressure.write(); cPressure.write();
} }

View file

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

View file

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

View file

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

View file

@ -23,138 +23,131 @@ philipc
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
//***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****// //***** FIX INCORRECT POINT ON PATCHES WITH FACEZONE *****//
contactPatchPairList& contacts = contact; contactPatchPairList& contacts = contact;
forAll(contacts, contactI) forAll(contacts, contactI)
{ {
label masterID = contacts[contactI].masterPatch().index(); label masterID = contacts[contactI].masterPatch().index();
label slaveID = contacts[contactI].slavePatch().index(); label slaveID = contacts[contactI].slavePatch().index();
primitivePatchInterpolation masterInterpolator primitivePatchInterpolation masterInterpolator
(mesh.boundaryMesh()[masterID]); (
primitivePatchInterpolation slaveInterpolator mesh.boundaryMesh()[masterID]
(mesh.boundaryMesh()[slaveID]); );
primitivePatchInterpolation slaveInterpolator
(
mesh.boundaryMesh()[slaveID]
);
//- U must be interpolated to the vertices, this ignores the faceZone //- U must be interpolated to the vertices, this ignores the faceZone
//- points with no U (unlike volPointInterpolation) //- points with no U (unlike volPointInterpolation)
vectorField correctMasterPointU = vectorField correctMasterPointU =
masterInterpolator.faceToPointInterpolate<vector> masterInterpolator.faceToPointInterpolate<vector>
( (
U.boundaryField()[masterID] U.boundaryField()[masterID]
); );
vectorField correctSlavePointU = vectorField correctSlavePointU =
slaveInterpolator.faceToPointInterpolate<vector> slaveInterpolator.faceToPointInterpolate<vector>
( (
U.boundaryField()[slaveID] U.boundaryField()[slaveID]
); );
vectorField oldMasterPoints = vectorField oldMasterPoints =
mesh.boundaryMesh()[masterID].localPoints(); mesh.boundaryMesh()[masterID].localPoints();
vectorField oldSlavePoints = vectorField oldSlavePoints =
mesh.boundaryMesh()[slaveID].localPoints(); 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];
}
}
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 *****// //***** NOW FIX AND SYNCHRONISE ALL THE FACEZONE POINTS *****//
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
//- find the patch corresponding to this faceZone //- find the patch corresponding to this faceZone
//- assuming that the FZ is called <patch_name>FaceZone //- assuming that the FZ is called <patch_name>FaceZone
string faceZoneName = mesh.faceZones().names()[faceZoneI]; string faceZoneName = mesh.faceZones().names()[faceZoneI];
//- remove the string FaceZone from the end of the face zone name to get the patch name //- 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)); string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
label patchID = mesh.boundaryMesh().findPatchID(patchName); label patchID = mesh.boundaryMesh().findPatchID(patchName);
if(patchID == -1) if(patchID == -1)
{ {
FatalError << "Patch " << patchName << " not found corresponding for faceZone" FatalError << "Patch " << patchName << " not found corresponding for faceZone"
<< faceZoneName << exit(FatalError); << faceZoneName << exit(FatalError);
} }
vectorField globalFZpoints = vectorField globalFZpoints =
mesh.faceZones()[faceZoneI]().localPoints(); mesh.faceZones()[faceZoneI]().localPoints();
//- new points for the face zone //- new points for the face zone
vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero); vectorField globalFZnewPoints(globalFZpoints.size(), vector::zero);
//- inter-proc points are shared by multiple procs //- inter-proc points are shared by multiple procs
//- pointNumProc is the number of procs which a point lies on //- pointNumProc is the number of procs which a point lies on
scalarField pointNumProcs(globalFZpoints.size(), 0.0); scalarField pointNumProcs(globalFZpoints.size(), 0.0);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
//if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size()) //if(localPoint < mesh.boundaryMesh()[patchID].localPoints().size())
if(pointOnLocalProcPatch[faceZoneI][localPoint]) if(pointOnLocalProcPatch[faceZoneI][localPoint])
{ {
label procPoint = label procPoint =
mesh.faceZones()[faceZoneI]().meshPoints()[localPoint]; mesh.faceZones()[faceZoneI]().meshPoints()[localPoint];
globalFZnewPoints[globalPointI] = globalFZnewPoints[globalPointI] = newPoints[procPoint];
newPoints[procPoint]; pointNumProcs[globalPointI] = 1;
pointNumProcs[globalPointI] = 1; }
} }
}
reduce(globalFZnewPoints, sumOp<vectorField>()); reduce(globalFZnewPoints, sumOp<vectorField>());
reduce(pointNumProcs, sumOp<scalarField>()); reduce(pointNumProcs, sumOp<scalarField>());
//- now average the newPoints between all procs //- now average the newPoints between all procs
if(min(pointNumProcs) < 1) if(min(pointNumProcs) < 1)
{ {
FatalError << "pointNumProc has not been set for all points" << exit(FatalError); FatalError << "pointNumProc has not been set for all points" << exit(FatalError);
} }
globalFZnewPoints /= pointNumProcs; globalFZnewPoints /= pointNumProcs;
//- the globalFZnewPoints now contains the correct FZ new points in //- the globalFZnewPoints now contains the correct FZ new points in
//- a global order, now convert them back into the local proc order //- a global order, now convert them back into the local proc order
vectorField procFZnewPoints(globalFZpoints.size(), vector::zero); vectorField procFZnewPoints(globalFZpoints.size(), vector::zero);
forAll(globalFZnewPoints, globalPointI) forAll(globalFZnewPoints, globalPointI)
{ {
label localPoint = procToGlobalFZmap[faceZoneI][globalPointI]; label localPoint = procToGlobalFZmap[faceZoneI][globalPointI];
procFZnewPoints[localPoint] = globalFZnewPoints[globalPointI];
}
procFZnewPoints[localPoint] = //- now fix the newPoints points on the globalFaceZones
globalFZnewPoints[globalPointI]; labelList procFZmeshPoints = mesh.faceZones()[faceZoneI]().meshPoints();
}
//- now fix the newPoints points on the globalFaceZones forAll(procFZmeshPoints, pointI)
labelList procFZmeshPoints = {
mesh.faceZones()[faceZoneI]().meshPoints(); label procPoint = procFZmeshPoints[pointI];
newPoints[procPoint] = procFZnewPoints[pointI];
forAll(procFZmeshPoints, pointI) }
{ }
label procPoint = procFZmeshPoints[pointI]; }
newPoints[procPoint] =
procFZnewPoints[pointI];
}
}
}

View file

@ -59,22 +59,22 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divSigmaExp volVectorField divSigmaExp
( (
IOobject IOobject
( (
"divSigmaExp", "divSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
// read rheology properties // read rheology properties
rheologyModel rheology(sigma); rheologyModel rheology(sigma);
volScalarField rho = rheology.rho(); volScalarField rho = rheology.rho();

View file

@ -19,117 +19,117 @@ philipc
//- these are read if present to allow restarting of contact cases //- these are read if present to allow restarting of contact cases
IOList<labelList> procToGlobalFZmap IOList<labelList> procToGlobalFZmap
( (
IOobject IOobject
( (
"procToGlobalFZmap", "procToGlobalFZmap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
IOList<labelList> pointOnLocalProcPatch IOList<labelList> pointOnLocalProcPatch
( (
IOobject IOobject
( (
"pointOnLocalProcPatch", "pointOnLocalProcPatch",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh.faceZones().size() mesh.faceZones().size()
); );
//- if they have been read then don't recalculate it //- if they have been read then don't recalculate it
bool globalFaceZoneMappingSet = false; bool globalFaceZoneMappingSet = false;
if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0) if(gMax(procToGlobalFZmap[0]) > 0 && gMax(pointOnLocalProcPatch[0]) > 0)
{ {
Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases" Info << "Reading procToGlobalFZmap and pointOnLocalProcPatch allowing restart of contact cases"
<< endl; << endl;
globalFaceZoneMappingSet = true; globalFaceZoneMappingSet = true;
} }
else else
{ {
Info << "procToGlobalFZmap and pointOnLocalProcPatch will be calculated as it has not been found" << nl 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 message should only appear starting a new analysis" << endl;
} }
//- this is only needed in a parallel runs //- this is only needed in a parallel runs
if(Pstream::parRun()) if(Pstream::parRun())
{ {
if(!globalFaceZoneMappingSet) if(!globalFaceZoneMappingSet)
{ {
forAll(mesh.faceZones(), faceZoneI) forAll(mesh.faceZones(), faceZoneI)
{ {
vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField globalFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0); procToGlobalFZmap[faceZoneI].setSize(globalFZpoints.size(), 0);
//- set all slave points to zero because only the master order is used //- set all slave points to zero because only the master order is used
if(!Pstream::master()) if(!Pstream::master())
globalFZpoints *= 0.0; {
globalFZpoints *= 0.0;
}
//- pass points to all procs //- pass points to all procs
reduce(globalFZpoints, sumOp<vectorField>()); reduce(globalFZpoints, sumOp<vectorField>());
//- now every proc has the master's list of FZ points //- now every proc has the master's list of FZ points
//- every proc must now find the mapping from their local FZpoints to //- every proc must now find the mapping from their local FZpoints to
//- the globalFZpoints //- the globalFZpoints
vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints(); vectorField procFZpoints = mesh.faceZones()[faceZoneI]().localPoints();
forAll(globalFZpoints, globalPointI) forAll(globalFZpoints, globalPointI)
{ {
forAll(procFZpoints, procPointI) forAll(procFZpoints, procPointI)
{ {
if(procFZpoints[procPointI] == globalFZpoints[globalPointI]) if(procFZpoints[procPointI] == globalFZpoints[globalPointI])
{ {
procToGlobalFZmap[faceZoneI][globalPointI] = procPointI; procToGlobalFZmap[faceZoneI][globalPointI] = procPointI;
break; 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 //- check what points are on the current proc patch
pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0); pointOnLocalProcPatch[faceZoneI].setSize(globalFZpoints.size(), 0);
//- find corresponding patch //- find corresponding patch
string faceZoneName = mesh.faceZones().names()[faceZoneI]; string faceZoneName = mesh.faceZones().names()[faceZoneI];
//- remove the string FaceZone from the end of the face zone name to get the patch name //- 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)); string patchName = faceZoneName.substr(0, (faceZoneName.size()-8));
label patchID = mesh.boundaryMesh().findPatchID(patchName); label patchID = mesh.boundaryMesh().findPatchID(patchName);
if(patchID == -1) if(patchID == -1)
{ {
FatalError << "Patch " << patchName << " not found corresponding for faceZone" FatalError << "Patch " << patchName << " not found corresponding for faceZone"
<< faceZoneName << exit(FatalError); << faceZoneName << exit(FatalError);
} }
forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi) forAll(mesh.faceZones()[faceZoneI]().localPoints(), fzpi)
{ {
forAll(mesh.boundaryMesh()[patchID].localPoints(), pi) forAll(mesh.boundaryMesh()[patchID].localPoints(), pi)
{ {
if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi]) if(mesh.faceZones()[faceZoneI]().localPoints()[fzpi] == mesh.boundaryMesh()[patchID].localPoints()[pi])
{ {
pointOnLocalProcPatch[faceZoneI][fzpi] = 1; pointOnLocalProcPatch[faceZoneI][fzpi] = 1;
break; break;
} }
} }
} }
} }
} //- end if(!globalFaceZoneMappingSet) } //- end if(!globalFaceZoneMappingSet)
} }
//- write to disk to allow restart of cases //- write to disk to allow restart of cases
//- because it is not possible to calculate the //- because it is not possible to calculate the
//- mapping after the meshes have moved //- mapping after the meshes have moved
if(!globalFaceZoneMappingSet) if(!globalFaceZoneMappingSet)
{ {
procToGlobalFZmap.write(); procToGlobalFZmap.write();
pointOnLocalProcPatch.write(); pointOnLocalProcPatch.write();
} }

View file

@ -66,140 +66,140 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivSigmaExpMethod.H" # include "readDivSigmaExpMethod.H"
# include "createGlobalToLocalFaceZonePointMap.H" # include "createGlobalToLocalFaceZonePointMap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl; Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++) for (runTime++; !runTime.end(); runTime++)
{ {
Info<< "Time: " << runTime.timeName() << endl; Info<< "Time: " << runTime.timeName() << endl;
# include "readContactControls.H" # include "readContactControls.H"
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
//-- for moving the mesh and then back again //-- for moving the mesh and then back again
vectorField oldMeshPoints = mesh.allPoints(); vectorField oldMeshPoints = mesh.allPoints();
int iCorr = 0; int iCorr = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
word solverName; word solverName;
lduMatrix::debug = 0; lduMatrix::debug = 0;
scalar residual = GREAT; scalar residual = GREAT;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
//- Predictor step //- Predictor step
if (predictor) if (predictor)
{ {
Info << "\nPredicting U, gradU and snGradU based on V, gradV and snGradV\n" << endl; Info << "\nPredicting U, gradU and snGradU based on V, gradV and snGradV\n" << endl;
U += V*runTime.deltaT(); U += V*runTime.deltaT();
gradU += gradV*runTime.deltaT(); gradU += gradV*runTime.deltaT();
snGradU += snGradV*runTime.deltaT(); snGradU += snGradV*runTime.deltaT();
} }
do //- start of momentum loop do //- start of momentum loop
{ {
U.storePrevIter(); U.storePrevIter();
//- correct the contact boundaries //- correct the contact boundaries
if(iCorr % uEqnContactCorrFreq == 0) if(iCorr % uEqnContactCorrFreq == 0)
{ {
Info << "\t\tCorrecting contact in the momentum loop " Info << "\t\tCorrecting contact in the momentum loop "
<< "iteration: " << iCorr << "iteration: " << iCorr
<< ", residual: " << residual << ", residual: " << residual
<< endl; << endl;
//# include "moveMeshLeastSquares.H" //# include "moveMeshLeastSquares.H"
# include "moveSolidMesh.H" # include "moveSolidMesh.H"
contact.correct(); contact.correct();
mesh.movePoints(oldMeshPoints); mesh.movePoints(oldMeshPoints);
} }
# include "calculateDivSigmaExp.H" # include "calculateDivSigmaExp.H"
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
fvm::d2dt2(rho, U) fvm::d2dt2(rho, U)
== ==
fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)") fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)")
+ divSigmaExp + divSigmaExp
); );
solverPerf = UEqn.solve(); solverPerf = UEqn.solve();
U.relax(); U.relax();
solverName = solverPerf.solverName(); solverName = solverPerf.solverName();
gradU = fvc::grad(U); gradU = fvc::grad(U);
snGradU = fvc::snGrad(U); snGradU = fvc::snGrad(U);
residual = solverPerf.initialResidual(); residual = solverPerf.initialResidual();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << U.name() << ", Solving for " << U.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} //- end of momentum loop } //- end of momentum loop
while while
( (
//relativeResidual > convergenceTolerance //relativeResidual > convergenceTolerance
residual > convergenceTolerance residual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
// Print out info per contact iteration // Print out info per contact iteration
Info << "\t\tSolving for " << U.name() Info << "\t\tSolving for " << U.name()
<< " using " << solverName << " using " << solverName
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl; << ", No outer iterations " << iCorr << endl;
lduMatrix::debug = 1; lduMatrix::debug = 1;
V = fvc::ddt(U); V = fvc::ddt(U);
gradV = fvc::ddt(gradU); gradV = fvc::ddt(gradU);
snGradV = (snGradU - snGradU.oldTime())/runTime.deltaT(); snGradV = (snGradU - snGradU.oldTime())/runTime.deltaT();
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"
# include "writeFields.H" # include "writeFields.H"
//# include "moveMeshLeastSquares.H" //# include "moveMeshLeastSquares.H"
//# include "moveSolidMesh.H" //# include "moveSolidMesh.H"
//# include "printContactResults.H" //# include "printContactResults.H"
//mesh.movePoints(oldMeshPoints); // mesh.movePoints(oldMeshPoints);
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl << endl; << endl << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

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

View file

@ -1,5 +1,5 @@
if (runTime.outputTime()) if (runTime.outputTime())
{ {
// FAILS IN PARALLEL - FIX // FAILS IN PARALLEL - FIX
// Info << "Print contact area" << endl; // Info << "Print contact area" << endl;
//volScalarField ca = contact.contactArea(); //volScalarField ca = contact.contactArea();
@ -12,39 +12,39 @@ if (runTime.outputTime())
//-------------------------------------------------------------// //-------------------------------------------------------------//
//This is the point distance for each contact vertex //This is the point distance for each contact vertex
pointScalarField cGapPoints pointScalarField cGapPoints
( (
IOobject IOobject
( (
"pointContactGap", "pointContactGap",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedScalar("scalar", dimLength, 0.0), dimensionedScalar("scalar", dimLength, 0.0),
"calculated" "calculated"
); );
contact.contactGapPoints(cGapPoints); contact.contactGapPoints(cGapPoints);
cGapPoints.write(); cGapPoints.write();
//- This is the point distance for each contact vertex //- This is the point distance for each contact vertex
pointVectorField cPointForce pointVectorField cPointForce
( (
IOobject IOobject
( (
"pointContactForce", "pointContactForce",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimForce, vector::zero), dimensionedVector("vector", dimForce, vector::zero),
"calculated" "calculated"
); );
contact.contactPointForce(cPointForce); contact.contactPointForce(cPointForce);
cPointForce.write(); cPointForce.write();
@ -52,4 +52,4 @@ if (runTime.outputTime())
//- SHOULD THIS BE A REF TO A TMP...? //- SHOULD THIS BE A REF TO A TMP...?
volScalarField cPressure = contact.contactPressure(); volScalarField cPressure = contact.contactPressure();
cPressure.write(); cPressure.write();
} }

View file

@ -4,6 +4,6 @@ Info << divSigmaExpMethod << " method chosen for calculation of sigmaExp" << end
if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian") if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian")
{ {
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian" << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

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

View file

@ -42,19 +42,19 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divSigmaExp volVectorField divSigmaExp
( (
IOobject IOobject
( (
"divSigmaExp", "divSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
rheologyModel rheology(sigma); rheologyModel rheology(sigma);

View file

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

View file

@ -52,116 +52,116 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivSigmaExpMethod.H" # include "readDivSigmaExpMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
# include "readGravity.H" # include "readGravity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
while(runTime.loop()) while(runTime.loop())
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = 0; scalar initialResidual = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::debug=0; lduMatrix::debug=0;
do do
{ {
U.storePrevIter(); U.storePrevIter();
# include "calculateDivSigmaExp.H" # include "calculateDivSigmaExp.H"
//- linear momentum equation //- linear momentum equation
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
fvm::d2dt2(rho, U) fvm::d2dt2(rho, U)
== ==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)") fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+ divSigmaExp + divSigmaExp
+ rho*gravity + rho*gravity
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->correct(UEqn); solidInterfacePtr->correct(UEqn);
} }
solverPerf = UEqn.solve(); solverPerf = UEqn.solve();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
U.relax(); U.relax();
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
gradU = solidInterfacePtr->grad(U); gradU = solidInterfacePtr->grad(U);
} }
else else
{ {
gradU = fvc::grad(U); gradU = fvc::grad(U);
} }
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << U.name() << ", Solving for " << U.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} }
while while
( (
//solverPerf.initialResidual() > convergenceTolerance //solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance relativeResidual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() << ", Solving for " << U.name() Info << nl << "Time " << runTime.value() << ", Solving for " << U.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", Relative residual = " << relativeResidual << ", Relative residual = " << relativeResidual
<< ", No outer iterations " << iCorr << ", No outer iterations " << iCorr
<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
lduMatrix::debug=0; lduMatrix::debug=0;
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"
# include "writeFields.H" # include "writeFields.H"
Info<< "ExecutionTime = " Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime() << runTime.elapsedCpuTime()
<< " s\n\n" << endl; << " s\n\n" << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

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

View file

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

View file

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

View file

@ -83,19 +83,19 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
rheologyModel rheology(sigma); rheologyModel rheology(sigma);

View file

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

View file

@ -51,118 +51,118 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
while(runTime.loop()) while(runTime.loop())
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
lduMatrix::debug = 0; lduMatrix::debug = 0;
do do
{ {
DU.storePrevIter(); DU.storePrevIter();
# include "calculateDivDSigmaExp.H" # include "calculateDivDSigmaExp.H"
//- linear momentum equation //- linear momentum equation
fvVectorMatrix DUEqn fvVectorMatrix DUEqn
( (
fvm::d2dt2(rho, DU) fvm::d2dt2(rho, DU)
== ==
fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)") fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+ divDSigmaExp + divDSigmaExp
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->correct(DUEqn); solidInterfacePtr->correct(DUEqn);
} }
solverPerf = DUEqn.solve(); solverPerf = DUEqn.solve();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
DU.relax(); DU.relax();
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
gradDU = solidInterfacePtr->grad(DU); gradDU = solidInterfacePtr->grad(DU);
} }
else else
{ {
gradDU = fvc::grad(DU); gradDU = fvc::grad(DU);
} }
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << DU.name() << ", Solving for " << DU.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} }
while while
( (
//solverPerf.initialResidual() > convergenceTolerance //solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance relativeResidual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name() Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", Relative residual = " << relativeResidual << ", Relative residual = " << relativeResidual
<< ", No outer iterations " << iCorr << ", No outer iterations " << iCorr
<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
lduMatrix::debug=0; lduMatrix::debug=0;
# include "calculateDEpsilonDSigma.H" # include "calculateDEpsilonDSigma.H"
U += DU; U += DU;
sigma += DSigma; sigma += DSigma;
epsilon += DEpsilon; epsilon += DEpsilon;
# include "writeFields.H" # include "writeFields.H"
# include "calculateNetForces.H" # include "calculateNetForces.H"
Info<< "ExecutionTime = " Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime() << runTime.elapsedCpuTime()
<< " s\n\n" << endl; << " s\n\n" << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

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

View file

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

View file

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

View file

@ -1,33 +1,35 @@
{ {
forAll(mesh.boundary(), patchID) forAll(mesh.boundary(), patchID)
{ {
if(U.boundaryField()[patchID].type() if
== solidDirectionMixedFvPatchVectorField::typeName (
) U.boundaryField()[patchID].type()
{ == solidDirectionMixedFvPatchVectorField::typeName
solidDirectionMixedFvPatchVectorField& loadingPatch = )
refCast<solidDirectionMixedFvPatchVectorField> {
( solidDirectionMixedFvPatchVectorField& loadingPatch =
U.boundaryField()[patchID] refCast<solidDirectionMixedFvPatchVectorField>
); (
U.boundaryField()[patchID]
);
tensorField Finv = inv(I + gradU); tensorField Finv = inv(I + gradU);
vectorField newN = Finv & n.boundaryField()[patchID]; vectorField newN = Finv & n.boundaryField()[patchID];
newN /= mag(newN); newN /= mag(newN);
loadingPatch.valueFraction() = sqr(newN); loadingPatch.valueFraction() = sqr(newN);
//- set gradient //- set gradient
loadingPatch.refGrad() = loadingPatch.refGrad() =
( (
//Traction //Traction
( (mu.boundaryField()[patchID] + lambda.boundaryField()[patchID]) * (n.boundaryField()[patchID] & gradU.boundaryField()[patchID]) ) ( (mu.boundaryField()[patchID] + lambda.boundaryField()[patchID]) * (n.boundaryField()[patchID] & gradU.boundaryField()[patchID]) )
- ( mu.boundaryField()[patchID] * (n.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) ) - ( mu.boundaryField()[patchID] * (n.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) )
- ( mu.boundaryField()[patchID] * ( n.boundaryField()[patchID] & (gradU.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) ) ) - ( mu.boundaryField()[patchID] * ( n.boundaryField()[patchID] & (gradU.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) ) )
- ( lambda.boundaryField()[patchID] * tr(gradU.boundaryField()[patchID]) * n.boundaryField()[patchID] ) - ( lambda.boundaryField()[patchID] * tr(gradU.boundaryField()[patchID]) * n.boundaryField()[patchID] )
- ( 0.5 * lambda.boundaryField()[patchID] * tr(gradU.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) * n.boundaryField()[patchID] ) - ( 0.5 * lambda.boundaryField()[patchID] * tr(gradU.boundaryField()[patchID] & gradU.boundaryField()[patchID].T()) * n.boundaryField()[patchID] )
) )
/ /
(2.0*mu.boundaryField()[patchID] + lambda.boundaryField()[patchID]); (2.0*mu.boundaryField()[patchID] + lambda.boundaryField()[patchID]);
} }
} }
} }

View file

@ -25,7 +25,7 @@
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
symm(gradU) + 0.5*symm(gradU.T() & gradU) symm(gradU) + 0.5*symm(gradU.T() & gradU)
); );
//- second Piloa-Kirchhoff stress tensor //- second Piloa-Kirchhoff stress tensor

View file

@ -47,104 +47,105 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
while(runTime.loop()) while(runTime.loop())
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = 0; scalar initialResidual = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::debug=0; lduMatrix::debug=0;
do do
{ {
U.storePrevIter(); U.storePrevIter();
# include "correctDirectionMixedTL.H" # include "correctDirectionMixedTL.H"
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
fvm::d2dt2(rho, U) fvm::d2dt2(rho, U)
== ==
fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)") fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)")
+ fvc::div( + fvc::div
-( (mu + lambda) * gradU ) (
+ ( mu * gradU.T() ) - ( (mu + lambda) * gradU )
+ ( mu * (gradU & gradU.T()) ) + ( mu * gradU.T() )
+ ( lambda * tr(gradU) * I ) + ( mu * (gradU & gradU.T()) )
+ ( 0.5 * lambda * tr(gradU & gradU.T()) * I ) + ( lambda * tr(gradU) * I )
+ ( sigma & gradU ), + ( 0.5 * lambda * tr(gradU & gradU.T()) * I )
"div(sigma)" + ( sigma & gradU ),
) "div(sigma)"
); )
);
solverPerf = UEqn.solve(); solverPerf = UEqn.solve();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
U.relax(); U.relax();
gradU = fvc::grad(U); gradU = fvc::grad(U);
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << U.name() << ", Solving for " << U.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} }
while while
( (
solverPerf.initialResidual() > convergenceTolerance solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance //relativeResidual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() << ", Solving for " << U.name() Info << nl << "Time " << runTime.value() << ", Solving for " << U.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", Relative residual = " << relativeResidual << ", Relative residual = " << relativeResidual
<< ", No outer iterations " << iCorr << ", No outer iterations " << iCorr
<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
lduMatrix::debug=0; lduMatrix::debug=0;
# include "writeFields.H" # include "writeFields.H"
Info<< "ExecutionTime = " Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime() << runTime.elapsedCpuTime()
<< " s\n\n" << endl; << " s\n\n" << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

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

View file

@ -1,36 +1,36 @@
if (runTime.outputTime()) if (runTime.outputTime())
{ {
volScalarField epsilonEq volScalarField epsilonEq
( (
IOobject IOobject
( (
"epsilonEq", "epsilonEq",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
sqrt((2.0/3.0)*magSqr(dev(epsilon))) sqrt((2.0/3.0)*magSqr(dev(epsilon)))
); );
Info<< "Max epsilonEq = " << max(epsilonEq).value() Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl; << endl;
volScalarField sigmaEq volScalarField sigmaEq
( (
IOobject IOobject
( (
"sigmaEq", "sigmaEq",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
sqrt((3.0/2.0)*magSqr(dev(sigma))) sqrt((3.0/2.0)*magSqr(dev(sigma)))
); );
Info<< "Max sigmaEq = " << max(sigmaEq).value() Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl; << endl;
//- Calculate Cauchy stress //- Calculate Cauchy stress
volTensorField F = I + gradU; volTensorField F = I + gradU;
@ -40,34 +40,34 @@ if (runTime.outputTime())
rho = rho/J; rho = rho/J;
volSymmTensorField sigmaCauchy volSymmTensorField sigmaCauchy
( (
IOobject IOobject
( (
"sigmaCauchy", "sigmaCauchy",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
(1/J) * symm(F.T() & sigma & F) (1/J) * symm(F.T() & sigma & F)
); );
//- Cauchy von Mises stress //- Cauchy von Mises stress
volScalarField sigmaCauchyEq volScalarField sigmaCauchyEq
( (
IOobject IOobject
( (
"sigmaCauchyEq", "sigmaCauchyEq",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
sqrt((3.0/2.0)*magSqr(dev(sigmaCauchy))) sqrt((3.0/2.0)*magSqr(dev(sigmaCauchy)))
); );
Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value() Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value()
<< endl; << endl;
//- write boundary forces //- write boundary forces
//- integrate (sigma2PK & F) over reference area //- integrate (sigma2PK & F) over reference area
@ -75,30 +75,30 @@ if (runTime.outputTime())
//- over the deformed area //- over the deformed area
Info << nl; Info << nl;
forAll(mesh.boundary(), patchi) forAll(mesh.boundary(), patchi)
{ {
Info << "Patch " << mesh.boundary()[patchi].name() << endl; Info << "Patch " << mesh.boundary()[patchi].name() << endl;
tensorField F = I + gradU.boundaryField()[patchi]; tensorField F = I + gradU.boundaryField()[patchi];
vectorField totalForce = mesh.Sf().boundaryField()[patchi] & (sigma.boundaryField()[patchi] & F); vectorField totalForce = mesh.Sf().boundaryField()[patchi] & (sigma.boundaryField()[patchi] & F);
vector force = sum( totalForce ); vector force = sum( totalForce );
Info << "\ttotal force is " << force << " N" << endl; Info << "\ttotal force is " << force << " N" << endl;
tensorField Finv = inv(F); tensorField Finv = inv(F);
vectorField nCurrent = Finv & n.boundaryField()[patchi]; vectorField nCurrent = Finv & n.boundaryField()[patchi];
nCurrent /= mag(nCurrent); nCurrent /= mag(nCurrent);
scalar normalForce = sum( nCurrent & totalForce ); scalar normalForce = sum( nCurrent & totalForce );
Info << "\tnormal force is " << normalForce << " N" << endl; Info << "\tnormal force is " << normalForce << " N" << endl;
scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce )); scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce ));
Info << "\tshear force is " << shearForce << " N" << endl << endl;; Info << "\tshear force is " << shearForce << " N" << endl << endl;;
} }
//- move mesh for visualisation and move it back after writing //- move mesh for visualisation and move it back after writing
vectorField oldPoints = mesh.allPoints(); vectorField oldPoints = mesh.allPoints();
#include "moveMeshLeastSquares.H" # include "moveMeshLeastSquares.H"
runTime.write(); runTime.write();
//- move mesh back //- move mesh back
mesh.movePoints(oldPoints); mesh.movePoints(oldPoints);
} }

View file

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

View file

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

View file

@ -15,59 +15,58 @@
FieldField<Field, vector> extraVecs(ptc.size()); FieldField<Field, vector> extraVecs(ptc.size());
{ {
const labelListList& pfaces = mesh.pointFaces(); const labelListList& pfaces = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = pfaces[curPoint]; const labelList& curFaces = pfaces[curPoint];
// extraVecs.hook(new vectorField(curFaces.size())); //- no hook function // extraVecs.hook(new vectorField(curFaces.size())); //- no hook function
extraVecs.set extraVecs.set
( (
pointI, pointI,
new vectorField(curFaces.size()) new vectorField(curFaces.size())
); );
vectorField& curExtraVectors = extraVecs[pointI]; vectorField& curExtraVectors = extraVecs[pointI];
label nFacesAroundPoint = 0; label nFacesAroundPoint = 0;
const vector& pointLoc = mesh.points()[curPoint]; const vector& pointLoc = mesh.points()[curPoint];
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(bm[patchID]) !isA<emptyFvPatch>(bm[patchID])
&& !bm[patchID].coupled() && !bm[patchID].coupled()
) )
{ {
// Found a face for extrapolation // Found a face for extrapolation
curExtraVectors[nFacesAroundPoint] = curExtraVectors[nFacesAroundPoint] =
pointLoc pointLoc
- centres.boundaryField()[patchID] - centres.boundaryField()[patchID]
[bm[patchID].patch().whichFace(curFaces[faceI])]; [bm[patchID].patch().whichFace(curFaces[faceI])];
nFacesAroundPoint++; nFacesAroundPoint++;
} }
} }
} }
curExtraVectors.setSize(nFacesAroundPoint); curExtraVectors.setSize(nFacesAroundPoint);
} }
} }

View file

@ -8,114 +8,116 @@
FieldField<Field, scalar> w(ptc.size()); FieldField<Field, scalar> w(ptc.size());
{ {
const labelListList& pf = mesh.pointFaces(); const labelListList& pf = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
pointScalarField volPointSumWeights pointScalarField volPointSumWeights
( (
IOobject IOobject
( (
"volPointSumWeights", "volPointSumWeights",
mesh.polyMesh::instance(), mesh.polyMesh::instance(),
mesh mesh
), ),
pMesh, pMesh,
dimensionedScalar("zero", dimless, 0) 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)
forAll (ptc, pointI)
{ {
w[pointI] /= volPointSumWeights[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]];
} }
} }

View file

@ -91,32 +91,32 @@
//- explicit terms in the momentum equation //- explicit terms in the momentum equation
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero) dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero)
); );
volVectorField divDSigmaLargeStrainExp volVectorField divDSigmaLargeStrainExp
( (
IOobject IOobject
( (
"divDSigmaLargeStrainExp", "divDSigmaLargeStrainExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero) dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero)
); );
rheologyModel rheology(sigma); rheologyModel rheology(sigma);
volScalarField mu = rheology.mu(); volScalarField mu = rheology.mu();

View file

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

View file

@ -53,131 +53,131 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "readDivDSigmaLargeStrainExpMethod.H" # include "readDivDSigmaLargeStrainExpMethod.H"
# include "readMoveMeshMethod.H" # include "readMoveMeshMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info << "\nStarting time loop\n" << endl; Info << "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++) for (runTime++; !runTime.end(); runTime++)
{ {
Info<< "Time = " << runTime.timeName() << nl << endl; Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::debug = 0; lduMatrix::debug = 0;
do do
{
DU.storePrevIter();
divDSigmaLargeStrainExp.storePrevIter();
# include "calculateDivDSigmaExp.H"
# include "calculateDivDSigmaLargeStrainExp.H"
//----------------------------------------------------//
//- updated lagrangian large strain momentum equation
//----------------------------------------------------//
fvVectorMatrix DUEqn
(
fvm::d2dt2(rho,DU)
==
fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+ divDSigmaExp
+ divDSigmaLargeStrainExp
);
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 "calculateDEpsilonDSigma.H"
# include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << DU.name()
<< " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual()
<< ", residualDU = " << relativeResidual
<< ", inner iterations " << solverPerf.nIterations() << endl;
}
while
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&& ++iCorr < nCorr
);
lduMatrix::debug = 1;
Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl;
# include "rotateFields.H"
# include "moveMesh.H"
# include "writeFields.H"
//- total force
forAll(mesh.boundary(), patchi)
{ {
DU.storePrevIter();
divDSigmaLargeStrainExp.storePrevIter();
# include "calculateDivDSigmaExp.H"
# include "calculateDivDSigmaLargeStrainExp.H"
//----------------------------------------------------//
//- updated lagrangian large strain momentum equation
//----------------------------------------------------//
fvVectorMatrix DUEqn
(
fvm::d2dt2(rho,DU)
==
fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+ divDSigmaExp
+ divDSigmaLargeStrainExp
);
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 "calculateDEpsilonDSigma.H"
# include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << DU.name()
<< " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual()
<< ", residualDU = " << relativeResidual
<< ", inner iterations " << solverPerf.nIterations() << endl;
}
while
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&& ++iCorr < nCorr
);
lduMatrix::debug = 1;
Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl;
# include "rotateFields.H"
# include "moveMesh.H"
# include "writeFields.H"
//- total force
forAll(mesh.boundary(), patchi)
{
vector force = sum(mesh.Sf().boundaryField()[patchi] & sigma.boundaryField()[patchi]); vector force = sum(mesh.Sf().boundaryField()[patchi] & sigma.boundaryField()[patchi]);
Info << "force on " << mesh.boundary()[patchi].name() Info << "force on " << mesh.boundary()[patchi].name()
<< " is " << force << endl; << " is " << force << endl;
} }
Info << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" Info << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }
// ************************************************************************* // // ************************************************************************* //

View file

@ -14,10 +14,10 @@ forAll (bm, patchI)
( (
!isA<emptyFvPatch>(bm[patchI]) !isA<emptyFvPatch>(bm[patchI])
&& !( && !(
bm[patchI].coupled() bm[patchI].coupled()
//&& Pstream::parRun() //&& Pstream::parRun()
//&& !mesh.parallelData().cyclicParallel() //&& !mesh.parallelData().cyclicParallel()
) )
) )
{ {
const labelList& bp = bm[patchI].patch().boundaryPoints(); const labelList& bp = bm[patchI].patch().boundaryPoints();
@ -25,9 +25,9 @@ forAll (bm, patchI)
const labelList& meshPoints = bm[patchI].patch().meshPoints(); const labelList& meshPoints = bm[patchI].patch().meshPoints();
forAll (bp, pointI) forAll (bp, pointI)
{ {
pointsCorrectionMap.insert(meshPoints[bp[pointI]]); pointsCorrectionMap.insert(meshPoints[bp[pointI]]);
} }
} }
} }

View file

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

View file

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

View file

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

View file

@ -7,38 +7,40 @@ pointVectorField& pf = pointDU;
// Do the correction // Do the correction
//GeometricField<Type, pointPatchField, pointMesh> pfCorr //GeometricField<Type, pointPatchField, pointMesh> pfCorr
/*pointVectorField pfCorr /*
pointVectorField pfCorr
( (
IOobject IOobject
( (
// "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr", // "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
"edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr", "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
//vf.instance(), //vf.instance(),
DU, DU,
pMesh, pMesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
//dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero), //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
dimensionedVector("zero", pf.dimensions(), vector::zero), dimensionedVector("zero", pf.dimensions(), vector::zero),
pf.boundaryField().types() pf.boundaryField().types()
);*/ );
*/
pointVectorField pfCorr pointVectorField pfCorr
( (
IOobject IOobject
( (
"pointDUcorr", "pointDUcorr",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimLength, vector::zero), dimensionedVector("vector", dimLength, vector::zero),
"calculated" "calculated"
); );
//const labelList& ptc = boundaryPoints(); //const labelList& ptc = boundaryPoints();
#include "findBoundaryPoints.H" #include "findBoundaryPoints.H"
@ -57,64 +59,68 @@ const labelListList& PointFaces = mesh.pointFaces();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = PointFaces[curPoint]; const labelList& curFaces = PointFaces[curPoint];
label fI = 0; label fI = 0;
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(mesh.boundary()[patchID]) !isA<emptyFvPatch>(mesh.boundary()[patchID])
&& !mesh.boundary()[patchID].coupled() && !mesh.boundary()[patchID].coupled()
) )
{ {
label faceInPatchID = label faceInPatchID =
bm[patchID].patch().whichFace(curFaces[faceI]); bm[patchID].patch().whichFace(curFaces[faceI]);
pfCorr[curPoint] += pfCorr[curPoint] +=
w[pointI][fI]* w[pointI][fI]*
( (
extraVecs[pointI][fI] extraVecs[pointI][fI]
& gradDU.boundaryField()[patchID][faceInPatchID] & gradDU.boundaryField()[patchID][faceInPatchID]
); );
fI++; fI++;
} }
} }
} }
} }
// Update coupled boundaries // Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].initAddField(); pfCorr.boundaryField()[patchI].initAddField();
} }
}*/ }
*/
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].addField(pfCorr.internalField()); pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
} }
}*/ }
*/
//Info << "pfCorr: " << pfCorr << endl; //Info << "pfCorr: " << pfCorr << endl;
pfCorr.correctBoundaryConditions(); pfCorr.correctBoundaryConditions();
//pfCorr.write(); //pfCorr.write();

View file

@ -1,12 +1,15 @@
//- the method used to calculate the explicit component of sigma //- the method used to calculate the explicit component of sigma
word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp")); word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp"));
Info << "Calculation of divDSigmaExp method: " << divDSigmaExpMethod << endl; Info << "Calculation of divDSigmaExp method: " << divDSigmaExpMethod << endl;
if(divDSigmaExpMethod != "standard" if
&& divDSigmaExpMethod != "surface" (
&& divDSigmaExpMethod != "decompose" divDSigmaExpMethod != "standard"
&& divDSigmaExpMethod != "laplacian") && divDSigmaExpMethod != "surface"
{ && divDSigmaExpMethod != "decompose"
&& divDSigmaExpMethod != "laplacian"
)
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian" << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

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

View file

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

View file

@ -15,59 +15,58 @@
FieldField<Field, vector> extraVecs(ptc.size()); FieldField<Field, vector> extraVecs(ptc.size());
{ {
const labelListList& pfaces = mesh.pointFaces(); const labelListList& pfaces = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = pfaces[curPoint]; const labelList& curFaces = pfaces[curPoint];
// extraVecs.hook(new vectorField(curFaces.size())); //- no hook function // extraVecs.hook(new vectorField(curFaces.size())); //- no hook function
extraVecs.set extraVecs.set
( (
pointI, pointI,
new vectorField(curFaces.size()) new vectorField(curFaces.size())
); );
vectorField& curExtraVectors = extraVecs[pointI]; vectorField& curExtraVectors = extraVecs[pointI];
label nFacesAroundPoint = 0; label nFacesAroundPoint = 0;
const vector& pointLoc = mesh.points()[curPoint]; const vector& pointLoc = mesh.points()[curPoint];
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(bm[patchID]) !isA<emptyFvPatch>(bm[patchID])
&& !bm[patchID].coupled() && !bm[patchID].coupled()
) )
{ {
// Found a face for extrapolation // Found a face for extrapolation
curExtraVectors[nFacesAroundPoint] = curExtraVectors[nFacesAroundPoint] =
pointLoc pointLoc
- centres.boundaryField()[patchID] - centres.boundaryField()[patchID]
[bm[patchID].patch().whichFace(curFaces[faceI])]; [bm[patchID].patch().whichFace(curFaces[faceI])];
nFacesAroundPoint++; nFacesAroundPoint++;
} }
} }
} }
curExtraVectors.setSize(nFacesAroundPoint); curExtraVectors.setSize(nFacesAroundPoint);
} }
} }

View file

@ -8,114 +8,116 @@
FieldField<Field, scalar> w(ptc.size()); FieldField<Field, scalar> w(ptc.size());
{ {
const labelListList& pf = mesh.pointFaces(); const labelListList& pf = mesh.pointFaces();
const volVectorField& centres = mesh.C(); const volVectorField& centres = mesh.C();
const fvBoundaryMesh& bm = mesh.boundary(); const fvBoundaryMesh& bm = mesh.boundary();
pointScalarField volPointSumWeights pointScalarField volPointSumWeights
( (
IOobject IOobject
( (
"volPointSumWeights", "volPointSumWeights",
mesh.polyMesh::instance(), mesh.polyMesh::instance(),
mesh mesh
), ),
pMesh, pMesh,
dimensionedScalar("zero", dimless, 0) 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)
forAll (ptc, pointI)
{ {
w[pointI] /= volPointSumWeights[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]];
} }
} }

View file

@ -57,33 +57,33 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero) dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero)
); );
volVectorField divDSigmaLargeStrainExp volVectorField divDSigmaLargeStrainExp
( (
IOobject IOobject
( (
"divDSigmaLargeStrainExp", "divDSigmaLargeStrainExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero) dimensionedVector("zero", dimensionSet(1, -2, -2, 0, 0, 0, 0), vector::zero)
); );
volSymmTensorField epsilon volSymmTensorField epsilon
( (

View file

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

View file

@ -61,151 +61,151 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "readDivDSigmaLargeStrainExpMethod.H" # include "readDivDSigmaLargeStrainExpMethod.H"
# include "readMoveMeshMethod.H" # include "readMoveMeshMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info << "\nStarting time loop\n" << endl; Info << "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++) for (runTime++; !runTime.end(); runTime++)
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::debug = 0; lduMatrix::debug = 0;
const volSymmTensorField& DEpsilonP = rheology.DEpsilonP(); const volSymmTensorField& DEpsilonP = rheology.DEpsilonP();
do do
{ {
DU.storePrevIter(); DU.storePrevIter();
divDSigmaLargeStrainExp.storePrevIter(); divDSigmaLargeStrainExp.storePrevIter();
# include "calculateDivDSigmaExp.H" # include "calculateDivDSigmaExp.H"
# include "calculateDivDSigmaLargeStrainExp.H" # include "calculateDivDSigmaLargeStrainExp.H"
//----------------------------------------------------// //----------------------------------------------------//
//- updated lagrangian large strain momentum equation //- updated lagrangian large strain momentum equation
//----------------------------------------------------// //----------------------------------------------------//
fvVectorMatrix DUEqn fvVectorMatrix DUEqn
( (
fvm::d2dt2(rho, DU) fvm::d2dt2(rho, DU)
== ==
fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)") fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+ divDSigmaExp + divDSigmaExp
+ divDSigmaLargeStrainExp + divDSigmaLargeStrainExp
- fvc::div(2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP))) - fvc::div(2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP)))
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->correct(DUEqn); solidInterfacePtr->correct(DUEqn);
} }
solverPerf = DUEqn.solve(); solverPerf = DUEqn.solve();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
DU.relax(); DU.relax();
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
gradDU = solidInterfacePtr->grad(DU); gradDU = solidInterfacePtr->grad(DU);
} }
else else
{ {
gradDU = fvc::grad(DU); gradDU = fvc::grad(DU);
} }
DF = gradDU.T(); DF = gradDU.T();
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
rheology.correct(); rheology.correct();
mu = rheology.newMu(); mu = rheology.newMu();
lambda = rheology.newLambda(); lambda = rheology.newLambda();
muf = fvc::interpolate(rheology.newMu()); muf = fvc::interpolate(rheology.newMu());
lambdaf = fvc::interpolate(rheology.newLambda()); lambdaf = fvc::interpolate(rheology.newLambda());
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->modifyProperties(muf, lambdaf); solidInterfacePtr->modifyProperties(muf, lambdaf);
} }
# include "calculateDEpsilonDSigma.H" # include "calculateDEpsilonDSigma.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << DU.name() << ", Solving for " << DU.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} }
while while
( (
//relativeResidual //relativeResidual
solverPerf.initialResidual() > convergenceTolerance solverPerf.initialResidual() > convergenceTolerance
&& ++iCorr < nCorr && ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name() Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl; << ", No outer iterations " << iCorr << endl;
lduMatrix::debug = 1; lduMatrix::debug = 1;
U += DU; U += DU;
epsilon += DEpsilon; epsilon += DEpsilon;
epsilonP += DEpsilonP; epsilonP += DEpsilonP;
volSymmTensorField DEpsilonE = DEpsilon - DEpsilonP; volSymmTensorField DEpsilonE = DEpsilon - DEpsilonP;
epsilonE += DEpsilonE; epsilonE += DEpsilonE;
sigma += DSigma; sigma += DSigma;
rheology.updateYieldStress(); rheology.updateYieldStress();
# include "rotateFields.H" # include "rotateFields.H"
# include "moveMesh.H" # include "moveMesh.H"
# include "writeFields.H" # include "writeFields.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl; << nl << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

@ -8,26 +8,26 @@ const fvBoundaryMesh& bm = mesh.boundary();
forAll (bm, patchI) forAll (bm, patchI)
{ {
// If the patch is empty, skip it // If the patch is empty, skip it
// If the patch is coupled, and there are no cyclic parallels, skip it // If the patch is coupled, and there are no cyclic parallels, skip it
if if
( (
!isA<emptyFvPatch>(bm[patchI]) !isA<emptyFvPatch>(bm[patchI])
&& !( && !(
bm[patchI].coupled() bm[patchI].coupled()
//&& Pstream::parRun() //&& Pstream::parRun()
//&& !mesh.parallelData().cyclicParallel() //&& !mesh.parallelData().cyclicParallel()
) )
) )
{ {
const labelList& bp = bm[patchI].patch().boundaryPoints(); const labelList& bp = bm[patchI].patch().boundaryPoints();
const labelList& meshPoints = bm[patchI].patch().meshPoints(); const labelList& meshPoints = bm[patchI].patch().meshPoints();
forAll (bp, pointI) forAll (bp, pointI)
{ {
pointsCorrectionMap.insert(meshPoints[bp[pointI]]); pointsCorrectionMap.insert(meshPoints[bp[pointI]]);
} }
} }
} }

View file

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

View file

@ -2,7 +2,7 @@
//- move mesh //- move mesh
//--------------------------------------------------// //--------------------------------------------------//
if(min(J.internalField()) > 0) if(min(J.internalField()) > 0)
{ {
Info << "Move solid mesh using inverse distance interpolation" << endl; Info << "Move solid mesh using inverse distance interpolation" << endl;
// Create point mesh // Create point mesh
@ -12,24 +12,24 @@ if(min(J.internalField()) > 0)
volPointInterpolation pointInterpolation(mesh); volPointInterpolation pointInterpolation(mesh);
wordList types wordList types
( (
pMesh.boundary().size(), pMesh.boundary().size(),
//fixedValueFvPatchVectorField::typeName //fixedValueFvPatchVectorField::typeName
calculatedFvPatchVectorField::typeName calculatedFvPatchVectorField::typeName
); );
pointVectorField pointDU pointVectorField pointDU
( (
IOobject IOobject
( (
"pointDU", "pointDU",
runTime.timeName(), runTime.timeName(),
mesh mesh
), ),
pMesh, pMesh,
dimensionedVector("zero", dimLength, vector::zero), dimensionedVector("zero", dimLength, vector::zero),
types types
); );
// Calculate mesh points displacement // Calculate mesh points displacement
pointInterpolation.interpolate(DU, pointDU); pointInterpolation.interpolate(DU, pointDU);
@ -41,26 +41,25 @@ if(min(J.internalField()) > 0)
//pointDU.write(); //pointDU.write();
const vectorField& pointDUI = const vectorField& pointDUI = pointDU.internalField();
pointDU.internalField();
// Move mesh // Move mesh
vectorField newPoints = mesh.allPoints(); vectorField newPoints = mesh.allPoints();
forAll (pointDUI, pointI) forAll (pointDUI, pointI)
{ {
newPoints[pointI] += pointDUI[pointI]; newPoints[pointI] += pointDUI[pointI];
} }
twoDPointCorrector twoDCorrector(mesh); twoDPointCorrector twoDCorrector(mesh);
twoDCorrector.correctPoints(newPoints); twoDCorrector.correctPoints(newPoints);
mesh.movePoints(newPoints); mesh.movePoints(newPoints);
mesh.V00(); mesh.V00();
mesh.moving(false); mesh.moving(false);
} }
else else
{ {
FatalErrorIn(args.executable()) FatalErrorIn(args.executable())
<< "Negative Jacobian" << "Negative Jacobian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

@ -7,38 +7,40 @@ pointVectorField& pf = pointDU;
// Do the correction // Do the correction
//GeometricField<Type, pointPatchField, pointMesh> pfCorr //GeometricField<Type, pointPatchField, pointMesh> pfCorr
/*pointVectorField pfCorr /*
pointVectorField pfCorr
( (
IOobject IOobject
( (
// "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr", // "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
"edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr", "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
//vf.instance(), //vf.instance(),
DU, DU,
pMesh, pMesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
//dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero), //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
dimensionedVector("zero", pf.dimensions(), vector::zero), dimensionedVector("zero", pf.dimensions(), vector::zero),
pf.boundaryField().types() pf.boundaryField().types()
);*/ );
*/
pointVectorField pfCorr pointVectorField pfCorr
( (
IOobject IOobject
( (
"pointDUcorr", "pointDUcorr",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
pMesh, pMesh,
dimensionedVector("vector", dimLength, vector::zero), dimensionedVector("vector", dimLength, vector::zero),
"calculated" "calculated"
); );
//const labelList& ptc = boundaryPoints(); //const labelList& ptc = boundaryPoints();
#include "findBoundaryPoints.H" #include "findBoundaryPoints.H"
@ -57,64 +59,68 @@ const labelListList& PointFaces = mesh.pointFaces();
forAll (ptc, pointI) forAll (ptc, pointI)
{ {
const label curPoint = ptc[pointI]; const label curPoint = ptc[pointI];
const labelList& curFaces = PointFaces[curPoint]; const labelList& curFaces = PointFaces[curPoint];
label fI = 0; label fI = 0;
// Go through all the faces // Go through all the faces
forAll (curFaces, faceI) forAll (curFaces, faceI)
{ {
if (!mesh.isInternalFace(curFaces[faceI])) if (!mesh.isInternalFace(curFaces[faceI]))
{ {
// This is a boundary face. If not in the empty patch // This is a boundary face. If not in the empty patch
// or coupled calculate the extrapolation vector // or coupled calculate the extrapolation vector
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(curFaces[faceI]); mesh.boundaryMesh().whichPatch(curFaces[faceI]);
if if
( (
!isA<emptyFvPatch>(mesh.boundary()[patchID]) !isA<emptyFvPatch>(mesh.boundary()[patchID])
&& !mesh.boundary()[patchID].coupled() && !mesh.boundary()[patchID].coupled()
) )
{ {
label faceInPatchID = label faceInPatchID =
bm[patchID].patch().whichFace(curFaces[faceI]); bm[patchID].patch().whichFace(curFaces[faceI]);
pfCorr[curPoint] += pfCorr[curPoint] +=
w[pointI][fI]* w[pointI][fI]*
( (
extraVecs[pointI][fI] extraVecs[pointI][fI]
& gradDU.boundaryField()[patchID][faceInPatchID] & gradDU.boundaryField()[patchID][faceInPatchID]
); );
fI++; fI++;
} }
} }
} }
} }
// Update coupled boundaries // Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].initAddField(); pfCorr.boundaryField()[patchI].initAddField();
} }
}*/ }
*/
/*forAll (pfCorr.boundaryField(), patchI) /*
forAll (pfCorr.boundaryField(), patchI)
{ {
if (pfCorr.boundaryField()[patchI].coupled()) if (pfCorr.boundaryField()[patchI].coupled())
{ {
pfCorr.boundaryField()[patchI].addField(pfCorr.internalField()); pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
} }
}*/ }
*/
//Info << "pfCorr: " << pfCorr << endl; //Info << "pfCorr: " << pfCorr << endl;
pfCorr.correctBoundaryConditions(); pfCorr.correctBoundaryConditions();
//pfCorr.write(); //pfCorr.write();

View file

@ -1,12 +1,15 @@
//- the method used to calculate the explicit component of sigma //- the method used to calculate the explicit component of sigma
word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp")); word divDSigmaExpMethod(mesh.solutionDict().subDict("stressedFoam").lookup("divDSigmaExp"));
Info << "Calculation of divDSigmaExp method: " << divDSigmaExpMethod << endl; Info << "Calculation of divDSigmaExp method: " << divDSigmaExpMethod << endl;
if(divDSigmaExpMethod != "standard" if
&& divDSigmaExpMethod != "surface" (
&& divDSigmaExpMethod != "decompose" divDSigmaExpMethod != "standard"
&& divDSigmaExpMethod != "laplacian") && divDSigmaExpMethod != "surface"
{ && divDSigmaExpMethod != "decompose"
&& divDSigmaExpMethod != "laplacian"
)
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian" << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

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

View file

@ -33,15 +33,15 @@
( (
IOobject IOobject
( (
"DEpsilon", "DEpsilon",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::READ_IF_PRESENT, IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
mesh, mesh,
dimensionedSymmTensor("zero", dimless, symmTensor::zero) dimensionedSymmTensor("zero", dimless, symmTensor::zero)
); );
volSymmTensorField DSigma volSymmTensorField DSigma
( (
@ -100,19 +100,19 @@
); );
volVectorField divDSigmaExp volVectorField divDSigmaExp
( (
IOobject IOobject
( (
"divDSigmaExp", "divDSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
plasticityModel rheology(gradDU, epsilon, sigma); plasticityModel rheology(gradDU, epsilon, sigma);

View file

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

View file

@ -50,127 +50,127 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivDSigmaExpMethod.H" # include "readDivDSigmaExpMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
for (runTime++; !runTime.end(); runTime++) for (runTime++; !runTime.end(); runTime++)
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = 0; scalar initialResidual = 0;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
lduMatrix::debug = 0; lduMatrix::debug = 0;
const volSymmTensorField& DEpsilonP = rheology.DEpsilonP(); const volSymmTensorField& DEpsilonP = rheology.DEpsilonP();
do do
{ {
DU.storePrevIter(); DU.storePrevIter();
# include "calculateDivDSigmaExp.H" # include "calculateDivDSigmaExp.H"
fvVectorMatrix DUEqn fvVectorMatrix DUEqn
( (
fvm::d2dt2(rho, DU) fvm::d2dt2(rho, DU)
== ==
fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)") fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
+ divDSigmaExp + divDSigmaExp
- fvc::div(2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP))) - fvc::div(2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP)))
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{
solidInterfacePtr->correct(DUEqn);
}
solverPerf = DUEqn.solve();
if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); solidInterfacePtr->correct(DUEqn);
} }
DU.relax(); solverPerf = DUEqn.solve();
if(solidInterfaceCorr) if(iCorr == 0)
{ {
gradDU = solidInterfacePtr->grad(DU); initialResidual = solverPerf.initialResidual();
} }
else
{
gradDU = fvc::grad(DU);
}
# include "calculateRelativeResidual.H" DU.relax();
rheology.correct(); if(solidInterfaceCorr)
mu = rheology.newMu(); {
lambda = rheology.newLambda(); gradDU = solidInterfacePtr->grad(DU);
muf = fvc::interpolate(rheology.newMu()); }
lambdaf = fvc::interpolate(rheology.newLambda()); else
if(solidInterfaceCorr) {
{ gradDU = fvc::grad(DU);
solidInterfacePtr->modifyProperties(muf, lambdaf); }
}
Info << "\tTime " << runTime.value() # include "calculateRelativeResidual.H"
<< ", Corrector " << iCorr
<< ", Solving for " << DU.name()
<< " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << endl;
}
while
(
solverPerf.initialResidual() > convergenceTolerance
&& ++iCorr < nCorr
);
Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name() rheology.correct();
<< ", Initial residual = " << initialResidual mu = rheology.newMu();
<< ", Final residual = " << solverPerf.initialResidual() lambda = rheology.newLambda();
<< ", No outer iterations " << iCorr << endl; muf = fvc::interpolate(rheology.newMu());
lambdaf = fvc::interpolate(rheology.newLambda());
if(solidInterfaceCorr)
{
solidInterfacePtr->modifyProperties(muf, lambdaf);
}
lduMatrix::debug = 1; Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << DU.name()
<< " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << endl;
}
while
(
solverPerf.initialResidual() > convergenceTolerance
&& ++iCorr < nCorr
);
# include "calculateDEpsilonDSigma.H" Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl;
U += DU; lduMatrix::debug = 1;
epsilon += DEpsilon; # include "calculateDEpsilonDSigma.H"
epsilonP += rheology.DEpsilonP(); U += DU;
sigma += DSigma; epsilon += DEpsilon;
rheology.updateYieldStress(); epsilonP += rheology.DEpsilonP();
# include "writeFields.H" sigma += DSigma;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" rheology.updateYieldStress();
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl; # include "writeFields.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

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

View file

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

View file

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

View file

@ -42,19 +42,19 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
); );
volVectorField divSigmaExp volVectorField divSigmaExp
( (
IOobject IOobject
( (
"divSigmaExp", "divSigmaExp",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
); );
rheologyModel rheology(sigma); rheologyModel rheology(sigma);

View file

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

View file

@ -52,111 +52,111 @@ int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivSigmaExpMethod.H" # include "readDivSigmaExpMethod.H"
# include "createSolidInterface.H" # include "createSolidInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
while(runTime.loop()) while(runTime.loop())
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = 0; scalar initialResidual = 0;
lduMatrix::solverPerformance solverPerf; lduMatrix::solverPerformance solverPerf;
scalar relativeResidual = GREAT; scalar relativeResidual = GREAT;
lduMatrix::debug=0; lduMatrix::debug=0;
do do
{ {
U.storePrevIter(); U.storePrevIter();
# include "calculateDivSigmaExp.H" # include "calculateDivSigmaExp.H"
//- linear momentum equation //- linear momentum equation
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
fvm::d2dt2(rho, U) fvm::d2dt2(rho, U)
== ==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)") fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+ divSigmaExp + divSigmaExp
); );
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
solidInterfacePtr->correct(UEqn); solidInterfacePtr->correct(UEqn);
} }
solverPerf = UEqn.solve(); solverPerf = UEqn.solve();
if(iCorr == 0) if(iCorr == 0)
{ {
initialResidual = solverPerf.initialResidual(); initialResidual = solverPerf.initialResidual();
} }
U.relax(); U.relax();
if(solidInterfaceCorr) if(solidInterfaceCorr)
{ {
gradU = solidInterfacePtr->grad(U); gradU = solidInterfacePtr->grad(U);
} }
else else
{ {
gradU = fvc::grad(U); gradU = fvc::grad(U);
} }
# include "calculateRelativeResidual.H" # include "calculateRelativeResidual.H"
Info << "\tTime " << runTime.value() Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << ", Corrector " << iCorr
<< ", Solving for " << U.name() << ", Solving for " << U.name()
<< " using " << solverPerf.solverName() << " using " << solverPerf.solverName()
<< ", residual = " << solverPerf.initialResidual() << ", residual = " << solverPerf.initialResidual()
<< ", relative residual = " << relativeResidual << endl; << ", relative residual = " << relativeResidual << endl;
} }
while while
( (
//solverPerf.initialResidual() > convergenceTolerance //solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance relativeResidual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() << ", Solving for " << U.name() Info << nl << "Time " << runTime.value() << ", Solving for " << U.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual() << ", Final residual = " << solverPerf.initialResidual()
<< ", Relative residual = " << relativeResidual << ", Relative residual = " << relativeResidual
<< ", No outer iterations " << iCorr << ", No outer iterations " << iCorr
<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
lduMatrix::debug=0; lduMatrix::debug=0;
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"
# include "writeFields.H" # include "writeFields.H"
Info<< "ExecutionTime = " Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime() << runTime.elapsedCpuTime()
<< " s\n\n" << endl; << " s\n\n" << endl;
} }
Info<< "End\n" << endl; Info<< "End\n" << endl;
return(0); return(0);
} }

View file

@ -4,6 +4,6 @@ Info << "Selecting divSigmaExp calculation method " << divSigmaExpMethod << end
if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian") if(divSigmaExpMethod != "standard" && divSigmaExpMethod != "surface" && divSigmaExpMethod != "decompose" && divSigmaExpMethod != "laplacian")
{ {
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian" << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError); << exit(FatalError);
} }

View file

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

View file

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

View file

@ -45,110 +45,110 @@ Description
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readSigmaExpMethod.H" # include "readSigmaExpMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating displacement field\n" << endl; Info<< "\nCalculating displacement field\n" << endl;
while(runTime.loop()) while(runTime.loop())
{ {
Info<< "Time: " << runTime.timeName() << nl << endl; Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readStressedFoamControls.H" # include "readStressedFoamControls.H"
int iCorr = 0; int iCorr = 0;
scalar initialResidual = GREAT; scalar initialResidual = GREAT;
scalar residual = GREAT; scalar residual = GREAT;
lduMatrix::solverPerformance solverPerfU; lduMatrix::solverPerformance solverPerfU;
lduMatrix::solverPerformance solverPerfT; lduMatrix::solverPerformance solverPerfT;
lduMatrix::debug=0; lduMatrix::debug=0;
do do
{ {
U.storePrevIter(); U.storePrevIter();
# include "calculateSigmaExp.H" # include "calculateSigmaExp.H"
//- energy equation //- energy equation
fvScalarMatrix TEqn fvScalarMatrix TEqn
(
fvm::ddt(rhoC, T) == fvm::laplacian(k, T, "laplacian(k,T)")
);
solverPerfT = TEqn.solve();
T.relax();
Info << "\tTime " << runTime.value()
<< ", Corrector " << iCorr << nl
<< "\t\tSolving for " << T.name()
<< " using " << solverPerfT.solverName()
<< ", residual = " << solverPerfT.initialResidual() << endl;
//- linear momentum equaiton
fvVectorMatrix UEqn
( (
fvm::d2dt2(rho, U) fvm::ddt(rhoC, T) == fvm::laplacian(k, T, "laplacian(k,T)")
== );
fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)")
+ sigmaExp
- fvc::grad(threeKalpha*(T-T0),"grad(threeKalphaDeltaT)")
);
solverPerfU = UEqn.solve(); solverPerfT = TEqn.solve();
if(iCorr == 0) T.relax();
{
initialResidual = max
(
solverPerfU.initialResidual(),
solverPerfT.initialResidual()
);
}
residual = max Info << "\tTime " << runTime.value()
( << ", Corrector " << iCorr << nl
solverPerfU.initialResidual(), << "\t\tSolving for " << T.name()
solverPerfT.initialResidual() << " using " << solverPerfT.solverName()
); << ", residual = " << solverPerfT.initialResidual() << endl;
U.relax(); //- linear momentum equaiton
fvVectorMatrix UEqn
(
fvm::d2dt2(rho, U)
==
fvm::laplacian(2*mu + lambda, U, "laplacian(DU,U)")
+ sigmaExp
- fvc::grad(threeKalpha*(T-T0),"grad(threeKalphaDeltaT)")
);
gradU = fvc::grad(U); solverPerfU = UEqn.solve();
Info << "\t\tSolving for " << U.name() if(iCorr == 0)
<< " using " << solverPerfU.solverName() {
<< ", residual = " << solverPerfU.initialResidual() << endl; initialResidual = max
(
solverPerfU.initialResidual(),
solverPerfT.initialResidual()
);
}
residual = max
(
solverPerfU.initialResidual(),
solverPerfT.initialResidual()
);
U.relax();
gradU = fvc::grad(U);
Info << "\t\tSolving for " << U.name()
<< " using " << solverPerfU.solverName()
<< ", residual = " << solverPerfU.initialResidual() << endl;
} }
while while
( (
residual > convergenceTolerance residual > convergenceTolerance
&& &&
++iCorr < nCorr ++iCorr < nCorr
); );
Info << nl << "Time " << runTime.value() Info << nl << "Time " << runTime.value()
<< ", Solving for " << U.name() << ", Solving for " << U.name()
<< ", Solving for " << T.name() << ", Solving for " << T.name()
<< ", Initial residual = " << initialResidual << ", Initial residual = " << initialResidual
<< ", Final U residual = " << solverPerfU.initialResidual() << ", Final U residual = " << solverPerfU.initialResidual()
<< ", Final T residual = " << solverPerfT.initialResidual() << ", Final T residual = " << solverPerfT.initialResidual()
<< ", No outer iterations " << iCorr << ", No outer iterations " << iCorr
<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
lduMatrix::debug=0; lduMatrix::debug=0;
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"

View file

@ -4,6 +4,6 @@ Info << sigmaExpMethod << " method chosen for calculation of sigmaExp" << endl;
if(sigmaExpMethod != "standard" && sigmaExpMethod != "surface" && sigmaExpMethod != "decompose") if(sigmaExpMethod != "standard" && sigmaExpMethod != "surface" && sigmaExpMethod != "decompose")
{ {
FatalError << "sigmaExp method " << sigmaExpMethod << " not found!" << nl FatalError << "sigmaExp method " << sigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose" << "valid methods are:\nstandard\nsurface\ndecompose"
<< exit(FatalError); << exit(FatalError);
} }

Some files were not shown because too many files have changed in this diff Show more