From f644d9d27719b7db3ed2640c1483fcdedc847707 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Thu, 26 Feb 2009 09:02:13 +0000
Subject: [PATCH] handling empty
---
.../sampledSurface/isoSurface/isoSurface.C | 141 ++++++++----------
.../isoSurface/isoSurfaceTemplates.C | 50 ++++++-
2 files changed, 107 insertions(+), 84 deletions(-)
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 939d77cf1b8..78081050ab4 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -199,81 +199,47 @@ void Foam::isoSurface::calcCutTypes
}
}
}
+
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
label faceI = pp.start();
- if (isA<emptyPolyPatch>(pp))
+ forAll(pp, i)
{
- // Still needs special treatment?
-
- forAll(pp, i)
- {
- bool ownLower = (cVals[own[faceI]] < iso_);
+ bool ownLower = (cVals[own[faceI]] < iso_);
- scalar nbrValue;
- point nbrPoint;
- getNeighbour
- (
- boundaryRegion,
- cVals,
- own[faceI],
- faceI,
- nbrValue,
- nbrPoint
- );
+ scalar nbrValue;
+ point nbrPoint;
+ getNeighbour
+ (
+ boundaryRegion,
+ cVals,
+ own[faceI],
+ faceI,
+ nbrValue,
+ nbrPoint
+ );
- bool neiLower = (nbrValue < iso_);
+ bool neiLower = (nbrValue < iso_);
+ if (ownLower != neiLower)
+ {
+ faceCutType_[faceI] = CUT;
+ }
+ else
+ {
+ // Mesh edge.
const face f = mesh_.faces()[faceI];
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
faceCutType_[faceI] = CUT;
}
-
- faceI++;
}
- }
- else
- {
- forAll(pp, i)
- {
- bool ownLower = (cVals[own[faceI]] < iso_);
- scalar nbrValue;
- point nbrPoint;
- getNeighbour
- (
- boundaryRegion,
- cVals,
- own[faceI],
- faceI,
- nbrValue,
- nbrPoint
- );
-
- bool neiLower = (nbrValue < iso_);
-
- if (ownLower != neiLower)
- {
- faceCutType_[faceI] = CUT;
- }
- else
- {
- // Mesh edge.
- const face f = mesh_.faces()[faceI];
-
- if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
- {
- faceCutType_[faceI] = CUT;
- }
- }
-
- faceI++;
- }
+ faceI++;
}
}
@@ -1589,26 +1555,6 @@ Foam::isoSurface::isoSurface
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
- // Check
- forAll(patches, patchI)
- {
- if (isA<emptyPolyPatch>(patches[patchI]))
- {
- FatalErrorIn
- (
- "isoSurface::isoSurface\n"
- "(\n"
- " const volScalarField& cVals,\n"
- " const scalarField& pVals,\n"
- " const scalar iso,\n"
- " const bool regularise,\n"
- " const scalar mergeTol\n"
- ")\n"
- ) << "Iso surfaces not supported on case with empty patches."
- << exit(FatalError);
- }
- }
-
// Pre-calculate patch-per-face to avoid whichPatch call.
labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
@@ -1724,7 +1670,7 @@ Foam::isoSurface::isoSurface
// Generate field to interpolate. This is identical to the mesh.C()
- // except on separated coupled patches.
+ // except on separated coupled patches and on empty patches.
slicedVolVectorField meshC
(
IOobject
@@ -1746,10 +1692,12 @@ Foam::isoSurface::isoSurface
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
+ const polyPatch& pp = patches[patchI];
+
if
(
- patches[patchI].coupled()
- && refCast<const coupledPolyPatch>(patches[patchI]).separated()
+ pp.coupled()
+ && refCast<const coupledPolyPatch>(pp).separated()
)
{
fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
@@ -1758,9 +1706,32 @@ Foam::isoSurface::isoSurface
);
pfld.operator==
(
- patches[patchI].patchSlice(mesh_.faceCentres())
+ pp.patchSlice(mesh_.faceCentres())
);
}
+ else if (isA<emptyPolyPatch>(pp))
+ {
+ typedef slicedVolVectorField::GeometricBoundaryField bType;
+
+ bType& bfld = const_cast<bType&>(meshC.boundaryField());
+
+ // Clear old value. Cannot resize it since slice.
+ bfld.set(patchI, NULL);
+
+ // Set new value we can change
+ bfld.set
+ (
+ patchI,
+ new calculatedFvPatchField<vector>
+ (
+ mesh_.boundary()[patchI],
+ meshC
+ )
+ );
+
+ // Change to face centres
+ bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
+ }
}
}
@@ -1885,6 +1856,14 @@ Foam::isoSurface::isoSurface
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
//}
+
+
+ if (debug)
+ {
+ fileName stlFile = mesh_.time().path() + ".stl";
+ Pout<< "Dumping surface to " << stlFile << endl;
+ triSurface::write(stlFile);
+ }
}
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 7e7206edcc7..0760b06a525 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -514,9 +514,53 @@ void Foam::isoSurface::generateTriPoints
}
else if (isA<emptyPolyPatch>(pp))
{
- // Assume zero-gradient. But what about coordinates?
+ // Check if field is there (when generating geometry the
+ // empty patches have been rewritten to be the face centres),
+ // otherwise use zero-gradient.
+
label faceI = pp.start();
+
+ const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
+
+ // Owner value of cVals
+ scalarField internalVals;
+ if (fvp.size() == 0)
+ {
+ internalVals.setSize(pp.size());
+ forAll(pp, i)
+ {
+ internalVals[i] = cVals[own[pp.start()+i]];
+ }
+ }
+ const scalarField& bVals =
+ (
+ fvp.size() > 0
+ ? static_cast<const scalarField&>(fvp)
+ : internalVals
+ );
+
+
+ const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
+
+ // Owner value of cCoords
+ Field<Type> internalCoords;
+ if (pc.size() == 0)
+ {
+ internalCoords.setSize(pp.size());
+ forAll(pp, i)
+ {
+ internalCoords[i] = cCoords[own[pp.start()+i]];
+ }
+ }
+ const Field<Type>& bCoords =
+ (
+ pc.size() > 0
+ ? static_cast<const Field<Type>&>(pc)
+ : internalCoords
+ );
+
+
forAll(pp, i)
{
if (faceCutType_[faceI] != NOTCUT)
@@ -534,8 +578,8 @@ void Foam::isoSurface::generateTriPoints
snappedPoint,
faceI,
- cVals[own[faceI]],
- cCoords.boundaryField()[patchI][i],
+ bVals[i],
+ bCoords[i],
false, // fc not snapped
pTraits<Type>::zero,
--
GitLab