Merge commit 'OpenCFD/master' into olesenm

This commit is contained in:
Mark Olesen
2009-01-27 23:13:37 +01:00
7 changed files with 554 additions and 166 deletions

View File

@ -30,12 +30,17 @@ Description
- any face inbetween differing cellZones (-cellZones)
Output is:
- mesh with multiple regions
- mesh with multiple regions or
- mesh with cells put into cellZones (-makeCellZones)
Should work in parallel but cellZone interfaces cannot align with
Note:
- Should work in parallel but cellZone interfaces cannot align with
processor boundaries so use the correct option in decomposition to
preserve those interfaces.
- If a cell zone gets split into more than one region it can detect
the largest matching region (-sloppyCellZones). This will accept any
region that covers more than 50% of the zone. It has to be a subset
so cannot have any cells in any other zone.
\*---------------------------------------------------------------------------*/
@ -1036,78 +1041,181 @@ EdgeMap<label> addRegionPatches
}
// Checks if regionI in cellRegion corresponds to existing zone.
label findCorrespondingZone
//// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
//// if no zone found, zone otherwise
//label findCorrespondingSubZone
//(
// const cellZoneMesh& cellZones,
// const labelList& existingZoneID,
// const labelList& cellRegion,
// const label regionI
//)
//{
// // Zone corresponding to region. No corresponding zone.
// label zoneI = labelMax;
//
// labelList regionCells = findIndices(cellRegion, regionI);
//
// if (regionCells.empty())
// {
// // My local portion is empty. Maps to any empty cellZone. Mark with
// // special value which can get overwritten by other processors.
// zoneI = -1;
// }
// else
// {
// // Get zone for first element.
// zoneI = existingZoneID[regionCells[0]];
//
// if (zoneI == -1)
// {
// zoneI = labelMax;
// }
// else
// {
// // 1. All regionCells in zoneI?
// forAll(regionCells, i)
// {
// if (existingZoneID[regionCells[i]] != zoneI)
// {
// zoneI = labelMax;
// break;
// }
// }
// }
// }
//
// // Determine same zone over all processors.
// reduce(zoneI, maxOp<label>());
//
// if (zoneI == labelMax)
// {
// // Cells in region that are not in zoneI
// zoneI = -1;
// }
//
// return zoneI;
//}
//XXXXXXXXX
// Find region that covers most of cell zone
label findCorrespondingRegion
(
const cellZoneMesh& cellZones,
const labelList& existingZoneID,
const labelList& cellRegion,
const label regionI
const labelList& existingZoneID, // per cell the (unique) zoneID
const regionSplit& cellRegion,
const label zoneI,
const label minOverlapSize
)
{
// Zone corresponding to region. No corresponding zone.
label zoneI = labelMax;
// Per region the number of cells in zoneI
labelList cellsInZone(cellRegion.nRegions(), 0);
labelList regionCells = findIndices(cellRegion, regionI);
if (regionCells.empty())
forAll(cellRegion, cellI)
{
// My local portion is empty. Maps to any empty cellZone. Mark with
// special value which can get overwritten by other processors.
zoneI = -1;
if (existingZoneID[cellI] == zoneI)
{
cellsInZone[cellRegion[cellI]]++;
}
}
Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
Pstream::listCombineScatter(cellsInZone);
// Pick region with largest overlap of zoneI
label regionI = findMax(cellsInZone);
if (cellsInZone[regionI] < minOverlapSize)
{
// Region covers too little of zone. Not good enough.
regionI = -1;
}
else
{
// Get zone for first element.
zoneI = existingZoneID[regionCells[0]];
if (zoneI == -1)
// Check that region contains no cells that aren't in cellZone.
forAll(cellRegion, cellI)
{
zoneI = labelMax;
}
else
{
// 1. All regionCells in zoneI?
forAll(regionCells, i)
if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
{
if (existingZoneID[regionCells[i]] != zoneI)
{
zoneI = labelMax;
break;
}
}
}
}
// Determine same zone over all processors.
reduce(zoneI, maxOp<label>());
// 2. All of cellZone present?
if (zoneI == labelMax)
{
zoneI = -1;
}
else if (zoneI != -1)
{
const cellZone& cz = cellZones[zoneI];
forAll(cz, i)
{
if (cellRegion[cz[i]] != regionI)
{
zoneI = -1;
// cellI in regionI but not in zoneI
regionI = -1;
break;
}
}
// If one in error, all should be in error. Note that branch gets taken
// on all procs.
reduce(zoneI, minOp<label>());
reduce(regionI, minOp<label>());
}
return zoneI;
return regionI;
}
//XXXXXXXXX
//// Checks if cellZone has corresponding cellRegion.
//label findCorrespondingRegion
//(
// const cellZoneMesh& cellZones,
// const labelList& existingZoneID, // per cell the (unique) zoneID
// const regionSplit& cellRegion,
// const label zoneI
//)
//{
// // Region corresponding to zone. Start off with special value: no
// // corresponding region.
// label regionI = labelMax;
//
// const cellZone& cz = cellZones[zoneI];
//
// if (cz.empty())
// {
// // My local portion is empty. Maps to any empty cellZone. Mark with
// // special value which can get overwritten by other processors.
// regionI = -1;
// }
// else
// {
// regionI = cellRegion[cz[0]];
//
// forAll(cz, i)
// {
// if (cellRegion[cz[i]] != regionI)
// {
// regionI = labelMax;
// break;
// }
// }
// }
//
// // Determine same zone over all processors.
// reduce(regionI, maxOp<label>());
//
//
// // 2. All of region present?
//
// if (regionI == labelMax)
// {
// regionI = -1;
// }
// else if (regionI != -1)
// {
// forAll(cellRegion, cellI)
// {
// if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
// {
// // cellI in regionI but not in zoneI
// regionI = -1;
// break;
// }
// }
// // If one in error, all should be in error. Note that branch gets taken
// // on all procs.
// reduce(regionI, minOp<label>());
// }
//
// return regionI;
//}
// Main program:
@ -1120,6 +1228,8 @@ int main(int argc, char *argv[])
argList::validOptions.insert("largestOnly", "");
argList::validOptions.insert("insidePoint", "point");
argList::validOptions.insert("overwrite", "");
argList::validOptions.insert("detectOnly", "");
argList::validOptions.insert("sloppyCellZones", "");
# include "setRootCase.H"
# include "createTime.H"
@ -1134,10 +1244,13 @@ int main(int argc, char *argv[])
<< blockedFacesName << nl << endl;
}
bool makeCellZones = args.options().found("makeCellZones");
bool largestOnly = args.options().found("largestOnly");
bool insidePoint = args.options().found("insidePoint");
bool useCellZones = args.options().found("cellZones");
bool overwrite = args.options().found("overwrite");
bool detectOnly = args.options().found("detectOnly");
bool sloppyCellZones = args.options().found("sloppyCellZones");
if (insidePoint && largestOnly)
{
@ -1281,6 +1394,31 @@ int main(int argc, char *argv[])
Info<< "Writing region per cell file (for manual decomposition) to "
<< cellToRegion.objectPath() << nl << endl;
}
// Write for postprocessing
{
volScalarField cellToRegion
(
IOobject
(
"cellToRegion",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0)
);
forAll(cellRegion, cellI)
{
cellToRegion[cellI] = cellRegion[cellI];
}
cellToRegion.write();
Info<< "Writing region per cell as volScalarField to "
<< cellToRegion.objectPath() << nl << endl;
}
// Sizes per region
@ -1307,40 +1445,131 @@ int main(int argc, char *argv[])
Info<< endl;
// Sizes per cellzone
// ~~~~~~~~~~~~~~~~~~
labelList zoneSizes(cellZones.size(), 0);
if (useCellZones || makeCellZones || sloppyCellZones)
{
List<wordList> zoneNames(Pstream::nProcs());
zoneNames[Pstream::myProcNo()] = cellZones.names();
Pstream::gatherList(zoneNames);
Pstream::scatterList(zoneNames);
forAll(zoneNames, procI)
{
if (zoneNames[procI] != zoneNames[0])
{
FatalErrorIn(args.executable())
<< "cellZones not synchronised across processors." << endl
<< "Master has cellZones " << zoneNames[0] << endl
<< "Processor " << procI
<< " has cellZones " << zoneNames[procI]
<< exit(FatalError);
}
}
forAll(cellZones, zoneI)
{
zoneSizes[zoneI] = returnReduce
(
cellZones[zoneI].size(),
sumOp<label>()
);
}
}
// Whether region corresponds to a cellzone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Info<< "Region\tZone\tName" << nl
<< "------\t----\t----" << endl;
// Region per zone
labelList regionToZone(cellRegion.nRegions());
labelList regionToZone(cellRegion.nRegions(), -1);
// Name of region
wordList regionNames(cellRegion.nRegions());
// Zone to region
labelList zoneToRegion(cellZones.size(), -1);
if (sloppyCellZones)
{
Info<< "Trying to match regions to existing cell zones;"
<< " region can be subset of cell zone." << nl << endl;
forAll(cellZones, zoneI)
{
label regionI = findCorrespondingRegion
(
zoneID,
cellRegion,
zoneI,
label(0.5*zoneSizes[zoneI]) // minimum overlap
);
if (regionI != -1)
{
Info<< "Sloppily matched region " << regionI
<< " size " << regionSizes[regionI]
<< " to zone " << zoneI << " size " << zoneSizes[zoneI]
<< endl;
zoneToRegion[zoneI] = regionI;
regionToZone[regionI] = zoneI;
regionNames[regionI] = cellZones[zoneI].name();
}
}
}
else
{
Info<< "Trying to match regions to existing cell zones." << nl << endl;
forAll(cellZones, zoneI)
{
label regionI = findCorrespondingRegion
(
zoneID,
cellRegion,
zoneI,
1 // minimum overlap
);
if (regionI != -1)
{
zoneToRegion[zoneI] = regionI;
regionToZone[regionI] = zoneI;
regionNames[regionI] = cellZones[zoneI].name();
}
}
}
// Allocate region names for unmatched regions.
forAll(regionToZone, regionI)
{
regionToZone[regionI] = findCorrespondingZone
(
cellZones,
zoneID,
cellRegion,
regionI
);
if (regionToZone[regionI] != -1)
{
regionNames[regionI] = cellZones[regionToZone[regionI]].name();
}
else
if (regionToZone[regionI] == -1)
{
regionNames[regionI] = "domain" + Foam::name(regionI);
}
}
// Print region to zone
Info<< "Region\tZone\tName" << nl
<< "------\t----\t----" << endl;
forAll(regionToZone, regionI)
{
Info<< regionI << '\t' << regionToZone[regionI] << '\t'
<< regionNames[regionI] << nl;
}
Info<< endl;
//// Print zone to region
//Info<< "Zone\tName\tRegion" << nl
// << "----\t----\t------" << endl;
//forAll(zoneToRegion, zoneI)
//{
// Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
// << zoneToRegion[zoneI] << nl;
//}
//Info<< endl;
// Since we're going to mess with patches make sure all non-processor ones
// are on all processors.
@ -1361,7 +1590,8 @@ int main(int argc, char *argv[])
interfaceSizes
);
Info<< "Region\tRegion\tFaces" << nl
Info<< "Sizes inbetween regions:" << nl << nl
<< "Region\tRegion\tFaces" << nl
<< "------\t------\t-----" << endl;
forAll(interfaces, interI)
@ -1373,7 +1603,10 @@ int main(int argc, char *argv[])
Info<< endl;
if (detectOnly)
{
return 0;
}
// Read objects in time directory
@ -1423,7 +1656,7 @@ int main(int argc, char *argv[])
{
Info<< "Only one region. Doing nothing." << endl;
}
else if (args.options().found("makeCellZones"))
else if (makeCellZones)
{
Info<< "Putting cells into cellZones instead of splitting mesh."
<< endl;

View File

@ -247,6 +247,7 @@ public:
inline const_iterator begin() const;
inline const const_iterator& end() const;
private:
//- iterator returned by end()
@ -254,7 +255,6 @@ private:
//- const_iterator returned by end()
static const_iterator endConstIter_;
};

View File

@ -94,7 +94,7 @@ void Foam::transform
}
else
{
rtf = vector::zero;
rtf = tf;
}
}
}

View File

@ -538,7 +538,7 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
forAll (bm, patchI)
{
if (bm[patchI].start() != nextPatchStart)
if (bm[patchI].start() != nextPatchStart && !boundaryError)
{
boundaryError = true;
@ -547,7 +547,9 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
<< " of type " << bm[patchI].type()
<< ". The patch should start on face no " << nextPatchStart
<< " and the patch specifies " << bm[patchI].start()
<< "." << endl;
<< "." << endl
<< "Possibly consecutive patches have this same problem."
<< " Suppressing future warnings." << endl;
}
nextPatchStart += bm[patchI].size();

View File

@ -413,6 +413,16 @@ private:
labelList& cellToZone
) const;
//- Determines cell zone from cell region information.
bool calcRegionToZone
(
const label surfZoneI,
const label ownRegion,
const label neiRegion,
labelList& regionToCellZone
) const;
//- Finds zone per cell. Uses topological walk with all faces
// marked in namedSurfaceIndex regarded as blocked.
void findCellZoneTopo
@ -423,6 +433,12 @@ private:
labelList& cellToZone
) const;
void makeConsistentFaceIndex
(
const labelList& cellToZone,
labelList& namedSurfaceIndex
) const;
//- Disallow default bitwise copy construct
meshRefinement(const meshRefinement&);

View File

@ -1011,20 +1011,24 @@ void Foam::meshRefinement::findCellZoneGeometric
}
labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
for
(
label faceI = mesh_.nInternalFaces();
faceI < mesh_.nFaces();
faceI++
)
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
label own = mesh_.faceOwner()[faceI];
neiCellZone[faceI-mesh_.nInternalFaces()] = cellToZone[own];
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
@ -1061,6 +1065,64 @@ void Foam::meshRefinement::findCellZoneGeometric
}
bool Foam::meshRefinement::calcRegionToZone
(
const label surfZoneI,
const label ownRegion,
const label neiRegion,
labelList& regionToCellZone
) const
{
bool changed = false;
// Check whether inbetween different regions
if (ownRegion != neiRegion)
{
// Jump. Change one of the sides to my type.
// 1. Interface between my type and unset region.
// Set region to keepRegion
if (regionToCellZone[ownRegion] == -2)
{
if (regionToCellZone[neiRegion] == surfZoneI)
{
// Face between unset and my region. Put unset
// region into keepRegion
regionToCellZone[ownRegion] = -1;
changed = true;
}
else if (regionToCellZone[neiRegion] != -2)
{
// Face between unset and other region.
// Put unset region into my region
regionToCellZone[ownRegion] = surfZoneI;
changed = true;
}
}
else if (regionToCellZone[neiRegion] == -2)
{
if (regionToCellZone[ownRegion] == surfZoneI)
{
// Face between unset and my region. Put unset
// region into keepRegion
regionToCellZone[neiRegion] = -1;
changed = true;
}
else if (regionToCellZone[ownRegion] != -2)
{
// Face between unset and other region.
// Put unset region into my region
regionToCellZone[neiRegion] = surfZoneI;
changed = true;
}
}
}
return changed;
}
// Finds region per cell. Assumes:
// - region containing keepPoint does not go into a cellZone
// - all other regions can be found by crossing faces marked in
@ -1152,59 +1214,75 @@ void Foam::meshRefinement::findCellZoneTopo
{
bool changed = false;
// Internal faces
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
{
// Get cell zone that surface cells are in
label surfZoneI = surfaceToCellZone[surfI];
// Calculate region to zone from cellRegions on either side
// of internal face.
bool changedCell = calcRegionToZone
(
surfaceToCellZone[surfI],
cellRegion[mesh_.faceOwner()[faceI]],
cellRegion[mesh_.faceNeighbour()[faceI]],
regionToCellZone
);
// Check whether inbetween different regions
label ownRegion = cellRegion[mesh_.faceOwner()[faceI]];
label neiRegion = cellRegion[mesh_.faceNeighbour()[faceI]];
changed = changed | changedCell;
}
}
if (ownRegion != neiRegion)
// Boundary faces
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Get coupled neighbour cellRegion
labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
// Jump. Change one of the sides to my type.
label faceI = pp.start()+i;
neiCellRegion[faceI-mesh_.nInternalFaces()] =
cellRegion[mesh_.faceOwner()[faceI]];
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiCellRegion, false);
// 1. Interface between my type and unset region.
// Set region to keepRegion
// Calculate region to zone from cellRegions on either side of coupled
// face.
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (regionToCellZone[ownRegion] == -2)
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
{
if (regionToCellZone[neiRegion] == surfZoneI)
{
// Face between unset and my region. Put unset
// region into keepRegion
regionToCellZone[ownRegion] = -1;
changed = true;
}
else if (regionToCellZone[neiRegion] != -2)
{
// Face between unset and other region.
// Put unset region into my region
regionToCellZone[ownRegion] = surfZoneI;
changed = true;
}
}
else if (regionToCellZone[neiRegion] == -2)
{
if (regionToCellZone[ownRegion] == surfZoneI)
{
// Face between unset and my region. Put unset
// region into keepRegion
regionToCellZone[neiRegion] = -1;
changed = true;
}
else if (regionToCellZone[ownRegion] != -2)
{
// Face between unset and other region.
// Put unset region into my region
regionToCellZone[neiRegion] = surfZoneI;
changed = true;
}
bool changedCell = calcRegionToZone
(
surfaceToCellZone[surfI],
cellRegion[mesh_.faceOwner()[faceI]],
neiCellRegion[faceI-mesh_.nInternalFaces()],
regionToCellZone
);
changed = changed | changedCell;
}
}
}
@ -1258,6 +1336,88 @@ void Foam::meshRefinement::findCellZoneTopo
}
// Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
// faces inbetween same cell zone.
void Foam::meshRefinement::makeConsistentFaceIndex
(
const labelList& cellToZone,
labelList& namedSurfaceIndex
) const
{
const labelList& faceOwner = mesh_.faceOwner();
const labelList& faceNeighbour = mesh_.faceNeighbour();
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = cellToZone[faceNeighbour[faceI]];
if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
{
namedSurfaceIndex[faceI] = -1;
}
else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
{
FatalErrorIn("meshRefinement::zonify()")
<< "Different cell zones on either side of face " << faceI
<< " at " << mesh_.faceCentres()[faceI]
<< " but face not marked with a surface."
<< abort(FatalError);
}
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Get coupled neighbour cellZone
labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
neiCellZone[faceI-mesh_.nInternalFaces()] =
cellToZone[mesh_.faceOwner()[faceI]];
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
// Use coupled cellZone to do check
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
{
namedSurfaceIndex[faceI] = -1;
}
else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
{
FatalErrorIn("meshRefinement::zonify()")
<< "Different cell zones on either side of face "
<< faceI << " at " << mesh_.faceCentres()[faceI]
<< " but face not marked with a surface."
<< abort(FatalError);
}
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Split off unreachable areas of mesh.
@ -2166,36 +2326,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
);
}
//{
// Pout<< "** finding out blocked faces." << endl;
//
// cellSet zonedCellsGeom(mesh_, "zonedCellsGeom", 100);
// forAll(cellToZone, cellI)
// {
// if (cellToZone[cellI] >= 0)
// {
// zonedCellsGeom.insert(cellI);
// }
// }
// Pout<< "Writing zoned cells to " << zonedCellsGeom.objectPath()
// << endl;
// zonedCellsGeom.write();
//
//
// faceSet zonedFaces(mesh_, "zonedFaces", 100);
// forAll(namedSurfaceIndex, faceI)
// {
// label surfI = namedSurfaceIndex[faceI];
//
// if (surfI != -1)
// {
// zonedFaces.insert(faceI);
// }
// }
// Pout<< "Writing zoned faces to " << zonedFaces.objectPath() << endl;
// zonedFaces.write();
//}
// Set using walking
// ~~~~~~~~~~~~~~~~~
@ -2212,6 +2342,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
}
// Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
// Topochange container
polyTopoChange meshMod(mesh_);
@ -2314,6 +2448,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
mesh_.clearOut();
}
// None of the faces has changed, only the zones. Still...
updateMesh(map, labelList());
return map;
}

View File

@ -108,13 +108,13 @@ spectEddyVisc::spectEddyVisc
tmp<volScalarField> spectEddyVisc::k() const
{
volScalarField Eps = 2*nuEff()*magSqr(symm(fvc::grad(U())));
volScalarField eps = 2*nuEff()*magSqr(symm(fvc::grad(U())));
return
cK1_*pow(delta(), 2.0/3.0)*pow(Eps, 2.0/3.0)
*exp(-cK2_*pow(delta(), -4.0/3.0)*nu()/pow(Eps, 1.0/3.0))
- cK3_*pow(Eps*nu(), 1.0/2.0)
*erfc(cK4_*pow(delta(), -2.0/3.0)*pow(Eps, -1.0/6.0));
cK1_*pow(delta()*eps, 2.0/3.0)
*exp(-cK2_*pow(delta(), -4.0/3.0)*nu()/pow(eps, 1.0/3.0))
- cK3_*sqrt(eps*nu())
*erfc(cK4_*pow(delta(), -2.0/3.0)*sqrt(nu())*pow(eps, -1.0/6.0));
}