diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index a993ecf92ae31cf47dcd94125b5e2cee77ebfbb8..7bd992fdc1f2149bc14443e62487f3be459186d6 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -225,6 +225,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
// Dump halves
{
OFstream str(prefix+cycPatch.name()+"_half0.obj");
+ Pout<< "Dumping cycPatch.name() half0 faces to " << str.name()
+ << endl;
meshTools::writeOBJ
(
str,
@@ -241,6 +243,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
}
{
OFstream str(prefix+cycPatch.name()+"_half1.obj");
+ Pout<< "Dumping cycPatch.name() half1 faces to " << str.name()
+ << endl;
meshTools::writeOBJ
(
str,
@@ -262,6 +266,9 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
OFstream str(prefix+cycPatch.name()+"_match.obj");
label vertI = 0;
+ Pout<< "Dumping cyclic match as lines between face centres to "
+ << str.name() << endl;
+
for (label faceI = 0; faceI < halfSize; faceI++)
{
const point& fc0 = mesh.faceCentres()[cycPatch.start()+faceI];
@@ -773,6 +780,8 @@ int main(int argc, char *argv[])
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, true);
mesh.movePoints(map().preMotionPoints());
+ dumpCyclicMatch("coupled_", mesh);
+
// Synchronise points.
if (!pointSync)
{
@@ -890,6 +899,8 @@ int main(int argc, char *argv[])
filterPatches(mesh);
+ dumpCyclicMatch("final_", mesh);
+
// Set the precision of the points data to 10
IOstream::defaultPrecision(10);
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
index 5f3597f21af7ebb740e5dac56e5076c26fee4953..9052daa0a8d56c2ab9c82be8d3ccc42511089532 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
@@ -1,77 +1,96 @@
-/*---------------------------------------------------------------------------*\
+/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
-| \\ / O peration | Version: 1.0 |
-| \\ / A nd | Web: http://www.openfoam.org |
+| \\ / O peration | Version: 1.5 |
+| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
-
FoamFile
{
- version 2.0;
- format ascii;
-
- root "";
- case "";
- instance "system";
- local "";
-
- class dictionary;
- object createPatcheDict;
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object createPatchDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// This application/dictionary controls:
+// - optional: create new patches from boundary faces (either given as
+// a set of patches or as a faceSet)
+// - always: order faces on coupled patches such that they are opposite. This
+// is done for all coupled faces, not just for any patches created.
+// - optional: synchronise points on coupled patches.
+
+// 1. Create cyclic:
+// - specify where the faces should come from
+// - specify the type of cyclic. If a rotational specify the rotationAxis
+// and centre to make matching easier
+// - pointSync true to guarantee points to line up.
+
+// 2. Correct incorrect cyclic:
+// This will usually fail upon loading:
+// "face 0 area does not match neighbour 2 by 0.0100005%"
+// " -- possible face ordering problem."
+// - change patch type from 'cyclic' to 'patch' in the polyMesh/boundary file.
+// - loosen match tolerance to get case to load
+// - regenerate cyclic as above
+
+
// Tolerance used in matching faces. Absolute tolerance is span of
-// face times this factor.
-matchTolerance 1E-6;
+// face times this factor. To load incorrectly matches meshes set this
+// to a higher value.
+matchTolerance 1E-3;
-// Do a synchronisation of coupled points.
+// Do a synchronisation of coupled points after creation of any patches.
pointSync true;
-
// Patches to create.
-// If no patches does a coupled point and face synchronisation anyway.
patches
(
{
// Name of new patch
name sidePatches;
- // Dictionary for new patch
- dictionary
- {
+ // Type of new patch
+ dictionary
+ {
type cyclic;
- // Optional: used when matching and synchronising points.
+
+ // Optional: explicitly set transformation tensor.
+ // Used when matching and synchronising points.
//transform translational;
//separationVector (-2289 0 0);
- }
+ transform rotational;
+ rotationAxis (1 0 0);
+ rotationCentre (0 0 0);
+ }
- // How to construct: either 'patches' or 'set'
+ // How to construct: either from 'patches' or 'set'
constructFrom patches;
// If constructFrom = patches : names of patches
- //patches (periodic-1 periodic-2);
- patches (outlet-side1 outlet-side2);
+ patches (periodic-1 periodic-2);
// If constructFrom = set : name of faceSet
set f0;
}
- //{
- // name bottom;
- // // Dictionary for new patch
- // dictionary
- // {
- // type patch;
- // }
- //
- // constructFrom set;
- //
- // patches (half0 half1);
- //
- // set bottomFaces;
- //}
+ {
+ name bottom;
+
+ // Type of new patch
+ dictionary
+ {
+ type wall;
+ }
+
+ constructFrom set;
+
+ patches ();
+
+ set bottomFaces;
+ }
);
diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.H b/src/OpenFOAM/meshes/MeshObject/MeshObject.H
index 3bf11b5f6a3f588596ce5d6a65295135db4ee250..eac6096c793cabeaca7c8a4826ffd736e9a6122c 100644
--- a/src/OpenFOAM/meshes/MeshObject/MeshObject.H
+++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.H
@@ -104,12 +104,12 @@ public:
);
- // Destructor
-
- static bool Delete(const Mesh& mesh);
+ // Destructors
virtual ~MeshObject();
+ static bool Delete(const Mesh& mesh);
+
// Member Functions
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index 73aff8a1fd2c711402e5ddd3a201583166df211b..5fbaa3fa95f10fa529c2b933ba840f2579aed1b4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -320,9 +320,9 @@ void Foam::coupledPolyPatch::calcTransformTensors
if (debug)
{
- Pout<< " rotation " << sum(mag(forwardT_ - forwardT_[0]))
- << " more than local tolerance " << error
- << ". Assuming uniform rotation." << endl;
+ Pout<< " difference in rotation less than"
+ << " local tolerance "
+ << error << ". Assuming uniform rotation." << endl;
}
}
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 171fb099f1a08ce774a19e2421c94ee83249cdb4..c2458833dac92a1c02c2039354303d95458216f8 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -34,6 +34,7 @@ License
#include "matchPoints.H"
#include "EdgeMap.H"
#include "Time.H"
+#include "transformList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -89,6 +90,9 @@ void Foam::cyclicPolyPatch::calcTransforms()
{
const pointField& points = this->points();
+ // Determine geometric quantities on the two halves
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
primitivePatch half0
(
SubList<face>
@@ -199,15 +203,69 @@ void Foam::cyclicPolyPatch::calcTransforms()
}
- // Calculate transformation tensors
- calcTransformTensors
- (
- half0Ctrs,
- half1Ctrs,
- half0Normals,
- half1Normals,
- half0Tols
- );
+ // See if transformation is prescribed
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ switch (transform_)
+ {
+ case ROTATIONAL:
+ {
+ // Specified single rotation tensor.
+
+ // Get best fitting face and its opposite number
+ label face0 = getConsistentRotationFace(half0Ctrs);
+ label face1 = face0;
+
+ vector n0 =
+ (
+ (half0Ctrs[face0] - rotationCentre_)
+ ^ rotationAxis_
+ );
+ vector n1 =
+ (
+ (half1Ctrs[face1] - rotationCentre_)
+ ^ -rotationAxis_
+ );
+ n0 /= mag(n0) + VSMALL;
+ n1 /= mag(n1) + VSMALL;
+
+ if (debug)
+ {
+ Pout<< "cyclicPolyPatch::calcTransforms :"
+ << " Specified rotation :"
+ << " n0:" << n0 << " n1:" << n1 << endl;
+ }
+
+ // Calculate transformation tensors from face0,1 only.
+ // Note: can use tight tolerance now.
+ calcTransformTensors
+ (
+ pointField(1, half0Ctrs[face0]),
+ pointField(1, half1Ctrs[face1]),
+ vectorField(1, n0),
+ vectorField(1, n1),
+ scalarField(1, half0Tols[face0]),
+ 1E-4
+ );
+
+ break;
+ }
+
+ default:
+ {
+ // Calculate transformation tensors from all faces.
+ calcTransformTensors
+ (
+ half0Ctrs,
+ half1Ctrs,
+ half0Normals,
+ half1Normals,
+ half0Tols
+ );
+
+ break;
+ }
+ }
}
}
@@ -402,6 +460,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
const faceList& half0Faces,
const faceList& half1Faces,
+ pointField& ppPoints,
pointField& half0Ctrs,
pointField& half1Ctrs,
pointField& anchors0,
@@ -442,6 +501,8 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
}
+ ppPoints = Foam::transform(reverseT, pp.points());
+
break;
}
//- Problem: usually specified translation is not accurate enough
@@ -501,6 +562,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
anchors0[faceI]
);
}
+ ppPoints = Foam::transform(reverseT, pp.points());
}
else
{
@@ -524,6 +586,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
half0Ctrs += ctr1 - ctr0;
anchors0 += ctr1 - ctr0;
+ ppPoints = pp.points() + ctr1 - ctr0;
}
break;
}
@@ -1079,7 +1142,7 @@ bool Foam::cyclicPolyPatch::order
faceList half1Faces(IndirectList<face>(pp, half1ToPatch));
// Get geometric quantities
- pointField half0Ctrs, half1Ctrs, anchors0;
+ pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
scalarField tols;
getCentresAndAnchors
(
@@ -1087,6 +1150,7 @@ bool Foam::cyclicPolyPatch::order
half0Faces,
half1Faces,
+ ppPoints,
half0Ctrs,
half1Ctrs,
anchors0,
@@ -1108,6 +1172,44 @@ bool Foam::cyclicPolyPatch::order
{
Pout<< "cyclicPolyPatch::order : test if already ordered:"
<< matchedAll << endl;
+
+ // Dump halves
+ fileName nm0("match1_"+name()+"_half0_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half0"
+ << " faces to OBJ file " << nm0 << endl;
+ writeOBJ(nm0, half0Faces, ppPoints);
+
+ fileName nm1("match1_"+name()+"_half1_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half1"
+ << " faces to OBJ file " << nm1 << endl;
+ writeOBJ(nm1, half1Faces, pp.points());
+
+ OFstream ccStr
+ (
+ boundaryMesh().mesh().time().path()
+ /"match1_"+ name() + "_faceCentres.obj"
+ );
+ Pout<< "cyclicPolyPatch::order : "
+ << "Dumping currently found cyclic match as lines between"
+ << " corresponding face centres to file " << ccStr.name()
+ << endl;
+
+ // Recalculate untransformed face centres
+ //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+ label vertI = 0;
+
+ forAll(half1Ctrs, i)
+ {
+ //if (from1To0[i] != -1)
+ {
+ // Write edge between c1 and c0
+ //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+ //const point& c0 = half0Ctrs[from1To0[i]];
+ const point& c0 = half0Ctrs[i];
+ const point& c1 = half1Ctrs[i];
+ writeOBJ(ccStr, c0, c1, vertI);
+ }
+ }
}
@@ -1133,6 +1235,7 @@ bool Foam::cyclicPolyPatch::order
half0Faces,
half1Faces,
+ ppPoints,
half0Ctrs,
half1Ctrs,
anchors0,
@@ -1153,6 +1256,42 @@ bool Foam::cyclicPolyPatch::order
{
Pout<< "cyclicPolyPatch::order : test if pairwise ordered:"
<< matchedAll << endl;
+ // Dump halves
+ fileName nm0("match2_"+name()+"_half0_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half0"
+ << " faces to OBJ file " << nm0 << endl;
+ writeOBJ(nm0, half0Faces, ppPoints);
+
+ fileName nm1("match2_"+name()+"_half1_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half1"
+ << " faces to OBJ file " << nm1 << endl;
+ writeOBJ(nm1, half1Faces, pp.points());
+
+ OFstream ccStr
+ (
+ boundaryMesh().mesh().time().path()
+ /"match2_"+name()+"_faceCentres.obj"
+ );
+ Pout<< "cyclicPolyPatch::order : "
+ << "Dumping currently found cyclic match as lines between"
+ << " corresponding face centres to file " << ccStr.name()
+ << endl;
+
+ // Recalculate untransformed face centres
+ //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+ label vertI = 0;
+
+ forAll(half1Ctrs, i)
+ {
+ if (from1To0[i] != -1)
+ {
+ // Write edge between c1 and c0
+ //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+ const point& c0 = half0Ctrs[from1To0[i]];
+ const point& c1 = half1Ctrs[i];
+ writeOBJ(ccStr, c0, c1, vertI);
+ }
+ }
}
}
@@ -1209,6 +1348,7 @@ bool Foam::cyclicPolyPatch::order
half0Faces,
half1Faces,
+ ppPoints,
half0Ctrs,
half1Ctrs,
anchors0,
@@ -1229,8 +1369,43 @@ bool Foam::cyclicPolyPatch::order
{
Pout<< "cyclicPolyPatch::order : test if baffles:"
<< matchedAll << endl;
- }
+ // Dump halves
+ fileName nm0("match3_"+name()+"_half0_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half0"
+ << " faces to OBJ file " << nm0 << endl;
+ writeOBJ(nm0, half0Faces, ppPoints);
+
+ fileName nm1("match3_"+name()+"_half1_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half1"
+ << " faces to OBJ file " << nm1 << endl;
+ writeOBJ(nm1, half1Faces, pp.points());
+
+ OFstream ccStr
+ (
+ boundaryMesh().mesh().time().path()
+ /"match3_"+ name() + "_faceCentres.obj"
+ );
+ Pout<< "cyclicPolyPatch::order : "
+ << "Dumping currently found cyclic match as lines between"
+ << " corresponding face centres to file " << ccStr.name()
+ << endl;
+
+ // Recalculate untransformed face centres
+ //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+ label vertI = 0;
+ forAll(half1Ctrs, i)
+ {
+ if (from1To0[i] != -1)
+ {
+ // Write edge between c1 and c0
+ //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+ const point& c0 = half0Ctrs[from1To0[i]];
+ const point& c1 = half1Ctrs[i];
+ writeOBJ(ccStr, c0, c1, vertI);
+ }
+ }
+ }
}
}
@@ -1259,6 +1434,7 @@ bool Foam::cyclicPolyPatch::order
half0Faces,
half1Faces,
+ ppPoints,
half0Ctrs,
half1Ctrs,
anchors0,
@@ -1279,6 +1455,42 @@ bool Foam::cyclicPolyPatch::order
{
Pout<< "cyclicPolyPatch::order : automatic ordering result:"
<< matchedAll << endl;
+ // Dump halves
+ fileName nm0("match4_"+name()+"_half0_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half0"
+ << " faces to OBJ file " << nm0 << endl;
+ writeOBJ(nm0, half0Faces, ppPoints);
+
+ fileName nm1("match4_"+name()+"_half1_faces.obj");
+ Pout<< "cyclicPolyPatch::order : Writing half1"
+ << " faces to OBJ file " << nm1 << endl;
+ writeOBJ(nm1, half1Faces, pp.points());
+
+ OFstream ccStr
+ (
+ boundaryMesh().mesh().time().path()
+ /"match4_"+ name() + "_faceCentres.obj"
+ );
+ Pout<< "cyclicPolyPatch::order : "
+ << "Dumping currently found cyclic match as lines between"
+ << " corresponding face centres to file " << ccStr.name()
+ << endl;
+
+ // Recalculate untransformed face centres
+ //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+ label vertI = 0;
+
+ forAll(half1Ctrs, i)
+ {
+ if (from1To0[i] != -1)
+ {
+ // Write edge between c1 and c0
+ //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+ const point& c0 = half0Ctrs[from1To0[i]];
+ const point& c1 = half1Ctrs[i];
+ writeOBJ(ccStr, c0, c1, vertI);
+ }
+ }
}
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index 7db828985fb8e37dcb78820e0d1dec0989484893..af363e8a36cfe84ff636000370b60c84620c9abf 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -141,6 +141,7 @@ private:
const faceList& half0Faces,
const faceList& half1Faces,
+ pointField& ppPoints,
pointField& half0Ctrs,
pointField& half1Ctrs,
pointField& anchors0,
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 71707195f37e79af5bd4eeb3d9edc539ed4ca70c..d981c4eda52c9a24ff5864edf79b383aed0d3a1d 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -152,6 +152,11 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
{
mesh.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh.clearOut();
+ }
faceCombiner.updateMesh(map);
@@ -301,6 +306,11 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
{
mesh.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh.clearOut();
+ }
faceCombiner.updateMesh(map);
@@ -363,6 +373,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
{
mesh.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh.clearOut();
+ }
pointRemover.updateMesh(map);
meshRefiner_.updateMesh(map, labelList(0));
@@ -411,6 +426,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
{
mesh.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh.clearOut();
+ }
pointRemover.updateMesh(map);
meshRefiner_.updateMesh(map, labelList(0));
@@ -2782,6 +2802,10 @@ void Foam::autoLayerDriver::addLayers
const_cast<Time&>(mesh.time())++;
Info<< "Writing shrunk mesh to " << mesh.time().timeName() << endl;
+
+ // See comment in autoSnapDriver why we should not remove meshPhi
+ // using mesh.clearPout().
+
mesh.write();
}
@@ -2970,6 +2994,11 @@ void Foam::autoLayerDriver::addLayers
{
mesh.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh.clearOut();
+ }
meshRefiner_.updateMesh(map, labelList(0));
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 3a6c79987c2ae84937b402759e40458616255be0..fdf89f4123b23c22d5e092eedb01e906ace6c6aa 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -930,6 +930,7 @@ void Foam::autoSnapDriver::preSmoothPatch
const_cast<Time&>(mesh.time())++;
Pout<< "Writing patch smoothed mesh to time " << mesh.time().timeName()
<< endl;
+
mesh.write();
}
@@ -1220,6 +1221,11 @@ void Foam::autoSnapDriver::smoothDisplacement
const_cast<Time&>(mesh.time())++;
Pout<< "Writing smoothed mesh to time " << mesh.time().timeName()
<< endl;
+
+ // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
+ // but this will also delete all pointMesh but not pointFields which
+ // gives an illegal situation.
+
mesh.write();
Pout<< "Writing displacement field ..." << endl;
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 15b78e96536f4396fda7f05e8c6adfd133a5ac8f..7b2ecfc395226c9d680f63fdf5301ab98afaa148 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -458,6 +458,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
{
mesh_.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh_.clearOut();
+ }
//- Redo the intersections on the newly create baffle faces. Note that
// this changes also the cell centre positions.
@@ -1448,6 +1453,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
{
mesh_.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh_.clearOut();
+ }
// Update intersections. Recalculate intersections on merged faces since
// this seems to give problems? Note: should not be nessecary since
@@ -2405,6 +2415,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
{
mesh_.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh_.clearOut();
+ }
// Update intersections. Is mapping only (no faces created, positions stay
// same) so no need to recalculate intersections.
@@ -2828,6 +2843,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
{
mesh_.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh_.clearOut();
+ }
return map;
}
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index b457f0f5c7f5f919e9cbf0807b07fa46155a6af5..72e752f40b9e9ddeb9d9b56b4ba43e6535caced8 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -835,7 +835,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// minLevel) and cache per cell the max surface level and the local normal
// on that surface.
labelList cellMaxLevel(mesh_.nCells(), -1);
- vectorField cellMaxNormal(mesh_.nCells());
+ vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
{
// Per segment the normals of the surfaces
@@ -1226,6 +1226,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
{
mesh_.movePoints(map().preMotionPoints());
}
+ else
+ {
+ // Delete mesh volumes.
+ mesh_.clearOut();
+ }
// Update intersection info
updateMesh(map, getChangedFaces(map, cellsToRefine));
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index e82103bdb935f31242bf5e0a191131ea8d1de836..efdd09d359a59f6de923ba4154e7ac1d1eb2bf9e 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.C
@@ -780,9 +780,6 @@ Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
- //!!! Workaround for movePoints bug
- const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh()).movePoints(newPoints);
-
pp_.movePoints(mesh_.points());
return tsweptVol;
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 7bf27115a57beb33d3ce99c1b17b5204e9b60038..2f301c040687b8da341e9e392dcbaa202a1c2c5f 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -186,10 +186,8 @@ $(schemes)/harmonic/harmonic.C
$(schemes)/localBlended/localBlended.C
$(schemes)/localMax/localMax.C
$(schemes)/localMin/localMin.C
-//$(schemes)/linearFit/linearFit.C
-//$(schemes)/linearFit/linearFitData.C
-$(schemes)/quadraticFit/quadraticFit.C
-$(schemes)/quadraticFit/quadraticFitData.C
+$(schemes)/linearFit/linearFit.C
+$(schemes)/quadraticLinearFit/quadraticLinearFit.C
limitedSchemes = $(surfaceInterpolation)/limitedSchemes
$(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
index 12d3179e62fc1f2a7ffdca439e0164887a952816..25ba241f1ec200ccf228de732c4b1642c88ade18 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
@@ -131,7 +131,6 @@ public:
const GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<scalar> >& stencilWeights
);
-
};
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H
index 704f9a92d459394366abf605d3f301276495c0b0..e36be07993bc097c65eb82bb6e0785ea2d59f73e 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H
@@ -70,10 +70,9 @@ public:
{}
- // Destructor
-
- virtual ~centredCFCStencilObject()
- {}
+ //- Destructor
+ virtual ~centredCFCStencilObject()
+ {}
};
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
index 480cf15a2a9412224c5e465dafde0fb9fc766b31..ec024775d357d0417d9c0af09e4097593b245bdd 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
@@ -135,19 +135,22 @@ Foam::extendedStencil::weightedSum
{
fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
- label faceI = pSfCorr.patch().patch().start();
-
- forAll(pSfCorr, i)
+ if (pSfCorr.coupled())
{
- const List<Type>& stField = stencilFld[faceI];
- const List<scalar>& stWeight = stencilWeights[faceI];
+ label faceI = pSfCorr.patch().patch().start();
- forAll(stField, j)
+ forAll(pSfCorr, i)
{
- pSfCorr[i] += stField[j]*stWeight[j];
- }
+ const List<Type>& stField = stencilFld[faceI];
+ const List<scalar>& stWeight = stencilWeights[faceI];
- faceI++;
+ forAll(stField, j)
+ {
+ pSfCorr[i] += stField[j]*stWeight[j];
+ }
+
+ faceI++;
+ }
}
}
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C
index 1cfa79eaa55ff6154579a2db7cbc7f200d6ae20d..00cd58e17d793060bc99dd7149c0e1f84c2ceffe 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C
@@ -102,32 +102,35 @@ Foam::extendedUpwindStencil::weightedSum
{
fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
- label faceI = pSfCorr.patch().patch().start();
-
- forAll(pSfCorr, i)
+ if (pSfCorr.coupled())
{
- if (phi[faceI] > 0)
- {
- // Flux out of owner. Use upwind (= owner side) stencil.
- const List<Type>& stField = ownFld[faceI];
- const List<scalar>& stWeight = ownWeights[faceI];
+ label faceI = pSfCorr.patch().patch().start();
- forAll(stField, j)
+ forAll(pSfCorr, i)
+ {
+ if (phi[faceI] > 0)
{
- pSfCorr[i] += stField[j]*stWeight[j];
+ // Flux out of owner. Use upwind (= owner side) stencil.
+ const List<Type>& stField = ownFld[faceI];
+ const List<scalar>& stWeight = ownWeights[faceI];
+
+ forAll(stField, j)
+ {
+ pSfCorr[i] += stField[j]*stWeight[j];
+ }
}
- }
- else
- {
- const List<Type>& stField = neiFld[faceI];
- const List<scalar>& stWeight = neiWeights[faceI];
-
- forAll(stField, j)
+ else
{
- pSfCorr[i] += stField[j]*stWeight[j];
+ const List<Type>& stField = neiFld[faceI];
+ const List<scalar>& stWeight = neiWeights[faceI];
+
+ forAll(stField, j)
+ {
+ pSfCorr[i] += stField[j]*stWeight[j];
+ }
}
+ faceI++;
}
- faceI++;
}
}
diff --git a/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H b/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H
index 93e9d9117aba8dd73803079c0683c05751290d1f..67179e46079ca9af3c9edea03f6b80b8391777a0 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H
+++ b/src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H
@@ -68,7 +68,6 @@ public:
//- Construct from mesh
explicit cellFaceCellStencil(const polyMesh& mesh);
-
};
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 3c2e6756ac4c65f454aee953bdf4604391f8a5f9..ff714d00eb6ad2eae99b4fba1f98ca659e07d33d 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -42,9 +42,9 @@ License
#include "extendedLeastSquaresVectors.H"
#include "extendedLeastSquaresVectors.H"
#include "leastSquaresVectors.H"
-//#include "linearFitData.H"
-#include "quadraticFitData.H"
-//#include "quadraticFitSnGradData.H"
+#include "CentredFitData.H"
+#include "linearFitPolynomial.H"
+#include "quadraticLinearFitPolynomial.H"
#include "skewCorrectionVectors.H"
#include "centredCECStencilObject.H"
@@ -89,15 +89,13 @@ void Foam::fvMesh::clearGeom()
// Mesh motion flux cannot be deleted here because the old-time flux
// needs to be saved.
-
// Things geometry dependent that are not updated.
volPointInterpolation::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
- //linearFitData::Delete(*this);
- quadraticFitData::Delete(*this);
- //quadraticFitSnGradData::Delete(*this);
+ CentredFitData<linearFitPolynomial>::Delete(*this);
+ CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
}
@@ -112,9 +110,8 @@ void Foam::fvMesh::clearAddressing()
extendedLeastSquaresVectors::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
- //linearFitData::Delete(*this);
- quadraticFitData::Delete(*this);
- //quadraticFitSnGradData::Delete(*this);
+ CentredFitData<linearFitPolynomial>::Delete(*this);
+ CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
centredCECStencilObject::Delete(*this);
@@ -299,7 +296,6 @@ Foam::fvMesh::~fvMesh()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-// Helper function for construction from pieces
void Foam::fvMesh::addFvPatches
(
const List<polyPatch*> & p,
@@ -487,6 +483,30 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
}
+// Temporary helper function to call move points on
+// MeshObjects
+template<class Type>
+void MeshObjectMovePoints(const Foam::fvMesh& mesh)
+{
+ if
+ (
+ mesh.db().objectRegistry::foundObject<Type>
+ (
+ Type::typeName
+ )
+ )
+ {
+ const_cast<Type&>
+ (
+ mesh.db().objectRegistry::lookupObject<Type>
+ (
+ Type::typeName
+ )
+ ).movePoints();
+ }
+}
+
+
Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
{
// Grab old time volumes if the time has been incremented
@@ -577,133 +597,12 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
// Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
// movePoints function.
-
- // volPointInterpolation
- if
- (
- db().objectRegistry::foundObject<volPointInterpolation>
- (
- volPointInterpolation::typeName
- )
- )
- {
- const_cast<volPointInterpolation&>
- (
- db().objectRegistry::lookupObject<volPointInterpolation>
- (
- volPointInterpolation::typeName
- )
- ).movePoints();
- }
-
- // extendedLeastSquaresVectors
- if
- (
- db().objectRegistry::foundObject<extendedLeastSquaresVectors>
- (
- extendedLeastSquaresVectors::typeName
- )
- )
- {
- const_cast<extendedLeastSquaresVectors&>
- (
- db().objectRegistry::lookupObject<extendedLeastSquaresVectors>
- (
- extendedLeastSquaresVectors::typeName
- )
- ).movePoints();
- }
-
- // leastSquaresVectors
- if
- (
- db().objectRegistry::foundObject<leastSquaresVectors>
- (
- leastSquaresVectors::typeName
- )
- )
- {
- const_cast<leastSquaresVectors&>
- (
- db().objectRegistry::lookupObject<leastSquaresVectors>
- (
- leastSquaresVectors::typeName
- )
- ).movePoints();
- }
-
- //// linearFitData
- //if
- //(
- // db().objectRegistry::foundObject<linearFitData>
- // (
- // linearFitData::typeName
- // )
- //)
- //{
- // const_cast<linearFitData&>
- // (
- // db().objectRegistry::lookupObject<linearFitData>
- // (
- // linearFitData::typeName
- // )
- // ).movePoints();
- //}
-
- // quadraticFitData
- if
- (
- db().objectRegistry::foundObject<quadraticFitData>
- (
- quadraticFitData::typeName
- )
- )
- {
- const_cast<quadraticFitData&>
- (
- db().objectRegistry::lookupObject<quadraticFitData>
- (
- quadraticFitData::typeName
- )
- ).movePoints();
- }
-
- //// quadraticFitSnGradData
- //if
- //(
- // db().objectRegistry::foundObject<quadraticFitSnGradData>
- // (
- // quadraticFitSnGradData::typeName
- // )
- //)
- //{
- // const_cast<quadraticFitSnGradData&>
- // (
- // db().objectRegistry::lookupObject<quadraticFitSnGradData>
- // (
- // quadraticFitSnGradData::typeName
- // )
- // ).movePoints();
- //}
-
- // skewCorrectionVectors
- if
- (
- db().objectRegistry::foundObject<skewCorrectionVectors>
- (
- skewCorrectionVectors::typeName
- )
- )
- {
- const_cast<skewCorrectionVectors&>
- (
- db().objectRegistry::lookupObject<skewCorrectionVectors>
- (
- skewCorrectionVectors::typeName
- )
- ).movePoints();
- }
-
+ MeshObjectMovePoints<volPointInterpolation>(*this);
+ MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
+ MeshObjectMovePoints<leastSquaresVectors>(*this);
+ MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
+ MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
+ MeshObjectMovePoints<skewCorrectionVectors>(*this);
return tsweptVols;
}
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C
similarity index 55%
rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C
index 868d33b6f3f132b872d08d4c11a8464810f70a21..dea1b41143792cd55dc44e47ae0455d04eed511c 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.C
@@ -24,125 +24,78 @@ License
\*---------------------------------------------------------------------------*/
-#include "quadraticFitData.H"
+#include "CentredFitData.H"
#include "surfaceFields.H"
#include "volFields.H"
#include "SVD.H"
#include "syncTools.H"
#include "extendedCentredStencil.H"
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
- defineTypeNameAndDebug(quadraticFitData, 0);
-}
-
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-Foam::quadraticFitData::quadraticFitData
+template<class Polynomial>
+Foam::CentredFitData<Polynomial>::CentredFitData
(
const fvMesh& mesh,
const extendedCentredStencil& stencil,
const scalar linearLimitFactor,
- const scalar cWeight
+ const scalar centralWeight
)
:
- MeshObject<fvMesh, quadraticFitData>(mesh),
+ MeshObject<fvMesh, CentredFitData<Polynomial> >(mesh),
+ stencil_(stencil),
linearLimitFactor_(linearLimitFactor),
- centralWeight_(cWeight),
+ centralWeight_(centralWeight),
# ifdef SPHERICAL_GEOMETRY
dim_(2),
# else
dim_(mesh.nGeometricD()),
# endif
- minSize_
- (
- dim_ == 1 ? 3 :
- dim_ == 2 ? 6 :
- dim_ == 3 ? 7 : 0
- ),
- fit_(mesh.nInternalFaces())
+ minSize_(Polynomial::nTerms(dim_)),
+ coeffs_(mesh.nFaces())
{
if (debug)
{
- Info << "Contructing quadraticFitData" << endl;
+ Info<< "Contructing CentredFitData<Polynomial>" << endl;
}
- // check input
- if (centralWeight_ < 1 - SMALL)
+ // Check input
+ if (linearLimitFactor > 1)
{
- FatalErrorIn("quadraticFitData::quadraticFitData")
- << "centralWeight requested = " << centralWeight_
+ FatalErrorIn("CentredFitData<Polynomial>::CentredFitData")
+ << "linearLimitFactor requested = " << linearLimitFactor
<< " should not be less than one"
<< exit(FatalError);
}
- if (minSize_ == 0)
- {
- FatalErrorIn("quadraticFitSnGradData")
- << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
- }
-
- // store the polynomial size for each cell to write out
- surfaceScalarField interpPolySize
- (
- IOobject
- (
- "quadraticFitInterpPolySize",
- "constant",
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- mesh,
- dimensionedScalar("quadraticFitInterpPolySize", dimless, scalar(0))
- );
-
- // Get the cell/face centres in stencil order.
- // Centred face stencils no good for triangles of tets. Need bigger stencils
- List<List<point> > stencilPoints(mesh.nFaces());
- stencil.collectData
- (
- mesh.C(),
- stencilPoints
- );
-
- // find the fit coefficients for every face in the mesh
-
- for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
- {
- interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
- }
+ calcFit();
if (debug)
{
- Info<< "quadraticFitData::quadraticFitData() :"
+ Info<< "CentredFitData<Polynomial>::CentredFitData() :"
<< "Finished constructing polynomialFit data"
<< endl;
-
- interpPolySize.write();
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-void Foam::quadraticFitData::findFaceDirs
+template<class Polynomial>
+void Foam::CentredFitData<Polynomial>::findFaceDirs
(
vector& idir, // value changed in return
vector& jdir, // value changed in return
vector& kdir, // value changed in return
const fvMesh& mesh,
- const label faci
+ const label facei
)
{
- idir = mesh.Sf()[faci];
+ idir = mesh.faceAreas()[facei];
idir /= mag(idir);
# ifndef SPHERICAL_GEOMETRY
- if (mesh.nGeometricD() <= 2) // find the normal direcion
+ if (mesh.nGeometricD() <= 2) // find the normal direction
{
if (mesh.directions()[0] == -1)
{
@@ -157,14 +110,14 @@ void Foam::quadraticFitData::findFaceDirs
kdir = vector(0, 0, 1);
}
}
- else // 3D so find a direction in the place of the face
+ else // 3D so find a direction in the plane of the face
{
- const face& f = mesh.faces()[faci];
- kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
+ const face& f = mesh.faces()[facei];
+ kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
}
# else
// Spherical geometry so kdir is the radial direction
- kdir = mesh.Cf()[faci];
+ kdir = mesh.faceCentres()[facei];
# endif
if (mesh.nGeometricD() == 3)
@@ -189,94 +142,125 @@ void Foam::quadraticFitData::findFaceDirs
}
-Foam::label Foam::quadraticFitData::calcFit
+template<class Polynomial>
+void Foam::CentredFitData<Polynomial>::calcFit()
+{
+ const fvMesh& mesh = this->mesh();
+
+ // Get the cell/face centres in stencil order.
+ // Centred face stencils no good for triangles of tets.
+ // Need bigger stencils
+ List<List<point> > stencilPoints(mesh.nFaces());
+ stencil_.collectData(mesh.C(), stencilPoints);
+
+ // find the fit coefficients for every face in the mesh
+
+ const surfaceScalarField& w = this->mesh().surfaceInterpolation::weights();
+
+ for(label facei = 0; facei < mesh.nInternalFaces(); facei++)
+ {
+ calcFit(stencilPoints[facei], w[facei], facei);
+ }
+
+ const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField();
+
+ forAll(bw, patchi)
+ {
+ const fvsPatchScalarField& pw = bw[patchi];
+
+ if (pw.coupled())
+ {
+ label facei = pw.patch().patch().start();
+
+ forAll(pw, i)
+ {
+ calcFit(stencilPoints[facei], pw[i], facei);
+ facei++;
+ }
+ }
+ }
+}
+
+
+template<class Polynomial>
+Foam::label Foam::CentredFitData<Polynomial>::calcFit
(
const List<point>& C,
- const label faci
+ const scalar wLin,
+ const label facei
)
{
vector idir(1,0,0);
vector jdir(0,1,0);
vector kdir(0,0,1);
- findFaceDirs(idir, jdir, kdir, mesh(), faci);
+ findFaceDirs(idir, jdir, kdir, this->mesh(), facei);
+ // Setup the point weights
scalarList wts(C.size(), scalar(1));
wts[0] = centralWeight_;
wts[1] = centralWeight_;
- point p0 = mesh().faceCentres()[faci];
- scalar scale = 0;
+ // Reference point
+ point p0 = this->mesh().faceCentres()[facei];
+
+ // p0 -> p vector in the face-local coordinate system
+ vector d;
- // calculate the matrix of the polynomial components
+ // Local coordinate scaling
+ scalar scale = 1;
+
+ // Matrix of the polynomial components
scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
for(label ip = 0; ip < C.size(); ip++)
{
const point& p = C[ip];
- scalar px = (p - p0)&idir;
- scalar py = (p - p0)&jdir;
+ d.x() = (p - p0)&idir;
+ d.y() = (p - p0)&jdir;
# ifndef SPHERICAL_GEOMETRY
- scalar pz = (p - p0)&kdir;
+ d.z() = (p - p0)&kdir;
# else
- scalar pz = mag(p) - mag(p0);
+ d.z() = mag(p) - mag(p0);
# endif
if (ip == 0)
{
- scale = max(max(mag(px), mag(py)), mag(pz));
+ scale = cmptMax(cmptMag((d)));
}
- px /= scale;
- py /= scale;
- pz /= scale;
-
- label is = 0;
-
- B[ip][is++] = wts[0]*wts[ip];
- B[ip][is++] = wts[0]*wts[ip]*px;
- B[ip][is++] = wts[ip]*sqr(px);
+ // Scale the radius vector
+ d /= scale;
- if (dim_ >= 2)
- {
- B[ip][is++] = wts[ip]*py;
- B[ip][is++] = wts[ip]*px*py;
- //B[ip][is++] = wts[ip]*sqr(py);
- }
- if (dim_ == 3)
- {
- B[ip][is++] = wts[ip]*pz;
- B[ip][is++] = wts[ip]*px*pz;
- //B[ip][is++] = wts[ip]*py*pz;
- //B[ip][is++] = wts[ip]*sqr(pz);
- }
+ Polynomial::addCoeffs
+ (
+ B[ip],
+ d,
+ wts[ip],
+ dim_
+ );
}
// Set the fit
label stencilSize = C.size();
- fit_[faci].setSize(stencilSize);
+ coeffs_[facei].setSize(stencilSize);
scalarList singVals(minSize_);
label nSVDzeros = 0;
- const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
- mesh().surfaceInterpolation::weights();
-
bool goodFit = false;
for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
{
SVD svd(B, SMALL);
- scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
- scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
-
- //goodFit = (fit0 > 0 && fit1 > 0);
+ scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
+ scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
goodFit =
- (mag(fit0 - w[faci])/w[faci] < linearLimitFactor_)
- && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < linearLimitFactor_);
+ (mag(fit0 - wLin) < linearLimitFactor_*wLin)
+ && (mag(fit1 - (1 - wLin)) < linearLimitFactor_*(1 - wLin));
- //scalar w0Err = fit0/w[faci];
- //scalar w1Err = fit1/(1 - w[faci]);
+ //scalar w0Err = fit0/wLin;
+ //scalar w1Err = fit1/(1 - wLin);
//goodFit =
// (w0Err > linearLimitFactor_ && w0Err < (1 + linearLimitFactor_))
@@ -284,12 +268,12 @@ Foam::label Foam::quadraticFitData::calcFit
if (goodFit)
{
- fit_[faci][0] = fit0;
- fit_[faci][1] = fit1;
+ coeffs_[facei][0] = fit0;
+ coeffs_[facei][1] = fit1;
for(label i=2; i<stencilSize; i++)
{
- fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
+ coeffs_[facei][i] = wts[i]*svd.VSinvUt()[0][i];
}
singVals = svd.S();
@@ -300,12 +284,6 @@ Foam::label Foam::quadraticFitData::calcFit
wts[0] *= 10;
wts[1] *= 10;
- for(label i = 0; i < B.n(); i++)
- {
- B[i][0] *= 10;
- B[i][1] *= 10;
- }
-
for(label j = 0; j < B.m(); j++)
{
B[0][j] *= 10;
@@ -327,40 +305,54 @@ Foam::label Foam::quadraticFitData::calcFit
// (
// min
// (
- // min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
- // min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
+ // min(alpha - beta*mag(coeffs_[facei][0] - wLin)/wLin, 1),
+ // min(alpha - beta*mag(coeffs_[facei][1] - (1 - wLin))
+ // /(1 - wLin), 1)
// ), 0
// );
- // Remove the uncorrected linear coefficients
- fit_[faci][0] -= w[faci];
- fit_[faci][1] -= 1 - w[faci];
+ //Info<< wLin << " " << coeffs_[facei][0]
+ // << " " << (1 - wLin) << " " << coeffs_[facei][1] << endl;
+
+ // Remove the uncorrected linear ocefficients
+ coeffs_[facei][0] -= wLin;
+ coeffs_[facei][1] -= 1 - wLin;
// if (limiter < 0.99)
// {
// for(label i = 0; i < stencilSize; i++)
// {
- // fit_[faci][i] *= limiter;
+ // coeffs_[facei][i] *= limiter;
// }
// }
}
else
{
- Pout<< "Could not fit face " << faci
- << " " << fit_[faci][0] << " " << w[faci]
- << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
+ if (debug)
+ {
+ WarningIn
+ (
+ "CentredFitData<Polynomial>::calcFit"
+ "(const List<point>& C, const label facei"
+ ) << "Could not fit face " << facei
+ << ", reverting to linear." << nl
+ << " Weights "
+ << coeffs_[facei][0] << " " << wLin << nl
+ << " Linear weights "
+ << coeffs_[facei][1] << " " << 1 - wLin << endl;
+ }
- fit_[faci] = 0;
+ coeffs_[facei] = 0;
}
return minSize_ - nSVDzeros;
}
-bool Foam::quadraticFitData::movePoints()
+template<class Polynomial>
+bool Foam::CentredFitData<Polynomial>::movePoints()
{
- notImplemented("quadraticFitData::movePoints()");
-
+ calcFit();
return true;
}
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H
similarity index 75%
rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H
rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H
index f0da60901daf96db6f4ab3bcbe623eec8e09ad07..1630e06a927419c0747d33988a7d48ab48ec8483 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitData.H
@@ -23,18 +23,18 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
- Foam::quadraticFitData
+ Foam::CentredFitData
Description
Data for the quadratic fit correction interpolation scheme
SourceFiles
- quadraticFitData.C
+ CentredFitData.C
\*---------------------------------------------------------------------------*/
-#ifndef quadraticFitData_H
-#define quadraticFitData_H
+#ifndef CentredFitData_H
+#define CentredFitData_H
#include "MeshObject.H"
#include "fvMesh.H"
@@ -47,15 +47,19 @@ namespace Foam
class extendedCentredStencil;
/*---------------------------------------------------------------------------*\
- Class quadraticFitData Declaration
+ Class CentredFitData Declaration
\*---------------------------------------------------------------------------*/
-class quadraticFitData
+template<class Polynomial>
+class CentredFitData
:
- public MeshObject<fvMesh, quadraticFitData>
+ public MeshObject<fvMesh, CentredFitData<Polynomial> >
{
// Private data
+ //- The stencil the fit is based on
+ const extendedCentredStencil& stencil_;
+
//- Factor the fit is allowed to deviate from linear.
// This limits the amount of high-order correction and increases
// stability on bad meshes
@@ -72,7 +76,7 @@ class quadraticFitData
//- For each cell in the mesh store the values which multiply the
// values of the stencil to obtain the gradient for each direction
- List<scalarList> fit_;
+ List<scalarList> coeffs_;
// Private member functions
@@ -87,17 +91,28 @@ class quadraticFitData
const label faci
);
- label calcFit(const List<point>&, const label faci);
+ //- Calculate the fit for the all the mesh faces
+ // and set the coefficients
+ void calcFit();
+
+ //- Calculate the fit for the specified face and set the coefficients
+ label calcFit
+ (
+ const List<point>&, // Stencil points
+ const scalar wLin, // Linear weight
+ const label faci // Current face index
+ );
public:
- TypeName("quadraticFitData");
+ TypeName("CentredFitData");
// Constructors
- explicit quadraticFitData
+ //- Construct from components
+ CentredFitData
(
const fvMesh& mesh,
const extendedCentredStencil& stencil,
@@ -107,16 +122,16 @@ public:
//- Destructor
- virtual ~quadraticFitData()
+ virtual ~CentredFitData()
{}
// Member functions
//- Return reference to fit coefficients
- const List<scalarList>& fit() const
+ const List<scalarList>& coeffs() const
{
- return fit_;
+ return coeffs_;
}
//- Delete the data when the mesh moves not implemented
@@ -130,6 +145,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+# include "CentredFitData.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
#endif
// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H
similarity index 56%
rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H
index c1308c6ba3eb1c42bac70877201126aaa630e268..6c5ad7aa022eabf2767fed22e6208bb6e1fc3685 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CentredFitScheme/CentredFitScheme.H
@@ -23,23 +23,19 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
- Foam::quadraticFit
+ Foam::CentredFitScheme
Description
- Quadratic fit interpolation scheme which applies an explicit correction to
- linear.
-
-SourceFiles
- quadraticFit.C
+ Centred fit surface interpolation scheme which applies an explicit
+ correction to linear.
\*---------------------------------------------------------------------------*/
-#ifndef quadraticFit_H
-#define quadraticFit_H
+#ifndef CentredFitScheme_H
+#define CentredFitScheme_H
-#include "quadraticFitData.H"
+#include "CentredFitData.H"
#include "linear.H"
-#include "centredCFCStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -47,11 +43,11 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class quadraticFit Declaration
+ Class CentredFitScheme Declaration
\*---------------------------------------------------------------------------*/
-template<class Type>
-class quadraticFit
+template<class Type, class Polynomial, class Stencil>
+class CentredFitScheme
:
public linear<Type>
{
@@ -69,31 +65,31 @@ class quadraticFit
// Private Member Functions
//- Disallow default bitwise copy construct
- quadraticFit(const quadraticFit&);
+ CentredFitScheme(const CentredFitScheme&);
//- Disallow default bitwise assignment
- void operator=(const quadraticFit&);
+ void operator=(const CentredFitScheme&);
public:
//- Runtime type information
- TypeName("quadraticFit");
+ TypeName("CentredFitScheme");
// Constructors
//- Construct from mesh and Istream
- quadraticFit(const fvMesh& mesh, Istream& is)
+ CentredFitScheme(const fvMesh& mesh, Istream& is)
:
linear<Type>(mesh),
linearLimitFactor_(readScalar(is)),
- centralWeight_(readScalar(is))
+ centralWeight_(1000)
{}
//- Construct from mesh, faceFlux and Istream
- quadraticFit
+ CentredFitScheme
(
const fvMesh& mesh,
const surfaceScalarField& faceFlux,
@@ -102,7 +98,7 @@ public:
:
linear<Type>(mesh),
linearLimitFactor_(readScalar(is)),
- centralWeight_(readScalar(is))
+ centralWeight_(1000)
{}
@@ -123,13 +119,13 @@ public:
{
const fvMesh& mesh = this->mesh();
- const extendedCentredStencil& stencil =
- centredCFCStencilObject::New
+ const extendedCentredStencil& stencil = Stencil::New
(
mesh
);
- const quadraticFitData& cfd = quadraticFitData::New
+ const CentredFitData<Polynomial>& cfd =
+ CentredFitData<Polynomial>::New
(
mesh,
stencil,
@@ -137,7 +133,7 @@ public:
centralWeight_
);
- const List<scalarList>& f = cfd.fit();
+ const List<scalarList>& f = cfd.coeffs();
return stencil.weightedSum(vf, f);
}
@@ -148,6 +144,34 @@ public:
} // End namespace Foam
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Add the patch constructor functions to the hash tables
+
+#define makeCentredFitSurfaceInterpolationTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \
+ \
+typedef CentredFitScheme<TYPE, POLYNOMIAL, STENCIL> \
+ CentredFitScheme##TYPE##POLYNOMIAL##STENCIL##_; \
+defineTemplateTypeNameAndDebugWithName \
+ (CentredFitScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \
+ \
+surfaceInterpolationScheme<TYPE>::addMeshConstructorToTable \
+<CentredFitScheme<TYPE, POLYNOMIAL, STENCIL> > \
+ add##SS##STENCIL##TYPE##MeshConstructorToTable_; \
+ \
+surfaceInterpolationScheme<TYPE>::addMeshFluxConstructorToTable \
+<CentredFitScheme<TYPE, POLYNOMIAL, STENCIL> > \
+ add##SS##STENCIL##TYPE##MeshFluxConstructorToTable_;
+
+#define makeCentredFitSurfaceInterpolationScheme(SS, POLYNOMIAL, STENCIL) \
+ \
+makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \
+makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \
+makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \
+makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor)\
+makeCentredFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,tensor)
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C
index a7dead621f9dfa472a2ad370487a990511a72e2c..beff7881b270bbda13da67a64ac03e14ca6395bb 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C
@@ -24,13 +24,26 @@ License
\*---------------------------------------------------------------------------*/
-#include "linearFit.H"
+#include "CentredFitScheme.H"
+#include "linearFitPolynomial.H"
+#include "centredCFCStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
- makeSurfaceInterpolationScheme(linearFit);
+ defineTemplateTypeNameAndDebug
+ (
+ CentredFitData<linearFitPolynomial>,
+ 0
+ );
+
+ makeCentredFitSurfaceInterpolationScheme
+ (
+ linearFit,
+ linearFitPolynomial,
+ centredCFCStencilObject
+ );
}
// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H
deleted file mode 100644
index 882f3e94bcc2dc046c226062820b1144a1a4e5ff..0000000000000000000000000000000000000000
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H
+++ /dev/null
@@ -1,138 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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
-
-Class
- linearFit
-
-Description
- Linear fit interpolation scheme which applies an explicit correction to
- linear.
-
-SourceFiles
- linearFit.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef linearFit_H
-#define linearFit_H
-
-#include "linear.H"
-#include "linearFitData.H"
-#include "extendedStencil.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
- Class linearFit Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class linearFit
-:
- public linear<Type>
-{
- // Private Data
- const scalar centralWeight_;
-
- // Private Member Functions
-
- //- Disallow default bitwise copy construct
- linearFit(const linearFit&);
-
- //- Disallow default bitwise assignment
- void operator=(const linearFit&);
-
-
-public:
-
- //- Runtime type information
- TypeName("linearFit");
-
-
- // Constructors
-
- //- Construct from mesh and Istream
- linearFit(const fvMesh& mesh, Istream& is)
- :
- linear<Type>(mesh),
- centralWeight_(readScalar(is))
- {}
-
-
- //- Construct from mesh, faceFlux and Istream
- linearFit
- (
- const fvMesh& mesh,
- const surfaceScalarField& faceFlux,
- Istream& is
- )
- :
- linear<Type>(mesh),
- centralWeight_(readScalar(is))
- {}
-
-
- // Member Functions
-
- //- Return true if this scheme uses an explicit correction
- virtual bool corrected() const
- {
- return true;
- }
-
- //- Return the explicit correction to the face-interpolate
- virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
- correction
- (
- const GeometricField<Type, fvPatchField, volMesh>& vf
- ) const
- {
- const fvMesh& mesh = this->mesh();
-
- const linearFitData& cfd = linearFitData::New
- (
- mesh,
- centralWeight_
- );
-
- const extendedStencil& stencil = cfd.stencil();
- const List<scalarList>& f = cfd.fit();
-
- return stencil.weightedSum(vf, f);
- }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C
deleted file mode 100644
index 996d915e1cc56aad3734772bbf7be4e5667f0e90..0000000000000000000000000000000000000000
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C
+++ /dev/null
@@ -1,371 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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
-
-\*---------------------------------------------------------------------------*/
-
-#include "linearFitData.H"
-#include "surfaceFields.H"
-#include "volFields.H"
-#include "SVD.H"
-#include "syncTools.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
- defineTypeNameAndDebug(linearFitData, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-static int count = 0;
-
-Foam::linearFitData::linearFitData
-(
- const fvMesh& mesh,
- const scalar cWeight
-)
-:
- MeshObject<fvMesh, linearFitData>(mesh),
- centralWeight_(cWeight),
-# ifdef SPHERICAL_GEOMETRY
- dim_(2),
-# else
- dim_(mesh.nGeometricD()),
-# endif
- minSize_
- (
- dim_ == 1 ? 2 :
- dim_ == 2 ? 3 :
- dim_ == 3 ? 4 : 0
- ),
- stencil_(mesh),
- fit_(mesh.nInternalFaces())
-{
- if (debug)
- {
- Info << "Contructing linearFitData" << endl;
- }
-
- // check input
- if (centralWeight_ < 1 - SMALL)
- {
- FatalErrorIn("linearFitData::linearFitData")
- << "centralWeight requested = " << centralWeight_
- << " should not be less than one"
- << exit(FatalError);
- }
-
- if (minSize_ == 0)
- {
- FatalErrorIn("linearFitSnGradData")
- << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
- }
-
- // store the polynomial size for each cell to write out
- surfaceScalarField interpPolySize
- (
- IOobject
- (
- "linearFitInterpPolySize",
- "constant",
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- mesh,
- dimensionedScalar("linearFitInterpPolySize", dimless, scalar(0))
- );
-
- // Get the cell/face centres in stencil order.
- // Centred face stencils no good for triangles of tets. Need bigger stencils
- List<List<point> > stencilPoints(stencil_.stencil().size());
- stencil_.collectData
- (
- mesh.C(),
- stencilPoints
- );
-
- // find the fit coefficients for every face in the mesh
-
- for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
- {
- interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
- }
-
- Pout<< "count = " << count << endl;
-
- if (debug)
- {
- Info<< "linearFitData::linearFitData() :"
- << "Finished constructing polynomialFit data"
- << endl;
-
- interpPolySize.write();
- }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-void Foam::linearFitData::findFaceDirs
-(
- vector& idir, // value changed in return
- vector& jdir, // value changed in return
- vector& kdir, // value changed in return
- const fvMesh& mesh,
- const label faci
-)
-{
- idir = mesh.Sf()[faci];
- //idir = mesh.C()[mesh.neighbour()[faci]] - mesh.C()[mesh.owner()[faci]];
- idir /= mag(idir);
-
-# ifndef SPHERICAL_GEOMETRY
- if (mesh.nGeometricD() <= 2) // find the normal direcion
- {
- if (mesh.directions()[0] == -1)
- {
- kdir = vector(1, 0, 0);
- }
- else if (mesh.directions()[1] == -1)
- {
- kdir = vector(0, 1, 0);
- }
- else
- {
- kdir = vector(0, 0, 1);
- }
- }
- else // 3D so find a direction in the place of the face
- {
- const face& f = mesh.faces()[faci];
- kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
- }
-# else
- // Spherical geometry so kdir is the radial direction
- kdir = mesh.Cf()[faci];
-# endif
-
- if (mesh.nGeometricD() == 3)
- {
- // Remove the idir component from kdir and normalise
- kdir -= (idir & kdir)*idir;
-
- scalar magk = mag(kdir);
-
- if (magk < SMALL)
- {
- FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
- << exit(FatalError);
- }
- else
- {
- kdir /= magk;
- }
- }
-
- jdir = kdir ^ idir;
-}
-
-
-Foam::label Foam::linearFitData::calcFit
-(
- const List<point>& C,
- const label faci
-)
-{
- vector idir(1,0,0);
- vector jdir(0,1,0);
- vector kdir(0,0,1);
- findFaceDirs(idir, jdir, kdir, mesh(), faci);
-
- scalarList wts(C.size(), scalar(1));
- wts[0] = centralWeight_;
- wts[1] = centralWeight_;
-
- point p0 = mesh().faceCentres()[faci];
- scalar scale = 0;
-
- // calculate the matrix of the polynomial components
- scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
-
- for(label ip = 0; ip < C.size(); ip++)
- {
- const point& p = C[ip];
-
- scalar px = (p - p0)&idir;
- scalar py = (p - p0)&jdir;
-# ifndef SPHERICAL_GEOMETRY
- scalar pz = (p - p0)&kdir;
-# else
- scalar pz = mag(p) - mag(p0);
-# endif
-
- if (ip == 0)
- {
- scale = max(max(mag(px), mag(py)), mag(pz));
- }
-
- px /= scale;
- py /= scale;
- pz /= scale;
-
- label is = 0;
-
- B[ip][is++] = wts[0]*wts[ip];
- B[ip][is++] = wts[0]*wts[ip]*px;
-
- if (dim_ >= 2)
- {
- B[ip][is++] = wts[ip]*py;
- }
- if (dim_ == 3)
- {
- B[ip][is++] = wts[ip]*pz;
- }
- }
-
- // Set the fit
- label stencilSize = C.size();
- fit_[faci].setSize(stencilSize);
- scalarList singVals(minSize_);
- label nSVDzeros = 0;
-
- const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
- mesh().surfaceInterpolation::weights();
-
- bool goodFit = false;
- for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
- {
- SVD svd(B, SMALL);
-
- scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
- scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
-
- //goodFit = (fit0 > 0 && fit1 > 0);
-
- goodFit =
- (mag(fit0 - w[faci])/w[faci] < 0.5)
- && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5);
-
- //scalar w0Err = fit0/w[faci];
- //scalar w1Err = fit1/(1 - w[faci]);
-
- //goodFit =
- // (w0Err > 0.5 && w0Err < 1.5)
- // && (w1Err > 0.5 && w1Err < 1.5);
-
- if (goodFit)
- {
- fit_[faci][0] = fit0;
- fit_[faci][1] = fit1;
-
- for(label i=2; i<stencilSize; i++)
- {
- fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
- }
-
- singVals = svd.S();
- nSVDzeros = svd.nZeros();
- }
- else // (not good fit so increase weight in the centre and for linear)
- {
- wts[0] *= 10;
- wts[1] *= 10;
-
- for(label i = 0; i < B.n(); i++)
- {
- B[i][0] *= 10;
- B[i][1] *= 10;
- }
-
- for(label j = 0; j < B.m(); j++)
- {
- B[0][j] *= 10;
- B[1][j] *= 10;
- }
- }
- }
-
- // static const scalar L = 0.1;
- // static const scalar R = 0.2;
-
- // static const scalar beta = 1.0/(R - L);
- // static const scalar alpha = R*beta;
-
- if (goodFit)
- {
- if ((mag(fit_[faci][0] - w[faci])/w[faci] < 0.15)
- && (mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) < 0.15))
- {
- count++;
- //Pout<< "fit " << mag(fit_[faci][0] - w[faci])/w[faci] << " " << mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) << endl;
- }
-
- // scalar limiter =
- // max
- // (
- // min
- // (
- // min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
- // min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
- // ), 0
- // );
-
- // Remove the uncorrected linear coefficients
- fit_[faci][0] -= w[faci];
- fit_[faci][1] -= 1 - w[faci];
-
- // if (limiter < 0.99)
- // {
- // for(label i = 0; i < stencilSize; i++)
- // {
- // fit_[faci][i] *= limiter;
- // }
- // }
- }
- else
- {
- Pout<< "Could not fit face " << faci
- << " " << fit_[faci][0] << " " << w[faci]
- << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
- fit_[faci] = 0;
- }
-
- return minSize_ - nSVDzeros;
-}
-
-
-bool Foam::linearFitData::movePoints()
-{
- notImplemented("linearFitData::movePoints()");
-
- return true;
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.H
deleted file mode 100644
index 7a2bc7d28e31888f059dec2e35e862daf72fc95c..0000000000000000000000000000000000000000
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.H
+++ /dev/null
@@ -1,140 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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
-
-Class
- linearFitData
-
-Description
- Data for the linear fit correction interpolation scheme
-
-SourceFiles
- linearFitData.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef linearFitData_H
-#define linearFitData_H
-
-#include "MeshObject.H"
-#include "fvMesh.H"
-#include "extendedStencil.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class globalIndex;
-
-/*---------------------------------------------------------------------------*\
- Class linearFitData Declaration
-\*---------------------------------------------------------------------------*/
-
-class linearFitData
-:
- public MeshObject<fvMesh, linearFitData>
-{
- // Private data
-
- //- weights for central stencil
- const scalar centralWeight_;
-
- //- dimensionality of the geometry
- const label dim_;
-
- //- minimum stencil size
- const label minSize_;
-
- //- Extended stencil addressing
- extendedStencil stencil_;
-
- //- For each cell in the mesh store the values which multiply the
- // values of the stencil to obtain the gradient for each direction
- List<scalarList> fit_;
-
-
- // Private member functions
-
- //- Find the normal direction and i, j and k directions for face faci
- static void findFaceDirs
- (
- vector& idir, // value changed in return
- vector& jdir, // value changed in return
- vector& kdir, // value changed in return
- const fvMesh& mesh,
- const label faci
- );
-
- label calcFit(const List<point>&, const label faci);
-
-
-public:
-
- TypeName("linearFitData");
-
-
- // Constructors
-
- explicit linearFitData
- (
- const fvMesh& mesh,
- scalar cWeightDim
- );
-
-
- // Destructor
-
- virtual ~linearFitData()
- {}
-
-
- // Member functions
-
-
- //- Return reference to the stencil
- const extendedStencil& stencil() const
- {
- return stencil_;
- }
-
- //- Return reference to fit coefficients
- const List<scalarList>& fit() const
- {
- return fit_;
- }
-
- //- Delete the data when the mesh moves not implemented
- virtual bool movePoints();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H
new file mode 100644
index 0000000000000000000000000000000000000000..70275f9ee17fe008de2b6296c7842c8566fd9198
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitPolynomial.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2008 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
+
+Class
+ Foam::linearFitPolynomial
+
+Description
+ Linear polynomial for interpolation fitting.
+ Can be used with the CentredFit scheme to crate a linear surface
+ interpolation scheme
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearFitPolynomial_H
+#define linearFitPolynomial_H
+
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class linearFitPolynomial Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearFitPolynomial
+{
+public:
+
+ // Member functions
+
+ static label nTerms(const direction dim)
+ {
+ return
+ (
+ dim == 1 ? 2 :
+ dim == 2 ? 3 :
+ dim == 3 ? 4 : 0
+ );
+ }
+
+ static void addCoeffs
+ (
+ scalar* coeffs,
+ const vector& d,
+ const scalar weight,
+ const direction dim
+ )
+ {
+ register label curIdx = 0;
+
+ coeffs[curIdx++] = weight;
+ coeffs[curIdx++] = weight*d.x();
+
+ if (dim >= 2)
+ {
+ coeffs[curIdx++] = weight*d.y();
+ }
+ if (dim == 3)
+ {
+ coeffs[curIdx++] = weight*d.z();
+ }
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C
similarity index 78%
rename from src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
rename to src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C
index a215576330b96524c911ddbcc13a622c877fc3d2..0ecd5ece93bb10b829d2d4a3cc6da821e61592c6 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFit.C
@@ -24,13 +24,26 @@ License
\*---------------------------------------------------------------------------*/
-#include "quadraticFit.H"
+#include "CentredFitScheme.H"
+#include "quadraticLinearFitPolynomial.H"
+#include "centredCFCStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
- makeSurfaceInterpolationScheme(quadraticFit);
+ defineTemplateTypeNameAndDebug
+ (
+ CentredFitData<quadraticLinearFitPolynomial>,
+ 0
+ );
+
+ makeCentredFitSurfaceInterpolationScheme
+ (
+ quadraticLinearFit,
+ quadraticLinearFitPolynomial,
+ centredCFCStencilObject
+ );
}
// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H
new file mode 100644
index 0000000000000000000000000000000000000000..4250e71c1e211771aafc912a30bfcd156faa4cd6
--- /dev/null
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticLinearFit/quadraticLinearFitPolynomial.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2008 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
+
+Class
+ Foam::quadraticLinearFitPolynomial
+
+Description
+ Quadratic/linear polynomial for interpolation fitting:
+ quadratic normal to the face,
+ linear in the plane of the face for consistency with 2nd-order Gauss.
+
+ Can be used with the CentredFit scheme to crate a quadratic surface
+ interpolation scheme
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticLinearFitPolynomial_H
+#define quadraticLinearFitPolynomial_H
+
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class quadraticLinearFitPolynomial Declaration
+\*---------------------------------------------------------------------------*/
+
+class quadraticLinearFitPolynomial
+{
+public:
+
+ // Member functions
+
+ static label nTerms(const direction dim)
+ {
+ return
+ (
+ dim == 1 ? 3 :
+ dim == 2 ? 5 :
+ dim == 3 ? 7 : 0
+ );
+ }
+
+ static void addCoeffs
+ (
+ scalar* coeffs,
+ const vector& d,
+ const scalar weight,
+ const direction dim
+ )
+ {
+ register label curIdx = 0;
+
+ coeffs[curIdx++] = weight;
+ coeffs[curIdx++] = weight*d.x();
+ coeffs[curIdx++] = weight*sqr(d.x());
+
+ if (dim >= 2)
+ {
+ coeffs[curIdx++] = weight*d.y();
+ coeffs[curIdx++] = weight*d.x()*d.y();
+ }
+ if (dim == 3)
+ {
+ coeffs[curIdx++] = weight*d.z();
+ coeffs[curIdx++] = weight*d.x()*d.z();
+ }
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index a42408e9bcbdab58611442c28a92ad2f05a1d01c..adfe29401e4a1e08219e506538121b7a6a559838 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -2192,9 +2192,10 @@ bool Foam::distributedTriSurfaceMesh::writeObject
// Make sure dictionary goes to same directory as surface
const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
+ // Dictionary needs to be written in ascii - binary output not supported.
return
triSurfaceMesh::writeObject(fmt, ver, cmp)
- && dict_.writeObject(fmt, ver, cmp);
+ && dict_.writeObject(IOstream::ASCII, ver, cmp);
}