Commit 13bddac8 authored by mattijs's avatar mattijs
Browse files

singleCellFvMesh and application

parent 5a2eff8e
singleCellMesh.C
EXE = $(FOAM_APPBIN)/singleCellMesh
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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 2 of the License, or (at your
option) any later version.
OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
singleCellMesh
Description
Removes all but one cells of the mesh. Used to generate mesh and fields
that can be used for boundary-only data.
Might easily result in illegal mesh though so only look at boundaries
in paraview.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "Time.H"
#include "ReadFields.H"
#include "singleCellFvMesh.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class GeoField>
void interpolateFields
(
const singleCellFvMesh& scMesh,
const PtrList<GeoField>& flds
)
{
forAll(flds, i)
{
tmp<GeoField> scFld = scMesh.interpolate(flds[i]);
GeoField* scFldPtr = scFld.ptr();
scFldPtr->writeOpt() = IOobject::AUTO_WRITE;
scFldPtr->store();
}
}
// Main program:
int main(int argc, char *argv[])
{
Foam::argList::validOptions.insert("overwrite", "");
# include "addTimeOptions.H"
# include "setRootCase.H"
# include "createTime.H"
// Get times list
instantList Times = runTime.times();
# include "checkTimeOptions.H"
runTime.setTime(Times[startTime], startTime);
# include "createMesh.H"
const word oldInstance = mesh.pointsInstance();
bool overwrite = args.optionFound("overwrite");
// Read objects in time directory
IOobjectList objects(mesh, runTime.timeName());
// Read vol fields.
PtrList<volScalarField> vsFlds;
ReadFields(mesh, objects, vsFlds);
PtrList<volVectorField> vvFlds;
ReadFields(mesh, objects, vvFlds);
PtrList<volSphericalTensorField> vstFlds;
ReadFields(mesh, objects, vstFlds);
PtrList<volSymmTensorField> vsymtFlds;
ReadFields(mesh, objects, vsymtFlds);
PtrList<volTensorField> vtFlds;
ReadFields(mesh, objects, vtFlds);
if (!overwrite)
{
runTime++;
}
// Create the mesh
singleCellFvMesh scMesh
(
IOobject
(
mesh.name(),
mesh.polyMesh::instance(),
runTime,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Map and store the fields on the scMesh.
interpolateFields(scMesh, vsFlds);
interpolateFields(scMesh, vvFlds);
interpolateFields(scMesh, vstFlds);
interpolateFields(scMesh, vsymtFlds);
interpolateFields(scMesh, vtFlds);
// Write
Info<< "Writing mesh to time " << runTime.timeName() << endl;
scMesh.write();
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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 2 of the License, or (at your
option) any later version.
OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Typedef
Foam::uindirectPrimitivePatch
Description
Foam::uindirectPrimitivePatch
\*---------------------------------------------------------------------------*/
#ifndef uindirectPrimitivePatch_H
#define uindirectPrimitivePatch_H
#include "PrimitivePatch.H"
#include "face.H"
#include "UIndirectList.H"
#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef PrimitivePatch<face, UIndirectList, const pointField&>
uindirectPrimitivePatch;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
fvMesh/fvMeshGeometry.C
fvMesh/fvMesh.C
fvMesh/singleCellFvMesh/singleCellFvMesh.C
fvMesh/fvMeshSubset/fvMeshSubset.C
fvBoundaryMesh = fvMesh/fvBoundaryMesh
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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 2 of the License, or (at your
option) any later version.
OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "singleCellFvMesh.H"
#include "syncTools.H"
#include "uindirectPrimitivePatch.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Conversion is a two step process:
// - from original (fine) patch faces to agglomerations (aggloms might not
// be in correct patch order)
// - from agglomerations to coarse patch faces
void Foam::singleCellFvMesh::agglomerateMesh
(
const fvMesh& mesh,
const labelListList& agglom
)
{
const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
// Check agglomeration within patch face range and continuous
labelList nAgglom(oldPatches.size());
forAll(oldPatches, patchI)
{
const polyPatch& pp = oldPatches[patchI];
nAgglom[patchI] = max(agglom[patchI])+1;
forAll(pp, i)
{
if (agglom[patchI][i] < 0 || agglom[patchI][i] >= pp.size())
{
FatalErrorIn
(
"singleCellFvMesh::agglomerateMesh(..)"
) << "agglomeration on patch " << patchI
<< " is out of range 0.." << pp.size()-1
<< exit(FatalError);
}
}
}
// Check agglomeration is sync
{
// Get neighbouring agglomeration
labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces());
forAll(oldPatches, patchI)
{
const polyPatch& pp = oldPatches[patchI];
if (pp.coupled())
{
label offset = pp.start()-mesh.nInternalFaces();
forAll(pp, i)
{
nbrAgglom[offset+i] = agglom[patchI][i];
}
}
}
syncTools::swapBoundaryFaceList(mesh, nbrAgglom, false);
// Get correspondence between this agglomeration and remote one
Map<label> localToNbr(nbrAgglom.size()/10);
forAll(oldPatches, patchI)
{
const polyPatch& pp = oldPatches[patchI];
if (pp.coupled())
{
label offset = pp.start()-mesh.nInternalFaces();
forAll(pp, i)
{
label bFaceI = offset+i;
label myZone = agglom[patchI][i];
label nbrZone = nbrAgglom[bFaceI];
Map<label>::const_iterator iter = localToNbr.find(myZone);
if (iter == localToNbr.end())
{
// First occurence of this zone. Store correspondence
// to remote zone number.
localToNbr.insert(myZone, nbrZone);
}
else
{
// Check that zone numbers are still the same.
if (iter() != nbrZone)
{
FatalErrorIn
(
"singleCellFvMesh::agglomerateMesh(..)"
) << "agglomeration is not synchronised across"
<< " coupled patch " << pp.name()
<< endl
<< "Local agglomeration " << myZone
<< ". Remote agglomeration " << nbrZone
<< exit(FatalError);
}
}
}
}
}
}
label coarseI = 0;
forAll(nAgglom, patchI)
{
coarseI += nAgglom[patchI];
}
// New faces
faceList patchFaces(coarseI);
// New patch start and size
labelList patchStarts(oldPatches.size());
labelList patchSizes(oldPatches.size());
// From new patch face back to agglomeration
patchFaceMap_.setSize(oldPatches.size());
// Face counter
coarseI = 0;
forAll(oldPatches, patchI)
{
const polyPatch& pp = oldPatches[patchI];
if (pp.size() > 0)
{
patchFaceMap_[patchI].setSize(nAgglom[patchI]);
// Patchfaces per agglomeration
labelListList agglomToPatch
(
invertOneToMany(nAgglom[patchI], agglom[patchI])
);
// From agglomeration to compact patch face
labelList agglomToFace(nAgglom[patchI], -1);
patchStarts[patchI] = coarseI;
forAll(pp, i)
{
label myAgglom = agglom[patchI][i];
if (agglomToFace[myAgglom] == -1)
{
// Agglomeration not yet done. We now have:
// - coarseI : current coarse mesh face
// - patchStarts[patchI] : coarse mesh patch start
// - myAgglom : agglomeration
// - agglomToPatch[myAgglom] : fine mesh faces for zone
label coarsePatchFaceI = coarseI - patchStarts[patchI];
patchFaceMap_[patchI][coarsePatchFaceI] = myAgglom;
agglomToFace[myAgglom] = coarsePatchFaceI;
const labelList& fineFaces = agglomToPatch[myAgglom];
// Create overall map from fine mesh faces to coarseI.
forAll(fineFaces, fineI)
{
reverseFaceMap_[pp.start()+fineFaces[fineI]] = coarseI;
}
// Construct single face
uindirectPrimitivePatch upp
(
UIndirectList<face>(pp, fineFaces),
pp.points()
);
if (upp.edgeLoops().size() != 1)
{
FatalErrorIn
(
"singleCellFvMesh::agglomerateMesh(..)"
) << "agglomeration does not create a"
<< " single, non-manifold"
<< " face for agglomeration " << coarseI
<< exit(FatalError);
}
patchFaces[coarseI++] = face
(
renumber
(
upp.meshPoints(),
upp.edgeLoops()[0]
)
);
}
}
patchSizes[patchI] = coarseI-patchStarts[patchI];
}
}
//Pout<< "patchStarts:" << patchStarts << endl;
//Pout<< "patchSizes:" << patchSizes << endl;
// Compact numbering for points
reversePointMap_.setSize(mesh.nPoints());
reversePointMap_.labelList::operator=(-1);
label newI = 0;
forAll(patchFaces, coarseI)
{
face& f = patchFaces[coarseI];
forAll(f, fp)
{
if (reversePointMap_[f[fp]] == -1)
{
reversePointMap_[f[fp]] = newI++;
}
f[fp] = reversePointMap_[f[fp]];
}
}
pointMap_ = invert(newI, reversePointMap_);
// Subset used points
pointField boundaryPoints(mesh.points(), pointMap_);
// Add patches (on still zero sized mesh)
List<polyPatch*> newPatches(oldPatches.size());
forAll(oldPatches, patchI)
{
newPatches[patchI] = oldPatches[patchI].clone
(
boundaryMesh(),
patchI,
0,
0
).ptr();
}
addFvPatches(newPatches);
// Owner, neighbour is trivial
labelList owner(patchFaces.size(), 0);
labelList neighbour(0);
// actually change the mesh
resetPrimitives
(
xferMove(boundaryPoints),
xferMove(patchFaces),
xferMove(owner),
xferMove(neighbour),
patchSizes,
patchStarts,
true //syncPar
);
// Adapt the zones
cellZones().clear();
cellZones().setSize(mesh.cellZones().size());
{
forAll(mesh.cellZones(), zoneI)
{
const cellZone& oldFz = mesh.cellZones()[zoneI];
DynamicList<label> newAddressing;
//Note: uncomment if you think it makes sense. Note that value
// of cell0 is the average.
//// Was old cell0 in this cellZone?
//if (oldFz.localID(0) != -1)
//{
// newAddressing.append(0);
//}
cellZones().set
(
zoneI,
oldFz.clone
(
newAddressing,
zoneI,
cellZones()
)
);
}
}
faceZones().clear();
faceZones().setSize(mesh.faceZones().size());
{
forAll(mesh.faceZones(), zoneI)
{
const faceZone& oldFz = mesh.faceZones()[zoneI];
DynamicList<label> newAddressing(oldFz.size());
DynamicList<bool> newFlipMap(oldFz.size());
forAll(oldFz, i)
{
label newFaceI = reverseFaceMap_[oldFz[i]];
if (newFaceI != -1)
{
newAddressing.append(newFaceI);
newFlipMap.append(oldFz.flipMap()[i]);
}
}
faceZones().set
(
zoneI,
oldFz.clone
(
newAddressing,
newFlipMap,
zoneI,
faceZones()
)
);
}
}
pointZones().clear();
pointZones().setSize(mesh.pointZones().size());
{
forAll(mesh.pointZones(), zoneI)
{
const pointZone& oldFz = mesh.pointZones()[zoneI];
DynamicList<label> newAddressing(oldFz.size());
forAll(oldFz, i)
{
label newPointI = reversePointMap_[oldFz[i]];
if (newPointI != -1)
{
newAddressing.append(newPointI);
}
}
pointZones().set
(
zoneI,
oldFz.clone
(
pointZones(),
zoneI,
newAddressing
)
);
}
}
}