From 139c22716d7897e8c0f651048260ff0e45ddc8ab Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 22 Feb 2012 16:34:58 +0000 Subject: [PATCH 1/5] ENH: renumberMesh: make frontwidth calculation an option --- .../manipulation/renumberMesh/Make/options | 2 +- .../manipulation/renumberMesh/renumberMesh.C | 35 +++++++++++++------ .../renumberMesh/renumberMeshDict | 1 + 3 files changed, 26 insertions(+), 12 deletions(-) diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Make/options b/applications/utilities/mesh/manipulation/renumberMesh/Make/options index 72eeafc9d1..23661b3e95 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options +++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -DFULLDEBUG -g -O0 \ + /* -DFULLDEBUG -g -O0 */ \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index 4120f8c67a..9c5d0b8de2 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -108,6 +108,7 @@ label getBand(const labelList& owner, const labelList& neighbour) // Calculate band of matrix void getBand ( + const bool calculateIntersect, const label nCells, const labelList& owner, const labelList& neighbour, @@ -129,17 +130,22 @@ void getBand cellBandwidth[nei] = max(cellBandwidth[nei], diff); } - forAll(nIntersect, cellI) - { - for (label colI = cellI-cellBandwidth[cellI]; colI <= cellI; colI++) - { - nIntersect[colI]++; - } - } - bandwidth = max(cellBandwidth); - profile = sum(cellBandwidth); - sumSqrIntersect = sum(Foam::sqr(nIntersect)); + profile = sum(1.0*cellBandwidth); + + sumSqrIntersect = 0.0; + if (calculateIntersect) + { + forAll(nIntersect, cellI) + { + for (label colI = cellI-cellBandwidth[cellI]; colI <= cellI; colI++) + { + nIntersect[colI] += 1.0; + } + } + + sumSqrIntersect = sum(Foam::sqr(nIntersect)); + } } @@ -565,6 +571,11 @@ int main(int argc, char *argv[]) "dict", "renumber according to system/renumberMeshDict" ); + argList::addBoolOption + ( + "frontWidth", + "calculate the rms of the frontwidth" + ); # include "setRootCase.H" # include "createTime.H" @@ -582,7 +593,7 @@ int main(int argc, char *argv[]) const word oldInstance = mesh.pointsInstance(); const bool readDict = args.optionFound("dict"); - + const bool doFrontWidth = args.optionFound("frontWidth"); const bool overwrite = args.optionFound("overwrite"); label band; @@ -590,6 +601,7 @@ int main(int argc, char *argv[]) scalar sumSqrIntersect; getBand ( + doFrontWidth, mesh.nCells(), mesh.faceOwner(), mesh.faceNeighbour(), @@ -1028,6 +1040,7 @@ int main(int argc, char *argv[]) scalar sumSqrIntersect; getBand ( + doFrontWidth, mesh.nCells(), mesh.faceOwner(), mesh.faceNeighbour(), diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict index 8878469010..7430ea3441 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict @@ -41,6 +41,7 @@ method CuthillMcKee; //method manual; //method random; //method spring; +//method boundaryFirst; //CuthillMcKeeCoeffs //{ From f9262f751f30ec47e7242ceeae66cbaa903ae3b1 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 23 Feb 2012 14:50:23 +0000 Subject: [PATCH 2/5] STYLE: randomRenumber: updated comment --- src/renumber/renumberMethods/randomRenumber/randomRenumber.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/renumber/renumberMethods/randomRenumber/randomRenumber.H b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H index f2860eac2a..571efc2833 100644 --- a/src/renumber/renumberMethods/randomRenumber/randomRenumber.H +++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H @@ -76,7 +76,7 @@ public: //- Return the order in which cells need to be visited, i.e. // from ordered back to original cell label. - // Use the mesh connectivity (if needed) + // This is only defined for geometric renumberMethods. virtual labelList renumber(const pointField&) const; //- Return the order in which cells need to be visited, i.e. From f3590edfa0737e33090f45c76e8f3a510ed2bc72 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 23 Feb 2012 14:55:28 +0000 Subject: [PATCH 3/5] BUG: Time: restart with non-standard Time (e.g. engineTime) compared values instead of names. --- src/OpenFOAM/db/Time/Time.C | 49 ++++++++++++++++++++++-------- src/OpenFOAM/db/Time/TimeIO.C | 9 ++++-- src/engine/engineTime/engineTime.C | 9 +++--- 3 files changed, 47 insertions(+), 20 deletions(-) diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C index 7ba18b3cc8..19a6c2b321 100644 --- a/src/OpenFOAM/db/Time/Time.C +++ b/src/OpenFOAM/db/Time/Time.C @@ -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-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -221,20 +221,43 @@ void Foam::Time::setControls() timeIndex_ = startTimeIndex_; } - scalar timeValue; - if (timeDict.readIfPresent("value", timeValue)) - { - word storedTimeName(timeName(timeValue)); - if (storedTimeName != timeName()) + // Check if values stored in time dictionary are consistent + + // 1. Based on time name + bool checkValue = true; + + string storedTimeName; + if (timeDict.readIfPresent("name", storedTimeName)) + { + if (storedTimeName == timeName()) { - IOWarningIn("Time::setControls()", timeDict) - << "Time read from time dictionary " << storedTimeName - << " differs from actual time " << timeName() << '.' << nl - << " This may cause unexpected database behaviour." - << " If you are not interested" << nl - << " in preserving time state delete the time dictionary." - << endl; + // Same time. No need to check stored value + checkValue = false; + } + } + + // 2. Based on time value + // (consistent up to the current time writing precision so it won't + // trigger if we just change the write precision) + if (checkValue) + { + scalar storedTimeValue; + if (timeDict.readIfPresent("value", storedTimeValue)) + { + word storedTimeName(timeName(storedTimeValue)); + + if (storedTimeName != timeName()) + { + IOWarningIn("Time::setControls()", timeDict) + << "Time read from time dictionary " << storedTimeName + << " differs from actual time " << timeName() << '.' << nl + << " This may cause unexpected database behaviour." + << " If you are not interested" << nl + << " in preserving time state delete" + << " the time dictionary." + << endl; + } } } } diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C index fcf0a969e2..557ebdddcb 100644 --- a/src/OpenFOAM/db/Time/TimeIO.C +++ b/src/OpenFOAM/db/Time/TimeIO.C @@ -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-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -293,12 +293,14 @@ bool Foam::Time::writeObject { if (outputTime()) { + const word tmName(timeName()); + IOdictionary timeDict ( IOobject ( "time", - timeName(), + tmName, "uniform", *this, IOobject::NO_READ, @@ -308,6 +310,7 @@ bool Foam::Time::writeObject ); timeDict.add("value", value()); + timeDict.add("name", string(tmName)); timeDict.add("index", timeIndex_); timeDict.add("deltaT", deltaT_); timeDict.add("deltaT0", deltaT0_); @@ -317,7 +320,7 @@ bool Foam::Time::writeObject if (writeOK && purgeWrite_) { - previousOutputTimes_.push(timeName()); + previousOutputTimes_.push(tmName); while (previousOutputTimes_.size() > purgeWrite_) { diff --git a/src/engine/engineTime/engineTime.C b/src/engine/engineTime/engineTime.C index f662f58596..0933218bd7 100644 --- a/src/engine/engineTime/engineTime.C +++ b/src/engine/engineTime/engineTime.C @@ -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-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,9 +91,10 @@ Foam::engineTime::engineTime timeAdjustment(); - startTime_ = degToTime(startTime_); - value() = degToTime(value()); - deltaT0_ = deltaT_; + startTime_ = degToTime(startTime_); + value() = degToTime(value()); + deltaTSave_ = deltaT_; + deltaT0_ = deltaT_; } From 0d4c02c839ada050764c4259da78b735f4e50854 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 23 Feb 2012 15:25:03 +0000 Subject: [PATCH 4/5] BUG: mergePolyMesh: zone indexing bug --- .../utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C index 8c63dcaf9f..df572a3619 100644 --- a/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C +++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C @@ -203,7 +203,7 @@ void Foam::mergePolyMesh::addMesh(const polyMesh& m) // Grab zone ID. If a point is not in a zone, it will return -1 zoneID = pz.whichZone(pointI); - if (zoneID > 0) + if (zoneID >= 0) { // Translate zone ID into the new index zoneID = pointZoneIndices[zoneID]; @@ -240,7 +240,7 @@ void Foam::mergePolyMesh::addMesh(const polyMesh& m) // Grab zone ID. If a cell is not in a zone, it will return -1 zoneID = cz.whichZone(cellI); - if (zoneID > 0) + if (zoneID >= 0) { // Translate zone ID into the new index zoneID = cellZoneIndices[zoneID]; @@ -343,7 +343,7 @@ void Foam::mergePolyMesh::addMesh(const polyMesh& m) newZone = fz.whichZone(faceI); newZoneFlip = false; - if (newZone > -1) + if (newZone >= 0) { newZoneFlip = fz[newZone].flipMap()[fz[newZone].whichFace(faceI)]; From 7e286757f703c67cd9061e7e744fd3f6bce66026 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 23 Feb 2012 15:47:58 +0000 Subject: [PATCH 5/5] ENH: tetrahedron: face access function --- .../primitiveShapes/tetrahedron/tetrahedron.H | 6 +++- .../tetrahedron/tetrahedronI.H | 35 +++++++++++++++++++ 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index 1a76022205..5ff25795b2 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H @@ -46,6 +46,7 @@ SourceFiles #include "Random.H" #include "FixedList.H" #include "UList.H" +#include "triPointRef.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -196,6 +197,8 @@ public: inline const Point& d() const; + //- Return i-th face + inline triPointRef tri(const label faceI) const; // Properties @@ -242,7 +245,8 @@ public: List& bary ) const; - //- Return nearest point to p on tetrahedron + //- Return nearest point to p on tetrahedron. Is p itself + // if inside. inline pointHit nearestPoint(const point& p) const; //- Return true if point is inside tetrahedron diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 8a69e88e04..bc1791a946 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -97,6 +97,41 @@ inline const Point& Foam::tetrahedron::d() const } +template +inline Foam::triPointRef Foam::tetrahedron::tri +( + const label faceI +) const +{ + // Warning. Ordering of faces needs to be the same for a tetrahedron + // class, a tetrahedron cell shape model and a tetCell + + if (faceI == 0) + { + return triPointRef(b_, c_, d_); + } + else if (faceI == 1) + { + return triPointRef(a_, d_, c_); + } + else if (faceI == 2) + { + return triPointRef(a_, b_, d_); + } + else if (faceI == 3) + { + return triPointRef(a_, c_, b_); + } + else + { + FatalErrorIn("tetrahedron::tri(const label faceI) const") + << "index out of range 0 -> 3. faceI = " << faceI + << abort(FatalError); + return triPointRef(b_, c_, d_); + } +} + + template inline Foam::vector Foam::tetrahedron::Sa() const {