Fixing indentation in applications/solvers/solidMechanics

This commit is contained in:
Henrik Rusche 2015-05-17 16:25:35 +02:00
parent b46695ce1e
commit c1cd77a15f
138 changed files with 5164 additions and 5127 deletions

View file

@ -22,9 +22,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest U

View file

@ -19,13 +19,12 @@ if(divSigmaExpMethod == "standard")
{
snGradU = fvc::snGrad(U);
surfaceTensorField shearGradU =
((I - n*n)&fvc::interpolate(gradU));
surfaceTensorField shearGradU = ((I - n*n) & fvc::interpolate(gradU));
divSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(snGradU & (I - n*n))
+ lambdaf*tr(shearGradU & (I - n*n))*n
+ muf*(shearGradU&n)
@ -36,14 +35,10 @@ if(divSigmaExpMethod == "standard")
{
divSigmaExp =
- fvc::laplacian(mu + lambda, U, "laplacian(DU,U)")
+ fvc::div
(
mu*gradU.T()
+ lambda*(I*tr(gradU)),
"div(sigma)"
);
+ fvc::div(mu*gradU.T() + lambda*(I*tr(gradU)), "div(sigma)");
}
else
{
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
}

View file

@ -5,8 +5,8 @@
vector netForce = vector::zero;
forAll(mesh.boundary(), patchi)
{
netForce +=
sum(
netForce += sum
(
mesh.Sf().boundaryField()[patchi]
&
(

View file

@ -44,7 +44,6 @@
// );
// philipc: I have moved cohesive stuff to constitutiveModel
// cohesiveZone is an index field
// which allows the user to limit the crack to certain areas at runtime
// 1 for faces within cohesiveZone
@ -124,6 +123,7 @@
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
forAll(cohesiveZone.internalField(), facei)
{
if(cohesiveZone.internalField()[facei])
@ -132,6 +132,7 @@
cohesiveZoneVol.internalField()[mesh.neighbour()[facei]] = 1.0;
}
}
forAll(cohesiveZone.boundaryField(), patchi)
{
forAll(cohesiveZone.boundaryField()[patchi], facei)

View file

@ -86,7 +86,7 @@ int main(int argc, char *argv[])
runTime++;
Info<< "\nTime: " << runTime.timeName() << " s\n" << endl;
Info << "Time = " << runTime.timeName() << nl << endl;
volScalarField rho = rheology.rho();
volScalarField mu = rheology.mu();
@ -196,8 +196,7 @@ int main(int argc, char *argv[])
(
solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&&
iCorr < nCorr
&& iCorr < nCorr
)
);

View file

@ -41,6 +41,7 @@ nCoupledFacesToBreak = 0;
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+ (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
}
maxEffTractionFraction = gMax(effTractionFraction);
SLList<label> facesToBreakList;
@ -77,9 +78,14 @@ nCoupledFacesToBreak = 0;
scalar faceToBreakEffTractionFraction = 0;
forAll(facesToBreakEffTractionFraction, faceI)
{
if (facesToBreakEffTractionFraction[faceI] > faceToBreakEffTractionFraction)
if
(
facesToBreakEffTractionFraction[faceI]
> faceToBreakEffTractionFraction
)
{
faceToBreakEffTractionFraction = facesToBreakEffTractionFraction[faceI];
faceToBreakEffTractionFraction =
facesToBreakEffTractionFraction[faceI];
faceToBreakIndex = facesToBreak[faceI];
}
}
@ -92,7 +98,11 @@ nCoupledFacesToBreak = 0;
bool procHasFaceToBreak = false;
if (nFacesToBreak > 0)
{
if ( mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction) < SMALL )
if
(
mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
< SMALL
)
{
// philipc - Maximum traction fraction is on this processor
procHasFaceToBreak = true;
@ -132,7 +142,10 @@ nCoupledFacesToBreak = 0;
scalarField pNormalTraction =
cohesiveZone.boundaryField()[patchI]*
( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
pNormalTraction = max(pNormalTraction, scalar(0)); // only consider tensile tractions
// only consider tensile tractions
pNormalTraction = max(pNormalTraction, scalar(0));
scalarField pShearTraction =
cohesiveZone.boundaryField()[patchI]*
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
@ -146,13 +159,15 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchUPtr)
{
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+ (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+ (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
}
label start = mesh.boundaryMesh()[patchI].start();
@ -318,6 +333,7 @@ nCoupledFacesToBreak = 0;
vector faceToBreakNormal = vector::zero;
scalar faceToBreakSigmaMax = 0.0;
scalar faceToBreakTauMax = 0.0;
// Set faces to break
if (nFacesToBreak > 0)
{
@ -334,19 +350,27 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchUPtr)
{
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) );
)
);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) );
)
);
}
faceToBreakTraction *= scaleFactor;
@ -373,19 +397,27 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchUPtr)
{
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) );
)
);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) );
)
);
}
faceToBreakTraction *= scaleFactor;
@ -549,7 +581,6 @@ nCoupledFacesToBreak = 0;
{
traction.boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
traction.oldTime().boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
@ -560,8 +591,10 @@ nCoupledFacesToBreak = 0;
}
else
{
cohesivePatchUFixedModePtr->traction()[i] = -faceToBreakTraction;
cohesivePatchUFixedModePtr->initiationTraction()[i] = -faceToBreakTraction;
cohesivePatchUFixedModePtr->traction()[i] =
-faceToBreakTraction;
cohesivePatchUFixedModePtr->initiationTraction()[i] =
-faceToBreakTraction;
}
}
}

View file

@ -12,7 +12,8 @@
// so a quick fix is to not set a reference on regions
// with a processor boundary
//if (U.boundaryField()[patchI].fixesValue())
if (
if
(
U.boundaryField()[patchI].fixesValue()
||
mesh.boundaryMesh()[patchI].type()

View file

@ -98,6 +98,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
volScalarField GI
(
IOobject
@ -112,6 +113,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
volScalarField GII
(
IOobject
@ -126,6 +128,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
forAll(U.boundaryField(), patchi)
{
// if(U.boundaryField()[patchi].type() == cohesiveLawMultiMatFvPatchVectorField::typeName)

View file

@ -19,8 +19,8 @@ if(historyPatchID != -1)
//- be dotted with the surface normal to give the actual traction/force
//- you cannot just take the component of the sigma tensor
//scalar forcePatchIntegrateMethod = gSum(
// mesh.magSf().boundaryField()[historyPatchID]
// *sigma.boundaryField()[historyPatchID].component(symmTensor::XY)
// mesh.magSf().boundaryField()[historyPatchID]*
// sigma.boundaryField()[historyPatchID].component(symmTensor::XY)
//);
vector avDisp = gAverage(U.boundaryField()[historyPatchID]);

View file

@ -22,9 +22,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest DU

View file

@ -19,13 +19,12 @@ if(divDSigmaExpMethod == "standard")
{
snGradDU = fvc::snGrad(DU);
surfaceTensorField shearGradDU =
((I - n*n)&fvc::interpolate(gradDU));
surfaceTensorField shearGradDU = ((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(snGradDU & (I - n*n))
+ lambdaf*tr(shearGradDU & (I - n*n))*n
+ muf*(shearGradDU & n)
@ -36,14 +35,11 @@ if(divDSigmaExpMethod == "standard")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -5,8 +5,8 @@
vector netForce = vector::zero;
forAll(mesh.boundary(), patchi)
{
netForce +=
sum(
netForce += sum
(
mesh.Sf().boundaryField()[patchi]
&
(

View file

@ -59,7 +59,6 @@
IOobject::AUTO_WRITE
),
mesh,
//dimensionedScalar("one", dimless, 1.0)
dimensionedScalar("zero", dimless, 0.0)
);
@ -143,6 +142,7 @@
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
forAll(cohesiveZone.internalField(), facei)
{
if(cohesiveZone.internalField()[facei])
@ -151,11 +151,12 @@
cohesiveZoneVol.internalField()[mesh.neighbour()[facei]] = 1.0;
}
}
forAll(cohesiveZone.boundaryField(), patchi)
{
forAll(cohesiveZone.boundaryField()[patchi], facei)
{
if(cohesiveZone.boundaryField()[patchi][facei])
if(cohesiveZone.boundaryField()[patchi][facei] > 0.0)
{
cohesiveZoneVol.boundaryField()[patchi][facei] = 1.0;
}

View file

@ -157,6 +157,7 @@ constitutiveModel rheology(sigma, DU);
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
// aitken relaxation factor
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

View file

@ -34,6 +34,7 @@ Author
Philip Cardiff UCD
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "constitutiveModel.H"
//#include "componentReferenceList.H"
@ -73,7 +74,7 @@ int main(int argc, char *argv[])
runTime++;
Info<< "\nTime: " << runTime.timeName() << " s\n" << endl;
Info<< "\nTime = " << runTime.timeName() << " s\n" << endl;
volScalarField rho = rheology.rho();
volScalarField mu = rheology.mu();
@ -211,7 +212,8 @@ int main(int argc, char *argv[])
returnReduce
(
cohesivePatchDUFixedModePtr->size(),
sumOp<label>())
sumOp<label>()
)
)
{
Pout << "Number of faces in crack: "

View file

@ -12,13 +12,13 @@ nCoupledFacesToBreak = 0;
// only consider tensile tractions
normalTraction = max(normalTraction, scalar(0));
scalarField shearTraction =
cohesiveZone.internalField()*
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
// the traction fraction is monitored to decide which faces to break:
// ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
const surfaceScalarField sigmaMax = rheology.cohLaw().sigmaMax();
const surfaceScalarField tauMax = rheology.cohLaw().tauMax();
@ -91,6 +91,7 @@ nCoupledFacesToBreak = 0;
faceToBreakIndex = facesToBreak[faceI];
}
}
scalar gMaxEffTractionFraction =
returnReduce(faceToBreakEffTractionFraction, maxOp<scalar>());
@ -111,7 +112,6 @@ nCoupledFacesToBreak = 0;
}
// Check if maximum is present on more then one processors
label procID = Pstream::nProcs();
if (procHasFaceToBreak)
{
@ -132,7 +132,6 @@ nCoupledFacesToBreak = 0;
SLList<label> coupledFacesToBreakList;
SLList<scalar> coupledFacesToBreakEffTractionFractionList;
forAll(mesh.boundary(), patchI)
{
if (mesh.boundary()[patchI].coupled())
@ -164,13 +163,15 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchDUPtr)
{
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+ (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+ (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
}
label start = mesh.boundaryMesh()[patchI].start();
@ -227,7 +228,6 @@ nCoupledFacesToBreak = 0;
}
}
scalar gMaxCoupledEffTractionFraction =
returnReduce(coupledFaceToBreakEffTractionFraction, maxOp<scalar>());
@ -249,7 +249,6 @@ nCoupledFacesToBreak = 0;
}
// Check if maximum is present on more then one processors
label procID = Pstream::nProcs();
if (procHasCoupledFaceToBreak)
{
@ -283,6 +282,7 @@ nCoupledFacesToBreak = 0;
{
label patchID =
mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
label start = mesh.boundaryMesh()[patchID].start();
label localIndex = coupledFaceToBreakIndex - start;
@ -358,19 +358,27 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchDUPtr)
{
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) );
)
);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) );
)
);
}
faceToBreakTraction *= scaleFactor;
@ -391,7 +399,7 @@ nCoupledFacesToBreak = 0;
faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
normalTrac = max(normalTrac, scalar(0));
normalTrac = max(normalTrac, 0.0);
scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
scalar scaleFactor = 1;
if(cohesivePatchDUPtr)
@ -401,8 +409,7 @@ nCoupledFacesToBreak = 0;
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*
(normalTrac/faceToBreakSigmaMax)
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
)
);
@ -415,10 +422,8 @@ nCoupledFacesToBreak = 0;
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*
(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*
(shearTrac/faceToBreakSigmaMax)
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
)
);
}
@ -614,7 +619,6 @@ nCoupledFacesToBreak = 0;
{
traction.boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
traction.oldTime().boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
@ -625,8 +629,10 @@ nCoupledFacesToBreak = 0;
}
else
{
cohesivePatchDUFixedModePtr->traction()[i] = -faceToBreakTraction;
cohesivePatchDUFixedModePtr->initiationTraction()[i] = -faceToBreakTraction;
cohesivePatchDUFixedModePtr->traction()[i] =
-faceToBreakTraction;
cohesivePatchDUFixedModePtr->initiationTraction()[i] =
-faceToBreakTraction;
}
}
}

View file

@ -12,7 +12,8 @@
// so a quick fix is to not set a reference on regions
// with a processor boundary
//if (U.boundaryField()[patchI].fixesValue())
if (
if
(
U.boundaryField()[patchI].fixesValue()
||
mesh.boundaryMesh()[patchI].type()

View file

@ -123,6 +123,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
volScalarField GI
(
IOobject
@ -137,6 +138,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
volScalarField GII
(
IOobject
@ -151,6 +153,7 @@ if (runTime.outputTime() || topoChange)
dimensionedScalar("zero", dimless, 0.0),
calculatedFvPatchVectorField::typeName
);
forAll(DU.boundaryField(), patchi)
{
// if(DU.boundaryField()[patchi].type() == cohesiveLawMultiMatFvPatchVectorField::typeName)

View file

@ -8,7 +8,8 @@ if(historyPatchID != -1)
vector direction(0, 1, 0);
//- for small strain or moving mesh
scalar force = gSum(
scalar force = gSum
(
direction &
(mesh.boundary()[historyPatchID].Sf() & sigma.boundaryField()[historyPatchID])
);

View file

@ -17,8 +17,7 @@ else if(divDSigmaExpMethod == "surface")
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU =
((I - n*n) & fvc::interpolate(gradDU));
surfaceTensorField shearGradDU = ((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp = fvc::div
(
@ -32,11 +31,13 @@ else if(divDSigmaExpMethod == "decompose")
}
else if(divDSigmaExpMethod == "laplacian")
{
divDSigmaExp = - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod
<< " not found!" << endl;
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -63,7 +63,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"

View file

@ -7,7 +7,8 @@ if(leftPatchID == -1)
}
//- 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])
);
@ -15,9 +16,10 @@ scalar leftForce = gSum(
//- patchIntegrate utility integrates it this way but this is worng because the sigma tensor should
//- be dotted with the surface normal to give the actual traction/force
//- you cannot just take the component of the sigma tensor
//scalar leftForcePatchIntegrateMethod = gSum(
// mesh.magSf().boundaryField()[leftPatchID]
// *sigma.boundaryField()[leftPatchID].component(symmTensor::XY)
//scalar leftForcePatchIntegrateMethod = gSum
//(
// mesh.magSf().boundaryField()[leftPatchID]*
// sigma.boundaryField()[leftPatchID].component(symmTensor::XY)
//);
vector gaugeU1 = vector::zero;

View file

@ -54,7 +54,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -73,14 +73,18 @@ int main(int argc, char *argv[])
fvm::d2dt2(rho, DU)
==
fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div(
+ fvc::div
(
- ( (mu + lambda) * gradDU )
+ ( mu * (
+ (
mu *
(
gradDU.T()
+ (gradDU & gradU.T())
+ (gradU & gradDU.T())
+ (gradDU & gradDU.T())
) )
)
)
+ ( lambda * tr(DEpsilon) * I )
+ ( DSigma & gradU )
+ ( (sigma + DSigma) & gradDU ),
@ -114,8 +118,7 @@ int main(int argc, char *argv[])
(
solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&&
++iCorr < nCorr
&& ++iCorr < nCorr
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()

View file

@ -56,7 +56,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -137,8 +137,7 @@ int main(int argc, char *argv[])
(
solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&&
++iCorr < nCorr
&& ++iCorr < nCorr
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << U.name()
@ -152,9 +151,9 @@ int main(int argc, char *argv[])
# include "writeFields.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl;
}
Info<< "End\n" << endl;

View file

@ -31,8 +31,7 @@ if(min(J.internalField()) > 0)
pointInterpolation.interpolate(U, pointU);
const vectorField& pointUI =
pointU.internalField();
const vectorField& pointUI = pointU.internalField();
//- Move mesh
vectorField newPoints = mesh.allPoints();

View file

@ -17,13 +17,12 @@ if(divDSigmaExpMethod == "standard")
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU =
((I - n*n)&fvc::interpolate(gradDU));
surfaceTensorField shearGradDU = ((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(DU) & (I - n*n))
+ lambdaf*tr(shearGradDU & (I - n*n))*n
+ muf*(shearGradDU & n)
@ -34,14 +33,11 @@ if(divDSigmaExpMethod == "standard")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -3,8 +3,7 @@
//----------------------------------------------------//
if(divDSigmaLargeStrainExpMethod == "standard")
{
divDSigmaLargeStrainExp =
fvc::div
divDSigmaLargeStrainExp = fvc::div
(
mu*(gradDU & gradDU.T())
//+ 0.5*lambda*(gradDU && gradDU)*I //- equivalent to 0.5*lambda*(I*tr(gradDU & gradDU.T()))
@ -15,8 +14,7 @@ if(divDSigmaLargeStrainExpMethod == "standard")
}
else if(divDSigmaLargeStrainExpMethod == "surface")
{
divDSigmaLargeStrainExp =
fvc::div
divDSigmaLargeStrainExp = fvc::div
(
muf * (mesh.Sf() & fvc::interpolate(gradDU & gradDU.T()))
+ 0.5*lambdaf * (mesh.Sf() & (fvc::interpolate(gradDU && gradDU)*I))

View file

@ -52,8 +52,7 @@ FieldField<Field, vector> extraVecs(ptc.size());
if
(
!isA<emptyFvPatch>(bm[patchID])
&& !bm[patchID].coupled()
!isA<emptyFvPatch>(bm[patchID]) && !bm[patchID].coupled()
)
{
// Found a face for extrapolation
@ -69,5 +68,4 @@ FieldField<Field, vector> extraVecs(ptc.size());
curExtraVectors.setSize(nFacesAroundPoint);
}
}

View file

@ -19,4 +19,3 @@ if(solidInterfaceCorr)
<< exit(FatalError);
}
}

View file

@ -66,7 +66,7 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
@ -148,4 +148,5 @@ int main(int argc, char *argv[])
return(0);
}
// ************************************************************************* //

View file

@ -7,7 +7,8 @@ pointVectorField& pf = pointDU;
// Do the correction
//GeometricField<Type, pointPatchField, pointMesh> pfCorr
/*pointVectorField pfCorr
/*
pointVectorField pfCorr
(
IOobject
(
@ -23,7 +24,8 @@ pointVectorField& pf = pointDU;
//dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
dimensionedVector("zero", pf.dimensions(), vector::zero),
pf.boundaryField().types()
);*/
);
*/
pointVectorField pfCorr
(
@ -96,21 +98,25 @@ forAll (ptc, pointI)
}
// Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI)
/*
forAll (pfCorr.boundaryField(), patchI)
{
if (pfCorr.boundaryField()[patchI].coupled())
{
pfCorr.boundaryField()[patchI].initAddField();
}
}*/
}
*/
/*forAll (pfCorr.boundaryField(), patchI)
/*
forAll (pfCorr.boundaryField(), patchI)
{
if (pfCorr.boundaryField()[patchI].coupled())
{
pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
}
}*/
}
*/
//Info << "pfCorr: " << pfCorr << endl;

View file

@ -22,9 +22,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest U

View file

@ -6,8 +6,9 @@ if(divSigmaExpMethod == "standard")
else if(divSigmaExpMethod == "surface")
{
//- this form seems to have the best convergence
divSigmaExp =
fvc::div(mesh.magSf()*
divSigmaExp = fvc::div
(
mesh.magSf()*
(
(n & (Cf && fvc::interpolate(symm(gradU))))
- (n & (Kf & fvc::interpolate(gradU)))
@ -21,5 +22,6 @@ if(divSigmaExpMethod == "standard")
}
else
{
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
}

View file

@ -5,8 +5,8 @@
vector netForce = vector::zero;
forAll(mesh.boundary(), patchi)
{
netForce +=
sum(
netForce += sum
(
mesh.Sf().boundaryField()[patchi]
&
(

View file

@ -123,6 +123,7 @@
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
forAll(cohesiveZone.internalField(), facei)
{
if(cohesiveZone.internalField()[facei])
@ -131,6 +132,7 @@
cohesiveZoneVol.internalField()[mesh.neighbour()[facei]] = 1.0;
}
}
forAll(cohesiveZone.boundaryField(), patchi)
{
forAll(cohesiveZone.boundaryField()[patchi], facei)

View file

@ -124,6 +124,7 @@
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
// aitken relaxation factor
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

View file

@ -90,7 +90,7 @@ int main(int argc, char *argv[])
runTime++;
Info<< "\nTime: " << runTime.timeName() << " s\n" << endl;
Info<< "\nTime = " << runTime.timeName() << " s\n" << endl;
volScalarField rho = rheology.rho();
volDiagTensorField K = rheology.K();
@ -173,7 +173,6 @@ int main(int argc, char *argv[])
// use leastSquaresSolidInterface grad scheme
gradU = fvc::grad(U);
# include "calculateRelativeResidual.H"
if (iCorr % infoFrequency == 0)
@ -199,8 +198,7 @@ int main(int argc, char *argv[])
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&&
iCorr < nCorr
&& iCorr < nCorr
)
);
@ -236,7 +234,8 @@ int main(int argc, char *argv[])
returnReduce
(
cohesivePatchUFixedModePtr->size(),
sumOp<label>())
sumOp<label>()
)
)
{
Pout << "Number of faces in crack: "

View file

@ -41,6 +41,7 @@ nCoupledFacesToBreak = 0;
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+ (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
}
maxEffTractionFraction = gMax(effTractionFraction);
SLList<label> facesToBreakList;
@ -85,7 +86,6 @@ nCoupledFacesToBreak = 0;
{
faceToBreakEffTractionFraction =
facesToBreakEffTractionFraction[faceI];
faceToBreakIndex = facesToBreak[faceI];
}
}
@ -145,6 +145,7 @@ nCoupledFacesToBreak = 0;
// only consider tensile tractions
pNormalTraction = max(pNormalTraction, scalar(0));
scalarField pShearTraction =
cohesiveZone.boundaryField()[patchI]*
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
@ -154,8 +155,7 @@ nCoupledFacesToBreak = 0;
const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
const scalarField& pTauMax = tauMax.boundaryField()[patchI];
scalarField pEffTractionFraction(pNormalTraction.size(), 0);
scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
if(cohesivePatchUPtr)
{
pEffTractionFraction =
@ -333,6 +333,7 @@ nCoupledFacesToBreak = 0;
vector faceToBreakNormal = vector::zero;
scalar faceToBreakSigmaMax = 0.0;
scalar faceToBreakTauMax = 0.0;
// Set faces to break
if (nFacesToBreak > 0)
{
@ -349,19 +350,27 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchUPtr)
{
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) );
)
);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) );
)
);
}
faceToBreakTraction *= scaleFactor;
@ -388,19 +397,27 @@ nCoupledFacesToBreak = 0;
if(cohesivePatchUPtr)
{
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) );
)
);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor =
::sqrt(1 / (
Foam::sqrt
(
1 /
(
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) );
)
);
}
faceToBreakTraction *= scaleFactor;
@ -571,7 +588,6 @@ nCoupledFacesToBreak = 0;
{
traction.boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
traction.oldTime().boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction;
@ -584,7 +600,6 @@ nCoupledFacesToBreak = 0;
{
cohesivePatchUFixedModePtr->traction()[i] =
-faceToBreakTraction;
cohesivePatchUFixedModePtr->initiationTraction()[i] =
-faceToBreakTraction;
}

View file

@ -12,7 +12,8 @@
// so a quick fix is to not set a reference on regions
// with a processor boundary
//if (U.boundaryField()[patchI].fixesValue())
if (
if
(
U.boundaryField()[patchI].fixesValue()
||
mesh.boundaryMesh()[patchI].type()

View file

@ -6,8 +6,10 @@ if(rheology.planeStress())
forAll(gradDU.internalField(), celli)
{
gradDU.internalField()[celli][tensor::ZZ] =
((-C.internalField()[celli][symmTensor4thOrder::XXZZ]*DEpsilon.internalField()[celli][symmTensor::XX]
- C.internalField()[celli][symmTensor4thOrder::YYZZ]*DEpsilon.internalField()[celli][symmTensor::YY])
(
(-C.internalField()[celli][symmTensor4thOrder::XXZZ]*DEpsilon.internalField()[celli][symmTensor::XX]
- C.internalField()[celli][symmTensor4thOrder::YYZZ]*DEpsilon.internalField()[celli][symmTensor::YY]
)
/
C.internalField()[celli][symmTensor4thOrder::ZZZZ])
- higherTerms.internalField()[celli];
@ -17,10 +19,16 @@ if(rheology.planeStress())
forAll(gradDU.boundaryField()[patchi], facei)
{
gradDU.boundaryField()[patchi][facei][tensor::ZZ] =
((-C.boundaryField()[patchi][facei][symmTensor4thOrder::XXZZ]*DEpsilon.boundaryField()[patchi][facei][symmTensor::XX]
- C.boundaryField()[patchi][facei][symmTensor4thOrder::YYZZ]*DEpsilon.boundaryField()[patchi][facei][symmTensor::YY])
(
(
- C.boundaryField()[patchi][facei][symmTensor4thOrder::XXZZ]*
DEpsilon.boundaryField()[patchi][facei][symmTensor::XX]
- C.boundaryField()[patchi][facei][symmTensor4thOrder::YYZZ]*
DEpsilon.boundaryField()[patchi][facei][symmTensor::YY]
)
/
C.boundaryField()[patchi][facei][symmTensor4thOrder::ZZZZ])
C.boundaryField()[patchi][facei][symmTensor4thOrder::ZZZZ]
)
- higherTerms.boundaryField()[patchi][facei];
}
}

View file

@ -29,6 +29,7 @@
dimensionedVector("zero", dimLength, vector::zero)
);
Info << "Reading accumulated strain field epsilon\n" << endl;
volSymmTensorField epsilon
(
IOobject
@ -57,6 +58,7 @@
dimensionedSymmTensor("zero", dimless, symmTensor::zero)
);
Info << "Reading accumulated stress field sigma\n" << endl;
volSymmTensorField sigma
(
IOobject
@ -72,6 +74,7 @@
);
Info << "Reading incremental stress field DSigma\n" << endl;
volSymmTensorField DSigma
(
IOobject

View file

@ -71,9 +71,9 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -105,7 +105,8 @@ int main(int argc, char *argv[])
+ fvc::d2dt2(rho, U)
==
fvm::laplacian(K, DU, "laplacian(K,DU)")
+ fvc::div(
+ fvc::div
(
DSigma
- (K & gradDU)
+ ( (sigma + DSigma) & gradDU ),
@ -155,8 +156,7 @@ int main(int argc, char *argv[])
while
(
solverPerf.initialResidual() > convergenceTolerance
&&
++iCorr < nCorr
&& ++iCorr < nCorr
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
@ -171,9 +171,9 @@ int main(int argc, char *argv[])
# include "rotateFields.H"
# include "writeFields.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s\n\n"
<< endl;
}
Info<< "End\n" << endl;

View file

@ -31,8 +31,7 @@
pointInterpolation.interpolate(DU, pointDU);
const vectorField& pointDUI =
pointDU.internalField();
const vectorField& pointDUI = pointDU.internalField();
//- Move mesh
vectorField newPoints = mesh.allPoints();

View file

@ -15,8 +15,7 @@ if (runTime.outputTime())
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField sigmaEq
(
@ -31,8 +30,7 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
// volVectorField traction
// (

View file

@ -1,18 +1,14 @@
if(divSigmaExpMethod == "standard")
{
//- calculating the full gradient has good convergence and no high freq oscillations
divSigmaExp =
fvc::div(
(C && symm(gradU))
- (K & gradU),
"div(sigma)"
);
divSigmaExp = fvc::div((C && symm(gradU)) - (K & gradU), "div(sigma)");
}
else if(divSigmaExpMethod == "surface")
{
//- this form seems to have the best convergence
divSigmaExp =
fvc::div(mesh.magSf()*
divSigmaExp = fvc::div
(
mesh.magSf()*
(
(n & (Cf && fvc::interpolate(symm(gradU))))
- (n & (Kf & fvc::interpolate(gradU)))
@ -28,5 +24,6 @@ if(divSigmaExpMethod == "standard")
}
else
{
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
}

View file

@ -60,9 +60,9 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -121,8 +121,7 @@ int main(int argc, char *argv[])
while
(
solverPerf.initialResidual() > convergenceTolerance
&&
++iCorr < nCorr
&& ++iCorr < nCorr
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << U.name()
@ -137,9 +136,9 @@ int main(int argc, char *argv[])
# include "writeFields.H"
# include "writeHistory.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s\n\n"
<< endl;
}
Info<< "End\n" << endl;

View file

@ -4,8 +4,10 @@
forAll(gradU.internalField(), celli)
{
gradU.internalField()[celli].zz() =
(-C.internalField()[celli].xxzz()*epsilon.internalField()[celli].xx()
- C.internalField()[celli].yyzz()*epsilon.internalField()[celli].yy())
(
- C.internalField()[celli].xxzz()*epsilon.internalField()[celli].xx()
- C.internalField()[celli].yyzz()*epsilon.internalField()[celli].yy()
)
/
C.internalField()[celli].zzzz();
}

View file

@ -15,8 +15,7 @@ if (runTime.outputTime())
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField sigmaEq
(
@ -31,8 +30,7 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
runTime.write();
}

View file

@ -27,9 +27,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest DU

View file

@ -2,9 +2,7 @@ if(divDSigmaExpMethod == "standard")
{
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)"
);
}
@ -43,8 +41,8 @@ else if(divDSigmaExpMethod == "decompose")
// divDSigmaExp = fvc::div
// (
// mesh.magSf()
// *(
// mesh.magSf()*
// (
// - (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
// + lambdaf*tr(shearGradDU&(I - n*n))*n
// + muf*(shearGradDU&n)
@ -55,12 +53,7 @@ else if(divDSigmaExpMethod == "laplacian")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{

View file

@ -197,8 +197,8 @@ forAll(extrapGradDU.internalField(), facei)
scalar gradGradDUZZYNei = gradGradDUZZY.internalField()[nei];
scalar gradGradDUZZZNei = gradGradDUZZZ.internalField()[nei];
tensor deltaOwnDotgradGradDUOwn =
tensor(
tensor deltaOwnDotgradGradDUOwn = tensor
(
deltaOwn.x()*gradGradDUXXXOwn + deltaOwn.y()*gradGradDUYXXOwn + deltaOwn.z()*gradGradDUZXXOwn,
deltaOwn.x()*gradGradDUXXYOwn + deltaOwn.y()*gradGradDUYXYOwn + deltaOwn.z()*gradGradDUZXYOwn,
deltaOwn.x()*gradGradDUXXZOwn + deltaOwn.y()*gradGradDUYXZOwn + deltaOwn.z()*gradGradDUZXZOwn,
@ -212,8 +212,8 @@ forAll(extrapGradDU.internalField(), facei)
deltaOwn.x()*gradGradDUXZZOwn + deltaOwn.y()*gradGradDUYZZOwn + deltaOwn.z()*gradGradDUZZZOwn
);
tensor deltaNeiDotgradGradDUNei =
tensor(
tensor deltaNeiDotgradGradDUNei = tensor
(
deltaNei.x()*gradGradDUXXXNei + deltaNei.y()*gradGradDUYXXNei + deltaNei.z()*gradGradDUZXXNei,
deltaNei.x()*gradGradDUXXYNei + deltaNei.y()*gradGradDUYXYNei + deltaNei.z()*gradGradDUZXYNei,
deltaNei.x()*gradGradDUXXZNei + deltaNei.y()*gradGradDUYXZNei + deltaNei.z()*gradGradDUZXZNei,
@ -249,11 +249,12 @@ forAll(extrapGradDU.boundaryField(), patchi)
}
// calculate thirdOrderTerm
volVectorField divThirdOrderTerm (
volVectorField divThirdOrderTerm
(
"thirdOrderTerm",
fvc::div(
(2*muf+lambdaf)*mesh.Sf()
& (extrapGradDU - averageGradDU)
fvc::div
(
(2*muf+lambdaf)*mesh.Sf() & (extrapGradDU - averageGradDU)
)
);

View file

@ -14,6 +14,7 @@
volTensorField gradDU = fvc::grad(DU);
Info<< "Creating field U\n" << endl;
volVectorField U
(
IOobject
@ -148,7 +149,6 @@
constitutiveModel rheology(sigma, DU);
volScalarField rho = rheology.rho();
volScalarField mu = rheology.mu();
volScalarField lambda = rheology.lambda();
surfaceScalarField muf = fvc::interpolate(mu, "mu");
@ -173,6 +173,7 @@
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
// aitken relaxation factor
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

View file

@ -58,7 +58,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -198,9 +198,9 @@ int main(int argc, char *argv[])
# include "writeFields.H"
# include "writeHistory.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl;
}
Info<< "End\n" << endl;

View file

@ -19,7 +19,7 @@ if(historyPatchID != -1)
Info << "Writing strain-stress to file for patch " << historyPatchName
<< endl;
// avaerage stress strain
// average stress strain
symmTensor stress = gAverage(sigma.boundaryField()[historyPatchID]);
symmTensor strain = gAverage(epsilon.boundaryField()[historyPatchID]);

View file

@ -27,9 +27,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest DU

View file

@ -17,13 +17,12 @@ if(divDSigmaExpMethod == "standard")
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU =
((I - n*n)&fvc::interpolate(gradDU));
surfaceTensorField shearGradDU = ((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(DU) & (I - n*n))
+ lambdaf*tr(shearGradDU&(I - n*n))*n
+ muf*(shearGradDU&n)
@ -34,14 +33,11 @@ if(divDSigmaExpMethod == "standard")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -3,8 +3,7 @@
//----------------------------------------------------//
if(divDSigmaLargeStrainExpMethod == "standard")
{
divDSigmaLargeStrainExp =
fvc::div
divDSigmaLargeStrainExp = fvc::div
(
mu*(gradDU & gradDU.T())
+ 0.5*lambda*(gradDU && gradDU)*I //- equivalent to 0.5*lambda*(I*tr(gradDU & gradDU.T()))
@ -14,8 +13,7 @@ if(divDSigmaLargeStrainExpMethod == "standard")
}
else if(divDSigmaLargeStrainExpMethod == "surface")
{
divDSigmaLargeStrainExp =
fvc::div
divDSigmaLargeStrainExp = fvc::div
(
muf * (mesh.Sf() & fvc::interpolate(gradDU & gradDU.T()))
+ 0.5*lambdaf * (mesh.Sf() & (fvc::interpolate(gradDU && gradDU)*I))

View file

@ -10,10 +10,10 @@ if(divDSigmaNonLinExpMethod == "standard")
}
else if(divDSigmaNonLinExpMethod == "surface")
{
divDSigmaNonLinExp =
fvc::div(
mesh.magSf()
*(
divDSigmaNonLinExp = fvc::div
(
mesh.magSf()*
(
( muf * (n & fvc::interpolate( gradDU & gradDU.T() )) )
+ ( 0.5*lambdaf * (n * tr(fvc::interpolate( gradDU & gradDU.T() ))) )
+ (n & fvc::interpolate( (sigma + DSigma) & gradDU ))

View file

@ -52,8 +52,7 @@ FieldField<Field, vector> extraVecs(ptc.size());
if
(
!isA<emptyFvPatch>(bm[patchID])
&& !bm[patchID].coupled()
!isA<emptyFvPatch>(bm[patchID]) && !bm[patchID].coupled()
)
{
// Found a face for extrapolation
@ -69,5 +68,4 @@ FieldField<Field, vector> extraVecs(ptc.size());
curExtraVectors.setSize(nFacesAroundPoint);
}
}

View file

@ -136,7 +136,6 @@
constitutiveModel rheology(sigma, DU);
volScalarField rho = rheology.rho();
volScalarField mu = rheology.mu();
volScalarField lambda = rheology.lambda();
surfaceScalarField muf = fvc::interpolate(rheology.mu());
@ -161,6 +160,7 @@
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
// aitken relaxation factor
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

View file

@ -73,7 +73,7 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
@ -170,10 +170,11 @@ int main(int argc, char *argv[])
(
iCorr++ < 2
||
(//solverPerf.initialResidual() > convergenceTolerance
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&&
iCorr < nCorr)
&& iCorr < nCorr
)
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
@ -204,4 +205,5 @@ int main(int argc, char *argv[])
return(0);
}
// ************************************************************************* //

View file

@ -8,7 +8,8 @@ if(moveMeshMethod == "inverseDistance")
}
else
{
FatalError << "move mesh method " << moveMeshMethod << " not recognised" << nl
FatalErrorIn(args.executable())
<< "move mesh method " << moveMeshMethod << " not recognised" << nl
<< "available methods are:" << nl
<< "inverseDistance" << nl
<< "leastSquares" << exit(FatalError);

View file

@ -7,7 +7,8 @@ pointVectorField& pf = pointDU;
// Do the correction
//GeometricField<Type, pointPatchField, pointMesh> pfCorr
/*pointVectorField pfCorr
/*
pointVectorField pfCorr
(
IOobject
(
@ -23,7 +24,8 @@ pointVectorField& pf = pointDU;
//dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
dimensionedVector("zero", pf.dimensions(), vector::zero),
pf.boundaryField().types()
);*/
);
*/
pointVectorField pfCorr
(
@ -96,21 +98,25 @@ forAll (ptc, pointI)
}
// Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI)
/*
forAll (pfCorr.boundaryField(), patchI)
{
if (pfCorr.boundaryField()[patchI].coupled())
{
pfCorr.boundaryField()[patchI].initAddField();
}
}*/
}
*/
/*forAll (pfCorr.boundaryField(), patchI)
/*
forAll (pfCorr.boundaryField(), patchI)
{
if (pfCorr.boundaryField()[patchI].coupled())
{
pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
}
}*/
}
*/
//Info << "pfCorr: " << pfCorr << endl;

View file

@ -13,8 +13,7 @@ if (runTime.outputTime())
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField sigmaEq
(
@ -29,8 +28,7 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
pointMesh pMesh(mesh);
pointScalarField contactPointGap
@ -69,7 +67,8 @@ if (runTime.outputTime())
//- total force
/* forAll(mesh.boundary(), patchi)
/*
forAll(mesh.boundary(), patchi)
{
Info << "Patch " << mesh.boundary()[patchi].name() << endl;
vectorField totalForce = mesh.Sf().boundaryField()[patchi] & sigma.boundaryField()[patchi];
@ -108,7 +107,8 @@ if (runTime.outputTime())
// << "SfTLy*sigmayy " << (SfTL[vector::Y]*sigma0[symmTensor::YY]) << nl
// << endl;
// }
}*/
}
*/
runTime.write();
}

View file

@ -27,9 +27,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest DU

View file

@ -10,15 +10,16 @@ if(divDSigmaExpMethod == "standard")
{
surfaceTensorField gradDUf = fvc::interpolate(gradDU);
divDSigmaExp =
fvc::div(
mesh.magSf()
*(
divDSigmaExp = fvc::div
(
mesh.magSf()*
(
muf*(n & gradDUf.T())
+ lambdaf*tr(gradDUf)*n
- (muf + lambdaf)*(n & gradDUf)
)
);
// divDSigmaExp = fvc::div
// (
// muf*(mesh.Sf() & fvc::interpolate(gradDU.T()))
@ -28,13 +29,12 @@ if(divDSigmaExpMethod == "standard")
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU =
((I - n*n)&fvc::interpolate(gradDU));
surfaceTensorField shearGradDU = ((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(DU) & (I - n*n))
+ lambdaf*tr(shearGradDU & (I - n*n))*n
+ muf*(shearGradDU & n)
@ -45,14 +45,11 @@ if(divDSigmaExpMethod == "standard")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
+ fvc::div(mu*gradDU.T() + lambda*(I*tr(gradDU)), "div(sigma)");
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -1,4 +1,4 @@
Info<< "Reading field DU\n" << endl;
Info<< "Reading displacement increment field DU\n" << endl;
volVectorField DU
(
IOobject
@ -141,5 +141,6 @@
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

View file

@ -58,9 +58,9 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -132,10 +132,11 @@ int main(int argc, char *argv[])
(
iCorr++ < 2
||
(//solverPerf.initialResidual() > convergenceTolerance
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&&
iCorr < nCorr)
&& iCorr < nCorr
)
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()

View file

@ -12,8 +12,8 @@ if (runTime.outputTime())
),
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField epsilonPEq
(
@ -27,8 +27,8 @@ if (runTime.outputTime())
),
sqrt((2.0/3.0)*magSqr(dev(epsilonP)))
);
Info<< "Max epsilonPEq = " << max(epsilonPEq).value()
<< endl;
Info<< "Max epsilonPEq = " << max(epsilonPEq).value() << endl;
volScalarField sigmaEq
(
@ -42,8 +42,8 @@ if (runTime.outputTime())
),
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
volScalarField sigmaHyd
(
@ -61,8 +61,7 @@ if (runTime.outputTime())
+ sigma.component(symmTensor::ZZ)
)/3
);
Info<< "Max sigmaHyd = " << max(sigmaHyd).value()
<< endl;
Info<< "Max sigmaHyd = " << max(sigmaHyd).value() << endl;
runTime.write();
}

View file

@ -1,5 +1,4 @@
// aitken acceleration
aitkenDelta.storePrevIter();
// update delta
@ -12,8 +11,7 @@ if(iCorr == 0)
}
else
{
vectorField b = aitkenDelta.internalField()
- aitkenDelta.prevIter().internalField();
vectorField b = aitkenDelta.internalField() - aitkenDelta.prevIter().internalField();
scalar sumMagB = gSum(magSqr(b));
if(sumMagB < SMALL)
@ -24,8 +22,7 @@ else
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest U

View file

@ -19,13 +19,12 @@ if(divSigmaExpMethod == "standard")
{
snGradU = fvc::snGrad(U);
surfaceTensorField shearGradU =
((I - n*n)&fvc::interpolate(gradU));
surfaceTensorField shearGradU = ((I - n*n) & fvc::interpolate(gradU));
divSigmaExp = fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(snGradU & (I - n*n))
+ lambdaf*tr(shearGradU & (I - n*n))*n
+ muf*(shearGradU & n)
@ -36,14 +35,10 @@ if(divSigmaExpMethod == "standard")
{
divSigmaExp =
- fvc::laplacian(mu + lambda, U, "laplacian(DU,U)")
+ fvc::div
(
mu*gradU.T()
+ lambda*(I*tr(gradU)),
"div(sigma)"
);
+ fvc::div(mu*gradU.T() + lambda*(I*tr(gradU)), "div(sigma)");
}
else
{
FatalError << "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
FatalErrorIn(args.executable())
<< "divSigmaExp method " << divSigmaExpMethod << " not found!" << endl;
}

View file

@ -197,8 +197,8 @@ forAll(extrapGradU.internalField(), facei)
scalar gradGradUZZYNei = gradGradUZZY.internalField()[nei];
scalar gradGradUZZZNei = gradGradUZZZ.internalField()[nei];
tensor deltaOwnDotgradGradUOwn =
tensor(
tensor deltaOwnDotgradGradUOwn =tensor
(
deltaOwn.x()*gradGradUXXXOwn + deltaOwn.y()*gradGradUYXXOwn + deltaOwn.z()*gradGradUZXXOwn,
deltaOwn.x()*gradGradUXXYOwn + deltaOwn.y()*gradGradUYXYOwn + deltaOwn.z()*gradGradUZXYOwn,
deltaOwn.x()*gradGradUXXZOwn + deltaOwn.y()*gradGradUYXZOwn + deltaOwn.z()*gradGradUZXZOwn,
@ -212,8 +212,8 @@ forAll(extrapGradU.internalField(), facei)
deltaOwn.x()*gradGradUXZZOwn + deltaOwn.y()*gradGradUYZZOwn + deltaOwn.z()*gradGradUZZZOwn
);
tensor deltaNeiDotgradGradUNei =
tensor(
tensor deltaNeiDotgradGradUNei = tensor
(
deltaNei.x()*gradGradUXXXNei + deltaNei.y()*gradGradUYXXNei + deltaNei.z()*gradGradUZXXNei,
deltaNei.x()*gradGradUXXYNei + deltaNei.y()*gradGradUYXYNei + deltaNei.z()*gradGradUZXYNei,
deltaNei.x()*gradGradUXXZNei + deltaNei.y()*gradGradUYXZNei + deltaNei.z()*gradGradUZXZNei,
@ -249,11 +249,12 @@ forAll(extrapGradU.boundaryField(), patchi)
}
// calculate thirdOrderTerm
volVectorField divThirdOrderTerm (
volVectorField divThirdOrderTerm
(
"thirdOrderTerm",
fvc::div(
(2*muf+lambdaf)*mesh.Sf()
& (extrapGradU - averageGradU)
fvc::div
(
(2*muf+lambdaf)*mesh.Sf() & (extrapGradU - averageGradU)
)
);

View file

@ -63,7 +63,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -142,7 +142,8 @@ int main(int argc, char *argv[])
while
(
iCorr++ == 0
|| (
||
(
solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&& iCorr < nCorr
@ -169,9 +170,9 @@ int main(int argc, char *argv[])
# include "writeFields.H"
# include "writeHistory.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl << endl;
}
Info<< "End\n" << endl;

View file

@ -26,9 +26,7 @@ if(iCorr == 0)
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
gSum(aitkenDelta.prevIter().internalField() & b)/sumMagB;
}
// correction to the latest U

View file

@ -21,11 +21,9 @@ else if(divSigmaExpMethod == "decompose")
{
snGradU = fvc::snGrad(U);
surfaceTensorField shearGradU =
((I - n*n)&fvc::interpolate(gradU));
surfaceTensorField shearGradU = ((I - n*n) & fvc::interpolate(gradU));
divSigmaExp = fvc::div
(
(
mesh.magSf()*
(
@ -33,21 +31,17 @@ else if(divSigmaExpMethod == "decompose")
+ lambdaf*tr(shearGradU & (I - n*n))*n
+ muf*(shearGradU & n)
)
)
- threeKalphaDeltaTf
);
}
/* else if(divSigmaExpMethod == "expLaplacian")
/*
else if(divSigmaExpMethod == "expLaplacian")
{
divSigmaExp =
- fvc::laplacian(mu + lambda, U, "laplacian(DU,U)")
+ fvc::div
(
mu*gradU.T()
+ lambda*(I*tr(gradU)),
"div(sigma)"
);
}*/
+ fvc::div(mu*gradU.T() + lambda*(I*tr(gradU)), "div(sigma)");
}
*/
else
{
FatalErrorIn(args.executable())

View file

@ -58,7 +58,7 @@ int main(int argc, char *argv[])
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"

View file

@ -13,8 +13,7 @@ if (runTime.outputTime())
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField sigmaEq
(
@ -29,8 +28,7 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
runTime.write();
}

View file

@ -84,7 +84,7 @@ int main(int argc, char *argv[])
Info << "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while (runTime.loop())
{
Info << "Time = " << runTime.timeName() << nl << endl;
@ -138,4 +138,5 @@ int main(int argc, char *argv[])
return(0);
}
// ************************************************************************* //

View file

@ -60,7 +60,7 @@ int main(int argc, char *argv[])
Info << "\nCalculating displacement field\n" << endl;
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Iteration: " << runTime.timeName() << nl << endl;
@ -85,4 +85,5 @@ int main(int argc, char *argv[])
return(0);
}
// ************************************************************************* //

View file

@ -59,9 +59,9 @@ int main(int argc, char *argv[])
scalar m = 0.5;
surfaceVectorField n = mesh.Sf()/mesh.magSf();
for (runTime++; !runTime.end(); runTime++)
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -96,8 +96,8 @@ int main(int argc, char *argv[])
fvm::laplacian(2*muf+lambdaf, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mesh.magSf()
*(
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
+ lambdaf*tr(sGradDU&(I - n*n))*n
+ muf*(sGradDU&n)
@ -161,9 +161,9 @@ int main(int argc, char *argv[])
# include "writeFields.H"
# include "writeHistory.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;