twoDPointCorrector.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "twoDPointCorrector.H"
27 #include "polyMesh.H"
28 #include "wedgePolyPatch.H"
29 #include "emptyPolyPatch.H"
30 #include "SubField.H"
31 #include "meshTools.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
38 }
39 
40 const Foam::scalar Foam::twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 void Foam::twoDPointCorrector::calcAddressing() const
46 {
47  // Find geometry normal
48  planeNormalPtr_ = new vector(0, 0, 0);
49  vector& pn = *planeNormalPtr_;
50 
51  // Algorithm:
52  // Attempt to find wedge patch and work out the normal from it.
53  // If not found, find an empty patch with faces in it and use the
54  // normal of the first face. If neither is found, declare an
55  // error.
56 
57  // Try and find a wedge patch
58  const polyBoundaryMesh& patches = mesh().boundaryMesh();
59 
61  {
62  if (isA<wedgePolyPatch>(patches[patchi]))
63  {
64  isWedge_ = true;
65 
66  const wedgePolyPatch& wp =
67  refCast<const wedgePolyPatch>(patches[patchi]);
68 
69  pn = wp.centreNormal();
70 
71  wedgeAxis_ = wp.axis();
72  wedgeAngle_ = mag(acos(wp.cosAngle()));
73 
74  if (polyMesh::debug)
75  {
76  Pout<< "Found normal from wedge patch " << patchi;
77  }
78 
79  break;
80  }
81  }
82 
83  // Try to find an empty patch with faces
84  if (!isWedge_)
85  {
87  {
88  if (isA<emptyPolyPatch>(patches[patchi]) && patches[patchi].size())
89  {
90  pn = patches[patchi].faceAreas()[0];
91 
92  if (polyMesh::debug)
93  {
94  Pout<< "Found normal from empty patch " << patchi;
95  }
96 
97  break;
98  }
99  }
100  }
101 
102 
103  if (mag(pn) < vSmall)
104  {
106  << "Cannot determine normal vector from patches."
107  << abort(FatalError);
108  }
109  else
110  {
111  pn /= mag(pn);
112  }
113 
114  if (polyMesh::debug)
115  {
116  Pout<< " twoDPointCorrector normal: " << pn << endl;
117  }
118 
119  // Select edges to be included in check.
120  normalEdgeIndicesPtr_ = new labelList(mesh().nEdges());
121  labelList& neIndices = *normalEdgeIndicesPtr_;
122 
123  const edgeList& meshEdges = mesh().edges();
124 
125  const pointField& meshPoints = mesh().points();
126 
127  label nNormalEdges = 0;
128 
129  forAll(meshEdges, edgeI)
130  {
131  const edge& e = meshEdges[edgeI];
132 
133  vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + vSmall);
134 
135  if (mag(edgeVector & pn) > edgeOrthogonalityTol)
136  {
137  // this edge is normal to the plane. Add it to the list
138  neIndices[nNormalEdges++] = edgeI;
139  }
140  }
141 
142  neIndices.setSize(nNormalEdges);
143 
144  // Construction check: number of points in a read 2-D or wedge geometry
145  // should be odd and the number of edges normal to the plane should be
146  // exactly half the number of points
147  if (!isWedge_)
148  {
149  if (meshPoints.size() % 2 != 0)
150  {
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."
155  << endl;
156  }
157 
158  if (2*nNormalEdges != meshPoints.size())
159  {
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()
168  << endl;
169  }
170  }
171 }
172 
173 
174 void Foam::twoDPointCorrector::clearAddressing() const
175 {
176  deleteDemandDrivenData(planeNormalPtr_);
177  deleteDemandDrivenData(normalEdgeIndicesPtr_);
178 }
179 
180 
181 void Foam::twoDPointCorrector::snapToWedge
182 (
183  const vector& n,
184  const point& A,
185  point& p
186 ) const
187 {
188  scalar ADash = mag(A - wedgeAxis_*(wedgeAxis_ & A));
189  vector pDash = ADash*tan(wedgeAngle_)*planeNormal();
190 
191  p = A + sign(n & p)*pDash;
192 }
193 
194 
195 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
196 
198 :
200  <
201  polyMesh,
204  >(mesh),
205  required_(mesh.nGeometricD() == 2),
206  planeNormalPtr_(nullptr),
207  normalEdgeIndicesPtr_(nullptr),
208  isWedge_(false),
209  wedgeAxis_(Zero),
210  wedgeAngle_(0.0)
211 {}
212 
213 
214 
215 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
216 
218 {
219  clearAddressing();
220 }
221 
222 
223 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
224 
226 {
227  const vector& pn = planeNormal();
228 
229  if (mag(pn.x()) >= edgeOrthogonalityTol)
230  {
231  return vector::X;
232  }
233  else if (mag(pn.y()) >= edgeOrthogonalityTol)
234  {
235  return vector::Y;
236  }
237  else if (mag(pn.z()) >= edgeOrthogonalityTol)
238  {
239  return vector::Z;
240  }
241  else
242  {
244  << "Plane normal not aligned with the coordinate system" << nl
245  << " pn = " << pn
246  << abort(FatalError);
247 
248  return vector::Z;
249  }
250 }
251 
252 
254 {
255  if (!planeNormalPtr_)
256  {
257  calcAddressing();
258  }
259 
260  return *planeNormalPtr_;
261 }
262 
263 
265 {
266  if (!normalEdgeIndicesPtr_)
267  {
268  calcAddressing();
269  }
270 
271  return *normalEdgeIndicesPtr_;
272 }
273 
274 
276 {
277  if (!required_) return;
278 
279  // Algorithm:
280  // Loop through all edges. Calculate the average point position A for
281  // the front and the back. Correct the position of point P (in two planes)
282  // such that vectors AP and planeNormal are parallel
283 
284  // Get reference to edges
285  const edgeList& meshEdges = mesh().edges();
286 
287  const labelList& neIndices = normalEdgeIndices();
288  const vector& pn = planeNormal();
289 
290  forAll(neIndices, edgeI)
291  {
292  point& pStart = p[meshEdges[neIndices[edgeI]].start()];
293 
294  point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
295 
296  // calculate average point position
297  point A = 0.5*(pStart + pEnd);
299 
300  if (isWedge_)
301  {
302  snapToWedge(pn, A, pStart);
303  snapToWedge(pn, A, pEnd);
304  }
305  else
306  {
307  // correct point locations
308  pStart = A + pn*(pn & (pStart - A));
309  pEnd = A + pn*(pn & (pEnd - A));
310  }
311  }
312 }
313 
314 
316 (
317  const pointField& p,
318  vectorField& disp
319 ) const
320 {
321  if (!required_) return;
322 
323  // Algorithm:
324  // Loop through all edges. Calculate the average point position A for
325  // the front and the back. Correct the position of point P (in two planes)
326  // such that vectors AP and planeNormal are parallel
327 
328  // Get reference to edges
329  const edgeList& meshEdges = mesh().edges();
330 
331  const labelList& neIndices = normalEdgeIndices();
332  const vector& pn = planeNormal();
333 
334  forAll(neIndices, edgeI)
335  {
336  const edge& e = meshEdges[neIndices[edgeI]];
337 
338  label startPointi = e.start();
339  point pStart = p[startPointi] + disp[startPointi];
340 
341  label endPointi = e.end();
342  point pEnd = p[endPointi] + disp[endPointi];
343 
344  // calculate average point position
345  point A = 0.5*(pStart + pEnd);
347 
348  if (isWedge_)
349  {
350  snapToWedge(pn, A, pStart);
351  snapToWedge(pn, A, pEnd);
352  disp[startPointi] = pStart - p[startPointi];
353  disp[endPointi] = pEnd - p[endPointi];
354  }
355  else
356  {
357  // correct point locations
358  disp[startPointi] = (A + pn*(pn & (pStart - A))) - p[startPointi];
359  disp[endPointi] = (A + pn*(pn & (pEnd - A))) - p[endPointi];
360  }
361  }
362 }
363 
364 
366 {
367  clearAddressing();
368 }
369 
370 
372 {
373  clearAddressing();
374 }
375 
376 
378 {
379  clearAddressing();
380 }
381 
382 
384 {
385  return true;
386 }
387 
388 
389 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Templated abstract base-class for demand-driven mesh objects used to automate their allocation to the...
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:224
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
Class containing mesh-to-mesh mapping information.
Definition: polyMeshMap.H:51
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Class applies a two-dimensional correction to mesh motion point field.
virtual bool movePoints()
Correct weighting factors for moving mesh.
virtual void topoChange(const polyTopoChangeMap &)
Update topology.
virtual void distribute(const polyDistributionMap &)
Redistribute or update using the given distribution map.
void correctDisplacement(const pointField &p, vectorField &disp) const
Correct motion displacements.
virtual void mapMesh(const polyMeshMap &)
Update from another mesh using the given map.
twoDPointCorrector(const polyMesh &mesh)
Construct from mesh.
direction normalDir() const
Return direction normal to plane.
void correctPoints(pointField &p) const
Correct motion points.
const labelList & normalEdgeIndices() const
Return indices of normal edges.
const vector & planeNormal() const
Return plane normal.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
const fvPatchList & patches
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar e
Elementary charge.
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:613
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
const doubleScalar e
Definition: doubleScalar.H:106
dimensionedScalar tan(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensionedScalar sign(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
void deleteDemandDrivenData(DataType *&dataPtr)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
error FatalError
List< edge > edgeList
Definition: edgeList.H:38
static const char nl
Definition: Ostream.H:266
uint8_t direction
Definition: direction.H:45
dimensionedScalar acos(const dimensionedScalar &ds)
volScalarField & p