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-2018 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 {
37  defineTypeNameAndDebug(twoDPointCorrector, 0);
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 
60  forAll(patches, patchi)
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  {
86  forAll(patches, patchi)
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 
197 Foam::twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
198 :
200  required_(mesh_.nGeometricD() == 2),
201  planeNormalPtr_(nullptr),
202  normalEdgeIndicesPtr_(nullptr),
203  isWedge_(false),
204  wedgeAxis_(Zero),
205  wedgeAngle_(0.0)
206 {}
207 
208 
209 
210 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
211 
213 {
214  clearAddressing();
215 }
216 
217 
218 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
219 
221 {
222  const vector& pn = planeNormal();
223 
224  if (mag(pn.x()) >= edgeOrthogonalityTol)
225  {
226  return vector::X;
227  }
228  else if (mag(pn.y()) >= edgeOrthogonalityTol)
229  {
230  return vector::Y;
231  }
232  else if (mag(pn.z()) >= edgeOrthogonalityTol)
233  {
234  return vector::Z;
235  }
236  else
237  {
239  << "Plane normal not aligned with the coordinate system" << nl
240  << " pn = " << pn
241  << abort(FatalError);
242 
243  return vector::Z;
244  }
245 }
246 
247 
249 {
250  if (!planeNormalPtr_)
251  {
252  calcAddressing();
253  }
254 
255  return *planeNormalPtr_;
256 }
257 
258 
260 {
261  if (!normalEdgeIndicesPtr_)
262  {
263  calcAddressing();
264  }
265 
266  return *normalEdgeIndicesPtr_;
267 }
268 
269 
271 {
272  if (!required_) return;
273 
274  // Algorithm:
275  // Loop through all edges. Calculate the average point position A for
276  // the front and the back. Correct the position of point P (in two planes)
277  // such that vectors AP and planeNormal are parallel
278 
279  // Get reference to edges
280  const edgeList& meshEdges = mesh_.edges();
281 
282  const labelList& neIndices = normalEdgeIndices();
283  const vector& pn = planeNormal();
284 
285  forAll(neIndices, edgeI)
286  {
287  point& pStart = p[meshEdges[neIndices[edgeI]].start()];
288 
289  point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
290 
291  // calculate average point position
292  point A = 0.5*(pStart + pEnd);
294 
295  if (isWedge_)
296  {
297  snapToWedge(pn, A, pStart);
298  snapToWedge(pn, A, pEnd);
299  }
300  else
301  {
302  // correct point locations
303  pStart = A + pn*(pn & (pStart - A));
304  pEnd = A + pn*(pn & (pEnd - A));
305  }
306  }
307 }
308 
309 
311 (
312  const pointField& p,
313  vectorField& disp
314 ) const
315 {
316  if (!required_) return;
317 
318  // Algorithm:
319  // Loop through all edges. Calculate the average point position A for
320  // the front and the back. Correct the position of point P (in two planes)
321  // such that vectors AP and planeNormal are parallel
322 
323  // Get reference to edges
324  const edgeList& meshEdges = mesh_.edges();
325 
326  const labelList& neIndices = normalEdgeIndices();
327  const vector& pn = planeNormal();
328 
329  forAll(neIndices, edgeI)
330  {
331  const edge& e = meshEdges[neIndices[edgeI]];
332 
333  label startPointi = e.start();
334  point pStart = p[startPointi] + disp[startPointi];
335 
336  label endPointi = e.end();
337  point pEnd = p[endPointi] + disp[endPointi];
338 
339  // calculate average point position
340  point A = 0.5*(pStart + pEnd);
342 
343  if (isWedge_)
344  {
345  snapToWedge(pn, A, pStart);
346  snapToWedge(pn, A, pEnd);
347  disp[startPointi] = pStart - p[startPointi];
348  disp[endPointi] = pEnd - p[endPointi];
349  }
350  else
351  {
352  // correct point locations
353  disp[startPointi] = (A + pn*(pn & (pStart - A))) - p[startPointi];
354  disp[endPointi] = (A + pn*(pn & (pEnd - A))) - p[endPointi];
355  }
356  }
357 }
358 
359 
361 {
362  clearAddressing();
363 }
364 
365 
367 {
368  return true;
369 }
370 
371 
372 // ************************************************************************* //
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.
Definition: UList.H:428
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
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
uint8_t direction
Definition: direction.H:45
const labelList & normalEdgeIndices() const
Return indices of normal edges.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const Cmpt & z() const
Definition: VectorI.H:87
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:237
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
const Cmpt & y() const
Definition: VectorI.H:81
void correctPoints(pointField &p) const
Correct motion points.
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
Definition: MeshObject.H:85
void correctDisplacement(const pointField &p, vectorField &disp) const
Correct motion displacements.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
List< edge > edgeList
Definition: edgeList.H:38
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Class applies a two-dimensional correction to mesh motion point field.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
static const char nl
Definition: Ostream.H:265
const vector & planeNormal() const
Return plane normal.
defineTypeNameAndDebug(combustionModel, 0)
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
label patchi
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
bool movePoints()
Correct weighting factors for moving mesh.
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
direction normalDir() const
Return direction normal to plane.
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
dimensionedScalar tan(const dimensionedScalar &ds)
void deleteDemandDrivenData(DataPtr &dataPtr)
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:613
label start() const
Return start vertex label.
Definition: edgeI.H:81
Namespace for OpenFOAM.