Merge branch 'master' into feature/cvMesh

Conflicts:
	src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
	src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
	src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
	src/meshTools/indexedOctree/treeDataPrimitivePatch.C
	src/meshTools/indexedOctree/treeDataTriSurface.C
	src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
This commit is contained in:
laurence
2013-05-08 12:20:52 +01:00
459 changed files with 195180 additions and 4530 deletions

View File

@ -493,6 +493,7 @@ int main(int argc, char *argv[])
mesh.nFaces(), // start
patchI, // index
mesh.boundaryMesh(),// polyBoundaryMesh
Pstream::worldComm, // communicator
Pstream::myProcNo(),// myProcNo
nbrProcI // neighbProcNo
)

View File

@ -222,6 +222,11 @@ castellatedMeshControls
// are only on the boundary of corresponding cellZones or also allow
// free-standing zone faces. Not used if there are no faceZones.
allowFreeStandingZoneFaces true;
// Optional: whether all baffles get eroded away. WIP. Used for
// surface simplification.
//allowFreeStandingBaffles false;
}
// Settings for the snapping.

View File

@ -148,7 +148,7 @@ bool Foam::checkWedges
{
Info<< " ***Wedge patch " << pp.name() << " not planar."
<< " Point " << pt << " is not in patch plane by "
<< d << " meter."
<< d << " metre."
<< endl;
}
return true;
@ -854,7 +854,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
if (allGeometry)
{
faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
if (mesh.checkVolRatio(true, 0.05, &faces))
if (mesh.checkVolRatio(true, 0.01, &faces))
{
noFailedChecks++;

View File

@ -62,7 +62,7 @@ int main(int argc, char *argv[])
boundBox bb(points);
Info<< "bounding box: min = " << bb.min()
<< " max = " << bb.max() << " meters."
<< " max = " << bb.max() << " metres."
<< endl;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -54,7 +54,7 @@ Foam::label Foam::mergePolyMesh::patchIndex(const polyPatch& p)
{
if (patchNames_[patchI] == pName)
{
if (patchTypes_[patchI] == pType)
if (word(patchDicts_[patchI]["type"]) == pType)
{
// Found name and types match
return patchI;
@ -68,7 +68,11 @@ Foam::label Foam::mergePolyMesh::patchIndex(const polyPatch& p)
}
// Patch not found. Append to the list
patchTypes_.append(pType);
{
OStringStream os;
p.write(os);
patchDicts_.append(dictionary(IStringStream(os.str())()));
}
if (nameFound)
{
@ -121,20 +125,22 @@ Foam::mergePolyMesh::mergePolyMesh(const IOobject& io)
:
polyMesh(io),
meshMod_(*this),
patchTypes_(2*boundaryMesh().size()),
patchNames_(2*boundaryMesh().size()),
patchDicts_(2*boundaryMesh().size()),
pointZoneNames_(),
faceZoneNames_(),
cellZoneNames_()
{
// Insert the original patches into the list
wordList curPatchTypes = boundaryMesh().types();
wordList curPatchNames = boundaryMesh().names();
forAll(curPatchTypes, patchI)
forAll(boundaryMesh(), patchI)
{
patchTypes_.append(curPatchTypes[patchI]);
patchNames_.append(curPatchNames[patchI]);
patchNames_.append(boundaryMesh()[patchI].name());
OStringStream os;
boundaryMesh()[patchI].write(os);
patchDicts_.append(dictionary(IStringStream(os.str())()));
}
// Insert point, face and cell zones into the list
@ -379,7 +385,7 @@ void Foam::mergePolyMesh::addMesh(const polyMesh& m)
void Foam::mergePolyMesh::merge()
{
Info<< "patch names: " << patchNames_ << nl
<< "patch types: " << patchTypes_ << nl
<< "patch dicts: " << patchDicts_ << nl
<< "point zone names: " << pointZoneNames_ << nl
<< "face zone names: " << faceZoneNames_ << nl
<< "cell zone names: " << cellZoneNames_ << endl;
@ -409,14 +415,16 @@ void Foam::mergePolyMesh::merge()
for (; patchI < patchNames_.size(); patchI++)
{
// Add a patch
dictionary dict(patchDicts_[patchI]);
dict.set("nFaces", 0);
dict.set("startFace", endOfLastPatch);
newPatches[patchI] =
(
polyPatch::New
(
patchTypes_[patchI],
patchNames_[patchI],
0,
endOfLastPatch,
dict,
patchI,
oldPatches
).ptr()

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -58,12 +58,12 @@ class mergePolyMesh
//- Topological change to accumulated all mesh changes
polyTopoChange meshMod_;
//- Patch types
DynamicList<word> patchTypes_;
//- Patch names
DynamicList<word> patchNames_;
//- Patch dictionaries
DynamicList<dictionary> patchDicts_;
//- Point zone names
DynamicList<word> pointZoneNames_;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -448,6 +448,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI]
);
@ -475,6 +476,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI],
referPatch,

View File

@ -1043,7 +1043,6 @@ int main(int argc, char *argv[])
surfaceMeshWriter writer
(
vMesh,
binary,
pp,
fz.name(),

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,14 +30,12 @@ License
Foam::surfaceMeshWriter::surfaceMeshWriter
(
const vtkMesh& vMesh,
const bool binary,
const indirectPrimitivePatch& pp,
const word& name,
const fileName& fName
)
:
vMesh_(vMesh),
binary_(binary),
pp_(pp),
fName_(fName),
@ -78,8 +76,4 @@ Foam::surfaceMeshWriter::surfaceMeshWriter
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,13 +53,11 @@ namespace Foam
class volPointInterpolation;
/*---------------------------------------------------------------------------*\
Class surfaceMeshWriter Declaration
Class surfaceMeshWriter Declaration
\*---------------------------------------------------------------------------*/
class surfaceMeshWriter
{
const vtkMesh& vMesh_;
const bool binary_;
const indirectPrimitivePatch& pp_;
@ -68,9 +66,6 @@ class surfaceMeshWriter
std::ofstream os_;
// label nPoints_;
//
// label nFaces_;
public:
@ -79,7 +74,6 @@ public:
//- Construct from components
surfaceMeshWriter
(
const vtkMesh&,
const bool binary,
const indirectPrimitivePatch& pp,
const word& name,
@ -94,19 +88,6 @@ public:
return os_;
}
// label nPoints() const
// {
// return nPoints_;
// }
//
// label nFaces() const
// {
// return nFaces_;
// }
//
// //- Write cellIDs
// void writePatchIDs();
//- Extract face data
template<class Type>
tmp<Field<Type> > getFaceField

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -179,12 +179,13 @@ void Foam::vtkPV398Foam::convertMeshPatches
const word patchName = getPartName(partId);
labelHashSet patchIds(patches.patchSet(List<wordRe>(1, patchName)));
labelHashSet
patchIds(patches.patchSet(List<wordRe>(1, wordRe(patchName))));
if (debug)
{
Info<< "Creating VTK mesh for patches [" << patchIds <<"] "
<< patchName << endl;
<< patchName << endl;
}
vtkPolyData* vtkmesh = NULL;

View File

@ -97,6 +97,33 @@ int main(int argc, char *argv[])
cci.write();
}
volScalarField V
(
IOobject
(
"V",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("V", mesh.V().dimensions(), 0.0),
calculatedFvPatchField<scalar>::typeName
);
V.dimensionedInternalField() = mesh.V();
forAll(V.boundaryField(), patchI)
{
V.boundaryField()[patchI] =
V.boundaryField()[patchI].patch().magSf();
}
Info<< "Writing cellVolumes and patch faceAreas to " << V.name()
<< " in " << runTime.timeName() << endl;
V.write();
}
Info<< "\nEnd\n" << endl;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,6 +40,7 @@ Description
#include "pairPatchAgglomeration.H"
#include "labelListIOList.H"
#include "syncTools.H"
#include "globalIndex.H"
using namespace Foam;
@ -53,7 +54,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createNamedMesh.H"
const word dictName("faceAgglomerateDict");
const word dictName("viewFactorsDict");
#include "setConstantMeshDictionaryIO.H"
@ -79,16 +80,16 @@ int main(int argc, char *argv[])
boundary.size()
);
forAll(boundary, patchId)
label nCoarseFaces = 0;
forAllConstIter(dictionary, agglomDict, iter)
{
const polyPatch& pp = boundary[patchId];
label patchI = pp.index();
finalAgglom[patchI].setSize(pp.size(), 0);
if (!pp.coupled())
labelList patchIds = boundary.findIndices(iter().keyword());
forAll(patchIds, i)
{
if (agglomDict.found(pp.name()))
label patchI = patchIds[i];
const polyPatch& pp = boundary[patchI];
if (!pp.coupled())
{
Info << "\nAgglomerating patch : " << pp.name() << endl;
pairPatchAgglomeration agglomObject
@ -99,12 +100,11 @@ int main(int argc, char *argv[])
agglomObject.agglomerate();
finalAgglom[patchI] =
agglomObject.restrictTopBottomAddressing();
}
else
{
FatalErrorIn(args.executable())
<< "Patch " << pp.name() << " not found in dictionary: "
<< agglomDict.name() << exit(FatalError);
if (finalAgglom[patchI].size())
{
nCoarseFaces += max(finalAgglom[patchI] + 1);
}
}
}
}
@ -144,6 +144,7 @@ int main(int argc, char *argv[])
if (writeAgglom)
{
globalIndex index(nCoarseFaces);
volScalarField facesAgglomeration
(
IOobject
@ -158,14 +159,26 @@ int main(int argc, char *argv[])
dimensionedScalar("facesAgglomeration", dimless, 0)
);
label coarsePatchIndex = 0;
forAll(boundary, patchId)
{
fvPatchScalarField& bFacesAgglomeration =
facesAgglomeration.boundaryField()[patchId];
forAll(bFacesAgglomeration, j)
const polyPatch& pp = boundary[patchId];
if (pp.size() > 0)
{
bFacesAgglomeration[j] = finalAgglom[patchId][j];
fvPatchScalarField& bFacesAgglomeration =
facesAgglomeration.boundaryField()[patchId];
forAll(bFacesAgglomeration, j)
{
bFacesAgglomeration[j] =
index.toGlobal
(
Pstream::myProcNo(),
finalAgglom[patchId][j] + coarsePatchIndex
);
}
coarsePatchIndex += max(finalAgglom[patchId]) + 1;
}
}

View File

@ -17,6 +17,12 @@ FoamFile
// Write agglomeration as a volScalarField with calculated boundary values
writeFacesAgglomeration true;
//Debug option
debug 0;
//Dump connectivity rays
dumpRays false;
// Per patch (wildcard possible) the coarsening level
bottomAir_to_heater
{

View File

@ -0,0 +1,4 @@
mapLagrangian.C
mapFields.C
EXE = $(FOAM_APPBIN)/mapFieldsNew

View File

@ -0,0 +1,13 @@
EXE_INC = \
-DFULLDEBUG -g -O0 \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lsampling \
-lmeshTools \
-llagrangian \
-lfiniteVolume \
-lgenericPatchFields

View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
InNamespace
Foam
Description
Gets the indices of (source)particles that have been appended to the
target cloud and maps the lagrangian fields accordingly.
\*---------------------------------------------------------------------------*/
#ifndef MapLagrangianFields_H
#define MapLagrangianFields_H
#include "cloud.H"
#include "GeometricField.H"
#include "meshToMeshNew.H"
#include "IOobjectList.H"
#include "CompactIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Gets the indices of (source)particles that have been appended to the
// target cloud and maps the lagrangian fields accordingly.
template<class Type>
void MapLagrangianFields
(
const string& cloudName,
const IOobjectList& objects,
const polyMesh& meshTarget,
const labelList& addParticles
)
{
{
IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian field " << fieldName << endl;
// Read field (does not need mesh)
IOField<Type> fieldSource(*fieldIter());
// Map
IOField<Type> fieldTarget
(
IOobject
(
fieldName,
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(IOField<Field<Type> >::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
const word& fieldName = fieldIter()->name();
Info<< " mapping lagrangian fieldField " << fieldName << endl;
// Read field (does not need mesh)
IOField<Field<Type> > fieldSource(*fieldIter());
// Map - use CompactIOField to automatically write in
// compact form for binary format.
CompactIOField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldName,
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
{
IOobjectList fieldFields =
objects.lookupClass(CompactIOField<Field<Type>, Type>::typeName);
forAllIter(IOobjectList, fieldFields, fieldIter)
{
Info<< " mapping lagrangian fieldField "
<< fieldIter()->name() << endl;
// Read field (does not need mesh)
CompactIOField<Field<Type>, Type> fieldSource(*fieldIter());
// Map
CompactIOField<Field<Type>, Type> fieldTarget
(
IOobject
(
fieldIter()->name(),
meshTarget.time().timeName(),
cloud::prefix/cloudName,
meshTarget,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
addParticles.size()
);
forAll(addParticles, i)
{
fieldTarget[i] = fieldSource[addParticles[i]];
}
// Write field
fieldTarget.write();
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef MapMeshes_H
#define MapMeshes_H
#include "MapVolFields.H"
#include "mapLagrangian.H"
#include "UnMapped.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<template<class> class CombineOp>
void MapMesh
(
const meshToMeshNew& interp,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
{
const polyMesh& meshSource = interp.srcRegion();
// Search for list of objects for this time
IOobjectList objects(meshSource, meshSource.time().timeName());
// Map volFields
// ~~~~~~~~~~~~~
MapVolFields<scalar>
(
objects,
selectedFields,
interp,
CombineOp<scalar>()
);
MapVolFields<vector>
(
objects,
selectedFields,
interp,
CombineOp<vector>()
);
MapVolFields<sphericalTensor>
(
objects,
selectedFields,
interp,
CombineOp<sphericalTensor>()
);
MapVolFields<symmTensor>
(
objects,
selectedFields,
interp,
CombineOp<symmTensor>()
);
MapVolFields<tensor>
(
objects,
selectedFields,
interp,
CombineOp<tensor>()
);
}
{
const polyMesh& meshTarget = interp.tgtRegion();
// Search for list of target objects for this time
IOobjectList objects(meshTarget, meshTarget.time().timeName());
// Mark surfaceFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<surfaceScalarField>(objects);
UnMapped<surfaceVectorField>(objects);
UnMapped<surfaceSphericalTensorField>(objects);
UnMapped<surfaceSymmTensorField>(objects);
UnMapped<surfaceTensorField>(objects);
// Mark pointFields as unmapped
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
UnMapped<pointScalarField>(objects);
UnMapped<pointVectorField>(objects);
UnMapped<pointSphericalTensorField>(objects);
UnMapped<pointSymmTensorField>(objects);
UnMapped<pointTensorField>(objects);
}
if (!noLagrangian)
{
mapLagrangian(interp);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,104 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef MapConsistentVolFields_H
#define MapConsistentVolFields_H
#include "GeometricField.H"
#include "meshToMeshNew.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type, class CombineOp>
void MapVolFields
(
const IOobjectList& objects,
const HashSet<word>& selectedFields,
const meshToMeshNew& interp,
const CombineOp& cop
)
{
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const fvMesh& meshSource = static_cast<const fvMesh&>(interp.srcRegion());
const fvMesh& meshTarget = static_cast<const fvMesh&>(interp.tgtRegion());
IOobjectList fields = objects.lookupClass(fieldType::typeName);
forAllIter(IOobjectList, fields, fieldIter)
{
const word& fieldName = fieldIter()->name();
if (selectedFields.empty() || selectedFields.found(fieldName))
{
Info<< " interpolating " << fieldName << endl;
const fieldType fieldSource(*fieldIter(), meshSource);
IOobject targetIO
(
fieldName,
meshTarget.time().timeName(),
meshTarget,
IOobject::MUST_READ
);
if (targetIO.headerOk())
{
fieldType fieldTarget(targetIO, meshTarget);
interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
fieldTarget.write();
}
else
{
targetIO.readOpt() = IOobject::NO_READ;
tmp<fieldType>
tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
fieldType fieldTarget(targetIO, tfieldTarget);
fieldTarget.write();
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef UnMapped_H
#define UnMapped_H
#include "IOobjectList.H"
#include "OSspecific.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type>
void UnMapped(const IOobjectList& objects)
{
IOobjectList fields = objects.lookupClass(Type::typeName);
forAllConstIter(IOobjectList, fields, fieldIter)
{
mvBak(fieldIter()->objectPath(), "unmapped");
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,11 @@
Info<< "\nCreate databases as time" << endl;
HashTable<string> srcOptions(args.options());
srcOptions.erase("case");
srcOptions.insert("case", fileName(rootDirSource/caseDirSource));
argList argsSrc(args, srcOptions, false, false, false);
Time runTimeSource(Time::controlDictName, argsSrc);
Time runTimeTarget(Time::controlDictName, args);

View File

@ -0,0 +1,343 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Application
mapFields
Description
Maps volume fields from one mesh to another, reading and
interpolating all fields present in the time directory of both cases.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "meshToMeshNew.H"
#include "processorPolyPatch.H"
#include "MapMeshes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void mapConsistentMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const meshToMeshNew::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
Info<< nl << "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMeshNew interp(meshSource, meshTarget, mapMethod);
if (subtract)
{
MapMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
else
{
MapMesh<plusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
}
void mapSubMesh
(
const fvMesh& meshSource,
const fvMesh& meshTarget,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const meshToMeshNew::interpolationMethod& mapMethod,
const bool subtract,
const HashSet<word>& selectedFields,
const bool noLagrangian
)
{
Info<< nl << "Creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMeshNew interp
(
meshSource,
meshTarget,
mapMethod,
patchMap,
cuttingPatches
);
if (subtract)
{
MapMesh<minusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
else
{
MapMesh<plusEqOp>
(
interp,
selectedFields,
noLagrangian
);
}
}
wordList addProcessorPatches
(
const fvMesh& meshTarget,
const wordList& cuttingPatches
)
{
// Add the processor patches to the cutting list
HashSet<word> cuttingPatchTable;
forAll(cuttingPatches, i)
{
cuttingPatchTable.insert(cuttingPatches[i]);
}
const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();
forAll(pbm, patchI)
{
if (isA<processorPolyPatch>(pbm[patchI]))
{
const word& patchName = pbm[patchI].name();
cuttingPatchTable.insert(patchName);
}
}
return cuttingPatchTable.toc();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"map volume fields from one mesh to another"
);
argList::validArgs.append("sourceCase");
argList::addOption
(
"sourceTime",
"scalar|'latestTime'",
"specify the source time"
);
argList::addOption
(
"sourceRegion",
"word",
"specify the source region"
);
argList::addOption
(
"targetRegion",
"word",
"specify the target region"
);
argList::addBoolOption
(
"consistent",
"source and target geometry and boundary conditions identical"
);
argList::addOption
(
"mapMethod",
"word",
"specify the mapping method"
);
argList::addBoolOption
(
"subtract",
"subtract mapped source from target"
);
argList::addOption
(
"fields",
"list",
"specify a list of fields to be mapped. Eg, '(U T p)' - "
"regular expressions not currently supported"
);
argList::addBoolOption
(
"noLagrangian",
"skip mapping lagrangian positions and fields"
);
argList args(argc, argv);
fileName rootDirTarget(args.rootPath());
fileName caseDirTarget(args.globalCaseName());
const fileName casePath = args[1];
const fileName rootDirSource = casePath.path();
const fileName caseDirSource = casePath.name();
Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;
word sourceRegion = fvMesh::defaultRegion;
if (args.optionFound("sourceRegion"))
{
sourceRegion = args["sourceRegion"];
Info<< "Source region: " << sourceRegion << endl;
}
Info<< "Target: " << rootDirTarget << " " << caseDirTarget << endl;
word targetRegion = fvMesh::defaultRegion;
if (args.optionFound("targetRegion"))
{
targetRegion = args["targetRegion"];
Info<< "Target region: " << targetRegion << endl;
}
const bool consistent = args.optionFound("consistent");
meshToMeshNew::interpolationMethod mapMethod =
meshToMeshNew::imCellVolumeWeight;
if (args.optionFound("mapMethod"))
{
mapMethod = meshToMeshNew::interpolationMethodNames_[args["mapMethod"]];
Info<< "Mapping method: "
<< meshToMeshNew::interpolationMethodNames_[mapMethod] << endl;
}
const bool subtract = args.optionFound("subtract");
if (subtract)
{
Info<< "Subtracting mapped source field from target" << endl;
}
HashSet<word> selectedFields;
if (args.optionFound("fields"))
{
args.optionLookup("fields")() >> selectedFields;
}
const bool noLagrangian = args.optionFound("noLagrangian");
#include "createTimes.H"
HashTable<word> patchMap;
wordList cuttingPatches;
if (!consistent)
{
IOdictionary mapFieldsDict
(
IOobject
(
"mapFieldsDict",
runTimeTarget.system(),
runTimeTarget,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
);
mapFieldsDict.lookup("patchMap") >> patchMap;
mapFieldsDict.lookup("cuttingPatches") >> cuttingPatches;
}
#include "setTimeIndex.H"
Info<< "\nCreate meshes\n" << endl;
fvMesh meshSource
(
IOobject
(
sourceRegion,
runTimeSource.timeName(),
runTimeSource
)
);
fvMesh meshTarget
(
IOobject
(
targetRegion,
runTimeTarget.timeName(),
runTimeTarget
)
);
Info<< "Source mesh size: " << meshSource.nCells() << tab
<< "Target mesh size: " << meshTarget.nCells() << nl << endl;
if (consistent)
{
mapConsistentMesh
(
meshSource,
meshTarget,
mapMethod,
subtract,
selectedFields,
noLagrangian
);
}
else
{
mapSubMesh
(
meshSource,
meshTarget,
patchMap,
addProcessorPatches(meshTarget, cuttingPatches),
mapMethod,
subtract,
selectedFields,
noLagrangian
);
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,30 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// List of pairs of source/target patches for mapping
patchMap
(
lid movingWall
);
// List of target patches cutting the source domain (these need to be
// handled specially e.g. interpolated from internal values)
cuttingPatches
(
fixedWalls
);
// ************************************************************************* //

View File

@ -0,0 +1,303 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "MapLagrangianFields.H"
#include "passiveParticleCloud.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
static const scalar perturbFactor = 1e-6;
// Special version of findCell that generates a cell guaranteed to be
// compatible with tracking.
static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
{
label cellI = -1;
label tetFaceI = -1;
label tetPtI = -1;
const polyMesh& mesh = cloud.pMesh();
mesh.findCellFacePt(pt, cellI, tetFaceI, tetPtI);
if (cellI >= 0)
{
return cellI;
}
else
{
// See if particle on face by finding nearest face and shifting
// particle.
meshSearch meshSearcher
(
mesh,
polyMesh::FACEPLANES // no decomposition needed
);
label faceI = meshSearcher.findNearestBoundaryFace(pt);
if (faceI >= 0)
{
const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
mesh.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI);
return cellI;
}
}
return -1;
}
void mapLagrangian(const meshToMeshNew& interp)
{
// Determine which particles are in meshTarget
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const polyMesh& meshSource = interp.srcRegion();
const polyMesh& meshTarget = interp.tgtRegion();
const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
const pointField& targetCc = meshTarget.cellCentres();
fileNameList cloudDirs
(
readDir
(
meshSource.time().timePath()/cloud::prefix,
fileName::DIRECTORY
)
);
forAll(cloudDirs, cloudI)
{
// Search for list of lagrangian objects for this time
IOobjectList objects
(
meshSource,
meshSource.time().timeName(),
cloud::prefix/cloudDirs[cloudI]
);
IOobject* positionsPtr = objects.lookup(word("positions"));
if (positionsPtr)
{
Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
// Read positions & cell
passiveParticleCloud sourceParcels
(
meshSource,
cloudDirs[cloudI],
false
);
Info<< " read " << sourceParcels.size()
<< " parcels from source mesh." << endl;
// Construct empty target cloud
passiveParticleCloud targetParcels
(
meshTarget,
cloudDirs[cloudI],
IDLList<passiveParticle>()
);
particle::TrackingData<passiveParticleCloud> td(targetParcels);
label sourceParticleI = 0;
// Indices of source particles that get added to targetParcels
DynamicList<label> addParticles(sourceParcels.size());
// Unmapped particles
labelHashSet unmappedSource(sourceParcels.size());
// Initial: track from fine-mesh cell centre to particle position
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This requires there to be no boundary in the way.
forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
{
bool foundCell = false;
// Assume that cell from read parcel is the correct one...
if (iter().cell() >= 0)
{
const labelList& targetCells =
sourceToTarget[iter().cell()];
// Particle probably in one of the targetcells. Try
// all by tracking from their cell centre to the parcel
// position.
forAll(targetCells, i)
{
// Track from its cellcentre to position to make sure.
autoPtr<passiveParticle> newPtr
(
new passiveParticle
(
meshTarget,
targetCc[targetCells[i]],
targetCells[i]
)
);
passiveParticle& newP = newPtr();
label faceI = newP.track(iter().position(), td);
if (faceI < 0 && newP.cell() >= 0)
{
// Hit position.
foundCell = true;
addParticles.append(sourceParticleI);
targetParcels.addParticle(newPtr.ptr());
break;
}
}
}
if (!foundCell)
{
// Store for closer analysis
unmappedSource.insert(sourceParticleI);
}
sourceParticleI++;
}
Info<< " after meshToMesh addressing found "
<< targetParcels.size()
<< " parcels in target mesh." << endl;
// Do closer inspection for unmapped particles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (unmappedSource.size())
{
sourceParticleI = 0;
forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
{
if (unmappedSource.found(sourceParticleI))
{
label targetCell =
findCell(targetParcels, iter().position());
if (targetCell >= 0)
{
unmappedSource.erase(sourceParticleI);
addParticles.append(sourceParticleI);
iter().cell() = targetCell;
targetParcels.addParticle
(
sourceParcels.remove(&iter())
);
}
}
sourceParticleI++;
}
}
addParticles.shrink();
Info<< " after additional mesh searching found "
<< targetParcels.size() << " parcels in target mesh." << endl;
if (addParticles.size())
{
IOPosition<passiveParticleCloud>(targetParcels).write();
// addParticles now contains the indices of the sourceMesh
// particles that were appended to the target mesh.
// Map lagrangian fields
// ~~~~~~~~~~~~~~~~~~~~~
MapLagrangianFields<label>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<scalar>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<vector>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<sphericalTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<symmTensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
MapLagrangianFields<tensor>
(
cloudDirs[cloudI],
objects,
meshTarget,
addParticles
);
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
InNamespace
Foam
Description
Maps lagrangian positions and fields
SourceFiles
mapLagrangian.C
\*---------------------------------------------------------------------------*/
#ifndef mapLagrangian_H
#define mapLagrangian_H
#include "meshToMeshNew.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Maps lagrangian positions and fields
void mapLagrangian(const meshToMeshNew& interp);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,33 @@
{
instantList sourceTimes = runTimeSource.times();
label sourceTimeIndex = runTimeSource.timeIndex();
if (args.optionFound("sourceTime"))
{
if (args["sourceTime"] == "latestTime")
{
sourceTimeIndex = sourceTimes.size() - 1;
}
else
{
sourceTimeIndex = Time::findClosestTimeIndex
(
sourceTimes,
args.optionRead<scalar>("sourceTime")
);
}
}
else
{
sourceTimeIndex = Time::findClosestTimeIndex
(
sourceTimes,
runTimeTarget.time().value()
);
}
runTimeSource.setTime(sourceTimes[sourceTimeIndex], sourceTimeIndex);
Info<< "\nSource time: " << runTimeSource.value()
<< "\nTarget time: " << runTimeTarget.value()
<< endl;
}

View File

@ -32,6 +32,19 @@ forAll(patches, patchI)
}
}
labelList triSurfaceToAgglom(5*nFineFaces);
const triSurface localSurface = triangulate
(
patches,
includePatches,
finalAgglom,
triSurfaceToAgglom,
globalNumbering,
coarsePatches
);
distributedTriSurfaceMesh surfacesMesh
(
IOobject
@ -43,12 +56,13 @@ distributedTriSurfaceMesh surfacesMesh
IOobject::NO_READ,
IOobject::NO_WRITE
),
triSurfaceTools::triangulate
(
patches,
includePatches
),
localSurface,
dict
);
triSurfaceToAgglom.resize(surfacesMesh.size());
//surfacesMesh.searchableSurface::write();
surfacesMesh.setField(triSurfaceToAgglom);

View File

@ -2,7 +2,7 @@
// Pre-size by assuming a certain percentage is visible.
// Maximum lenght for dynamicList
const label maxDynListLength = 10000;
const label maxDynListLength = 100000;
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
@ -14,12 +14,17 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
DynamicList<label> startIndex(start.size());
DynamicList<label> endIndex(start.size());
DynamicList<label> startAgg(start.size());
DynamicList<label> endAgg(start.size());
const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
const pointField& remoteArea = remoteCoarseSf[procI];
const pointField& remoteFc = remoteCoarseCf[procI];
const labelField& remoteAgg = remoteCoarseAgg[procI];
label i = 0;
label j = 0;
@ -29,6 +34,7 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const point& fc = myFc[i];
const vector& fA = myArea[i];
const label& fAgg = myAgg[i];
for (; j < remoteFc.size(); j++)//
{
@ -36,26 +42,32 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const point& remFc = remoteFc[j];
const vector& remA = remoteArea[j];
const label& remAgg = remoteAgg[j];
const vector& d = remFc - fc;
if (((d & fA) < 0.) && ((d & remA) > 0))
{
start.append(fc + 0.0001*d);
start.append(fc + 0.001*d);
startIndex.append(i);
end.append(fc + 0.9999*d);
startAgg.append(globalNumbering.toGlobal(procI, fAgg));
end.append(fc + 0.999*d);
label globalI = globalNumbering.toGlobal(procI, j);
endIndex.append(globalI);
endAgg.append(globalNumbering.toGlobal(procI, remAgg));
if (startIndex.size() > maxDynListLength)
{
break;
FatalErrorIn
(
"shootRays"
) << "Dynamic list need from capacity."
<< "Actual size maxDynListLength : "
<< maxDynListLength
<< abort(FatalError);
}
}
}
}
if (startIndex.size() > maxDynListLength)
{
break;
}
if (j == remoteFc.size())
{
@ -63,23 +75,102 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
}
}
List<pointIndexHit> hitInfo(startIndex.size());
surfacesMesh.findLine(start, end, hitInfo);
}while (returnReduce(i < myFc.size(), orOp<bool>()));
forAll (hitInfo, rayI)
List<pointIndexHit> hitInfo(startIndex.size());
surfacesMesh.findLine(start, end, hitInfo);
// Return hit global agglo index
labelList aggHitIndex;
surfacesMesh.getField(hitInfo, aggHitIndex);
DynamicList<label> dRayIs;
// Collect the rays which has not abstacle in bettween in rayStartFace
// and rayEndFace. If the ray hit itself get stored in dRayIs
forAll (hitInfo, rayI)
{
if (!hitInfo[rayI].hit())
{
if (!hitInfo[rayI].hit())
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
}
else if (aggHitIndex[rayI] == startAgg[rayI])
{
dRayIs.append(rayI);
}
}
start.clear();
// Continue rays which hit themself. If they hit the target
// agglomeration are added to rayStartFace and rayEndFace
bool firstLoop = true;
DynamicField<point> startHitItself;
DynamicField<point> endHitItself;
label iter = 0;
do
{
labelField rayIs;
rayIs.transfer(dRayIs);
dRayIs.clear();
forAll (rayIs, rayI)
{
const label rayID = rayIs[rayI];
label hitIndex = -1;
if (firstLoop)
{
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
hitIndex = rayIs[rayI];
}
else
{
hitIndex = rayI;
}
if (hitInfo[hitIndex].hit())
{
if (aggHitIndex[hitIndex] == startAgg[rayID])
{
const vector& endP = end[rayID];
const vector& startP = hitInfo[hitIndex].hitPoint();
const vector& d = endP - startP;
startHitItself.append(startP + 0.01*d);
endHitItself.append(startP + 1.01*d);
dRayIs.append(rayID);
}
else if (aggHitIndex[hitIndex] == endAgg[rayID])
{
rayStartFace.append(startIndex[rayID]);
rayEndFace.append(endIndex[rayID]);
}
}
}
start.clear();
startIndex.clear();
end.clear();
endIndex.clear();
hitInfo.clear();
hitInfo.resize(dRayIs.size());
}while (returnReduce(i < myFc.size(), orOp<bool>()));
surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
surfacesMesh.getField(hitInfo, aggHitIndex);
endHitItself.clear();
startHitItself.clear();
firstLoop = false;
iter ++;
}while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
startIndex.clear();
end.clear();
endIndex.clear();
startAgg.clear();
endAgg.clear();
}

View File

@ -68,6 +68,101 @@ Description
using namespace Foam;
triSurface triangulate
(
const polyBoundaryMesh& bMesh,
const labelHashSet& includePatches,
const labelListIOList& finalAgglom,
labelList& triSurfaceToAgglom,
const globalIndex& globalNumbering,
const polyBoundaryMesh& coarsePatches
)
{
const polyMesh& mesh = bMesh.mesh();
// Storage for surfaceMesh. Size estimate.
DynamicList<labelledTri> triangles
(
mesh.nFaces() - mesh.nInternalFaces()
);
label newPatchI = 0;
label localTriFaceI = 0;
forAllConstIter(labelHashSet, includePatches, iter)
{
const label patchI = iter.key();
const polyPatch& patch = bMesh[patchI];
const pointField& points = patch.points();
label nTriTotal = 0;
forAll(patch, patchFaceI)
{
const face& f = patch[patchFaceI];
faceList triFaces(f.nTriangles(points));
label nTri = 0;
f.triangles(points, nTri, triFaces);
forAll(triFaces, triFaceI)
{
const face& f = triFaces[triFaceI];
triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
nTriTotal++;
triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
(
Pstream::myProcNo(),
finalAgglom[patchI][patchFaceI]
+ coarsePatches[patchI].start()
);
}
}
newPatchI++;
}
triSurfaceToAgglom.resize(localTriFaceI-1);
triangles.shrink();
// Create globally numbered tri surface
triSurface rawSurface(triangles, mesh.points());
// Create locally numbered tri surface
triSurface surface
(
rawSurface.localFaces(),
rawSurface.localPoints()
);
// Add patch names to surface
surface.patches().setSize(newPatchI);
newPatchI = 0;
forAllConstIter(labelHashSet, includePatches, iter)
{
const label patchI = iter.key();
const polyPatch& patch = bMesh[patchI];
surface.patches()[newPatchI].index() = patchI;
surface.patches()[newPatchI].name() = patch.name();
surface.patches()[newPatchI].geometricType() = patch.type();
newPatchI++;
}
return surface;
}
void writeRays
(
const fileName& fName,
@ -213,7 +308,7 @@ int main(int argc, char *argv[])
if (debug)
{
Info << "\nCreating single cell mesh..." << endl;
Pout << "\nCreating single cell mesh..." << endl;
}
singleCellFvMesh coarseMesh
@ -230,6 +325,11 @@ int main(int argc, char *argv[])
finalAgglom
);
if (debug)
{
Pout << "\nCreated single cell mesh..." << endl;
}
// Calculate total number of fine and coarse faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -283,6 +383,8 @@ int main(int argc, char *argv[])
DynamicList<point> localCoarseCf(nCoarseFaces);
DynamicList<point> localCoarseSf(nCoarseFaces);
DynamicList<label> localAgg(nCoarseFaces);
forAll (viewFactorsPatches, i)
{
const label patchID = viewFactorsPatches[i];
@ -296,11 +398,18 @@ int main(int argc, char *argv[])
const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
labelHashSet includePatches;
includePatches.insert(patchID);
forAll(coarseCf, faceI)
{
point cf = coarseCf[faceI];
const label coarseFaceI = coarsePatchFace[faceI];
const labelList& fineFaces = coarseToFine[coarseFaceI];
const label agglomI =
agglom[fineFaces[0]] + coarsePatches[patchID].start();
// Construct single face
uindirectPrimitivePatch upp
(
@ -308,6 +417,7 @@ int main(int argc, char *argv[])
pp.points()
);
List<point> availablePoints
(
upp.faceCentres().size()
@ -342,6 +452,7 @@ int main(int argc, char *argv[])
point sf = coarseSf[faceI];
localCoarseCf.append(cf);
localCoarseSf.append(sf);
localAgg.append(agglomI);
}
}
@ -350,9 +461,12 @@ int main(int argc, char *argv[])
List<pointField> remoteCoarseCf(Pstream::nProcs());
List<pointField> remoteCoarseSf(Pstream::nProcs());
List<labelField> remoteCoarseAgg(Pstream::nProcs());
remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
// Collect remote Cf and Sf on fine mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -370,8 +484,12 @@ int main(int argc, char *argv[])
Pstream::scatterList(remoteCoarseCf);
Pstream::gatherList(remoteCoarseSf);
Pstream::scatterList(remoteCoarseSf);
Pstream::gatherList(remoteCoarseAgg);
Pstream::scatterList(remoteCoarseAgg);
globalIndex globalNumbering(nCoarseFaces);
// Set up searching engine for obstacles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#include "searchingEngine.H"
@ -383,8 +501,6 @@ int main(int argc, char *argv[])
DynamicList<label> rayEndFace(rayStartFace.size());
globalIndex globalNumbering(nCoarseFaces);
// Return rayStartFace in local index andrayEndFace in global index
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

View File

@ -435,7 +435,7 @@ int main(int argc, char *argv[])
scalar smallDim = 1e-6 * bb.mag();
Info<< "Checking for points less than 1e-6 of bounding box ("
<< bb.span() << " meter) apart."
<< bb.span() << " metre) apart."
<< endl;
// Sort points

View File

@ -0,0 +1,3 @@
surfaceHookUp.C
EXE = $(FOAM_APPBIN)/surfaceHookUp

View File

@ -0,0 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools \
-lfileFormats \
-ltriSurface

View File

@ -0,0 +1,525 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Application
surfaceHookUp
Description
Find close open edges and stitches the surface along them
Usage
- surfaceHookUp hookDistance [OPTION]
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "triSurfaceMesh.H"
#include "indexedOctree.H"
#include "treeBoundBox.H"
#include "PackedBoolList.H"
#include "unitConversion.H"
#include "searchableSurfaces.H"
using namespace Foam;
// Split faceI along edgeI at position newPointI
void greenRefine
(
const triSurface& surf,
const label faceI,
const label edgeI,
const label newPointI,
DynamicList<labelledTri>& newFaces
)
{
const labelledTri& f = surf.localFaces()[faceI];
const edge& e = surf.edges()[edgeI];
// Find index of edge in face.
label fp0 = findIndex(f, e[0]);
label fp1 = f.fcIndex(fp0);
label fp2 = f.fcIndex(fp1);
if (f[fp1] == e[1])
{
// Edge oriented like face
newFaces.append
(
labelledTri
(
f[fp0],
newPointI,
f[fp2],
f.region()
)
);
newFaces.append
(
labelledTri
(
newPointI,
f[fp1],
f[fp2],
f.region()
)
);
}
else
{
newFaces.append
(
labelledTri
(
f[fp2],
newPointI,
f[fp1],
f.region()
)
);
newFaces.append
(
labelledTri
(
newPointI,
f[fp0],
f[fp1],
f.region()
)
);
}
}
//scalar checkEdgeAngle
//(
// const triSurface& surf,
// const label edgeIndex,
// const label pointIndex,
// const scalar& angle
//)
//{
// const edge& e = surf.edges()[edgeIndex];
// vector eVec = e.vec(surf.localPoints());
// eVec /= mag(eVec) + SMALL;
// const labelList& pEdges = surf.pointEdges()[pointIndex];
//
// forAll(pEdges, eI)
// {
// const edge& nearE = surf.edges()[pEdges[eI]];
// vector nearEVec = nearE.vec(surf.localPoints());
// nearEVec /= mag(nearEVec) + SMALL;
// const scalar dot = eVec & nearEVec;
// const scalar minCos = degToRad(angle);
// if (mag(dot) > minCos)
// {
// return false;
// }
// }
// return true;
//}
void createBoundaryEdgeTrees
(
const PtrList<triSurfaceMesh>& surfs,
PtrList<indexedOctree<treeDataEdge> >& bEdgeTrees,
labelListList& treeBoundaryEdges
)
{
forAll(surfs, surfI)
{
const triSurface& surf = surfs[surfI];
// Boundary edges
treeBoundaryEdges[surfI] =
labelList
(
identity(surf.nEdges() - surf.nInternalEdges())
+ surf.nInternalEdges()
);
Random rndGen(17301893);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are less face/edge aligned items.
treeBoundBox bb
(
treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bEdgeTrees.set
(
surfI,
new indexedOctree<treeDataEdge>
(
treeDataEdge
(
false, // cachebb
surf.edges(), // edges
surf.localPoints(), // points
treeBoundaryEdges[surfI] // selected edges
),
bb, // bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"hook surfaces to other surfaces by moving and retriangulating their"
"boundary edges to match other surface boundary edges"
);
argList::noParallel();
argList::validArgs.append("hookTolerance");
# include "addDictOption.H"
# include "setRootCase.H"
# include "createTime.H"
const word dictName("surfaceHookUpDict");
# include "setSystemRunTimeDictionaryIO.H"
Info<< "Reading " << dictName << nl << endl;
const IOdictionary dict(dictIO);
const scalar dist(args.argRead<scalar>(1));
const scalar matchTolerance(SMALL);
Info<< "Hooking distance = " << dist << endl;
searchableSurfaces surfs
(
IOobject
(
"surfacesToHook",
runTime.constant(),
"triSurface",
runTime
),
dict
);
Info<< nl << "Reading surfaces: " << endl;
forAll(surfs, surfI)
{
Info<< incrIndent;
Info<< nl << indent << "Surface = " << surfs.names()[surfI] << endl;
const wordList& regions = surfs[surfI].regions();
forAll(regions, surfRegionI)
{
Info<< incrIndent;
Info<< indent << "Regions = " << regions[surfRegionI] << endl;
Info<< decrIndent;
}
Info<< decrIndent;
}
PtrList<indexedOctree<treeDataEdge> > bEdgeTrees(surfs.size());
labelListList treeBoundaryEdges(surfs.size());
List<DynamicList<labelledTri> > newFaces(surfs.size());
List<DynamicList<point> > newPoints(surfs.size());
List<PackedBoolList> visitedFace(surfs.size());
PtrList<triSurfaceMesh> newSurfaces(surfs.size());
forAll(surfs, surfI)
{
const triSurfaceMesh& surf =
refCast<const triSurfaceMesh>(surfs[surfI]);
newSurfaces.set
(
surfI,
new triSurfaceMesh
(
IOobject
(
"hookedSurface_" + surfs.names()[surfI],
runTime.constant(),
"triSurface",
runTime
),
surf
)
);
}
label nChanged = 0;
label nIters = 0;
do
{
Info<< nl << "Iteration = " << nIters++ << endl;
nChanged = 0;
createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);
forAll(newSurfaces, surfI)
{
const triSurface& newSurf = newSurfaces[surfI];
newFaces[surfI] = newSurf.localFaces();
newPoints[surfI] = newSurf.localPoints();
visitedFace[surfI] = PackedBoolList(newSurf.size(), false);
}
forAll(newSurfaces, surfI)
{
const triSurface& surf = newSurfaces[surfI];
List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
labelList bPointsHitTree(surf.boundaryPoints().size(), -1);
const labelListList& pointEdges = surf.pointEdges();
forAll(bPointsTobEdges, bPointI)
{
pointIndexHit& nearestHit = bPointsTobEdges[bPointI];
const label pointI = surf.boundaryPoints()[bPointI];
const point& samplePt = surf.localPoints()[pointI];
const labelList& pEdges = pointEdges[pointI];
// Add edges connected to the edge to the shapeMask
DynamicList<label> shapeMask;
shapeMask.append(pEdges);
forAll(bEdgeTrees, treeI)
{
const indexedOctree<treeDataEdge>& bEdgeTree =
bEdgeTrees[treeI];
pointIndexHit currentHit =
bEdgeTree.findNearest
(
samplePt,
sqr(dist),
treeDataEdge::findNearestOpSubset
(
bEdgeTree,
shapeMask
)
);
if
(
currentHit.hit()
&&
(
!nearestHit.hit()
||
(
magSqr(currentHit.hitPoint() - samplePt)
< magSqr(nearestHit.hitPoint() - samplePt)
)
)
)
{
nearestHit = currentHit;
bPointsHitTree[bPointI] = treeI;
}
}
scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);
if (nearestHit.hit())
{
// bool rejectEdge =
// checkEdgeAngle
// (
// surf,
// nearestHit.index(),
// pointI,
// 30
// );
if (dist2 > Foam::sqr(dist))
{
nearestHit.setMiss();
}
}
}
forAll(bPointsTobEdges, bPointI)
{
const pointIndexHit& eHit = bPointsTobEdges[bPointI];
if (eHit.hit())
{
const label hitSurfI = bPointsHitTree[bPointI];
const triSurface& hitSurf = newSurfaces[hitSurfI];
const label eIndex =
treeBoundaryEdges[hitSurfI][eHit.index()];
const edge& e = hitSurf.edges()[eIndex];
const label pointI = surf.boundaryPoints()[bPointI];
const labelList& eFaces = hitSurf.edgeFaces()[eIndex];
if (eFaces.size() != 1)
{
WarningIn(args.executable())
<< "Edge is attached to " << eFaces.size()
<< " faces." << endl;
continue;
}
const label faceI = eFaces[0];
if (visitedFace[hitSurfI][faceI])
{
continue;
}
DynamicList<labelledTri> newFacesFromSplit(2);
const point& pt = surf.localPoints()[pointI];
if
(
(
magSqr(pt - hitSurf.localPoints()[e.start()])
< matchTolerance
)
|| (
magSqr(pt - hitSurf.localPoints()[e.end()])
< matchTolerance
)
)
{
continue;
}
nChanged++;
// Keep the points in the same place and move the edge
newPoints[hitSurfI].append(newPoints[surfI][pointI]);
// Move the points to the edges
//newPoints[pointI] = eHit.hitPoint();
//newPoints.append(eHit.hitPoint());
visitedFace[hitSurfI][faceI] = true;
// Split the other face.
greenRefine
(
hitSurf,
faceI,
eIndex,
newPoints[hitSurfI].size() - 1,
newFacesFromSplit
);
forAll(newFacesFromSplit, newFaceI)
{
if (newFaceI == 0)
{
newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
}
else
{
newFaces[hitSurfI].append
(
newFacesFromSplit[newFaceI]
);
}
}
}
}
}
Info<< " Number of edges split = " << nChanged << endl;
forAll(newSurfaces, surfI)
{
newSurfaces.set
(
surfI,
new triSurfaceMesh
(
IOobject
(
"hookedSurface_" + surfs.names()[surfI],
runTime.constant(),
"triSurface",
runTime
),
triSurface
(
newFaces[surfI],
newSurfaces[surfI].patches(),
pointField(newPoints[surfI])
)
)
);
}
} while (nChanged > 0);
Info<< endl;
forAll(newSurfaces, surfI)
{
const triSurfaceMesh& newSurf = newSurfaces[surfI];
Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
<< endl;
newSurf.searchableSurface::write();
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,22 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object surfaceHookUpDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
surface1.stl {type triSurfaceMesh;}
surface2.stl {type triSurfaceMesh;}
// ************************************************************************* //

View File

@ -54,7 +54,7 @@ int main(int argc, char *argv[])
const fileName outFileName = args[3];
Info<< "Reading surface from " << surfFileName << " ..." << endl;
Info<< "Merging points within " << mergeTol << " meter." << endl;
Info<< "Merging points within " << mergeTol << " metre." << endl;
triSurface surf1(surfFileName);