40 const Foam::scalar Foam::twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
45 void Foam::twoDPointCorrector::calcAddressing()
const 48 planeNormalPtr_ =
new vector(0, 0, 0);
49 vector& pn = *planeNormalPtr_;
58 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
62 if (isA<wedgePolyPatch>(patches[patchi]))
66 const wedgePolyPatch& wp =
67 refCast<const wedgePolyPatch>(patches[
patchi]);
69 pn = wp.centreNormal();
71 wedgeAxis_ = wp.axis();
72 wedgeAngle_ =
mag(
acos(wp.cosAngle()));
76 Pout<<
"Found normal from wedge patch " <<
patchi;
88 if (isA<emptyPolyPatch>(patches[patchi]) && patches[patchi].size())
90 pn = patches[
patchi].faceAreas()[0];
94 Pout<<
"Found normal from empty patch " <<
patchi;
103 if (
mag(pn) < vSmall)
106 <<
"Cannot determine normal vector from patches." 116 Pout<<
" twoDPointCorrector normal: " << pn <<
endl;
120 normalEdgeIndicesPtr_ =
new labelList(mesh_.nEdges());
121 labelList& neIndices = *normalEdgeIndicesPtr_;
123 const edgeList& meshEdges = mesh_.edges();
125 const pointField& meshPoints = mesh_.points();
127 label nNormalEdges = 0;
131 const edge& e = meshEdges[edgeI];
133 vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + vSmall);
135 if (
mag(edgeVector & pn) > edgeOrthogonalityTol)
138 neIndices[nNormalEdges++] = edgeI;
142 neIndices.setSize(nNormalEdges);
149 if (meshPoints.size() % 2 != 0)
152 <<
"the number of vertices in the geometry " 153 <<
"is odd - this should not be the case for a 2-D case. " 154 <<
"Please check the geometry." 158 if (2*nNormalEdges != meshPoints.size())
161 <<
"The number of points in the mesh is " 162 <<
"not equal to twice the number of edges normal to the plane " 163 <<
"- this may be OK only for wedge geometries.\n" 164 <<
" Please check the geometry or adjust " 165 <<
"the orthogonality tolerance.\n" <<
endl 166 <<
"Number of normal edges: " << nNormalEdges
167 <<
" number of points: " << meshPoints.size()
174 void Foam::twoDPointCorrector::clearAddressing()
const 181 void Foam::twoDPointCorrector::snapToWedge
188 scalar ADash =
mag(A - wedgeAxis_*(wedgeAxis_ & A));
189 vector pDash = ADash*
tan(wedgeAngle_)*planeNormal();
191 p = A +
sign(n & p)*pDash;
200 required_(mesh_.nGeometricD() == 2),
201 planeNormalPtr_(nullptr),
202 normalEdgeIndicesPtr_(nullptr),
224 if (
mag(pn.
x()) >= edgeOrthogonalityTol)
228 else if (
mag(pn.
y()) >= edgeOrthogonalityTol)
232 else if (
mag(pn.
z()) >= edgeOrthogonalityTol)
239 <<
"Plane normal not aligned with the coordinate system" <<
nl 250 if (!planeNormalPtr_)
255 return *planeNormalPtr_;
261 if (!normalEdgeIndicesPtr_)
266 return *normalEdgeIndicesPtr_;
272 if (!required_)
return;
287 point& pStart = p[meshEdges[neIndices[edgeI]].start()];
289 point& pEnd = p[meshEdges[neIndices[edgeI]].
end()];
292 point A = 0.5*(pStart + pEnd);
297 snapToWedge(pn, A, pStart);
298 snapToWedge(pn, A, pEnd);
303 pStart = A + pn*(pn & (pStart - A));
304 pEnd = A + pn*(pn & (pEnd - A));
316 if (!required_)
return;
331 const edge& e = meshEdges[neIndices[edgeI]];
334 point pStart = p[startPointi] + disp[startPointi];
337 point pEnd = p[endPointi] + disp[endPointi];
340 point A = 0.5*(pStart + pEnd);
345 snapToWedge(pn, A, pStart);
346 snapToWedge(pn, A, pEnd);
347 disp[startPointi] = pStart - p[startPointi];
348 disp[endPointi] = pEnd - p[endPointi];
353 disp[startPointi] = (A + pn*(pn & (pStart - A))) - p[startPointi];
354 disp[endPointi] = (A + pn*(pn & (pEnd - A))) - p[endPointi];
dimensionedScalar sign(const dimensionedScalar &ds)
void updateMesh(const mapPolyMesh &)
Update topology.
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList & normalEdgeIndices() const
Return indices of normal edges.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
iterator end()
Return an iterator to end traversing the UList.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
void correctPoints(pointField &p) const
Correct motion points.
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
void correctDisplacement(const pointField &p, vectorField &disp) const
Correct motion displacements.
vectorField pointField
pointField is a vectorField.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Class applies a two-dimensional correction to mesh motion point field.
List< label > labelList
A List of labels.
errorManip< error > abort(error &err)
const vector & planeNormal() const
Return plane normal.
defineTypeNameAndDebug(combustionModel, 0)
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
bool movePoints()
Correct weighting factors for moving mesh.
label end() const
Return end vertex label.
prefixOSstream Pout(cout, "Pout")
direction normalDir() const
Return direction normal to plane.
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
~twoDPointCorrector()
Destructor.
dimensionedScalar tan(const dimensionedScalar &ds)
void deleteDemandDrivenData(DataPtr &dataPtr)
label start() const
Return start vertex label.