Permalink
Browse files

ENH: PatchEdgeFaceWave: new wave method

  • Loading branch information...
1 parent 859ebf8 commit 8071d6bd29314fd0423082b81aa80acea1b96409 mattijs committed Oct 17, 2011
@@ -0,0 +1,3 @@
+Test-PatchEdgeFaceWave.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-PatchEdgeFaceWave
@@ -0,0 +1,7 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lmeshTools
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011 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/>.
+
+Description
+ Test PatchEdgeFaceWave.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "PatchEdgeFaceWave.H"
+#include "patchEdgeFaceInfo.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ argList::validArgs.append("patch");
+
+# include "setRootCase.H"
+# include "createTime.H"
+# include "createMesh.H"
+
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+ // Get name of patch
+ const word patchName = args[1];
+
+ const polyPatch& patch = patches[patchName];
+
+ // Data on all edges and faces
+ List<patchEdgeFaceInfo> allEdgeInfo(patch.nEdges());
+ List<patchEdgeFaceInfo> allFaceInfo(patch.size());
+
+ // Initial seed
+ DynamicList<label> initialEdges;
+ DynamicList<patchEdgeFaceInfo> initialEdgesInfo;
+
+
+ // Just set an edge on the master
+ if (Pstream::master())
+ {
+ label edgeI = 0;
+ Info<< "Starting walk on edge " << edgeI << endl;
+
+ initialEdges.append(edgeI);
+ const edge& e = patch.edges()[edgeI];
+ initialEdgesInfo.append
+ (
+ patchEdgeFaceInfo
+ (
+ e.centre(patch.localPoints()),
+ 0.0
+ )
+ );
+ }
+
+
+ // Walk
+ PatchEdgeFaceWave
+ <
+ primitivePatch,
+ patchEdgeFaceInfo
+ > calc
+ (
+ mesh,
+ patch,
+ initialEdges,
+ initialEdgesInfo,
+ allEdgeInfo,
+ allFaceInfo,
+ returnReduce(patch.nEdges(), sumOp<label>())
+ );
+
+
+ // Extract as patchField
+ volScalarField vsf
+ (
+ IOobject
+ (
+ "patchDist",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("patchDist", dimLength, 0.0)
+ );
+ scalarField pf(vsf.boundaryField()[patch.index()].size());
+ forAll(pf, faceI)
+ {
+ pf[faceI] = Foam::sqrt(allFaceInfo[faceI].distSqr());
+ }
+ vsf.boundaryField()[patch.index()] = pf;
+
+ Info<< "Writing patchDist volScalarField to " << runTime.value()
+ << endl;
+
+ vsf.write();
+
+
+ Info<< "\nEnd\n" << endl;
+ return 0;
+}
+
+
+// ************************************************************************* //
@@ -46,6 +46,20 @@ defineTypeNameAndDebug(Foam::globalMeshData, 0);
// Geometric matching tolerance. Factor of mesh bounding box.
const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8;
+namespace Foam
+{
+template<>
+class minEqOp<labelPair>
+{
+public:
+ void operator()(labelPair& x, const labelPair& y) const
+ {
+ x[0] = min(x[0], y[0]);
+ x[1] = min(x[1], y[1]);
+ }
+};
+}
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -1063,6 +1077,128 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
}
+void Foam::globalMeshData::calcGlobalEdgeOrientation() const
+{
+ if (debug)
+ {
+ Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
+ << " calculating edge orientation w.r.t. master edge." << endl;
+ }
+
+ const globalIndex& globalPoints = globalPointNumbering();
+
+ // 1. Determine master point
+ labelList masterPoint;
+ {
+ const mapDistribute& map = globalPointSlavesMap();
+
+ masterPoint.setSize(map.constructSize());
+ masterPoint = labelMax;
+
+ for (label pointI = 0; pointI < coupledPatch().nPoints(); pointI++)
+ {
+ masterPoint[pointI] = globalPoints.toGlobal(pointI);
+ }
+ syncData
+ (
+ masterPoint,
+ globalPointSlaves(),
+ globalPointTransformedSlaves(),
+ map,
+ minEqOp<label>()
+ );
+ }
+
+ // Now all points should know who is master by comparing their global
+ // pointID with the masterPointID. We now can use this information
+ // to find the orientation of the master edge.
+
+ {
+ const mapDistribute& map = globalEdgeSlavesMap();
+ const labelListList& slaves = globalEdgeSlaves();
+ const labelListList& transformedSlaves = globalEdgeTransformedSlaves();
+
+ // Distribute orientation of master edge (in masterPoint numbering)
+ labelPairList masterEdgeVerts(map.constructSize());
+ masterEdgeVerts = labelPair(labelMax, labelMax);
+
+ for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++)
+ {
+ if
+ (
+ (
+ slaves[edgeI].size()
+ + transformedSlaves[edgeI].size()
+ )
+ > 0
+ )
+ {
+ // I am master. Fill in my masterPoint equivalent.
+
+ const edge& e = coupledPatch().edges()[edgeI];
+ masterEdgeVerts[edgeI] = labelPair
+ (
+ masterPoint[e[0]],
+ masterPoint[e[1]]
+ );
+ }
+ }
+ syncData
+ (
+ masterEdgeVerts,
+ slaves,
+ transformedSlaves,
+ map,
+ minEqOp<labelPair>()
+ );
+
+ // Now check my edges on how they relate to the master's edgeVerts
+ globalEdgeOrientationPtr_.reset
+ (
+ new PackedBoolList(coupledPatch().nEdges())
+ );
+ PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_();
+
+ forAll(coupledPatch().edges(), edgeI)
+ {
+ const edge& e = coupledPatch().edges()[edgeI];
+ const labelPair masterE
+ (
+ masterPoint[e[0]],
+ masterPoint[e[1]]
+ );
+
+ label stat = labelPair::compare
+ (
+ masterE,
+ masterEdgeVerts[edgeI]
+ );
+ if (stat == 0)
+ {
+ FatalErrorIn
+ (
+ "globalMeshData::calcGlobalEdgeOrientation() const"
+ ) << "problem : my edge:" << e
+ << " in master points:" << masterE
+ << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI]
+ << exit(FatalError);
+ }
+ else
+ {
+ globalEdgeOrientation[edgeI] = (stat == 1);
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
+ << " finished calculating edge orientation."
+ << endl;
+ }
+}
+
+
// Calculate uncoupled boundary faces (without calculating
// primitiveMesh::pointFaces())
void Foam::globalMeshData::calcPointBoundaryFaces
@@ -1660,6 +1796,7 @@ void Foam::globalMeshData::clearOut()
globalEdgeNumberingPtr_.clear();
globalEdgeSlavesPtr_.clear();
globalEdgeTransformedSlavesPtr_.clear();
+ globalEdgeOrientationPtr_.clear();
globalEdgeSlavesMapPtr_.clear();
// Face
@@ -2095,6 +2232,16 @@ const
}
+const Foam::PackedBoolList& Foam::globalMeshData::globalEdgeOrientation() const
+{
+ if (!globalEdgeOrientationPtr_.valid())
+ {
+ calcGlobalEdgeOrientation();
+ }
+ return globalEdgeOrientationPtr_();
+}
+
+
const Foam::mapDistribute& Foam::globalMeshData::globalEdgeSlavesMap() const
{
if (!globalEdgeSlavesMapPtr_.valid())
@@ -95,6 +95,7 @@ class mapDistribute;
template<class T> class EdgeMap;
class globalIndex;
class globalIndexAndTransform;
+class PackedBoolList;
/*---------------------------------------------------------------------------*\
Class globalMeshData Declaration
@@ -191,6 +192,7 @@ class globalMeshData
mutable autoPtr<globalIndex> globalEdgeNumberingPtr_;
mutable autoPtr<labelListList> globalEdgeSlavesPtr_;
mutable autoPtr<labelListList> globalEdgeTransformedSlavesPtr_;
+ mutable autoPtr<PackedBoolList> globalEdgeOrientationPtr_;
mutable autoPtr<mapDistribute> globalEdgeSlavesMapPtr_;
@@ -297,6 +299,9 @@ class globalMeshData
//- Calculate global edge addressing.
void calcGlobalEdgeSlaves() const;
+ //- Calculate orientation w.r.t. edge master.
+ void calcGlobalEdgeOrientation() const;
+
// Global boundary face/cell addressing
@@ -539,6 +544,8 @@ public:
const labelListList& globalEdgeSlaves() const;
const labelListList& globalEdgeTransformedSlaves() const;
const mapDistribute& globalEdgeSlavesMap() const;
+ //- Is my edge same orientation master edge
+ const PackedBoolList& globalEdgeOrientation() const;
// Coupled point to boundary faces. These are uncoupled boundary
// faces only but include empty patches.
@@ -30,6 +30,7 @@ License
#include "PatchToolsSearch.C"
#include "PatchToolsSortEdges.C"
#include "PatchToolsNormals.C"
+#include "PatchToolsMatch.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Oops, something went wrong.

0 comments on commit 8071d6b

Please sign in to comment.