exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Feature: Single precision and long double precision port

Browse source
This commit is contained in:
Hrvoje Jasak 2015-05-01 13:21:04 +01:00 committed by Dominik Christ
parent 51cd34cc0e
commit 6b022758d1
48 changed files with 1129 additions and 600 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C
View file

@ -46,7 +46,7 @@ ersPointSource::ersPointSource
alpha_(readScalar(dict.lookup("alpha"))), alpha_(readScalar(dict.lookup("alpha"))),
direction_(dict.lookup("direction")) direction_(dict.lookup("direction"))
{ {
q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), 0.0); q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), scalar(0));
} }

12
applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C
View file

@ -30,18 +30,10 @@ Description
#include "contactPatchPair.H" #include "contactPatchPair.H"
#include "volFields.H" #include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components // Construct from components
contactPatchPair::contactPatchPair Foam::contactPatchPair::contactPatchPair
( (
const fvMesh& m, const fvMesh& m,
const label master, const label master,
@ -75,6 +67,4 @@ contactPatchPair::contactPatchPair
{} {}
} // End namespace Foam
// ************************************************************************* // // ************************************************************************* //

2
applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C
View file

@ -302,7 +302,7 @@ void Foam::contactPatchPair::correct
( (
slavePressure slavePressure
), ),
0.0 scalar(0)
); );
// Calculate master traction, using the positive part of // Calculate master traction, using the positive part of

50
applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H
View file

@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0;
// scalarField effTraction = // scalarField effTraction =
// cohesiveZone.internalField() * // cohesiveZone.internalField() *
// mag(traction.internalField()); // mag(traction.internalField());
scalarField normalTraction = scalarField normalTraction =
cohesiveZone.internalField() * cohesiveZone.internalField()*
( n.internalField() & traction.internalField() ); ( n.internalField() & traction.internalField() );
normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
scalarField shearTraction = // only consider tensile tractions
cohesiveZone.internalField() * normalTraction = max(normalTraction, scalar(0));
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); scalarField shearTraction =
cohesiveZone.internalField() *
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
// the traction fraction is monitored to decide which faces to break: // the traction fraction is monitored to decide which faces to break:
// ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@ -27,16 +29,18 @@ nCoupledFacesToBreak = 0;
scalarField effTractionFraction(normalTraction.size(), 0.0); scalarField effTractionFraction(normalTraction.size(), 0.0);
if(cohesivePatchUPtr) if(cohesivePatchUPtr)
{ {
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
}
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
}
maxEffTractionFraction = gMax(effTractionFraction); maxEffTractionFraction = gMax(effTractionFraction);
SLList<label> facesToBreakList; SLList<label> facesToBreakList;
@ -44,15 +48,21 @@ nCoupledFacesToBreak = 0;
forAll(effTractionFraction, faceI) forAll(effTractionFraction, faceI)
{ {
if (effTractionFraction[faceI] > 1.0) if (effTractionFraction[faceI] > 1.0)
{ {
facesToBreakList.insert(faceI); facesToBreakList.insert(faceI);
facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]); facesToBreakEffTractionFractionList.insert
(
effTractionFraction[faceI]
);
} }
} }
labelList facesToBreak(facesToBreakList); labelList facesToBreak(facesToBreakList);
List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList); List<scalar> facesToBreakEffTractionFraction
(
facesToBreakEffTractionFractionList
);
nFacesToBreak = facesToBreak.size(); nFacesToBreak = facesToBreak.size();
@ -122,7 +132,7 @@ nCoupledFacesToBreak = 0;
scalarField pNormalTraction = scalarField pNormalTraction =
cohesiveZone.boundaryField()[patchI] * cohesiveZone.boundaryField()[patchI] *
( n.boundaryField()[patchI] & traction.boundaryField()[patchI] ); ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions pNormalTraction = max(pNormalTraction, scalar(0)); // only consider tensile tractions
scalarField pShearTraction = scalarField pShearTraction =
cohesiveZone.boundaryField()[patchI] * cohesiveZone.boundaryField()[patchI] *
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] ); mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );

130
applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H
View file

@ -6,13 +6,16 @@ nCoupledFacesToBreak = 0;
// scalarField effTraction = // scalarField effTraction =
// cohesiveZone.internalField() * // cohesiveZone.internalField() *
// mag(traction.internalField()); // mag(traction.internalField());
scalarField normalTraction = scalarField normalTraction =
cohesiveZone.internalField() * cohesiveZone.internalField()*
( n.internalField() & traction.internalField() ); ( n.internalField() & traction.internalField() );
normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
scalarField shearTraction = // only consider tensile tractions
cohesiveZone.internalField() * normalTraction = max(normalTraction, scalar(0));
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
scalarField shearTraction =
cohesiveZone.internalField() *
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
// the traction fraction is monitored to decide which faces to break: // the traction fraction is monitored to decide which faces to break:
// ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@ -28,16 +31,18 @@ nCoupledFacesToBreak = 0;
scalarField effTractionFraction(normalTraction.size(), 0.0); scalarField effTractionFraction(normalTraction.size(), 0.0);
if(cohesivePatchDUPtr) if(cohesivePatchDUPtr)
{ {
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
}
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
}
maxEffTractionFraction = gMax(effTractionFraction); maxEffTractionFraction = gMax(effTractionFraction);
@ -46,15 +51,21 @@ nCoupledFacesToBreak = 0;
forAll(effTractionFraction, faceI) forAll(effTractionFraction, faceI)
{ {
if (effTractionFraction[faceI] > 1.0) if (effTractionFraction[faceI] > 1.0)
{ {
facesToBreakList.insert(faceI); facesToBreakList.insert(faceI);
facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]); facesToBreakEffTractionFractionList.insert
(
effTractionFraction[faceI]
);
} }
} }
labelList facesToBreak(facesToBreakList); labelList facesToBreak(facesToBreakList);
List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList); List<scalar> facesToBreakEffTractionFraction
(
facesToBreakEffTractionFractionList
);
nFacesToBreak = facesToBreak.size(); nFacesToBreak = facesToBreak.size();
@ -69,9 +80,14 @@ nCoupledFacesToBreak = 0;
scalar faceToBreakEffTractionFraction = 0; scalar faceToBreakEffTractionFraction = 0;
forAll(facesToBreakEffTractionFraction, faceI) forAll(facesToBreakEffTractionFraction, faceI)
{ {
if (facesToBreakEffTractionFraction[faceI] > faceToBreakEffTractionFraction) if
(
facesToBreakEffTractionFraction[faceI]
> faceToBreakEffTractionFraction
)
{ {
faceToBreakEffTractionFraction = facesToBreakEffTractionFraction[faceI]; faceToBreakEffTractionFraction =
facesToBreakEffTractionFraction[faceI];
faceToBreakIndex = facesToBreak[faceI]; faceToBreakIndex = facesToBreak[faceI];
} }
} }
@ -83,7 +99,11 @@ nCoupledFacesToBreak = 0;
bool procHasFaceToBreak = false; bool procHasFaceToBreak = false;
if (nFacesToBreak > 0) if (nFacesToBreak > 0)
{ {
if ( mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction) < SMALL ) if
(
mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
< SMALL
)
{ {
// philipc - Maximum traction fraction is on this processor // philipc - Maximum traction fraction is on this processor
procHasFaceToBreak = true; procHasFaceToBreak = true;
@ -123,9 +143,12 @@ nCoupledFacesToBreak = 0;
// scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI]; // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
scalarField pNormalTraction = scalarField pNormalTraction =
cohesiveZone.boundaryField()[patchI] * cohesiveZone.boundaryField()[patchI]*
( n.boundaryField()[patchI] & traction.boundaryField()[patchI] ); ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions
// only consider tensile tractions
pNormalTraction = max(pNormalTraction, scalar(0));
scalarField pShearTraction = scalarField pShearTraction =
cohesiveZone.boundaryField()[patchI] * cohesiveZone.boundaryField()[patchI] *
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] ); mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
@ -139,24 +162,24 @@ nCoupledFacesToBreak = 0;
// (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax); // (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
scalarField pEffTractionFraction(pNormalTraction.size(), 0.0); scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
if(cohesivePatchDUPtr) if(cohesivePatchDUPtr)
{ {
pEffTractionFraction = pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax); (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
} }
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
pEffTractionFraction = pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax); (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
} }
label start = mesh.boundaryMesh()[patchI].start(); label start = mesh.boundaryMesh()[patchI].start();
forAll(pEffTractionFraction, faceI) forAll(pEffTractionFraction, faceI)
{ {
if (pEffTractionFraction[faceI] > maxEffTractionFraction) if (pEffTractionFraction[faceI] > maxEffTractionFraction)
{ {
maxEffTractionFraction = pEffTractionFraction[faceI]; maxEffTractionFraction = pEffTractionFraction[faceI];
} }
if (pEffTractionFraction[faceI] > 1.0) if (pEffTractionFraction[faceI] > 1.0)
@ -368,30 +391,41 @@ nCoupledFacesToBreak = 0;
faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex]; faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex]; faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
scalar normalTrac = faceToBreakNormal & faceToBreakTraction; scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
normalTrac = max(normalTrac, 0.0); normalTrac = max(normalTrac, scalar(0));
scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction ); scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
scalar scaleFactor = 1; scalar scaleFactor = 1;
if(cohesivePatchDUPtr) if(cohesivePatchDUPtr)
{ {
scaleFactor = scaleFactor =
::sqrt(1 / ( Foam::sqrt
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax) 1/
) ); (
} (normalTrac/faceToBreakSigmaMax)*
(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
)
);
}
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor = scaleFactor =
::sqrt(1 / ( Foam::sqrt
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax) 1/
) ); (
} (normalTrac/faceToBreakSigmaMax)*
(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*
(shearTrac/faceToBreakSigmaMax)
)
);
}
faceToBreakTraction *= scaleFactor; faceToBreakTraction *= scaleFactor;
topoChange = true; topoChange = true;
} }
reduce(topoChange, orOp<bool>()); reduce(topoChange, orOp<bool>());

409
applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/updateCrack.H
View file

@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0;
// scalarField effTraction = // scalarField effTraction =
// cohesiveZone.internalField() * // cohesiveZone.internalField() *
// mag(traction.internalField()); // mag(traction.internalField());
scalarField normalTraction = scalarField normalTraction =
cohesiveZone.internalField() * cohesiveZone.internalField()*
( n.internalField() & traction.internalField() ); ( n.internalField() & traction.internalField() );
normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
scalarField shearTraction = // only consider tensile tractions
cohesiveZone.internalField() * normalTraction = max(normalTraction, scalar(0));
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); scalarField shearTraction =
cohesiveZone.internalField()*
mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
// the traction fraction is monitored to decide which faces to break: // the traction fraction is monitored to decide which faces to break:
// ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@ -26,17 +28,19 @@ nCoupledFacesToBreak = 0;
//scalarField effTractionFraction = effTraction/sigmaMax; //scalarField effTractionFraction = effTraction/sigmaMax;
scalarField effTractionFraction(normalTraction.size(), 0.0); scalarField effTractionFraction(normalTraction.size(), 0.0);
if(cohesivePatchUPtr) if (cohesivePatchUPtr)
{ {
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
}
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
effTractionFraction = effTractionFraction =
(normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
} + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
}
maxEffTractionFraction = gMax(effTractionFraction); maxEffTractionFraction = gMax(effTractionFraction);
SLList<label> facesToBreakList; SLList<label> facesToBreakList;
@ -44,15 +48,21 @@ nCoupledFacesToBreak = 0;
forAll(effTractionFraction, faceI) forAll(effTractionFraction, faceI)
{ {
if (effTractionFraction[faceI] > 1.0) if (effTractionFraction[faceI] > 1.0)
{ {
facesToBreakList.insert(faceI); facesToBreakList.insert(faceI);
facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]); facesToBreakEffTractionFractionList.insert
(
effTractionFraction[faceI]
);
} }
} }
labelList facesToBreak(facesToBreakList); labelList facesToBreak(facesToBreakList);
List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList); List<scalar> facesToBreakEffTractionFraction
(
facesToBreakEffTractionFractionList
);
nFacesToBreak = facesToBreak.size(); nFacesToBreak = facesToBreak.size();
@ -67,9 +77,15 @@ nCoupledFacesToBreak = 0;
scalar faceToBreakEffTractionFraction = 0; scalar faceToBreakEffTractionFraction = 0;
forAll(facesToBreakEffTractionFraction, faceI) forAll(facesToBreakEffTractionFraction, faceI)
{ {
if (facesToBreakEffTractionFraction[faceI] > faceToBreakEffTractionFraction) if
(
facesToBreakEffTractionFraction[faceI]
> faceToBreakEffTractionFraction
)
{ {
faceToBreakEffTractionFraction = facesToBreakEffTractionFraction[faceI]; faceToBreakEffTractionFraction =
facesToBreakEffTractionFraction[faceI];
faceToBreakIndex = facesToBreak[faceI]; faceToBreakIndex = facesToBreak[faceI];
} }
} }
@ -82,7 +98,11 @@ nCoupledFacesToBreak = 0;
bool procHasFaceToBreak = false; bool procHasFaceToBreak = false;
if (nFacesToBreak > 0) if (nFacesToBreak > 0)
{ {
if ( mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction) < SMALL ) if
(
mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
< SMALL
)
{ {
// philipc - Maximum traction fraction is on this processor // philipc - Maximum traction fraction is on this processor
procHasFaceToBreak = true; procHasFaceToBreak = true;
@ -90,7 +110,7 @@ nCoupledFacesToBreak = 0;
} }
// Check if maximum is present on more then one processors // Check if maximum is present on more then one processors
label procID = Pstream::nProcs(); label procID = Pstream::nProcs();
if (procHasFaceToBreak) if (procHasFaceToBreak)
{ {
procID = Pstream::myProcNo(); procID = Pstream::myProcNo();
@ -115,54 +135,59 @@ nCoupledFacesToBreak = 0;
if (mesh.boundary()[patchI].coupled()) if (mesh.boundary()[patchI].coupled())
{ {
// scalarField pEffTraction = // scalarField pEffTraction =
// cohesiveZone.boundaryField()[patchI] * // cohesiveZone.boundaryField()[patchI] *
// mag(traction.boundaryField()[patchI]); // mag(traction.boundaryField()[patchI]);
// scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI]; // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
scalarField pNormalTraction = scalarField pNormalTraction =
cohesiveZone.boundaryField()[patchI] * cohesiveZone.boundaryField()[patchI] *
( n.boundaryField()[patchI] & traction.boundaryField()[patchI] ); ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions
scalarField pShearTraction =
cohesiveZone.boundaryField()[patchI] *
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
// the traction fraction is monitored to decide which faces to break: // only consider tensile tractions
// ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face pNormalTraction = max(pNormalTraction, scalar(0));
const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI]; scalarField pShearTraction =
const scalarField& pTauMax = tauMax.boundaryField()[patchI]; cohesiveZone.boundaryField()[patchI] *
mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
scalarField pEffTractionFraction(pNormalTraction.size(), 0.0); // the traction fraction is monitored to decide which faces to break:
if(cohesivePatchUPtr) // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
{ const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
pEffTractionFraction = const scalarField& pTauMax = tauMax.boundaryField()[patchI];
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
}
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
}
label start = mesh.boundaryMesh()[patchI].start(); scalarField pEffTractionFraction(pNormalTraction.size(), 0);
forAll(pEffTractionFraction, faceI) if(cohesivePatchUPtr)
{ {
if (pEffTractionFraction[faceI] > maxEffTractionFraction) pEffTractionFraction =
{ (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
maxEffTractionFraction = pEffTractionFraction[faceI]; + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
} }
else
{
// solidCohesiveFixedModeMix only uses sigmaMax
pEffTractionFraction =
(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+ (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
}
if (pEffTractionFraction[faceI] > 1.0) label start = mesh.boundaryMesh()[patchI].start();
{
coupledFacesToBreakList.insert(start + faceI); forAll(pEffTractionFraction, faceI)
coupledFacesToBreakEffTractionFractionList.insert {
( if (pEffTractionFraction[faceI] > maxEffTractionFraction)
pEffTractionFraction[faceI] {
); maxEffTractionFraction = pEffTractionFraction[faceI];
} }
}
if (pEffTractionFraction[faceI] > 1.0)
{
coupledFacesToBreakList.insert(start + faceI);
coupledFacesToBreakEffTractionFractionList.insert
(
pEffTractionFraction[faceI]
);
}
}
} }
} }
@ -250,7 +275,7 @@ nCoupledFacesToBreak = 0;
if (nCoupledFacesToBreak) if (nCoupledFacesToBreak)
{ {
label patchID = label patchID =
mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex); mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
label start = mesh.boundaryMesh()[patchID].start(); label start = mesh.boundaryMesh()[patchID].start();
label localIndex = coupledFaceToBreakIndex - start; label localIndex = coupledFaceToBreakIndex - start;
@ -315,31 +340,31 @@ nCoupledFacesToBreak = 0;
faceToBreakNormal = n.internalField()[faceToBreakIndex]; faceToBreakNormal = n.internalField()[faceToBreakIndex];
// Scale broken face traction // Scale broken face traction
faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex]; faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex];
faceToBreakTauMax = tauMaxI[faceToBreakIndex]; faceToBreakTauMax = tauMaxI[faceToBreakIndex];
scalar normalTrac = faceToBreakNormal & faceToBreakTraction; scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
normalTrac = max(normalTrac, 0.0); normalTrac = max(normalTrac, 0.0);
scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction ); scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
scalar scaleFactor = 1; scalar scaleFactor = 1;
if(cohesivePatchUPtr) if(cohesivePatchUPtr)
{ {
scaleFactor = scaleFactor =
::sqrt(1 / ( ::sqrt(1 / (
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax) + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) ); ) );
} }
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor = scaleFactor =
::sqrt(1 / ( ::sqrt(1 / (
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax) + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) ); ) );
} }
faceToBreakTraction *= scaleFactor; faceToBreakTraction *= scaleFactor;
topoChange = true; topoChange = true;
} }
@ -354,29 +379,29 @@ nCoupledFacesToBreak = 0;
faceToBreakNormal = n.boundaryField()[patchID][localIndex]; faceToBreakNormal = n.boundaryField()[patchID][localIndex];
// Scale broken face traction // Scale broken face traction
faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex]; faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex]; faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
scalar normalTrac = faceToBreakNormal & faceToBreakTraction; scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
normalTrac = max(normalTrac, 0.0); normalTrac = max(normalTrac, 0.0);
scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction ); scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
scalar scaleFactor = 1; scalar scaleFactor = 1;
if(cohesivePatchUPtr) if(cohesivePatchUPtr)
{ {
scaleFactor = scaleFactor =
::sqrt(1 / ( ::sqrt(1 / (
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax) + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
) ); ) );
} }
else else
{ {
// solidCohesiveFixedModeMix only uses sigmaMax // solidCohesiveFixedModeMix only uses sigmaMax
scaleFactor = scaleFactor =
::sqrt(1 / ( ::sqrt(1 / (
(normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax) (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+ (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax) + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
) ); ) );
} }
faceToBreakTraction *= scaleFactor; faceToBreakTraction *= scaleFactor;
@ -416,20 +441,20 @@ nCoupledFacesToBreak = 0;
Cf = fvc::interpolate(C); Cf = fvc::interpolate(C);
Kf = fvc::interpolate(K); Kf = fvc::interpolate(K);
// we need to modify propertiess after cracking otherwise momentum equation is wrong // we need to modify propertiess after cracking otherwise momentum equation is wrong
// but solidInterface seems to hold some information about old mesh // but solidInterface seems to hold some information about old mesh
// so we will delete it and make another // so we will delete it and make another
// we could probably add a public clearout function // we could probably add a public clearout function
// create new solidInterface // create new solidInterface
//Pout << "Creating new solidInterface" << endl; //Pout << "Creating new solidInterface" << endl;
//delete solidInterfacePtr; //delete solidInterfacePtr;
//solidInterfacePtr = new solidInterface(mesh, rheology); //solidInterfacePtr = new solidInterface(mesh, rheology);
// delete demand driven data as the mesh has changed // delete demand driven data as the mesh has changed
if(rheology.solidInterfaceActive()) if(rheology.solidInterfaceActive())
{ {
rheology.solInterface().clearOut(); rheology.solInterface().clearOut();
solidInterfacePtr->modifyProperties(Cf, Kf); solidInterfacePtr->modifyProperties(Cf, Kf);
} }
// Local crack displacement // Local crack displacement
vectorField UpI = vectorField UpI =
@ -441,21 +466,21 @@ nCoupledFacesToBreak = 0;
vectorField globalUpI = mesh.globalCrackField(UpI); vectorField globalUpI = mesh.globalCrackField(UpI);
vectorField globalOldUpI = mesh.globalCrackField(oldUpI); vectorField globalOldUpI = mesh.globalCrackField(oldUpI);
// C and K field on new crack faces must be updated // C and K field on new crack faces must be updated
symmTensor4thOrderField CPI = C.boundaryField()[cohesivePatchID].patchInternalField(); symmTensor4thOrderField CPI = C.boundaryField()[cohesivePatchID].patchInternalField();
diagTensorField KPI = K.boundaryField()[cohesivePatchID].patchInternalField(); diagTensorField KPI = K.boundaryField()[cohesivePatchID].patchInternalField();
symmTensor4thOrderField globalCPI = mesh.globalCrackField(CPI); symmTensor4thOrderField globalCPI = mesh.globalCrackField(CPI);
diagTensorField globalKPI = mesh.globalCrackField(KPI); diagTensorField globalKPI = mesh.globalCrackField(KPI);
// cohesivePatchU.size() // cohesivePatchU.size()
int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size()); int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
// Initialise U for new cohesive face // Initialise U for new cohesive face
const labelList& gcfa = mesh.globalCrackFaceAddressing(); const labelList& gcfa = mesh.globalCrackFaceAddressing();
label globalIndex = mesh.localCrackStart(); label globalIndex = mesh.localCrackStart();
// for (label i=0; i<cohesivePatchU.size(); i++) // for (label i=0; i<cohesivePatchU.size(); i++)
for (label i=0; i<cohesivePatchSize; i++) for (label i=0; i<cohesivePatchSize; i++)
{ {
label oldFaceIndex = faceMap[start+i]; label oldFaceIndex = faceMap[start+i];
// If new face // If new face
@ -474,10 +499,10 @@ nCoupledFacesToBreak = 0;
+ globalOldUpI[gcfa[globalIndex]] + globalOldUpI[gcfa[globalIndex]]
); );
// initialise C and K on new faces // initialise C and K on new faces
// set new face value to value of internal cell // set new face value to value of internal cell
Cf.boundaryField()[cohesivePatchID][i] = globalCPI[globalIndex]; Cf.boundaryField()[cohesivePatchID][i] = globalCPI[globalIndex];
Kf.boundaryField()[cohesivePatchID][i] = globalKPI[globalIndex]; Kf.boundaryField()[cohesivePatchID][i] = globalKPI[globalIndex];
globalIndex++; globalIndex++;
} }
@ -488,81 +513,87 @@ nCoupledFacesToBreak = 0;
} }
// we must calculate grad using interface // we must calculate grad using interface
// U at the interface has not been calculated yet as interface.correct() // U at the interface has not been calculated yet as interface.correct()
// has not been called yet // has not been called yet
// not really a problem as gradU is correct in second outer iteration // not really a problem as gradU is correct in second outer iteration
// as long as this does not cause convergence problems for the first iterations. // as long as this does not cause convergence problems for the first iterations.
// we should be able to calculate the interface displacements without // we should be able to calculate the interface displacements without
// having to call interface.correct() // having to call interface.correct()
// todo: add calculateInterfaceU() function // todo: add calculateInterfaceU() function
// interface grad uses Gauss, we need least squares // interface grad uses Gauss, we need least squares
//gradU = solidInterfacePtr->grad(U); //gradU = solidInterfacePtr->grad(U);
gradU = fvc::grad(U); // leastSquaresSolidInterface grad scheme gradU = fvc::grad(U); // leastSquaresSolidInterface grad scheme
//snGradU = fvc::snGrad(U); //snGradU = fvc::snGrad(U);
# include "calculateTraction.H" # include "calculateTraction.H"
//if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write(); //if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
// Initialise initiation traction for new cohesive patch face // Initialise initiation traction for new cohesive patch face
// for (label i=0; i<cohesivePatchU.size(); i++) // for (label i=0; i<cohesivePatchU.size(); i++)
for (label i=0; i<cohesivePatchSize; i++) for (label i=0; i<cohesivePatchSize; i++)
{
label oldFaceIndex = faceMap[start+i];
// If new face
if
(
(oldFaceIndex == faceToBreakIndex)
|| (oldFaceIndex == coupledFaceToBreakIndex)
)
{ {
label oldFaceIndex = faceMap[start+i]; vector n0 =
mesh.Sf().boundaryField()[cohesivePatchID][i]
/mesh.magSf().boundaryField()[cohesivePatchID][i];
//vector n1 = -n0;
// If new face if ((n0&faceToBreakNormal) > SMALL)
if
(
(oldFaceIndex == faceToBreakIndex)
|| (oldFaceIndex == coupledFaceToBreakIndex)
)
{ {
vector n0 = traction.boundaryField()[cohesivePatchID][i] =
mesh.Sf().boundaryField()[cohesivePatchID][i]
/mesh.magSf().boundaryField()[cohesivePatchID][i];
//vector n1 = -n0;
if ((n0&faceToBreakNormal) > SMALL)
{
traction.boundaryField()[cohesivePatchID][i] =
faceToBreakTraction; faceToBreakTraction;
traction.oldTime().boundaryField()[cohesivePatchID][i] = traction.oldTime().boundaryField()[cohesivePatchID][i] =
faceToBreakTraction; faceToBreakTraction;
if(cohesivePatchUPtr) if(cohesivePatchUPtr)
{ {
cohesivePatchUPtr->traction()[i] = faceToBreakTraction; cohesivePatchUPtr->traction()[i] = faceToBreakTraction;
}
else
{
cohesivePatchUFixedModePtr->traction()[i] = faceToBreakTraction;
cohesivePatchUFixedModePtr->initiationTraction()[i] = faceToBreakTraction;
}
} }
else else
{ {
traction.boundaryField()[cohesivePatchID][i] = cohesivePatchUFixedModePtr->traction()[i] =
faceToBreakTraction;
cohesivePatchUFixedModePtr->initiationTraction()[i] =
faceToBreakTraction;
}
}
else
{
traction.boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction; -faceToBreakTraction;
traction.oldTime().boundaryField()[cohesivePatchID][i] = traction.oldTime().boundaryField()[cohesivePatchID][i] =
-faceToBreakTraction; -faceToBreakTraction;
//cohesivePatchU.traction()[i] = -faceToBreakTraction; //cohesivePatchU.traction()[i] = -faceToBreakTraction;
if(cohesivePatchUPtr) if(cohesivePatchUPtr)
{ {
cohesivePatchUPtr->traction()[i] = -faceToBreakTraction; cohesivePatchUPtr->traction()[i] = -faceToBreakTraction;
} }
else else
{ {
cohesivePatchUFixedModePtr->traction()[i] = -faceToBreakTraction; cohesivePatchUFixedModePtr->traction()[i] =
cohesivePatchUFixedModePtr->initiationTraction()[i] = -faceToBreakTraction; -faceToBreakTraction;
}
cohesivePatchUFixedModePtr->initiationTraction()[i] =
-faceToBreakTraction;
} }
} }
} }
}
// hmmnn we only need a reference for very small groups of cells // hmmnn we only need a reference for very small groups of cells
// turn off for now // turn off for now
//# include "updateReference.H" //# include "updateReference.H"
} }
} }

2
applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C
View file

@ -43,7 +43,7 @@ int main(int argc, char *argv[])
Info << "tf value: " << FadOneValue(tf) << endl; Info << "tf value: " << FadOneValue(tf) << endl;
Info << "tf deriv: " << FadOneDeriv(tf, 0) << endl; Info << "tf deriv: " << FadOneDeriv(tf, 0) << endl;
FadOneSetDeriv(tf, 0, Field<double>(3, 1)); FadOneSetDeriv(tf, 0, scalarField(3, 1));
fadScalar ten(10.0); fadScalar ten(10.0);
fadScalar bs = tf[0] - ten; fadScalar bs = tf[0] - ten;

2
etc/bashrc
View file

@ -138,7 +138,7 @@ export WM_COMPILER_LIB_ARCH=
# WM_ARCH_OPTION = 32 | 64 # WM_ARCH_OPTION = 32 | 64
: ${WM_ARCH_OPTION:=64}; export WM_ARCH_OPTION : ${WM_ARCH_OPTION:=64}; export WM_ARCH_OPTION
# WM_PRECISION_OPTION = DP | SP # WM_PRECISION_OPTION = LDP | DP | SP
: ${WM_PRECISION_OPTION:=DP}; export WM_PRECISION_OPTION : ${WM_PRECISION_OPTION:=DP}; export WM_PRECISION_OPTION
# WM_COMPILE_OPTION = Opt | Debug | Prof # WM_COMPILE_OPTION = Opt | Debug | Prof

2
src/Pstream/mpi/Pstream.C
View file

@ -37,6 +37,8 @@ License
# define MPI_SCALAR MPI_FLOAT # define MPI_SCALAR MPI_FLOAT
#elif defined(WM_DP) #elif defined(WM_DP)
# define MPI_SCALAR MPI_DOUBLE # define MPI_SCALAR MPI_DOUBLE
#elif defined(WM_LDP)
# define MPI_SCALAR MPI_LONG_DOUBLE
#endif #endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

14
src/cudaSolvers/cudaSolver/cudaTypes.H
View file

@ -38,15 +38,19 @@ Author
#ifndef cudaTypes_H #ifndef cudaTypes_H
#define cudaTypes_H #define cudaTypes_H
#if defined(WM_DP) #if defined(WM_SP)
#define MPI_SCALAR MPI_DOUBLE # define MPI_SCALAR MPI_FLOAT
typedef float ValueType;
#elif defined(WM_DP)
# define MPI_SCALAR MPI_DOUBLE
typedef double ValueType; typedef double ValueType;
#elif defined(WM_SP) #elif defined(WM_LDP)
#define MPI_SCALAR MPI_FLOAT # error Long double not supported in CUDA
typedef float ValueType;
#endif #endif

43
src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C
View file

@ -163,6 +163,8 @@ void Foam::mgMeshLevel::makeChild() const
child_.setSize(nCells()); child_.setSize(nCells());
int nMoves = -1; int nMoves = -1;
# ifdef WM_DP
MGridGen MGridGen
( (
nCells(), nCells(),
@ -175,10 +177,51 @@ void Foam::mgMeshLevel::makeChild() const
mgMaxClusterSize_, mgMaxClusterSize_,
options.begin(), options.begin(),
&nMoves, &nMoves,
&nCoarseCells,
child_.begin()
);
# else
// Conversion of type for MGridGen interface
const scalarField& vols = cellVolumes();
Field<double> dblVols(vols.size());
forAll (dblVols, i)
{
dblVols[i] = vols[i];
}
Field<double> dblAreas(boundaryAreas.size());
forAll (dblAreas, i)
{
dblAreas[i] = boundaryAreas[i];
}
Field<double> dblFaceWeights(faceWeights.size());
forAll (dblFaceWeights, i)
{
dblFaceWeights[i] = faceWeights[i];
}
MGridGen
(
nCells(),
cellCellOffsets.begin(),
dblVols.begin(),
dblAreas.begin(),
cellCells.begin(),
dblFaceWeights.begin(),
mgMinClusterSize_,
mgMaxClusterSize_,
options.begin(),
&nMoves,
&nCoarseCells_, &nCoarseCells_,
child_.begin() child_.begin()
); );
# endif
Info<< "Number of coarse cells = " << nCoarseCells_ << endl; Info<< "Number of coarse cells = " << nCoarseCells_ << endl;
} }

1
src/foam/Make/files
View file

@ -26,6 +26,7 @@ $(ints)/label/label.C
$(ints)/uLabel/uLabel.C $(ints)/uLabel/uLabel.C
primitives/Scalar/doubleScalar/doubleScalar.C primitives/Scalar/doubleScalar/doubleScalar.C
primitives/Scalar/longDoubleScalar/longDoubleScalar.C
primitives/Scalar/floatScalar/floatScalar.C primitives/Scalar/floatScalar/floatScalar.C
primitives/Scalar/scalar/scalar.C primitives/Scalar/scalar/scalar.C
primitives/DiagTensor/diagTensor/diagTensor.C primitives/DiagTensor/diagTensor/diagTensor.C

4
src/foam/db/IOstreams/Fstreams/IFstream.H
View file

@ -109,8 +109,8 @@ public:
IFstream IFstream
( (
const fileName& pathname, const fileName& pathname,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
); );

7
src/foam/db/IOstreams/IOstreams/Istream.H
View file

@ -75,8 +75,8 @@ public:
//- Set stream status //- Set stream status
Istream Istream
( (
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion, versionNumber version = currentVersion,
compressionType compression=UNCOMPRESSED compressionType compression=UNCOMPRESSED
) )
: :
@ -122,6 +122,9 @@ public:
//- Read a doubleScalar //- Read a doubleScalar
virtual Istream& read(doubleScalar&) = 0; virtual Istream& read(doubleScalar&) = 0;
//- Read a longDoubleScalar
virtual Istream& read(longDoubleScalar&) = 0;
//- Read binary block //- Read binary block
virtual Istream& read(char*, std::streamsize) = 0; virtual Istream& read(char*, std::streamsize) = 0;

7
src/foam/db/IOstreams/IOstreams/Ostream.H
View file

@ -77,8 +77,8 @@ public:
//- Set stream status //- Set stream status
Ostream Ostream
( (
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion, versionNumber version = currentVersion,
compressionType compression=UNCOMPRESSED compressionType compression=UNCOMPRESSED
) )
: :
@ -131,6 +131,9 @@ public:
//- Write doubleScalar //- Write doubleScalar
virtual Ostream& write(const doubleScalar) = 0; virtual Ostream& write(const doubleScalar) = 0;
//- Write doubleScalar
virtual Ostream& write(const longDoubleScalar) = 0;
//- Write binary block //- Write binary block
virtual Ostream& write(const char*, std::streamsize) = 0; virtual Ostream& write(const char*, std::streamsize) = 0;

22
src/foam/db/IOstreams/Pstreams/IPstream.C
View file

@ -211,6 +211,21 @@ Foam::Istream& Foam::IPstream::read(token& t)
return *this; return *this;
} }
// longDoubleScalar
case token::LONG_DOUBLE_SCALAR :
{
longDoubleScalar val;
if (read(val))
{
t = val;
}
else
{
t.setBad();
}
return *this;
}
// Character (returned as a single character word) or error // Character (returned as a single character word) or error
default: default:
{ {
@ -281,6 +296,13 @@ Foam::Istream& Foam::IPstream::read(doubleScalar& val)
} }
Foam::Istream& Foam::IPstream::read(longDoubleScalar& val)
{
readFromBuffer(val);
return *this;
}
Foam::Istream& Foam::IPstream::read(char* data, std::streamsize count) Foam::Istream& Foam::IPstream::read(char* data, std::streamsize count)
{ {
if (format() != BINARY) if (format() != BINARY)

5
src/foam/db/IOstreams/Pstreams/IPstream.H
View file

@ -84,7 +84,7 @@ public:
const int fromProcNo, const int fromProcNo,
const label bufSize = 0, const label bufSize = 0,
streamFormat format=BINARY, streamFormat format=BINARY,
versionNumber version=currentVersion versionNumber version = currentVersion
); );
@ -143,6 +143,9 @@ public:
//- Read a doubleScalar //- Read a doubleScalar
Istream& read(doubleScalar&); Istream& read(doubleScalar&);
//- Read a longDoubleScalar
Istream& read(longDoubleScalar&);
//- Read binary block //- Read binary block
Istream& read(char*, std::streamsize); Istream& read(char*, std::streamsize);

8
src/foam/db/IOstreams/Pstreams/OPstream.C
View file

@ -210,6 +210,14 @@ Foam::Ostream& Foam::OPstream::write(const doubleScalar val)
} }
Foam::Ostream& Foam::OPstream::write(const longDoubleScalar val)
{
write(char(token::LONG_DOUBLE_SCALAR));
writeToBuffer(val);
return *this;
}
Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count) Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count)
{ {
if (format() != BINARY) if (format() != BINARY)

5
src/foam/db/IOstreams/Pstreams/OPstream.H
View file

@ -83,7 +83,7 @@ public:
const int toProcNo, const int toProcNo,
const label bufSize = 0, const label bufSize = 0,
streamFormat format=BINARY, streamFormat format=BINARY,
versionNumber version=currentVersion versionNumber version = currentVersion
); );
@ -152,6 +152,9 @@ public:
//- Write doubleScalar //- Write doubleScalar
Ostream& write(const doubleScalar); Ostream& write(const doubleScalar);
//- Write longDoubleScalar
Ostream& write(const longDoubleScalar);
//- Write binary block //- Write binary block
Ostream& write(const char*, std::streamsize); Ostream& write(const char*, std::streamsize);

8
src/foam/db/IOstreams/Sstreams/ISstream.C
View file

@ -489,6 +489,14 @@ Foam::Istream& Foam::ISstream::read(doubleScalar& val)
} }
Foam::Istream& Foam::ISstream::read(longDoubleScalar& val)
{
is_ >> val;
setState(is_.rdstate());
return *this;
}
// read binary block // read binary block
Foam::Istream& Foam::ISstream::read(char* buf, std::streamsize count) Foam::Istream& Foam::ISstream::read(char* buf, std::streamsize count)
{ {

9
src/foam/db/IOstreams/Sstreams/ISstream.H
View file

@ -97,9 +97,9 @@ public:
( (
istream& is, istream& is,
const string& name, const string& name,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion, versionNumber version = currentVersion,
compressionType compression=UNCOMPRESSED compressionType compression = UNCOMPRESSED
); );
@ -165,6 +165,9 @@ public:
//- Read a doubleScalar //- Read a doubleScalar
virtual Istream& read(doubleScalar&); virtual Istream& read(doubleScalar&);
//- Read a longDoubleScalar
virtual Istream& read(longDoubleScalar&);
//- Read binary block //- Read binary block
virtual Istream& read(char*, std::streamsize); virtual Istream& read(char*, std::streamsize);

8
src/foam/db/IOstreams/Sstreams/OSstream.C
View file

@ -195,6 +195,14 @@ Foam::Ostream& Foam::OSstream::write(const doubleScalar val)
} }
Foam::Ostream& Foam::OSstream::write(const longDoubleScalar val)
{
os_ << val;
setState(os_.rdstate());
return *this;
}
Foam::Ostream& Foam::OSstream::write(const char* buf, std::streamsize count) Foam::Ostream& Foam::OSstream::write(const char* buf, std::streamsize count)
{ {
if (format() != BINARY) if (format() != BINARY)

7
src/foam/db/IOstreams/Sstreams/OSstream.H
View file

@ -90,8 +90,8 @@ public:
( (
ostream& os, ostream& os,
const string& name, const string& name,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion, versionNumber version = currentVersion,
compressionType compression=UNCOMPRESSED compressionType compression=UNCOMPRESSED
); );
@ -155,6 +155,9 @@ public:
//- Write doubleScalar //- Write doubleScalar
virtual Ostream& write(const doubleScalar); virtual Ostream& write(const doubleScalar);
//- Write longDoubleScalar
virtual Ostream& write(const longDoubleScalar);
//- Write binary block //- Write binary block
virtual Ostream& write(const char*, std::streamsize); virtual Ostream& write(const char*, std::streamsize);

7
src/foam/db/IOstreams/Sstreams/prefixOSstream.C
View file

@ -146,6 +146,13 @@ Foam::Ostream& Foam::prefixOSstream::write(const doubleScalar val)
} }
Foam::Ostream& Foam::prefixOSstream::write(const longDoubleScalar val)
{
checkWritePrefix();
return OSstream::write(val);
}
Foam::Ostream& Foam::prefixOSstream::write Foam::Ostream& Foam::prefixOSstream::write
( (
const char* buf, const char* buf,

7
src/foam/db/IOstreams/Sstreams/prefixOSstream.H
View file

@ -73,8 +73,8 @@ public:
( (
ostream& os, ostream& os,
const string& name, const string& name,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion, versionNumber version = currentVersion,
compressionType compression=UNCOMPRESSED compressionType compression=UNCOMPRESSED
); );
@ -130,6 +130,9 @@ public:
//- Write doubleScalar //- Write doubleScalar
virtual Ostream& write(const doubleScalar); virtual Ostream& write(const doubleScalar);
//- Write doubleScalar
virtual Ostream& write(const longDoubleScalar);
//- Write binary block //- Write binary block
virtual Ostream& write(const char*, std::streamsize); virtual Ostream& write(const char*, std::streamsize);

8
src/foam/db/IOstreams/StringStreams/IStringStream.H
View file

@ -60,8 +60,8 @@ public:
IStringStream IStringStream
( (
const string& buffer, const string& buffer,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
) )
: :
ISstream ISstream
@ -78,8 +78,8 @@ public:
IStringStream IStringStream
( (
const char* buffer, const char* buffer,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
) )
: :
ISstream ISstream

4
src/foam/db/IOstreams/StringStreams/OStringStream.H
View file

@ -59,8 +59,8 @@ public:
//- Construct and set stream status //- Construct and set stream status
OStringStream OStringStream
( (
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
) )
: :
OSstream OSstream

7
src/foam/db/IOstreams/Tstreams/ITstream.C
View file

@ -149,6 +149,13 @@ Foam::Istream& Foam::ITstream::read(doubleScalar&)
} }
Foam::Istream& Foam::ITstream::read(longDoubleScalar&)
{
notImplemented("Istream& ITstream::read(longDoubleScalar&)");
return *this;
}
Foam::Istream& Foam::ITstream::read(char*, std::streamsize) Foam::Istream& Foam::ITstream::read(char*, std::streamsize)
{ {
notImplemented("Istream& ITstream::read(char*, std::streamsize)"); notImplemented("Istream& ITstream::read(char*, std::streamsize)");

7
src/foam/db/IOstreams/Tstreams/ITstream.H
View file

@ -70,8 +70,8 @@ public:
( (
const string& name, const string& name,
const tokenList& tokens, const tokenList& tokens,
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
) )
: :
Istream(format, version), Istream(format, version),
@ -167,6 +167,9 @@ public:
//- Read a doubleScalar //- Read a doubleScalar
virtual Istream& read(doubleScalar&); virtual Istream& read(doubleScalar&);
//- Read a longDoubleScalar
virtual Istream& read(longDoubleScalar&);
//- Read binary block //- Read binary block
virtual Istream& read(char*, std::streamsize); virtual Istream& read(char*, std::streamsize);

4
src/foam/db/IOstreams/hashes/OSHA1stream.H
View file

@ -151,8 +151,8 @@ public:
//- Construct and set stream status //- Construct and set stream status
OSHA1stream OSHA1stream
( (
streamFormat format=ASCII, streamFormat format = ASCII,
versionNumber version=currentVersion versionNumber version = currentVersion
) )
: :
OSstream OSstream

11
src/foam/db/IOstreams/token/token.H
View file

@ -81,6 +81,7 @@ public:
LABEL, LABEL,
FLOAT_SCALAR, FLOAT_SCALAR,
DOUBLE_SCALAR, DOUBLE_SCALAR,
LONG_DOUBLE_SCALAR,
COMPOUND, COMPOUND,
ERROR ERROR
@ -258,6 +259,7 @@ private:
label labelToken_; label labelToken_;
floatScalar floatScalarToken_; floatScalar floatScalarToken_;
doubleScalar doubleScalarToken_; doubleScalar doubleScalarToken_;
longDoubleScalar longDoubleScalarToken_;
mutable compound* compoundTokenPtr_; mutable compound* compoundTokenPtr_;
}; };
@ -307,6 +309,9 @@ public:
//- Construct doubleScalar token //- Construct doubleScalar token
inline token(const doubleScalar, label lineNumber=0); inline token(const doubleScalar, label lineNumber=0);
//- Construct longDoubleScalar token
inline token(const longDoubleScalar, label lineNumber=0);
//- Construct from Istream //- Construct from Istream
token(Istream&); token(Istream&);
@ -344,6 +349,9 @@ public:
inline bool isDoubleScalar() const; inline bool isDoubleScalar() const;
inline doubleScalar doubleScalarToken() const; inline doubleScalar doubleScalarToken() const;
inline bool isLongDoubleScalar() const;
inline longDoubleScalar longDoubleScalarToken() const;
inline bool isScalar() const; inline bool isScalar() const;
inline scalar scalarToken() const; inline scalar scalarToken() const;
@ -391,6 +399,7 @@ public:
inline void operator=(const label); inline void operator=(const label);
inline void operator=(const floatScalar); inline void operator=(const floatScalar);
inline void operator=(const doubleScalar); inline void operator=(const doubleScalar);
inline void operator=(const longDoubleScalar);
inline void operator=(compound*); inline void operator=(compound*);
@ -404,6 +413,7 @@ public:
inline bool operator==(const label) const; inline bool operator==(const label) const;
inline bool operator==(const floatScalar) const; inline bool operator==(const floatScalar) const;
inline bool operator==(const doubleScalar) const; inline bool operator==(const doubleScalar) const;
inline bool operator==(const longDoubleScalar) const;
// Inequality // Inequality
@ -415,6 +425,7 @@ public:
inline bool operator!=(const label) const; inline bool operator!=(const label) const;
inline bool operator!=(const floatScalar) const; inline bool operator!=(const floatScalar) const;
inline bool operator!=(const doubleScalar) const; inline bool operator!=(const doubleScalar) const;
inline bool operator!=(const longDoubleScalar) const;
// IOstream operators // IOstream operators

70
src/foam/db/IOstreams/token/tokenI.H
View file

@ -103,6 +103,10 @@ inline token::token(const token& t)
doubleScalarToken_ = t.doubleScalarToken_; doubleScalarToken_ = t.doubleScalarToken_;
break; break;
case LONG_DOUBLE_SCALAR:
longDoubleScalarToken_ = t.longDoubleScalarToken_;
break;
case COMPOUND: case COMPOUND:
compoundTokenPtr_ = t.compoundTokenPtr_; compoundTokenPtr_ = t.compoundTokenPtr_;
compoundTokenPtr_->refCount::operator++(); compoundTokenPtr_->refCount::operator++();
@ -162,6 +166,15 @@ inline token::token(const doubleScalar s, label lineNumber)
{} {}
// Construct longDoubleScalar token
inline token::token(const longDoubleScalar s, label lineNumber)
:
type_(LONG_DOUBLE_SCALAR),
longDoubleScalarToken_(s),
lineNumber_(lineNumber)
{}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
// Delete token clearing the storage used by word or string // Delete token clearing the storage used by word or string
@ -289,6 +302,7 @@ inline bool token::isDoubleScalar() const
return (type_ == DOUBLE_SCALAR); return (type_ == DOUBLE_SCALAR);
} }
inline doubleScalar token::doubleScalarToken() const inline doubleScalar token::doubleScalarToken() const
{ {
if (type_ == DOUBLE_SCALAR) if (type_ == DOUBLE_SCALAR)
@ -303,9 +317,34 @@ inline doubleScalar token::doubleScalarToken() const
} }
inline bool token::isLongDoubleScalar() const
{
return (type_ == LONG_DOUBLE_SCALAR);
}
inline longDoubleScalar token::longDoubleScalarToken() const
{
if (type_ == LONG_DOUBLE_SCALAR)
{
return longDoubleScalarToken_;
}
else
{
parseError("longDoubleScalar");
return 0.0;
}
}
inline bool token::isScalar() const inline bool token::isScalar() const
{ {
return (type_ == FLOAT_SCALAR || type_ == DOUBLE_SCALAR); return
(
type_ == FLOAT_SCALAR
|| type_ == DOUBLE_SCALAR
|| type_ == LONG_DOUBLE_SCALAR
);
} }
inline scalar token::scalarToken() const inline scalar token::scalarToken() const
@ -318,6 +357,10 @@ inline scalar token::scalarToken() const
{ {
return doubleScalarToken_; return doubleScalarToken_;
} }
else if (type_ == LONG_DOUBLE_SCALAR)
{
return longDoubleScalarToken_;
}
else else
{ {
parseError("scalar"); parseError("scalar");
@ -420,6 +463,10 @@ inline void token::operator=(const token& t)
doubleScalarToken_ = t.doubleScalarToken_; doubleScalarToken_ = t.doubleScalarToken_;
break; break;
case LONG_DOUBLE_SCALAR:
longDoubleScalarToken_ = t.longDoubleScalarToken_;
break;
case COMPOUND: case COMPOUND:
compoundTokenPtr_ = t.compoundTokenPtr_; compoundTokenPtr_ = t.compoundTokenPtr_;
compoundTokenPtr_->refCount::operator++(); compoundTokenPtr_->refCount::operator++();
@ -432,6 +479,7 @@ inline void token::operator=(const token& t)
lineNumber_ = t.lineNumber_; lineNumber_ = t.lineNumber_;
} }
inline void token::operator=(const punctuationToken p) inline void token::operator=(const punctuationToken p)
{ {
clear(); clear();
@ -484,6 +532,13 @@ inline void token::operator=(const doubleScalar s)
doubleScalarToken_ = s; doubleScalarToken_ = s;
} }
inline void token::operator=(const longDoubleScalar s)
{
clear();
type_ = LONG_DOUBLE_SCALAR;
longDoubleScalarToken_ = s;
}
inline void token::operator=(token::compound* cPtr) inline void token::operator=(token::compound* cPtr)
{ {
clear(); clear();
@ -522,6 +577,9 @@ inline bool token::operator==(const token& t) const
case DOUBLE_SCALAR: case DOUBLE_SCALAR:
return equal(doubleScalarToken_, t.doubleScalarToken_); return equal(doubleScalarToken_, t.doubleScalarToken_);
case LONG_DOUBLE_SCALAR:
return equal(longDoubleScalarToken_, t.longDoubleScalarToken_);
case COMPOUND: case COMPOUND:
return compoundTokenPtr_ == t.compoundTokenPtr_; return compoundTokenPtr_ == t.compoundTokenPtr_;
@ -562,6 +620,11 @@ inline bool token::operator==(const doubleScalar s) const
return (type_ == DOUBLE_SCALAR && equal(doubleScalarToken_, s)); return (type_ == DOUBLE_SCALAR && equal(doubleScalarToken_, s));
} }
inline bool token::operator==(const longDoubleScalar s) const
{
return (type_ == LONG_DOUBLE_SCALAR && equal(longDoubleScalarToken_, s));
}
inline bool token::operator!=(const token& t) const inline bool token::operator!=(const token& t) const
{ {
return !operator==(t); return !operator==(t);
@ -592,6 +655,11 @@ inline bool token::operator!=(const doubleScalar s) const
return !operator==(s); return !operator==(s);
} }
inline bool token::operator!=(const longDoubleScalar s) const
{
return !operator==(s);
}
inline bool token::operator!=(const label l) const inline bool token::operator!=(const label l) const
{ {
return !operator==(l); return !operator==(l);

12
src/foam/db/IOstreams/token/tokenIO.C
View file

@ -85,6 +85,10 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const token& t)
os << t.doubleScalarToken_; os << t.doubleScalarToken_;
break; break;
case token::LONG_DOUBLE_SCALAR:
os << t.longDoubleScalarToken_;
break;
case token::COMPOUND: case token::COMPOUND:
os << *t.compoundTokenPtr_; os << *t.compoundTokenPtr_;
break; break;
@ -168,6 +172,10 @@ ostream& Foam::operator<<(ostream& os, const InfoProxy<token>& ip)
os << " the doubleScalar " << t.doubleScalarToken(); os << " the doubleScalar " << t.doubleScalarToken();
break; break;
case token::LONG_DOUBLE_SCALAR:
os << " the longDoubleScalar " << t.doubleScalarToken();
break;
case token::COMPOUND: case token::COMPOUND:
{ {
if (t.compoundToken().empty()) if (t.compoundToken().empty())
@ -238,6 +246,10 @@ Ostream& operator<<(Ostream& os, const InfoProxy<token>& ip)
os << " the doubleScalar " << t.doubleScalarToken(); os << " the doubleScalar " << t.doubleScalarToken();
break; break;
case token::LONG_DOUBLE_SCALAR:
os << " the longDoubleScalar " << t.longDoubleScalarToken();
break;
case token::COMPOUND: case token::COMPOUND:
{ {
if (t.compoundToken().empty()) if (t.compoundToken().empty())

2
src/foam/db/Time/Time.C
View file

@ -75,7 +75,7 @@ void Foam::Time::adjustDeltaT()
{ {
scalar timeToNextWrite = max scalar timeToNextWrite = max
( (
0.0, scalar(0),
(outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_) (outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_)
); );

4
src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
View file

@ -516,8 +516,8 @@ evaluate
// ZT, 22-07-2014 - point ordering is not same // ZT, 22-07-2014 - point ordering is not same
// at master and slave side after topology change // at master and slave side after topology change
const labelList& neiPoints = // const labelList& neiPoints =
procPatch_.procPolyPatch().neighbPoints(); // procPatch_.procPolyPatch().neighbPoints();
// VS, 2015-04-12 - This doesn't work in parallel! // VS, 2015-04-12 - This doesn't work in parallel!
// tpn = // tpn =

4
src/foam/meshes/meshShapes/face/face.C
View file

@ -225,7 +225,9 @@ Foam::label Foam::face::split
splitEdge /= Foam::mag(splitEdge) + VSMALL; splitEdge /= Foam::mag(splitEdge) + VSMALL;
const scalar splitCos = splitEdge & rightEdge; const scalar splitCos = splitEdge & rightEdge;
const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos))); const scalar splitAngle =
acos(max(scalar(-1), min(scalar(1), splitCos)));
const scalar angleDiff = fabs(splitAngle - bisectAngle); const scalar angleDiff = fabs(splitAngle - bisectAngle);
if (angleDiff < minDiff) if (angleDiff < minDiff)

10
src/foam/meshes/primitiveShapes/triangle/triangleI.H
View file

@ -400,7 +400,7 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
{ {
scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom; scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
return 0.5*Foam::sqrt(min(GREAT, max(0, a))); return 0.5*Foam::sqrt(min(GREAT, max(scalar(0), a)));
} }
} }
@ -411,9 +411,7 @@ inline scalar triangle<Point, PointRef>::quality() const
return return
mag() mag()
/( /(
mathematicalConstant::pi 0.413497*mathematicalConstant::pi*sqr(circumRadius())
*Foam::sqr(circumRadius())
*0.413497
+ VSMALL + VSMALL
); );
} }
@ -775,8 +773,8 @@ inline bool triangle<Point, PointRef>::classify
if (hit) if (hit)
{ {
// alpha,beta might get negative due to precision errors // alpha,beta might get negative due to precision errors
alpha = max(0.0, min(1.0, alpha)); alpha = max(scalar(0), min(scalar(1), alpha));
beta = max(0.0, min(1.0, beta)); beta = max(scalar(0), min(scalar(1), beta));
} }
nearType = NONE; nearType = NONE;

4
src/foam/primitives/Scalar/Scalar.H
View file

@ -161,11 +161,11 @@ inline Scalar sqrtSumSqr(const Scalar a, const Scalar b)
if (maga > magb) if (maga > magb)
{ {
return maga*sqrt(1.0 + sqr(magb/maga)); return maga*::sqrt(1.0 + sqr(magb/maga));
} }
else else
{ {
return magb < ScalarVSMALL ? 0.0 : magb*sqrt(1.0 + sqr(maga/magb)); return magb < ScalarVSMALL ? 0.0 : magb*::sqrt(1.0 + sqr(maga/magb));
} }
} }

7
src/foam/primitives/Scalar/doubleFloat.H
View file

@ -84,6 +84,13 @@ inline double pow(const type1 s, const type2 e) \
return ::pow(double(s), double(e)); \ return ::pow(double(s), double(e)); \
} }
MAXMINPOW(long double, long double, long double)
MAXMINPOW(long double, long double, double)
MAXMINPOW(long double, long double, float)
MAXMINPOW(long double, long double, int)
MAXMINPOW(long double, double, long double)
MAXMINPOW(long double, float, long double)
MAXMINPOW(long double, int, long double)
MAXMINPOW(double, double, double) MAXMINPOW(double, double, double)
MAXMINPOW(double, double, float) MAXMINPOW(double, double, float)

43
src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.C Normal file
View file

@ -0,0 +1,43 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "longDoubleScalar.H"
#include "IOstreams.H"
#include <sstream>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define Scalar longDoubleScalar
#define ScalarVGREAT longDoubleScalarVGREAT
#define ScalarVSMALL longDoubleScalarVSMALL
#define readScalar readLongDoubleScalar
#include "Scalar.C"
#undef Scalar
#undef ScalarVGREAT
#undef ScalarVSMALL
#undef readScalar
// ************************************************************************* //

115
src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.H Normal file
View file

@ -0,0 +1,115 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::doubleScalar
Description
Long double precision floating point scalar type.
SourceFiles
longDoubleScalar.C
\*---------------------------------------------------------------------------*/
#ifndef longDoubleScalar_H
#define longDoubleScalar_H
#include "doubleFloat.H"
#include "direction.H"
#include "word.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
typedef long double longDoubleScalar;
// Largest and smallest scalar values allowed in certain parts of the code.
// (15 is the number of significant figures in an
// IEEE double precision number. See limits.h or float.h)
static const longDoubleScalar longDoubleScalarGREAT = 1.0e+15;
static const longDoubleScalar longDoubleScalarVGREAT = 1.0e+300;
static const longDoubleScalar longDoubleScalarROOTVGREAT = 1.0e+150;
static const longDoubleScalar longDoubleScalarSMALL = 1.0e-15;
static const longDoubleScalar longDoubleScalarVSMALL = 1.0e-300;
static const longDoubleScalar longDoubleScalarROOTVSMALL = 1.0e-150;
#define Scalar longDoubleScalar
#define ScalarVGREAT longDoubleScalarVGREAT
#define ScalarVSMALL longDoubleScalarVSMALL
#define readScalar readLongDoubleScalar
inline Scalar mag(const Scalar s)
{
return ::fabs(s);
}
#define transFunc(func) \
inline Scalar func(const Scalar s) \
{ \
return ::func(s); \
}
#include "Scalar.H"
inline Scalar hypot(const Scalar x, const Scalar y)
{
return ::hypot(x, y);
}
inline Scalar atan2(const Scalar y, const Scalar x)
{
return ::atan2(y, x);
}
inline Scalar jn(const int n, const Scalar s)
{
return ::jn(n, s);
}
inline Scalar yn(const int n, const Scalar s)
{
return ::yn(n, s);
}
#undef Scalar
#undef ScalarVGREAT
#undef ScalarVSMALL
#undef readScalar
#undef transFunc
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

19
src/foam/primitives/Scalar/scalar/scalar.H
View file

@ -38,6 +38,7 @@ SourceFiles
#include "floatScalar.H" #include "floatScalar.H"
#include "doubleScalar.H" #include "doubleScalar.H"
#include "longDoubleScalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -77,6 +78,24 @@ namespace Foam
scalar readScalar(Istream& is); scalar readScalar(Istream& is);
} }
#elif defined(WM_LDP)
// Define scalar as a long double
namespace Foam
{
typedef longDoubleScalar scalar;
static const scalar GREAT = longDoubleScalarGREAT;
static const scalar VGREAT = longDoubleScalarVGREAT;
static const scalar ROOTVGREAT = longDoubleScalarROOTVGREAT;
static const scalar SMALL = longDoubleScalarSMALL;
static const scalar VSMALL = longDoubleScalarVSMALL;
static const scalar ROOTVSMALL = longDoubleScalarROOTVSMALL;
scalar readScalar(Istream& is);
}
#endif #endif

4
src/foam/primitives/Tensor/tensor/tensor.C
View file

@ -118,7 +118,7 @@ vector eigenValues(const tensor& t)
if (disc >= -SMALL) if (disc >= -SMALL)
{ {
scalar q = -0.5*sqrt(max(0.0, disc)); scalar q = -0.5*sqrt(max(scalar(0), disc));
i = 0; i = 0;
ii = -0.5*a + q; ii = -0.5*a + q;
@ -381,7 +381,7 @@ vector eigenValues(const symmTensor& t)
// If there is a zero root // If there is a zero root
if (mag(c) < SMALL) if (mag(c) < SMALL)
{ {
const scalar disc = Foam::max(sqr(a) - 4*b, 0.0); const scalar disc = Foam::max(sqr(a) - 4*b, scalar(0));
scalar q = -0.5*sqrt(max(scalar(0), disc)); scalar q = -0.5*sqrt(max(scalar(0), disc));

291
src/foam/primitives/transform/symmTransform.H
View file

@ -117,163 +117,163 @@ inline SymmTensor4thOrder<Cmpt> transform
const SymmTensor4thOrder<Cmpt>& C const SymmTensor4thOrder<Cmpt>& C
) )
{ {
//- represent fourth order tensors in 6x6 notation then rotation is given by //- represent fourth order tensors in 6x6 notation. Rotation is given by
//- C_rotated_af = A_ba * C_cd * A_ef //- C_rotated_af = A_ba * C_cd * A_ef
//- where A is a function of stt //- where A is a function of stt
const double s = ::sqrt(2); const scalar s = ::sqrt(2);
const double A[6][6] = const scalar A[6][6] =
{ {
{ stt.xx()*stt.xx(), stt.xy()*stt.xy(), stt.xz()*stt.xz(), s*stt.xx()*stt.xy(), s*stt.xy()*stt.xz(), s*stt.xz()*stt.xx() }, { stt.xx()*stt.xx(), stt.xy()*stt.xy(), stt.xz()*stt.xz(), s*stt.xx()*stt.xy(), s*stt.xy()*stt.xz(), s*stt.xz()*stt.xx() },
{ stt.xy()*stt.xy(), stt.yy()*stt.yy(), stt.yz()*stt.yz(), s*stt.xy()*stt.yy(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.xy() }, { stt.xy()*stt.xy(), stt.yy()*stt.yy(), stt.yz()*stt.yz(), s*stt.xy()*stt.yy(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.xy() },
{ stt.xz()*stt.xz(), stt.yz()*stt.yz(), stt.zz()*stt.zz(), s*stt.xz()*stt.yz(), s*stt.yz()*stt.zz(), s*stt.zz()*stt.xz() }, { stt.xz()*stt.xz(), stt.yz()*stt.yz(), stt.zz()*stt.zz(), s*stt.xz()*stt.yz(), s*stt.yz()*stt.zz(), s*stt.zz()*stt.xz() },
{ s*stt.xx()*stt.xy(), s*stt.xy()*stt.yy(), s*stt.xz()*stt.yz(), { s*stt.xx()*stt.xy(), s*stt.xy()*stt.yy(), s*stt.xz()*stt.yz(),
(stt.xy()*stt.xy()+stt.xx()*stt.yy()), (stt.xz()*stt.yy()+stt.xy()*stt.yz()), (stt.xx()*stt.yz()+stt.xz()*stt.xy()) }, (stt.xy()*stt.xy()+stt.xx()*stt.yy()), (stt.xz()*stt.yy()+stt.xy()*stt.yz()), (stt.xx()*stt.yz()+stt.xz()*stt.xy()) },
{ s*stt.xy()*stt.xz(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.zz(), { s*stt.xy()*stt.xz(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.zz(),
(stt.yy()*stt.xz()+stt.xy()*stt.yz()), (stt.yz()*stt.yz()+stt.yy()*stt.zz()), (stt.xy()*stt.zz()+stt.yz()*stt.xz()) }, (stt.yy()*stt.xz()+stt.xy()*stt.yz()), (stt.yz()*stt.yz()+stt.yy()*stt.zz()), (stt.xy()*stt.zz()+stt.yz()*stt.xz()) },
{ s*stt.xz()*stt.xx(), s*stt.yz()*stt.xy(), s*stt.zz()*stt.xz(), { s*stt.xz()*stt.xx(), s*stt.yz()*stt.xy(), s*stt.zz()*stt.xz(),
(stt.yz()*stt.xx()+stt.xz()*stt.xy()), (stt.zz()*stt.xy()+stt.yz()*stt.xz()), (stt.xz()*stt.xz()+stt.zz()*stt.xx()) } (stt.yz()*stt.xx()+stt.xz()*stt.xy()), (stt.zz()*stt.xy()+stt.yz()*stt.xz()), (stt.xz()*stt.xz()+stt.zz()*stt.xx()) }
}; };
return symmTensor4thOrder return symmTensor4thOrder
( (
// xxxx // xxxx
A[0][0] * C.xxxx() * A[0][0] + A[0][0] * C.xxxx() * A[0][0] +
A[1][0] * C.xxyy() * A[0][0] + A[1][0] * C.xxyy() * A[0][0] +
A[2][0] * C.xxzz() * A[0][0] + A[2][0] * C.xxzz() * A[0][0] +
A[0][0] * C.xxyy() * A[1][0] + A[0][0] * C.xxyy() * A[1][0] +
A[1][0] * C.yyyy() * A[1][0] + A[1][0] * C.yyyy() * A[1][0] +
A[2][0] * C.yyzz() * A[1][0] + A[2][0] * C.yyzz() * A[1][0] +
A[0][0] * C.xxzz() * A[2][0] + A[0][0] * C.xxzz() * A[2][0] +
A[1][0] * C.yyzz() * A[2][0] + A[1][0] * C.yyzz() * A[2][0] +
A[2][0] * C.zzzz() * A[2][0] + A[2][0] * C.zzzz() * A[2][0] +
A[3][0] * C.xyxy() * A[3][0] + A[3][0] * C.xyxy() * A[3][0] +
A[4][0] * C.yzyz() * A[4][0] + A[4][0] * C.yzyz() * A[4][0] +
A[5][0] * C.zxzx() * A[5][0], A[5][0] * C.zxzx() * A[5][0],
// xxyy // xxyy
A[0][0] * C.xxxx() * A[0][1] + A[0][0] * C.xxxx() * A[0][1] +
A[1][0] * C.xxyy() * A[0][1] + A[1][0] * C.xxyy() * A[0][1] +
A[2][0] * C.xxzz() * A[0][1] + A[2][0] * C.xxzz() * A[0][1] +
A[0][0] * C.xxyy() * A[1][1] + A[0][0] * C.xxyy() * A[1][1] +
A[1][0] * C.yyyy() * A[1][1] + A[1][0] * C.yyyy() * A[1][1] +
A[2][0] * C.yyzz() * A[1][1] + A[2][0] * C.yyzz() * A[1][1] +
A[0][0] * C.xxzz() * A[2][1] + A[0][0] * C.xxzz() * A[2][1] +
A[1][0] * C.yyzz() * A[2][1] + A[1][0] * C.yyzz() * A[2][1] +
A[2][0] * C.zzzz() * A[2][1] + A[2][0] * C.zzzz() * A[2][1] +
A[3][0] * C.xyxy() * A[3][1] + A[3][0] * C.xyxy() * A[3][1] +
A[4][0] * C.yzyz() * A[4][1] + A[4][0] * C.yzyz() * A[4][1] +
A[5][0] * C.zxzx() * A[5][1], A[5][0] * C.zxzx() * A[5][1],
// xzz // xzz
A[0][0] * C.xxxx() * A[0][2] + A[0][0] * C.xxxx() * A[0][2] +
A[1][0] * C.xxyy() * A[0][2] + A[1][0] * C.xxyy() * A[0][2] +
A[2][0] * C.xxzz() * A[0][2] + A[2][0] * C.xxzz() * A[0][2] +
A[0][0] * C.xxyy() * A[1][2] + A[0][0] * C.xxyy() * A[1][2] +
A[1][0] * C.yyyy() * A[1][2] + A[1][0] * C.yyyy() * A[1][2] +
A[2][0] * C.yyzz() * A[1][2] + A[2][0] * C.yyzz() * A[1][2] +
A[0][0] * C.xxzz() * A[2][2] + A[0][0] * C.xxzz() * A[2][2] +
A[1][0] * C.yyzz() * A[2][2] + A[1][0] * C.yyzz() * A[2][2] +
A[2][0] * C.zzzz() * A[2][2] + A[2][0] * C.zzzz() * A[2][2] +
A[3][0] * C.xyxy() * A[3][2] + A[3][0] * C.xyxy() * A[3][2] +
A[4][0] * C.yzyz() * A[4][2] + A[4][0] * C.yzyz() * A[4][2] +
A[5][0] * C.zxzx() * A[5][2], A[5][0] * C.zxzx() * A[5][2],
// yyyy // yyyy
A[0][1] * C.xxxx() * A[0][1] + A[0][1] * C.xxxx() * A[0][1] +
A[1][1] * C.xxyy() * A[0][1] + A[1][1] * C.xxyy() * A[0][1] +
A[2][1] * C.xxzz() * A[0][1] + A[2][1] * C.xxzz() * A[0][1] +
A[0][1] * C.xxyy() * A[1][1] + A[0][1] * C.xxyy() * A[1][1] +
A[1][1] * C.yyyy() * A[1][1] + A[1][1] * C.yyyy() * A[1][1] +
A[2][1] * C.yyzz() * A[1][1] + A[2][1] * C.yyzz() * A[1][1] +
A[0][1] * C.xxzz() * A[2][1] + A[0][1] * C.xxzz() * A[2][1] +
A[1][1] * C.yyzz() * A[2][1] + A[1][1] * C.yyzz() * A[2][1] +
A[2][1] * C.zzzz() * A[2][1] + A[2][1] * C.zzzz() * A[2][1] +
A[3][1] * C.xyxy() * A[3][1] + A[3][1] * C.xyxy() * A[3][1] +
A[4][1] * C.yzyz() * A[4][1] + A[4][1] * C.yzyz() * A[4][1] +
A[5][1] * C.zxzx() * A[5][1], A[5][1] * C.zxzx() * A[5][1],
// yyzz // yyzz
A[0][1] * C.xxxx() * A[0][2] + A[0][1] * C.xxxx() * A[0][2] +
A[1][1] * C.xxyy() * A[0][2] + A[1][1] * C.xxyy() * A[0][2] +
A[2][1] * C.xxzz() * A[0][2] + A[2][1] * C.xxzz() * A[0][2] +
A[0][1] * C.xxyy() * A[1][2] + A[0][1] * C.xxyy() * A[1][2] +
A[1][1] * C.yyyy() * A[1][2] + A[1][1] * C.yyyy() * A[1][2] +
A[2][1] * C.yyzz() * A[1][2] + A[2][1] * C.yyzz() * A[1][2] +
A[0][1] * C.xxzz() * A[2][2] + A[0][1] * C.xxzz() * A[2][2] +
A[1][1] * C.yyzz() * A[2][2] + A[1][1] * C.yyzz() * A[2][2] +
A[2][1] * C.zzzz() * A[2][2] + A[2][1] * C.zzzz() * A[2][2] +
A[3][1] * C.xyxy() * A[3][2] + A[3][1] * C.xyxy() * A[3][2] +
A[4][1] * C.yzyz() * A[4][2] + A[4][1] * C.yzyz() * A[4][2] +
A[5][1] * C.zxzx() * A[5][2], A[5][1] * C.zxzx() * A[5][2],
// zzzz // zzzz
A[0][2] * C.xxxx() * A[0][2] + A[0][2] * C.xxxx() * A[0][2] +
A[1][2] * C.xxyy() * A[0][2] + A[1][2] * C.xxyy() * A[0][2] +
A[2][2] * C.xxzz() * A[0][2] + A[2][2] * C.xxzz() * A[0][2] +
A[0][2] * C.xxyy() * A[1][2] + A[0][2] * C.xxyy() * A[1][2] +
A[1][2] * C.yyyy() * A[1][2] + A[1][2] * C.yyyy() * A[1][2] +
A[2][2] * C.yyzz() * A[1][2] + A[2][2] * C.yyzz() * A[1][2] +
A[0][2] * C.xxzz() * A[2][2] + A[0][2] * C.xxzz() * A[2][2] +
A[1][2] * C.yyzz() * A[2][2] + A[1][2] * C.yyzz() * A[2][2] +
A[2][2] * C.zzzz() * A[2][2] + A[2][2] * C.zzzz() * A[2][2] +
A[3][2] * C.xyxy() * A[3][2] + A[3][2] * C.xyxy() * A[3][2] +
A[4][2] * C.yzyz() * A[4][2] + A[4][2] * C.yzyz() * A[4][2] +
A[5][2] * C.zxzx() * A[5][2], A[5][2] * C.zxzx() * A[5][2],
// xyxy // xyxy
A[0][3] * C.xxxx() * A[0][3] + A[0][3] * C.xxxx() * A[0][3] +
A[1][3] * C.xxyy() * A[0][3] + A[1][3] * C.xxyy() * A[0][3] +
A[2][3] * C.xxzz() * A[0][3] + A[2][3] * C.xxzz() * A[0][3] +
A[0][3] * C.xxyy() * A[1][3] + A[0][3] * C.xxyy() * A[1][3] +
A[1][3] * C.yyyy() * A[1][3] + A[1][3] * C.yyyy() * A[1][3] +
A[2][3] * C.yyzz() * A[1][3] + A[2][3] * C.yyzz() * A[1][3] +
A[0][3] * C.xxzz() * A[2][3] + A[0][3] * C.xxzz() * A[2][3] +
A[1][3] * C.yyzz() * A[2][3] + A[1][3] * C.yyzz() * A[2][3] +
A[2][3] * C.zzzz() * A[2][3] + A[2][3] * C.zzzz() * A[2][3] +
A[3][3] * C.xyxy() * A[3][3] + A[3][3] * C.xyxy() * A[3][3] +
A[4][3] * C.yzyz() * A[4][3] + A[4][3] * C.yzyz() * A[4][3] +
A[5][3] * C.zxzx() * A[5][3], A[5][3] * C.zxzx() * A[5][3],
// yzyz // yzyz
A[0][4] * C.xxxx() * A[0][4] + A[0][4] * C.xxxx() * A[0][4] +
A[1][4] * C.xxyy() * A[0][4] + A[1][4] * C.xxyy() * A[0][4] +
A[2][4] * C.xxzz() * A[0][4] + A[2][4] * C.xxzz() * A[0][4] +
A[0][4] * C.xxyy() * A[1][4] + A[0][4] * C.xxyy() * A[1][4] +
A[1][4] * C.yyyy() * A[1][4] + A[1][4] * C.yyyy() * A[1][4] +
A[2][4] * C.yyzz() * A[1][4] + A[2][4] * C.yyzz() * A[1][4] +
A[0][4] * C.xxzz() * A[2][4] + A[0][4] * C.xxzz() * A[2][4] +
A[1][4] * C.yyzz() * A[2][4] + A[1][4] * C.yyzz() * A[2][4] +
A[2][4] * C.zzzz() * A[2][4] + A[2][4] * C.zzzz() * A[2][4] +
A[3][4] * C.xyxy() * A[3][4] + A[3][4] * C.xyxy() * A[3][4] +
A[4][4] * C.yzyz() * A[4][4] + A[4][4] * C.yzyz() * A[4][4] +
A[5][4] * C.zxzx() * A[5][4], A[5][4] * C.zxzx() * A[5][4],
// zxzx // zxzx
A[0][5] * C.xxxx() * A[0][5] + A[0][5] * C.xxxx() * A[0][5] +
A[1][5] * C.xxyy() * A[0][5] + A[1][5] * C.xxyy() * A[0][5] +
A[2][5] * C.xxzz() * A[0][5] + A[2][5] * C.xxzz() * A[0][5] +
A[0][5] * C.xxyy() * A[1][5] + A[0][5] * C.xxyy() * A[1][5] +
A[1][5] * C.yyyy() * A[1][5] + A[1][5] * C.yyyy() * A[1][5] +
A[2][5] * C.yyzz() * A[1][5] + A[2][5] * C.yyzz() * A[1][5] +
A[0][5] * C.xxzz() * A[2][5] + A[0][5] * C.xxzz() * A[2][5] +
A[1][5] * C.yyzz() * A[2][5] + A[1][5] * C.yyzz() * A[2][5] +
A[2][5] * C.zzzz() * A[2][5] + A[2][5] * C.zzzz() * A[2][5] +
A[3][5] * C.xyxy() * A[3][5] + A[3][5] * C.xyxy() * A[3][5] +
A[4][5] * C.yzyz() * A[4][5] + A[4][5] * C.yzyz() * A[4][5] +
A[5][5] * C.zxzx() * A[5][5] A[5][5] * C.zxzx() * A[5][5]
); );
} }
template<class Cmpt> template<class Cmpt>
inline DiagTensor<Cmpt> transform inline DiagTensor<Cmpt> transform
( (
const symmTensor& stt, const symmTensor& stt,
const DiagTensor<Cmpt>& dt const DiagTensor<Cmpt>& dt
) )
{ {
return dt; return dt;
} }
@ -328,21 +328,26 @@ inline symmTensor transformMask<symmTensor>(const symmTensor& st)
template<> template<>
inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const symmTensor& st) inline symmTensor4thOrder transformMask<symmTensor4thOrder>
(
const symmTensor& st
)
{ {
notImplemented("symmTransform.H\n" notImplemented
"template<>\n" (
"inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const symmTensor& st)\n" "template<>\n"
"not implemented"); "inline symmTensor4thOrder transformMask<symmTensor4thOrder>"
"(const symmTensor& st)"
);
return symmTensor4thOrder::zero; return symmTensor4thOrder::zero;
} }
template<> template<>
inline diagTensor transformMask<diagTensor>(const symmTensor& st) inline diagTensor transformMask<diagTensor>(const symmTensor& st)
{ {
return diagTensor(st.xx(), st.yy(), st.zz()); return diagTensor(st.xx(), st.yy(), st.zz());
} }

267
src/foam/primitives/transform/transform.H
View file

@ -146,151 +146,156 @@ inline SphericalTensor<Cmpt> transform
template<class Cmpt> template<class Cmpt>
inline SymmTensor4thOrder<Cmpt> transform(const tensor& L, const SymmTensor4thOrder<Cmpt>& C) inline SymmTensor4thOrder<Cmpt> transform
(
const tensor& L,
const SymmTensor4thOrder<Cmpt>& C
)
{ {
//- represent fourth order tensors in 6x6 notation then rotation is given by //- represent fourth order tensors in 6x6 notation. Rotation is given by
//- C_rotated_af = A_ba * C_cd * A_ef //- C_rotated_af = A_ba * C_cd * A_ef
//- where A is a function of L //- where A is a function of L
const double s = ::sqrt(2); const scalar s = ::sqrt(2);
const double A[6][6] = const scalar A[6][6] =
{ {
{ L.xx()*L.xx(), L.xy()*L.xy(), L.xz()*L.xz(), s*L.xx()*L.xy(), s*L.xy()*L.xz(), s*L.xz()*L.xx() }, {
{ L.yx()*L.yx(), L.yy()*L.yy(), L.yz()*L.yz(), s*L.yx()*L.yy(), s*L.yy()*L.yz(), s*L.yz()*L.yx() }, L.xx()*L.xx(), L.xy()*L.xy(), L.xz()*L.xz(), s*L.xx()*L.xy(), s*L.xy()*L.xz(), s*L.xz()*L.xx() },
{ L.zx()*L.zx(), L.zy()*L.zy(), L.zz()*L.zz(), s*L.zx()*L.zy(), s*L.zy()*L.zz(), s*L.zz()*L.zx() }, { L.yx()*L.yx(), L.yy()*L.yy(), L.yz()*L.yz(), s*L.yx()*L.yy(), s*L.yy()*L.yz(), s*L.yz()*L.yx() },
{ s*L.xx()*L.yx(), s*L.xy()*L.yy(), s*L.xz()*L.yz(), (L.xy()*L.yx()+L.xx()*L.yy()), (L.xz()*L.yy()+L.xy()*L.yz()), (L.xx()*L.yz()+L.xz()*L.yx()) }, { L.zx()*L.zx(), L.zy()*L.zy(), L.zz()*L.zz(), s*L.zx()*L.zy(), s*L.zy()*L.zz(), s*L.zz()*L.zx() },
{ s*L.yx()*L.zx(), s*L.yy()*L.zy(), s*L.yz()*L.zz(), (L.yy()*L.zx()+L.yx()*L.zy()), (L.yz()*L.zy()+L.yy()*L.zz()), (L.yx()*L.zz()+L.yz()*L.zx()) }, { s*L.xx()*L.yx(), s*L.xy()*L.yy(), s*L.xz()*L.yz(), (L.xy()*L.yx()+L.xx()*L.yy()), (L.xz()*L.yy()+L.xy()*L.yz()), (L.xx()*L.yz()+L.xz()*L.yx()) },
{ s*L.zx()*L.xx(), s*L.zy()*L.xy(), s*L.zz()*L.xz(), (L.zy()*L.xx()+L.zx()*L.xy()), (L.zz()*L.xy()+L.zy()*L.xz()), (L.zx()*L.xz()+L.zz()*L.xx()) } { s*L.yx()*L.zx(), s*L.yy()*L.zy(), s*L.yz()*L.zz(), (L.yy()*L.zx()+L.yx()*L.zy()), (L.yz()*L.zy()+L.yy()*L.zz()), (L.yx()*L.zz()+L.yz()*L.zx()) },
{ s*L.zx()*L.xx(), s*L.zy()*L.xy(), s*L.zz()*L.xz(), (L.zy()*L.xx()+L.zx()*L.xy()), (L.zz()*L.xy()+L.zy()*L.xz()), (L.zx()*L.xz()+L.zz()*L.xx()) }
}; };
return symmTensor4thOrder return symmTensor4thOrder
( (
// xxxx // xxxx
A[0][0] * C.xxxx() * A[0][0] + A[0][0] * C.xxxx() * A[0][0] +
A[1][0] * C.xxyy() * A[0][0] + A[1][0] * C.xxyy() * A[0][0] +
A[2][0] * C.xxzz() * A[0][0] + A[2][0] * C.xxzz() * A[0][0] +
A[0][0] * C.xxyy() * A[1][0] + A[0][0] * C.xxyy() * A[1][0] +
A[1][0] * C.yyyy() * A[1][0] + A[1][0] * C.yyyy() * A[1][0] +
A[2][0] * C.yyzz() * A[1][0] + A[2][0] * C.yyzz() * A[1][0] +
A[0][0] * C.xxzz() * A[2][0] + A[0][0] * C.xxzz() * A[2][0] +
A[1][0] * C.yyzz() * A[2][0] + A[1][0] * C.yyzz() * A[2][0] +
A[2][0] * C.zzzz() * A[2][0] + A[2][0] * C.zzzz() * A[2][0] +
A[3][0] * C.xyxy() * A[3][0] + A[3][0] * C.xyxy() * A[3][0] +
A[4][0] * C.yzyz() * A[4][0] + A[4][0] * C.yzyz() * A[4][0] +
A[5][0] * C.zxzx() * A[5][0], A[5][0] * C.zxzx() * A[5][0],
// xxyy // xxyy
A[0][0] * C.xxxx() * A[0][1] + A[0][0] * C.xxxx() * A[0][1] +
A[1][0] * C.xxyy() * A[0][1] + A[1][0] * C.xxyy() * A[0][1] +
A[2][0] * C.xxzz() * A[0][1] + A[2][0] * C.xxzz() * A[0][1] +
A[0][0] * C.xxyy() * A[1][1] + A[0][0] * C.xxyy() * A[1][1] +
A[1][0] * C.yyyy() * A[1][1] + A[1][0] * C.yyyy() * A[1][1] +
A[2][0] * C.yyzz() * A[1][1] + A[2][0] * C.yyzz() * A[1][1] +
A[0][0] * C.xxzz() * A[2][1] + A[0][0] * C.xxzz() * A[2][1] +
A[1][0] * C.yyzz() * A[2][1] + A[1][0] * C.yyzz() * A[2][1] +
A[2][0] * C.zzzz() * A[2][1] + A[2][0] * C.zzzz() * A[2][1] +
A[3][0] * C.xyxy() * A[3][1] + A[3][0] * C.xyxy() * A[3][1] +
A[4][0] * C.yzyz() * A[4][1] + A[4][0] * C.yzyz() * A[4][1] +
A[5][0] * C.zxzx() * A[5][1], A[5][0] * C.zxzx() * A[5][1],
// xzz // xzz
A[0][0] * C.xxxx() * A[0][2] + A[0][0] * C.xxxx() * A[0][2] +
A[1][0] * C.xxyy() * A[0][2] + A[1][0] * C.xxyy() * A[0][2] +
A[2][0] * C.xxzz() * A[0][2] + A[2][0] * C.xxzz() * A[0][2] +
A[0][0] * C.xxyy() * A[1][2] + A[0][0] * C.xxyy() * A[1][2] +
A[1][0] * C.yyyy() * A[1][2] + A[1][0] * C.yyyy() * A[1][2] +
A[2][0] * C.yyzz() * A[1][2] + A[2][0] * C.yyzz() * A[1][2] +
A[0][0] * C.xxzz() * A[2][2] + A[0][0] * C.xxzz() * A[2][2] +
A[1][0] * C.yyzz() * A[2][2] + A[1][0] * C.yyzz() * A[2][2] +
A[2][0] * C.zzzz() * A[2][2] + A[2][0] * C.zzzz() * A[2][2] +
A[3][0] * C.xyxy() * A[3][2] + A[3][0] * C.xyxy() * A[3][2] +
A[4][0] * C.yzyz() * A[4][2] + A[4][0] * C.yzyz() * A[4][2] +
A[5][0] * C.zxzx() * A[5][2], A[5][0] * C.zxzx() * A[5][2],
// yyyy // yyyy
A[0][1] * C.xxxx() * A[0][1] + A[0][1] * C.xxxx() * A[0][1] +
A[1][1] * C.xxyy() * A[0][1] + A[1][1] * C.xxyy() * A[0][1] +
A[2][1] * C.xxzz() * A[0][1] + A[2][1] * C.xxzz() * A[0][1] +
A[0][1] * C.xxyy() * A[1][1] + A[0][1] * C.xxyy() * A[1][1] +
A[1][1] * C.yyyy() * A[1][1] + A[1][1] * C.yyyy() * A[1][1] +
A[2][1] * C.yyzz() * A[1][1] + A[2][1] * C.yyzz() * A[1][1] +
A[0][1] * C.xxzz() * A[2][1] + A[0][1] * C.xxzz() * A[2][1] +
A[1][1] * C.yyzz() * A[2][1] + A[1][1] * C.yyzz() * A[2][1] +
A[2][1] * C.zzzz() * A[2][1] + A[2][1] * C.zzzz() * A[2][1] +
A[3][1] * C.xyxy() * A[3][1] + A[3][1] * C.xyxy() * A[3][1] +
A[4][1] * C.yzyz() * A[4][1] + A[4][1] * C.yzyz() * A[4][1] +
A[5][1] * C.zxzx() * A[5][1], A[5][1] * C.zxzx() * A[5][1],
// yyzz // yyzz
A[0][1] * C.xxxx() * A[0][2] + A[0][1] * C.xxxx() * A[0][2] +
A[1][1] * C.xxyy() * A[0][2] + A[1][1] * C.xxyy() * A[0][2] +
A[2][1] * C.xxzz() * A[0][2] + A[2][1] * C.xxzz() * A[0][2] +
A[0][1] * C.xxyy() * A[1][2] + A[0][1] * C.xxyy() * A[1][2] +
A[1][1] * C.yyyy() * A[1][2] + A[1][1] * C.yyyy() * A[1][2] +
A[2][1] * C.yyzz() * A[1][2] + A[2][1] * C.yyzz() * A[1][2] +
A[0][1] * C.xxzz() * A[2][2] + A[0][1] * C.xxzz() * A[2][2] +
A[1][1] * C.yyzz() * A[2][2] + A[1][1] * C.yyzz() * A[2][2] +
A[2][1] * C.zzzz() * A[2][2] + A[2][1] * C.zzzz() * A[2][2] +
A[3][1] * C.xyxy() * A[3][2] + A[3][1] * C.xyxy() * A[3][2] +
A[4][1] * C.yzyz() * A[4][2] + A[4][1] * C.yzyz() * A[4][2] +
A[5][1] * C.zxzx() * A[5][2], A[5][1] * C.zxzx() * A[5][2],
// zzzz // zzzz
A[0][2] * C.xxxx() * A[0][2] + A[0][2] * C.xxxx() * A[0][2] +
A[1][2] * C.xxyy() * A[0][2] + A[1][2] * C.xxyy() * A[0][2] +
A[2][2] * C.xxzz() * A[0][2] + A[2][2] * C.xxzz() * A[0][2] +
A[0][2] * C.xxyy() * A[1][2] + A[0][2] * C.xxyy() * A[1][2] +
A[1][2] * C.yyyy() * A[1][2] + A[1][2] * C.yyyy() * A[1][2] +
A[2][2] * C.yyzz() * A[1][2] + A[2][2] * C.yyzz() * A[1][2] +
A[0][2] * C.xxzz() * A[2][2] + A[0][2] * C.xxzz() * A[2][2] +
A[1][2] * C.yyzz() * A[2][2] + A[1][2] * C.yyzz() * A[2][2] +
A[2][2] * C.zzzz() * A[2][2] + A[2][2] * C.zzzz() * A[2][2] +
A[3][2] * C.xyxy() * A[3][2] + A[3][2] * C.xyxy() * A[3][2] +
A[4][2] * C.yzyz() * A[4][2] + A[4][2] * C.yzyz() * A[4][2] +
A[5][2] * C.zxzx() * A[5][2], A[5][2] * C.zxzx() * A[5][2],
// xyxy // xyxy
A[0][3] * C.xxxx() * A[0][3] + A[0][3] * C.xxxx() * A[0][3] +
A[1][3] * C.xxyy() * A[0][3] + A[1][3] * C.xxyy() * A[0][3] +
A[2][3] * C.xxzz() * A[0][3] + A[2][3] * C.xxzz() * A[0][3] +
A[0][3] * C.xxyy() * A[1][3] + A[0][3] * C.xxyy() * A[1][3] +
A[1][3] * C.yyyy() * A[1][3] + A[1][3] * C.yyyy() * A[1][3] +
A[2][3] * C.yyzz() * A[1][3] + A[2][3] * C.yyzz() * A[1][3] +
A[0][3] * C.xxzz() * A[2][3] + A[0][3] * C.xxzz() * A[2][3] +
A[1][3] * C.yyzz() * A[2][3] + A[1][3] * C.yyzz() * A[2][3] +
A[2][3] * C.zzzz() * A[2][3] + A[2][3] * C.zzzz() * A[2][3] +
A[3][3] * C.xyxy() * A[3][3] + A[3][3] * C.xyxy() * A[3][3] +
A[4][3] * C.yzyz() * A[4][3] + A[4][3] * C.yzyz() * A[4][3] +
A[5][3] * C.zxzx() * A[5][3], A[5][3] * C.zxzx() * A[5][3],
// yzyz // yzyz
A[0][4] * C.xxxx() * A[0][4] + A[0][4] * C.xxxx() * A[0][4] +
A[1][4] * C.xxyy() * A[0][4] + A[1][4] * C.xxyy() * A[0][4] +
A[2][4] * C.xxzz() * A[0][4] + A[2][4] * C.xxzz() * A[0][4] +
A[0][4] * C.xxyy() * A[1][4] + A[0][4] * C.xxyy() * A[1][4] +
A[1][4] * C.yyyy() * A[1][4] + A[1][4] * C.yyyy() * A[1][4] +
A[2][4] * C.yyzz() * A[1][4] + A[2][4] * C.yyzz() * A[1][4] +
A[0][4] * C.xxzz() * A[2][4] + A[0][4] * C.xxzz() * A[2][4] +
A[1][4] * C.yyzz() * A[2][4] + A[1][4] * C.yyzz() * A[2][4] +
A[2][4] * C.zzzz() * A[2][4] + A[2][4] * C.zzzz() * A[2][4] +
A[3][4] * C.xyxy() * A[3][4] + A[3][4] * C.xyxy() * A[3][4] +
A[4][4] * C.yzyz() * A[4][4] + A[4][4] * C.yzyz() * A[4][4] +
A[5][4] * C.zxzx() * A[5][4], A[5][4] * C.zxzx() * A[5][4],
// zxzx // zxzx
A[0][5] * C.xxxx() * A[0][5] + A[0][5] * C.xxxx() * A[0][5] +
A[1][5] * C.xxyy() * A[0][5] + A[1][5] * C.xxyy() * A[0][5] +
A[2][5] * C.xxzz() * A[0][5] + A[2][5] * C.xxzz() * A[0][5] +
A[0][5] * C.xxyy() * A[1][5] + A[0][5] * C.xxyy() * A[1][5] +
A[1][5] * C.yyyy() * A[1][5] + A[1][5] * C.yyyy() * A[1][5] +
A[2][5] * C.yyzz() * A[1][5] + A[2][5] * C.yyzz() * A[1][5] +
A[0][5] * C.xxzz() * A[2][5] + A[0][5] * C.xxzz() * A[2][5] +
A[1][5] * C.yyzz() * A[2][5] + A[1][5] * C.yyzz() * A[2][5] +
A[2][5] * C.zzzz() * A[2][5] + A[2][5] * C.zzzz() * A[2][5] +
A[3][5] * C.xyxy() * A[3][5] + A[3][5] * C.xyxy() * A[3][5] +
A[4][5] * C.yzyz() * A[4][5] + A[4][5] * C.yzyz() * A[4][5] +
A[5][5] * C.zxzx() * A[5][5] A[5][5] * C.zxzx() * A[5][5]
); );
} }

42
src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
View file

@ -146,6 +146,8 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
List<int> finalAgglom(nFineCells); List<int> finalAgglom(nFineCells);
int nMoves = -1; int nMoves = -1;
# ifdef WM_DP
MGridGen MGridGen
( (
nFineCells, nFineCells,
@ -162,6 +164,46 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
finalAgglom.begin() finalAgglom.begin()
); );
# else
// Conversion of type for MGridGen interface
Field<realtype> dblVols(V.size());
forAll (dblVols, i)
{
dblVols[i] = V[i];
}
Field<realtype> dblAreas(Sb.size());
forAll (dblAreas, i)
{
dblAreas[i] = Sb[i];
}
Field<realtype> dblFaceWeights(faceWeights.size());
forAll (dblFaceWeights, i)
{
dblFaceWeights[i] = faceWeights[i];
}
MGridGen
(
nFineCells,
cellCellOffsets.begin(),
dblVols.begin(),
dblAreas.begin(),
cellCells.begin(),
dblFaceWeights.begin(),
minSize,
maxSize,
options.begin(),
&nMoves,
&nCoarseCells,
finalAgglom.begin()
);
# endif
return tmp<labelField>(new labelField(finalAgglom)); return tmp<labelField>(new labelField(finalAgglom));
} }

2
src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C
View file

@ -729,7 +729,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs()
unloadingDeltaEff_, unloadingDeltaEff_,
Foam::sqrt Foam::sqrt
( (
max(deltaN_, 0.0)*max(deltaN_, 0.0) + deltaS_*deltaS_ sqr(max(deltaN_, scalar(0))) + sqr(deltaS_)
) )
); );

19
src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
View file

@ -261,8 +261,8 @@ void dirichletNeumann::correct
slaveDispMag = slaveDispMag =
localSlaveInterpolator.pointToFaceInterpolate<scalar> localSlaveInterpolator.pointToFaceInterpolate<scalar>
( (
min(slavePointPenetration, 0.0) min(slavePointPenetration, scalar(0))
); );
// for visualisation // for visualisation
slaveContactPointGap() = slavePointPenetration; slaveContactPointGap() = slavePointPenetration;
@ -277,11 +277,14 @@ void dirichletNeumann::correct
= mesh.boundaryMesh()[slavePatchIndex].start(); = mesh.boundaryMesh()[slavePatchIndex].start();
forAll(slaveDispMag, facei) forAll(slaveDispMag, facei)
{ {
slaveDispMag[facei] = slaveDispMag[facei] =
globalSlavePenetration globalSlavePenetration
[ [
mesh.faceZones()[slaveFaceZoneID()].whichFace(slavePatchStart + facei) mesh.faceZones()[slaveFaceZoneID()].whichFace
]; (
slavePatchStart + facei
)
];
//- when the master surface surrounds the slave (like the pelvis //- when the master surface surrounds the slave (like the pelvis
// and femur head) then // and femur head) then
@ -317,7 +320,7 @@ void dirichletNeumann::correct
slaveContactPointGap() = slavePointPenetration; slaveContactPointGap() = slavePointPenetration;
// set all positive penetrations to zero // set all positive penetrations to zero
slaveDispMag = min(slaveDispMag,0.0); slaveDispMag = min(slaveDispMag, scalar(0));
minSlavePointPenetration = min(slavePointPenetration); minSlavePointPenetration = min(slavePointPenetration);
reduce(minSlavePointPenetration, minOp<scalar>()); reduce(minSlavePointPenetration, minOp<scalar>());