diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
index a4bdba28e5e27253b5417f0b8b7f53de2f60ec2d..d7e0dc74ceb8d52673f3a55f8106a7eb175181b0 100644
--- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
@@ -79,9 +79,15 @@ protected:
//- Make patch weighting factors
virtual void makeWeights(scalarField&) const = 0;
+ //- Make patch geodesic distance between P and N
+ virtual void makeLPN(scalarField&) const = 0;
+
//- Make patch face - neighbour cell distances
virtual void makeDeltaCoeffs(scalarField&) const = 0;
+ //- Make non-orthogonality correction vectors
+ virtual void makeCorrectionVectors(vectorField&) const = 0;
+
//- Calculate the uniform transformation tensors
void calcTransformTensors
(
diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
index 5d0ffcb436fe3931b4b9406c2f6c9eca5d3aa56b..426cc281583f67b08d410502e401382eceeb3e9b 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
@@ -218,6 +218,13 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
}
+void Foam::cyclicFaPatch::makeLPN(scalarField& lPN) const
+{
+ makeDeltaCoeffs(lPN);
+ lPN = scalar(1)/lPN;
+}
+
+
void Foam::cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const
{
const scalarField deltas(edgeNormals() & coupledFaPatch::delta());
diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H
index 9e4059fdff423da2eaa7cc2685d6755653006cde..e9c559353ae06eb43d6de745d125e3ecbefec807 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H
@@ -78,9 +78,15 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
+ //- Make patch geodesic distance between P and N
+ void makeLPN(scalarField&) const;
+
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
+ //- Make non-orthogonality correction vectors
+ void makeCorrectionVectors(vectorField& cv) const { cv = Zero; }
+
public:
diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
index c60a08a0ba62a8f41d448e5badbf6e170a90161a..d38a88c165edb5a2281a9f27a4bc00e431c5688b 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
@@ -398,18 +398,21 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
{
if (Pstream::parRun())
{
- // The face normals point in the opposite direction on the other side
- scalarField neighbEdgeCentresCn
+ const vectorField& skew = skewCorrectionVectors();
+ const scalarField& PN = lPN();
+
+ const scalarField lEN
(
- neighbEdgeNormals()
- & (neighbEdgeCentres() - neighbEdgeFaceCentres())
+ mag(neighbEdgeFaceCentres() - edgeCentres() + skew)
);
- w = neighbEdgeCentresCn/
- (
- (edgeNormals() & coupledFaPatch::delta())
- + neighbEdgeCentresCn + VSMALL
- );
+ forAll(w, i)
+ {
+ if (mag(PN[i]) > SMALL)
+ {
+ w[i] = lEN[i]/PN[i];
+ }
+ }
}
else
{
@@ -418,12 +421,61 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
}
+void Foam::processorFaPatch::makeLPN(scalarField& lPN) const
+{
+ const vectorField& skew = skewCorrectionVectors();
+
+ if (Pstream::parRun())
+ {
+ lPN =
+ mag(edgeCentres() - skew - edgeFaceCentres()) // lPE
+ + mag(neighbEdgeFaceCentres() - edgeCentres() + skew); // lEN
+
+ for (scalar& lpn: lPN)
+ {
+ if (mag(lpn) < SMALL)
+ {
+ lpn = SMALL;
+ }
+ }
+ }
+ else
+ {
+ lPN =
+ mag(edgeCentres() - skew - edgeFaceCentres()) // lPE
+ + mag(edgeFaceCentres() - edgeCentres() + skew); // lEN
+ }
+}
+
+
void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const
{
if (Pstream::parRun())
{
- dc = (1.0 - weights())
- /((edgeNormals() & coupledFaPatch::delta()) + VSMALL);
+ const edgeList& edges = boundaryMesh().mesh().edges();
+ const pointField& points = boundaryMesh().mesh().points();
+ const vectorField& lengths = edgeLengths();
+ const scalarField& PN = lPN();
+
+ tmp<vectorField> tdelta = processorFaPatch::delta();
+ vectorField& unitDelta = tdelta.ref();
+
+ forAll(dc, i)
+ {
+ vector edgeNormal =
+ normalised(lengths[i] ^ edges[i + start()].vec(points));
+
+ unitDelta[i].removeCollinear(edgeNormal);
+ unitDelta[i].normalise();
+
+ edgeNormal = normalised(lengths[i]);
+
+ const scalar alpha = PN[i]*(edgeNormal & unitDelta[i]);
+ if (mag(alpha) > SMALL)
+ {
+ dc[i] = scalar(1)/max(alpha, 0.05*PN[i]);
+ }
+ }
}
else
{
@@ -433,6 +485,44 @@ void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const
}
+void Foam::processorFaPatch::makeCorrectionVectors(vectorField& cv) const
+{
+ if (Pstream::parRun())
+ {
+ const edgeList& edges = boundaryMesh().mesh().edges();
+ const pointField& points = boundaryMesh().mesh().points();
+ const vectorField& lengths = edgeLengths();
+
+ tmp<vectorField> tdelta = processorFaPatch::delta();
+ vectorField& unitDelta = tdelta.ref();
+
+ forAll(cv, i)
+ {
+ vector edgeNormal =
+ normalised(lengths[i] ^ edges[i + start()].vec(points));
+
+ unitDelta[i].removeCollinear(edgeNormal);
+ unitDelta.normalise();
+
+ edgeNormal = normalised(lengths[i]);
+
+ const scalar alpha = unitDelta[i] & edgeNormal;
+ scalar dc = SMALL;
+ if (mag(alpha) > SMALL)
+ {
+ dc = scalar(1)/alpha;
+ }
+
+ cv[i] = edgeNormal - dc*unitDelta[i];
+ }
+ }
+ else
+ {
+ cv = Zero;
+ }
+}
+
+
Foam::tmp<Foam::vectorField> Foam::processorFaPatch::delta() const
{
if (Pstream::parRun())
diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
index 6bcf0b05f9c24e0d2efc953e6c94d2c9b5e6dd7c..7eeb4eb3aa8cd9c2300584cf3ae0ba7d44f39cbd 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
@@ -105,9 +105,15 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
+ //- Make patch geodesic distance between P and N
+ void makeLPN(scalarField&) const;
+
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
+ //- Make non-orthogonality correction vectors
+ void makeCorrectionVectors(vectorField&) const;
+
//- Find non-globa patch points
void makeNonGlobalPatchPoints() const;
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index 735b9b2224ba3ad7bb9ed176939857ab71803488..9b103bd8728c07b56cfcb58d6a4f6db29da44738 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -500,12 +500,30 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
}
+void Foam::faPatch::makeLPN(scalarField& lPN) const
+{
+ lPN = (edgeNormals() & delta());
+}
+
+
+const Foam::scalarField& Foam::faPatch::lPN() const
+{
+ return boundaryMesh().mesh().lPN().boundaryField()[index()];
+}
+
+
void Foam::faPatch::makeDeltaCoeffs(scalarField& dc) const
{
dc = scalar(1)/(edgeNormals() & delta());
}
+const Foam::scalarField& Foam::faPatch::deltaCoeffs() const
+{
+ return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()];
+}
+
+
void Foam::faPatch::makeCorrectionVectors(vectorField& k) const
{
const vectorField unitDelta(delta()/mag(delta()));
@@ -514,9 +532,10 @@ void Foam::faPatch::makeCorrectionVectors(vectorField& k) const
}
-const Foam::scalarField& Foam::faPatch::deltaCoeffs() const
+const Foam::vectorField& Foam::faPatch::skewCorrectionVectors() const
{
- return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()];
+ return
+ boundaryMesh().mesh().skewCorrectionVectors().boundaryField()[index()];
}
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
index 448b82d1e0a1057a8108c8e19da612553a61aa1e..6bb87df10f1aa728d102834c86975d934fd33bd1 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
@@ -413,14 +413,24 @@ public:
//- Return patch weighting factors
const scalarField& weights() const;
+ //- Make patch geodesic distance between P and N
+ virtual void makeLPN(scalarField&) const;
+
+ //- Return patch geodesic distance between P and N
+ const scalarField& lPN() const;
+
//- Make patch edge - neighbour face distances
virtual void makeDeltaCoeffs(scalarField&) const;
- void makeCorrectionVectors(vectorField&) const;
-
//- Return patch edge - neighbour face distances
const scalarField& deltaCoeffs() const;
+ //- Make non-orthogonality correction vectors
+ virtual void makeCorrectionVectors(vectorField&) const;
+
+ //- Return patch skew-correction vectors
+ const vectorField& skewCorrectionVectors() const;
+
// Topological changes
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
index 5da45321ce37f7a182ec8c9334cdb16b2d9560ae..625a8a147111e06eb7b027a7767c44e974ec234b 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
@@ -98,9 +98,17 @@ const Foam::edgeVectorField& Foam::edgeInterpolation::correctionVectors() const
{
if (orthogonal())
{
- FatalErrorInFunction
- << "cannot return correctionVectors; mesh is orthogonal"
- << abort(FatalError);
+ return tmp<edgeVectorField>::New
+ (
+ IOobject
+ (
+ "correctionVectors",
+ mesh().pointsInstance(),
+ mesh().thisDb()
+ ),
+ mesh(),
+ dimensionedVector(dimless, Zero)
+ );
}
return (*correctionVectorsPtr_);
@@ -123,9 +131,17 @@ Foam::edgeInterpolation::skewCorrectionVectors() const
{
if (!skew())
{
- FatalErrorInFunction
- << "cannot return skewCorrectionVectors; mesh is now skewed"
- << abort(FatalError);
+ return tmp<edgeVectorField>::New
+ (
+ IOobject
+ (
+ "skewCorrectionVectors",
+ mesh().pointsInstance(),
+ mesh().thisDb()
+ ),
+ mesh(),
+ dimensionedVector(dimless, Zero)
+ );
}
return (*skewCorrectionVectorsPtr_);
@@ -233,12 +249,10 @@ void Foam::edgeInterpolation::makeLPN() const
forAll(lPN.boundaryField(), patchI)
{
- mesh().boundary()[patchI].makeDeltaCoeffs
+ mesh().boundary()[patchI].makeLPN
(
lPN.boundaryFieldRef()[patchI]
);
-
- lPN.boundaryFieldRef()[patchI] = 1.0/lPN.boundaryField()[patchI];
}
@@ -305,6 +319,7 @@ void Foam::edgeInterpolation::makeWeights() const
// weight = (0,1]
const scalar lPN = lPE + lEN;
+
if (mag(lPN) > SMALL)
{
weightingFactorsIn[edgeI] = lEN/lPN;
@@ -498,31 +513,6 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
mesh().boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]);
}
- scalar NonOrthogCoeff = 0.0;
-
- if (owner.size() > 0)
- {
- scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField()));
- sinAlpha.clamp_range(-1, 1);
-
- NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
- }
-
- reduce(NonOrthogCoeff, maxOp<scalar>());
-
- DebugInFunction
- << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
- << endl;
-
- if (NonOrthogCoeff < 0.1)
- {
- orthogonal_ = true;
- correctionVectorsPtr_.reset(nullptr);
- }
- else
- {
- orthogonal_ = false;
- }
DebugInFunction
<< "Finished constructing non-orthogonal correction vectors"
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes
index 8903be44f0c76a43ef323b0662892d634d32a298..1810f4b3ce6a461e45cac269df01f157f3aa0c49 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes
@@ -21,7 +21,7 @@ ddtSchemes
gradSchemes
{
- default none;
+ default Gauss linear;
grad(jouleHeatingSource:V_ceilingShell) Gauss linear;
}