60 weightingFactors_(nullptr),
61 differenceFactors_(nullptr),
62 correctionVectors_(nullptr),
63 skewCorrectionVectors_(nullptr),
92 if (!weightingFactors_)
97 return (*weightingFactors_);
103 if (!differenceFactors_)
108 return (*differenceFactors_);
114 if (orthogonal_ ==
false && !correctionVectors_)
116 makeCorrectionVectors();
128 <<
"cannot return correctionVectors; mesh is orthogonal"
132 return (*correctionVectors_);
138 if (skew_ ==
true && !skewCorrectionVectors_)
140 makeSkewCorrectionVectors();
153 <<
"cannot return skewCorrectionVectors; mesh is now skewed"
157 return (*skewCorrectionVectors_);
177void Foam::edgeInterpolation::makeLPN()
const
180 <<
"Constructing geodesic distance between points P and N"
189 faMesh_.time().constant(),
222 - faceCentres[owner[edgeI]]
228 faceCentres[neighbour[edgeI]]
233 lPNIn[edgeI] = (lPE + lEN);
249 <<
"Finished constructing geodesic distance PN"
254void Foam::edgeInterpolation::makeWeights()
const
257 <<
"Constructing weighting factors for edge interpolation"
266 mesh()().pointsInstance(),
284 scalarField& weightingFactorsIn = weightingFactors.primitiveFieldRef();
292 curSkewCorrVec = skewCorrectionVectors()[edgeI];
300 - faceCentres[owner[edgeI]]
306 faceCentres[neighbour[edgeI]]
311 weightingFactorsIn[edgeI] =
315# ifdef BAD_MESH_STABILISATION
326 weightingFactors.boundaryFieldRef()[patchI]
331 <<
"Finished constructing weighting factors for face interpolation"
336void Foam::edgeInterpolation::makeDeltaCoeffs()
const
339 <<
"Constructing differencing factors array for edge gradient"
350 "differenceFactors_",
351 mesh()().pointsInstance(),
382 faceCentres[neighbour[edgeI]]
383 - faceCentres[owner[edgeI]];
385 unitDelta.removeCollinear(edgeNormal);
386 unitDelta.normalise();
394 curSkewCorrVec = skewCorrectionVectors()[edgeI];
402 - faceCentres[owner[edgeI]]
408 faceCentres[neighbour[edgeI]]
413 scalar lPN = lPE + lEN;
423 dc[edgeI] = 1.0/
max((lPN*(unitDelta & edgeNormal)), 0.05*lPN);
427 forAll(DeltaCoeffs.boundaryField(), patchI)
431 DeltaCoeffs.boundaryFieldRef()[patchI]
437void Foam::edgeInterpolation::makeCorrectionVectors()
const
440 <<
"Constructing non-orthogonal correction vectors"
448 mesh()().pointsInstance(),
472 vectorField& CorrVecsIn = CorrVecs.primitiveFieldRef();
482 faceCentres[neighbour[edgeI]]
483 - faceCentres[owner[edgeI]];
485 unitDelta.removeCollinear(edgeNormal);
486 unitDelta.normalise();
492 deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);
497 - deltaCoeffs[edgeI]*unitDelta;
503 forAll(CorrVecs.boundaryField(), patchI)
505 mesh().
boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]);
508 scalar NonOrthogCoeff = 0.0;
510 if (owner.size() > 0)
516 sinAlpha[edgeI] =
max(-1,
min(sinAlpha[edgeI], 1));
522 reduce(NonOrthogCoeff, maxOp<scalar>());
525 <<
"non-orthogonality coefficient = " << NonOrthogCoeff <<
" deg."
528 if (NonOrthogCoeff < 0.1)
539 <<
"Finished constructing non-orthogonal correction vectors"
544void Foam::edgeInterpolation::makeSkewCorrectionVectors()
const
547 <<
"Constructing skew correction vectors"
554 "skewCorrectionVectors",
555 mesh()().pointsInstance(),
579 const vector& P =
C[owner[edgeI]];
580 const vector&
N =
C[neighbour[edgeI]];
585 -(((
N - P)^(S - P))&((
N - P)^
e ))/(((
N - P)^
e)&((
N - P)^
e));
589 SkewCorrVecs[edgeI] = Ce[edgeI] - E;
594 SkewCorrVecs.boundaryFieldRef();
596 forAll(SkewCorrVecs.boundaryField(), patchI)
600 if (patchSkewCorrVecs.coupled())
610 forAll(patchSkewCorrVecs, edgeI)
612 const vector& P =
C[edgeFaces[edgeI]];
618 - (((
N - P)^(S - P))&((
N - P)^
e))
619 /(((
N - P)^
e)&((
N - P)^
e));
623 patchSkewCorrVecs[edgeI] =
624 Ce.boundaryField()[patchI][edgeI] - E;
629 patchSkewCorrVecs = vector::zero;
634 scalar skewCoeff = 0.0;
645 - SkewCorrVecs[edgeI]
652 + SkewCorrVecs[edgeI]
658 skewCoeff =
max(
mag(SkewCorrVecs.internalField())/
mag(lPN));
662 <<
"skew coefficient = " << skewCoeff <<
endl;
675 <<
"Finished constructing skew correction vectors"
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of boundary fields.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
SubList< edge > subList
Declare type of subList.
Face to edge interpolation scheme. Included in faMesh.
bool movePoints() const
Do what is necessary if the mesh has moved.
const edgeScalarField & deltaCoeffs() const
Return reference to difference factors array.
const edgeVectorField & correctionVectors() const
Return reference to non-orthogonality correction vectors array.
bool orthogonal() const
Return whether mesh is orthogonal or not.
const edgeScalarField & lPN() const
Return reference to PN geodesic distance.
const edgeScalarField & weights() const
Return reference to weighting factors array.
~edgeInterpolation()
Destructor.
bool skew() const
Return whether mesh is skew or not.
void clearOut()
Clear all geometry and addressing.
const edgeVectorField & skewCorrectionVectors() const
Return reference to skew vectors array.
Finite area mesh (used for 2-D non-Euclidian finite area method) defined using a patch of faces on a ...
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
const labelUList & owner() const
Internal face owner. Note bypassing virtual mechanism so.
const labelUList & neighbour() const
Internal face neighbour.
virtual const pointField & points() const
Return raw points.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Skew-correction vectors for the skewness-corrected interpolation scheme.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define DebugInFunction
Report an information message using Foam::Info.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar asin(const dimensionedScalar &ds)
GeometricField< vector, faePatchField, edgeMesh > edgeVectorField
const dimensionSet dimless
Dimensionless.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
vectorField pointField
pointField is a vectorField.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
GeometricField< vector, faPatchField, areaMesh > areaVectorField
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
List< edge > edgeList
A List of edges.
faePatchField< vector > faePatchVectorField
UList< label > labelUList
A UList of labels.
void deleteDemandDrivenData(DataPtr &dataPtr)
dimensionedTensor skew(const dimensionedTensor &dt)
Calculate the matrix for the second temporal derivative.
#define forAll(list, i)
Loop across all elements in list.
const Vector< label > N(dict.get< Vector< label > >("N"))