freeSurfacePointDisplacement.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
30#include "emptyFaPatch.H"
31#include "processorFaPatch.H"
32#include "wedgeFaPatch.H"
33#include "wallFvPatch.H"
34#include "coordinateSystem.H"
35#include "unitConversion.H"
36#include "scalarMatrices.H"
37#include "tensor2D.H"
38
39// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40
42Foam::interfaceTrackingFvMesh::pointDisplacement()
43{
46
47 auto tdisplacement = tmp<vectorField>::New(points.size(), Zero);
48 auto& displacement = tdisplacement.ref();
49
50 // Calculate displacement of internal points
51 const vectorField& pointNormals = aMesh().pointAreaNormals();
52 const edgeList& edges = aMesh().patch().edges();
53 labelList internalPoints = aMesh().internalPoints();
54
55 for (const label curPoint : internalPoints)
56 {
57 const labelList& curPointFaces = pointFaces[curPoint];
58
59 vectorField lsPoints(curPointFaces.size(), Zero);
60
61 for (label i=0; i<curPointFaces.size(); i++)
62 {
63 label curFace = curPointFaces[i];
64
65 lsPoints[i] = controlPoints()[curFace];
66 }
67
68 vectorField pointAndNormal
69 (
70 lsPlanePointAndNormal
71 (
72 lsPoints,
73 points[curPoint],
74 pointNormals[curPoint]
75 )
76 );
77
78 vector& P = pointAndNormal[0];
79 vector& N = pointAndNormal[1];
80
81 displacement[curPoint] =
82 pointsDisplacementDir()[curPoint]
83 *((P - points[curPoint])&N)
84 /(pointsDisplacementDir()[curPoint]&N);
85 }
86
87 // Mirror control points
88 FieldField<Field, vector> patchMirrorPoints(aMesh().boundary().size());
89
90 // Old faMesh points
92 const labelList& meshPoints = aMesh().patch().meshPoints();
93 forAll(oldPoints, pI)
94 {
95 oldPoints[pI] =
96 mesh().oldPoints()[meshPoints[pI]];
97 }
98
99 forAll(patchMirrorPoints, patchI)
100 {
101 patchMirrorPoints.set
102 (
103 patchI,
104 new vectorField
105 (
106 aMesh().boundary()[patchI].faPatch::size(),
107 Zero
108 )
109 );
110
112 (
113 aMesh().boundary()[patchI].ngbPolyPatchFaceNormals()
114 );
115
116 const labelList& eFaces =
117 aMesh().boundary()[patchI].edgeFaces();
118
119 // Correct N according to specified contact angle
120 if (contactAnglePtr_)
121 {
122 label ngbPolyPatchID =
123 aMesh().boundary()[patchI].ngbPolyPatchIndex();
124
125 if (ngbPolyPatchID != -1)
126 {
127 if
128 (
129 mesh().boundary()[ngbPolyPatchID].type()
131 )
132 {
133 // Info<< aMesh().boundary()[patchI].name() << endl;
134
135 scalar rotAngle = degToRad
136 (
138 (
139 90
140 - contactAnglePtr_->boundaryField()[patchI]
141 )
142 );
143
144 const vectorField& pEdgN =
145 aMesh().edgeAreaNormals().boundaryField()[patchI];
146
147 vectorField rotationAxis( N^pEdgN );
148
149 const edgeList::subList patchEdges =
150 aMesh().boundary()[patchI].patchSlice(aMesh().edges());
151
152 forAll(rotationAxis, edgeI)
153 {
154 vector e = patchEdges[edgeI].vec(oldPoints);
155 // vector e = patchEdges[edgeI].vec(aMesh().points());
156
157 // Adjust direction
158 rotationAxis[edgeI] =
159 e*(e&rotationAxis[edgeI])
160 /mag((e&rotationAxis[edgeI]));
161 }
162 rotationAxis /= mag(rotationAxis) + SMALL;
163
164 vectorField rotationAxis2 = rotationAxis;
165 forAll(rotationAxis2, edgeI)
166 {
167 rotationAxis2[edgeI] =
168 (N[edgeI]^facesDisplacementDir()[eFaces[edgeI]]);
169
170 // Adjust direction
171 rotationAxis2[edgeI] =
172 rotationAxis2[edgeI]
173 *(rotationAxis2[edgeI]&rotationAxis[edgeI])
174 /(
175 mag((rotationAxis2[edgeI]&rotationAxis[edgeI]))
176 + SMALL
177 );
178 }
179 rotationAxis2 /= mag(rotationAxis2) + SMALL;
180
181 // Rodrigues' rotation formula
182 N = N*cos(rotAngle)
183 + rotationAxis*(rotationAxis & N)*(1 - cos(rotAngle))
184 + (rotationAxis^N)*sin(rotAngle);
185
186 N /= mag(N) + SMALL;
187
188 N = (rotationAxis^N);
189
190 N = (N^rotationAxis2);
191
192 N /= mag(N) + SMALL;
193
194 // Info<< N << endl;
195 }
196 }
197 }
198
199 const labelList peFaces =
201 (
202 aMesh().edgeOwner(),
203 aMesh().boundary()[patchI].faPatch::size(),
204 aMesh().boundary()[patchI].start()
205 );
206
207 const labelList& pEdges = aMesh().boundary()[patchI];
208
209 vectorField peCentres(pEdges.size(), Zero);
210 forAll(peCentres, edgeI)
211 {
212 peCentres[edgeI] =
213 edges[pEdges[edgeI]].centre(points);
214 }
215
217 (
218 vectorField(controlPoints(), peFaces)
219 - peCentres
220 );
221
222 // Info<< aMesh().boundary()[patchI].name() << endl;
223 // Info<< vectorField(controlPoints(), peFaces) << endl;
224
225 patchMirrorPoints[patchI] =
226 peCentres + ((I - 2*N*N)&delta);
227
228 // Info<< patchMirrorPoints[patchI] << endl;
229 }
230
231 // Calculate displacement of boundary points
232 labelList boundaryPoints = aMesh().boundaryPoints();
233
236
237 for (const label curPoint : boundaryPoints)
238 {
239 if (motionPointsMask()[curPoint] == 1)
240 {
241 // Calculating mirror points
242 const labelList& curPointEdges = pointEdges[curPoint];
243
244 vectorField mirrorPoints(2, Zero);
245
246 label counter = -1;
247
248 forAll(curPointEdges, edgeI)
249 {
250 label curEdge = curPointEdges[edgeI];
251
252 if(edgeFaces[curEdge].size() == 1)
253 {
254 label patchID = -1;
255 label edgeID = -1;
256 forAll(aMesh().boundary(), patchI)
257 {
258 const labelList& pEdges =
259 aMesh().boundary()[patchI];
260 label index = pEdges.find(curEdge);
261 if (index != -1)
262 {
263 patchID = patchI;
264 edgeID = index;
265 break;
266 }
267 }
268
269 mirrorPoints[++counter] =
270 patchMirrorPoints[patchID][edgeID];
271 }
272 }
273
274 // Calculating LS plane fit
275 const labelList& curPointFaces = pointFaces[curPoint];
276
277 vectorField lsPoints
278 (
279 curPointFaces.size() + mirrorPoints.size(),
280 Zero
281 );
282
283 counter = -1;
284
285 for (label i=0; i<curPointFaces.size(); i++)
286 {
287 label curFace = curPointFaces[i];
288
289 lsPoints[++counter] = controlPoints()[curFace];
290 }
291
292 for (label i=0; i<mirrorPoints.size(); i++)
293 {
294 lsPoints[++counter] = mirrorPoints[i];
295 }
296
297 vectorField pointAndNormal
298 (
299 lsPlanePointAndNormal
300 (
301 lsPoints,
302 points[curPoint],
303 pointNormals[curPoint]
304 )
305 );
306
307 vector& P = pointAndNormal[0];
308 vector& N = pointAndNormal[1];
309
310 displacement[curPoint] =
311 pointsDisplacementDir()[curPoint]
312 *((P - points[curPoint])&N)
313 /(pointsDisplacementDir()[curPoint]&N);
314 }
315 }
316
317 // Calculate displacement of processor patch points
318 forAll(aMesh().boundary(), patchI)
319 {
320 if
321 (
322 aMesh().boundary()[patchI].type()
324 )
325 {
326 const processorFaPatch& procPatch =
327 refCast<const processorFaPatch>(aMesh().boundary()[patchI]);
328
329 const labelList& patchPointLabels =
330 procPatch.pointLabels();
331
332 FieldField<Field, vector> lsPoints(patchPointLabels.size());
333 forAll(lsPoints, pointI)
334 {
335 lsPoints.set(pointI, new vectorField(0, Zero));
336 }
337
338 const labelList& nonGlobalPatchPoints =
339 procPatch.nonGlobalPatchPoints();
340
341 forAll(nonGlobalPatchPoints, pointI)
342 {
343 label curPatchPoint =
344 nonGlobalPatchPoints[pointI];
345
346 label curPoint =
347 patchPointLabels[curPatchPoint];
348
349 const labelList& curPointFaces = pointFaces[curPoint];
350
351 lsPoints[curPatchPoint].setSize(curPointFaces.size());
352
353 forAll(curPointFaces, faceI)
354 {
355 label curFace = curPointFaces[faceI];
356
357 lsPoints[curPatchPoint][faceI] = controlPoints()[curFace];
358 }
359
361 }
362
363 // Parallel data exchange
364 {
365 OPstream toNeighbProc
366 (
368 procPatch.neighbProcNo()
369 );
370
371 toNeighbProc << lsPoints;
372 }
373
374 FieldField<Field, vector> ngbLsPoints(patchPointLabels.size());
375 {
376 IPstream fromNeighbProc
377 (
379 procPatch.neighbProcNo()
380 );
381
382 fromNeighbProc >> ngbLsPoints;
383 }
384
385 forAll(nonGlobalPatchPoints, pointI)
386 {
387 label curPatchPoint =
388 nonGlobalPatchPoints[pointI];
389
390 label curPoint =
391 patchPointLabels[curPatchPoint];
392
393 label curNgbPoint = procPatch.neighbPoints()[curPatchPoint];
394
395 vectorField allLsPoints
396 (
397 lsPoints[curPatchPoint].size()
398 + ngbLsPoints[curNgbPoint].size(),
399 Zero
400 );
401
402 label counter = -1;
403 forAll(lsPoints[curPatchPoint], pointI)
404 {
405 allLsPoints[++counter] = lsPoints[curPatchPoint][pointI];
406 }
407 forAll(ngbLsPoints[curNgbPoint], pointI)
408 {
409 allLsPoints[++counter] = ngbLsPoints[curNgbPoint][pointI];
410 }
411
412 vectorField pointAndNormal
413 (
414 lsPlanePointAndNormal
415 (
416 allLsPoints,
417 points[curPoint],
418 pointNormals[curPoint]
419 )
420 );
421
422 vector& P = pointAndNormal[0];
423 vector& N = pointAndNormal[1];
424
425 if (motionPointsMask()[curPoint] != 0)
426 {
427 displacement[curPoint] =
428 pointsDisplacementDir()[curPoint]
429 *((P - points[curPoint])&N)
430 /(pointsDisplacementDir()[curPoint]&N);
431 }
432 }
433 }
434 }
435
436
437 // Calculate displacement of global processor patch points
438 if (aMesh().globalData().nGlobalPoints() > 0)
439 {
440 const labelList& spLabels =
442
443 const labelList& addr = aMesh().globalData().sharedPointAddr();
444
445 for (label k=0; k<aMesh().globalData().nGlobalPoints(); k++)
446 {
447 List<List<vector>> procLsPoints(Pstream::nProcs());
448
449 label curSharedPointIndex = addr.find(k);
450
451 if (curSharedPointIndex != -1)
452 {
453 label curPoint = spLabels[curSharedPointIndex];
454
455 const labelList& curPointFaces = pointFaces[curPoint];
456
457 procLsPoints[Pstream::myProcNo()] =
458 List<vector>(curPointFaces.size());
459
460 forAll(curPointFaces, faceI)
461 {
462 label curFace = curPointFaces[faceI];
463
464 procLsPoints[Pstream::myProcNo()][faceI] =
465 controlPoints()[curFace];
466 }
467 }
468
469 Pstream::allGatherList(procLsPoints);
470
471 if (curSharedPointIndex != -1)
472 {
473 label curPoint = spLabels[curSharedPointIndex];
474
475 label nAllPoints = 0;
476 forAll(procLsPoints, procI)
477 {
478 nAllPoints += procLsPoints[procI].size();
479 }
480
481 vectorField allPoints(nAllPoints, Zero);
482
483 label counter = 0;
484 forAll(procLsPoints, procI)
485 {
486 forAll(procLsPoints[procI], pointI)
487 {
488 allPoints[counter++] =
489 procLsPoints[procI][pointI];
490 }
491 }
492
493 vectorField pointAndNormal
494 (
495 lsPlanePointAndNormal
496 (
497 allPoints,
498 points[curPoint],
499 pointNormals[curPoint]
500 )
501 );
502
503 const vector& P = pointAndNormal[0];
504 const vector& N = pointAndNormal[1];
505
506 displacement[curPoint] =
507 pointsDisplacementDir()[curPoint]
508 *((P - points[curPoint])&N)
509 /(pointsDisplacementDir()[curPoint]&N);
510 }
511 }
512 }
513
514 return tdisplacement;
515}
516
517
519Foam::interfaceTrackingFvMesh::lsPlanePointAndNormal
520(
521 const vectorField& points,
522 const vector& origin,
523 const vector& axis
524) const
525{
526 // LS in local CS
527 vector dir = (points[0] - origin);
528 dir -= axis*(axis&dir);
529 dir /= mag(dir);
530 coordinateSystem cs("cs", origin, axis, dir);
531
532 vectorField localPoints(cs.localPosition(points));
533 vector avgLocalPoint = average(localPoints);
534
535 // scalarField W = 1.0/(mag(points - origin) + SMALL);
536 scalarField W(points.size(), scalar(1));
537
538 const label nCoeffs = 2;
540
541 scalar L = 2*max(mag(localPoints-avgLocalPoint));
542 for (label i=0; i<localPoints.size(); i++)
543 {
544 M[i][0] = (localPoints[i].x() - avgLocalPoint.x())/L;
545 M[i][1] = (localPoints[i].y() - avgLocalPoint.y())/L;
546 }
547
548 // Applying weights
549 for (label i=0; i<M.n(); i++)
550 {
551 for (label j=0; j<M.m(); j++)
552 {
553 M[i][j] *= W[i];
554 }
555 }
556
557 tensor2D lsM(Zero);
558 // scalarSquareMatrix lsM(nCoeffs, Zero);
559
560 for (label i=0; i<nCoeffs; i++)
561 {
562 for (label j=0; j<nCoeffs; j++)
563 {
564 for (label k=0; k<M.n(); k++)
565 {
566 lsM[i*nCoeffs+j] += M[k][i]*M[k][j];
567 // lsM(i,j) += M[k][i]*M[k][j];
568 }
569 }
570 }
571
572 // Calculate inverse
573 tensor2D invLsM = inv(lsM);
574
575 scalarRectangularMatrix curInvMatrix(nCoeffs, points.size(), Zero);
576
577 for (label i=0; i<nCoeffs; i++)
578 {
579 for (label j=0; j<M.n(); j++)
580 {
581 for (label k=0; k<nCoeffs; k++)
582 {
583 curInvMatrix[i][j] += invLsM[i*nCoeffs+k]*M[j][k]*W[j];
584 // curInvMatrix[i][j] += invLsM[i][k]*M[j][k]*W[j];
585 }
586 }
587 }
588
589 scalarField coeffs(nCoeffs, Zero);
590 scalarField source(points.size(), Zero);
591
592 for (label i=0; i<points.size(); i++)
593 {
594 source[i] = (localPoints[i].z() - avgLocalPoint.z())/L;
595 }
596
597 for (label i=0; i<nCoeffs; i++)
598 {
599 for (label j=0; j<source.size(); j++)
600 {
601 coeffs[i] += curInvMatrix[i][j]*source[j];
602 }
603 }
604
605 vector n0(-coeffs[0], -coeffs[1], 1.0);
606 n0 = cs.globalVector(n0);
607 n0 /= mag(n0);
608
609 vector p0 = avgLocalPoint;
610 p0 = cs.globalPosition(p0);
611
612 auto tpointAndNormal = tmp<vectorField>::New(2);
613 auto& pointAndNormal = tpointAndNormal.ref();
614
615 pointAndNormal[0] = p0;
616 pointAndNormal[1] = n0;
617
618 return tpointAndNormal;
619}
620
621
622// ************************************************************************* //
label k
#define M(I)
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:80
const Boundary & boundaryField() const
Return const-reference to the boundary field.
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Input inter-processor communications stream.
Definition: IPstream.H:57
SubList< label > subList
Declare type of subList.
Definition: List.H:111
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Output inter-processor communications stream.
Definition: OPstream.H:57
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const labelListList & pointFaces() const
Return point-face addressing.
const labelListList & edgeFaces() const
Return edge-face addressing.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
A List obtained as a section of another List.
Definition: SubList.H:70
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
Definition: Tensor2D.H:59
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Base class for coordinate system specification, the default coordinate system type is cartesian .
UPtrList< const labelUList > edgeFaces() const
Return a list of edgeFaces for each patch.
const labelList & sharedPointAddr() const noexcept
Return addressing into the complete globally shared points list.
const labelList & sharedPointLabels() const noexcept
Return indices of local points that are globally shared.
label nGlobalPoints() const noexcept
Return number of globally shared points.
labelList internalPoints() const
Return internal point labels.
const faBoundaryMesh & boundary() const noexcept
Return constant reference to boundary mesh.
Definition: faMeshI.H:38
const uindirectPrimitivePatch & patch() const
Return constant reference to primitive patch.
Definition: faMeshI.H:128
const faGlobalMeshData & globalData() const
Return parallel info.
Definition: faMesh.C:893
const edgeVectorField & edgeAreaNormals() const
Return edge area normals.
Definition: faMesh.C:841
const vectorField & pointAreaNormals() const
Return point area normals.
Definition: faMesh.C:852
labelList boundaryPoints() const
Return boundary point labels.
virtual label size() const
Patch size is the number of edge labels.
Definition: faPatch.H:311
const labelList & pointLabels() const
Return patch point labels.
Definition: faPatch.C:273
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:712
tmp< surfaceVectorField > delta() const
Return face deltas as surfaceVectorField.
vectorField & controlPoints()
Return control points.
vectorField & pointsDisplacementDir()
Return reference to point displacement direction field.
vectorField & facesDisplacementDir()
Return reference to control points displacement direction field.
faMesh & aMesh()
Return reference to finite area mesh.
labelList & motionPointsMask()
Return reference to motion points mask field.
virtual const pointField & oldPoints() const
Return old points (mesh motion)
Definition: polyMesh.C:1133
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & pointFaces() const
const labelListList & edgeFaces() const
const labelList & neighbPoints() const
Return neighbour point labels. This is for my local point the.
int neighbProcNo() const noexcept
Return neighbour processor number.
int myProcNo() const noexcept
Return processor number.
const labelList & nonGlobalPatchPoints() const
Return the set of labels of the processor patch points which are.
A class for managing temporary objects.
Definition: tmp.H:65
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const volScalarField & p0
Definition: EEqn.H:36
const pointField & points
label nPoints
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar sin(const dimensionedScalar &ds)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
static const Identity< scalar > I
Definition: Identity.H:94
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.
Unit conversion functions.
const vector L(dict.get< vector >("L"))
const Vector< label > N(dict.get< Vector< label > >("N"))