Merge branch 'master' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2
[foam-extend-3.2.git] / applications / solvers / surfaceTracking / freeSurface / freeSurfacePointDisplacement.C
bloba54cde36a02ae36fbba2dab368e474b0043302fe
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | foam-extend: Open Source CFD
4    \\    /   O peration     | Version:     3.2
5     \\  /    A nd           | Web:         http://www.foam-extend.org
6      \\/     M anipulation  | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9     This file is part of foam-extend.
11     foam-extend is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by the
13     Free Software Foundation, either version 3 of the License, or (at your
14     option) any later version.
16     foam-extend is distributed in the hope that it will be useful, but
17     WITHOUT ANY WARRANTY; without even the implied warranty of
18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19     General Public License for more details.
21     You should have received a copy of the GNU General Public License
22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
27 #include "freeSurface.H"
28 #include "primitivePatchInterpolation.H"
29 #include "emptyFaPatch.H"
30 #include "wedgeFaPatch.H"
31 #include "PstreamCombineReduceOps.H"
32 #include "coordinateSystem.H"
33 #include "scalarMatrices.H"
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 namespace Foam
41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
44 tmp<vectorField> freeSurface::pointDisplacement(const scalarField& deltaH)
46     const pointField& points = aMesh().patch().localPoints();
47     const labelListList& pointFaces = aMesh().patch().pointFaces();
49     controlPoints() += facesDisplacementDir()*deltaH;
51     tmp<vectorField> tdisplacement
52     (
53         new vectorField
54         (
55             points.size(),
56             vector::zero
57         )
58     );
60     vectorField& displacement = tdisplacement();
63     // Calculate displacement of internal points
64     const vectorField& pointNormals = aMesh().pointAreaNormals();
65     const edgeList& edges = aMesh().patch().edges();
66     labelList internalPoints = aMesh().internalPoints();
68     forAll (internalPoints, pointI)
69     {
70         label curPoint = internalPoints[pointI];
72         const labelList& curPointFaces = pointFaces[curPoint];
74         vectorField lsPoints(curPointFaces.size(), vector::zero);
76         for (label i=0; i<curPointFaces.size(); i++)
77         {
78             label curFace = curPointFaces[i];
80             lsPoints[i] = controlPoints()[curFace];
81         }
83         vectorField pointAndNormal =
84             lsPlanePointAndNormal
85             (
86                 lsPoints,
87                 points[curPoint],
88                 pointNormals[curPoint]
89             );
91         vector& P = pointAndNormal[0];
92         vector& N = pointAndNormal[1];
94         displacement[curPoint] =
95             pointsDisplacementDir()[curPoint]
96            *((P - points[curPoint])&N)
97            /(pointsDisplacementDir()[curPoint]&N);
98     }
101     // Mirror control points
102     FieldField<Field, vector> patchMirrorPoints(aMesh().boundary().size());
104     forAll(patchMirrorPoints, patchI)
105     {
106         patchMirrorPoints.set
107         (
108             patchI,
109             new vectorField
110             (
111                 aMesh().boundary()[patchI].faPatch::size(),
112                 vector::zero
113             )
114         );
116         vectorField N =
117             aMesh().boundary()[patchI].ngbPolyPatchFaceNormals();
119         const labelList peFaces =
120             labelList::subList
121             (
122                 aMesh().edgeOwner(),
123                 aMesh().boundary()[patchI].faPatch::size(),
124                 aMesh().boundary()[patchI].start()
125             );
127         const labelList& pEdges = aMesh().boundary()[patchI];
129         vectorField peCentres(pEdges.size(), vector::zero);
130         forAll(peCentres, edgeI)
131         {
132             peCentres[edgeI] =
133                 edges[pEdges[edgeI]].centre(points);
134         }
136         vectorField delta =
137             vectorField(controlPoints(), peFaces)
138           - peCentres;
140         patchMirrorPoints[patchI] =
141             peCentres + ((I - 2*N*N)&delta);
142     }
145     // Calculate displacement of boundary points
146     labelList boundaryPoints = aMesh().boundaryPoints();
148     const labelListList& edgeFaces = aMesh().patch().edgeFaces();
149     const labelListList& pointEdges = aMesh().patch().pointEdges();
151     forAll (boundaryPoints, pointI)
152     {
153         label curPoint = boundaryPoints[pointI];
155         if (motionPointsMask()[curPoint] == 1)
156         {
157             // Calculating mirror points
158             const labelList& curPointEdges = pointEdges[curPoint];
160             vectorField mirrorPoints(2, vector::zero);
162             label counter = -1;
164             forAll (curPointEdges, edgeI)
165             {
166                 label curEdge = curPointEdges[edgeI];
168                 if(edgeFaces[curEdge].size() == 1)
169                 {
170                     label patchID = -1;
171                     label edgeID = -1;
172                     forAll(aMesh().boundary(), patchI)
173                     {
174                         const labelList& pEdges =
175                             aMesh().boundary()[patchI];
176                         label index = findIndex(pEdges, curEdge);
177                         if (index != -1)
178                         {
179                             patchID = patchI;
180                             edgeID = index;
181                             break;
182                         }
183                     }
185                     mirrorPoints[++counter] =
186                         patchMirrorPoints[patchID][edgeID];
187                 }
188             }
190             // Calculating LS plane fit
191             const labelList& curPointFaces = pointFaces[curPoint];
193             vectorField lsPoints
194             (
195                 curPointFaces.size() + mirrorPoints.size(),
196                 vector::zero
197             );
199             counter = -1;
201             for (label i=0; i<curPointFaces.size(); i++)
202             {
203                 label curFace = curPointFaces[i];
205                 lsPoints[++counter] = controlPoints()[curFace];
206             }
208             for (label i=0; i<mirrorPoints.size(); i++)
209             {
210                 lsPoints[++counter] = mirrorPoints[i];
211             }
213             vectorField pointAndNormal =
214                 lsPlanePointAndNormal
215                 (
216                     lsPoints,
217                     points[curPoint],
218                     pointNormals[curPoint]
219                 );
221             vector& P = pointAndNormal[0];
222             vector& N = pointAndNormal[1];
224             displacement[curPoint] =
225                 pointsDisplacementDir()[curPoint]
226                *((P - points[curPoint])&N)
227                /(pointsDisplacementDir()[curPoint]&N);
228         }
229     }
232     // Calculate displacement of axis point
233     forAll (aMesh().boundary(), patchI)
234     {
235         if
236         (
237             aMesh().boundary()[patchI].type()
238          == wedgeFaPatch::typeName
239         )
240         {
241             const wedgeFaPatch& wedgePatch =
242                 refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
244             if(wedgePatch.axisPoint() > -1)
245             {
246                 label axisPoint = wedgePatch.axisPoint();
248                 displacement[axisPoint] =
249                     pointsDisplacementDir()[axisPoint]
250                    *(
251                         pointsDisplacementDir()[axisPoint]
252                        &(
253                             controlPoints()[pointFaces[axisPoint][0]]
254                           - points[axisPoint]
255                         )
256                     );
257             }
258         }
259     }
262     // Calculate displacement of processor patch points
263     forAll (aMesh().boundary(), patchI)
264     {
265         if
266         (
267             aMesh().boundary()[patchI].type()
268          == processorFaPatch::typeName
269         )
270         {
271             const processorFaPatch& procPatch =
272                 refCast<const processorFaPatch>(aMesh().boundary()[patchI]);
274             const labelList& patchPointLabels =
275                 procPatch.pointLabels();
277             FieldField<Field, vector> lsPoints(patchPointLabels.size());
278             forAll(lsPoints, pointI)
279             {
280                 lsPoints.set(pointI, new vectorField(0, vector::zero));
281             }
283             const labelList& nonGlobalPatchPoints =
284                 procPatch.nonGlobalPatchPoints();
286             forAll(nonGlobalPatchPoints, pointI)
287             {
288                 label curPatchPoint =
289                     nonGlobalPatchPoints[pointI];
291                 label curPoint =
292                     patchPointLabels[curPatchPoint];
294                 const labelList& curPointFaces = pointFaces[curPoint];
296                 lsPoints[curPatchPoint].setSize(curPointFaces.size());
298                 forAll(curPointFaces, faceI)
299                 {
300                     label curFace = curPointFaces[faceI];
302                     lsPoints[curPatchPoint][faceI] = controlPoints()[curFace];
303                 }
305 #               include "boundaryProcessorFaPatchPoints.H"
306             }
308             scalar lsPointsSize = 0;
309             forAll(lsPoints, pointI)
310             {
311                 lsPointsSize +=
312                     2*lsPoints[pointI].size()*sizeof(vector);
313             }
315             // Parallel data exchange
316             {
317                 OPstream toNeighbProc
318                 (
319                     Pstream::blocking,
320                     procPatch.neighbProcNo(),
321                     lsPointsSize
322                 );
324                 toNeighbProc << lsPoints;
325             }
327             FieldField<Field, vector> ngbLsPoints(patchPointLabels.size());
329             {
330                 IPstream fromNeighbProc
331                 (
332                     Pstream::blocking,
333                     procPatch.neighbProcNo(),
334                     lsPointsSize
335                 );
337                 fromNeighbProc >> ngbLsPoints;
338             }
340             forAll(nonGlobalPatchPoints, pointI)
341             {
342                 label curPatchPoint =
343                     nonGlobalPatchPoints[pointI];
345                 label curPoint =
346                     patchPointLabels[curPatchPoint];
348                 label curNgbPoint = procPatch.neighbPoints()[curPatchPoint];
350                 vectorField allLsPoints
351                 (
352                     lsPoints[curPatchPoint].size()
353                   + ngbLsPoints[curNgbPoint].size(),
354                     vector::zero
355                 );
357                 label counter = -1;
358                 forAll(lsPoints[curPatchPoint], pointI)
359                 {
360                     allLsPoints[++counter] = lsPoints[curPatchPoint][pointI];
361                 }
362                 forAll(ngbLsPoints[curNgbPoint], pointI)
363                 {
364                     allLsPoints[++counter] = ngbLsPoints[curNgbPoint][pointI];
365                 }
367                 vectorField pointAndNormal =
368                     lsPlanePointAndNormal
369                     (
370                         allLsPoints,
371                         points[curPoint],
372                         pointNormals[curPoint]
373                     );
375                 vector& P = pointAndNormal[0];
376                 vector& N = pointAndNormal[1];
378                 if (motionPointsMask()[curPoint] != 0)
379                 {
380                     displacement[curPoint] =
381                         pointsDisplacementDir()[curPoint]
382                        *((P - points[curPoint])&N)
383                        /(pointsDisplacementDir()[curPoint]&N);
384                 }
385             }
386         }
387     }
390     // Calculate displacement of global processor patch points
391     if (aMesh().globalData().nGlobalPoints() > 0)
392     {
393         const labelList& spLabels =
394             aMesh().globalData().sharedPointLabels();
396         const labelList& addr = aMesh().globalData().sharedPointAddr();
398         for (label k=0; k<aMesh().globalData().nGlobalPoints(); k++)
399         {
400             List<List<vector> > procLsPoints(Pstream::nProcs());
402             label curSharedPointIndex = findIndex(addr, k);
404             if (curSharedPointIndex != -1)
405             {
406                 label curPoint = spLabels[curSharedPointIndex];
408                 const labelList& curPointFaces = pointFaces[curPoint];
410                 procLsPoints[Pstream::myProcNo()] =
411                     List<vector>(curPointFaces.size());
413                 forAll (curPointFaces, faceI)
414                 {
415                     label curFace = curPointFaces[faceI];
417                     procLsPoints[Pstream::myProcNo()][faceI] =
418                         controlPoints()[curFace];
419                 }
420             }
422             Pstream::gatherList(procLsPoints);
423             Pstream::scatterList(procLsPoints);
425             if (curSharedPointIndex != -1)
426             {
427                 label curPoint = spLabels[curSharedPointIndex];
429                 label nAllPoints = 0;
430                 forAll(procLsPoints, procI)
431                 {
432                     nAllPoints += procLsPoints[procI].size();
433                 }
435                 vectorField allPoints(nAllPoints, vector::zero);
437                 label counter = 0;
438                 forAll(procLsPoints, procI)
439                 {
440                     forAll(procLsPoints[procI], pointI)
441                     {
442                         allPoints[counter++] =
443                             procLsPoints[procI][pointI];
444                     }
445                 }
447                 vectorField pointAndNormal =
448                     lsPlanePointAndNormal
449                     (
450                         allPoints,
451                         points[curPoint],
452                         pointNormals[curPoint]
453                     );
455                 const vector& P = pointAndNormal[0];
456                 const vector& N = pointAndNormal[1];
458                 displacement[curPoint] =
459                     pointsDisplacementDir()[curPoint]
460                    *((P - points[curPoint])&N)
461                    /(pointsDisplacementDir()[curPoint]&N);
462             }
463         }
464     }
466     return tdisplacement;
470 tmp<vectorField> freeSurface::lsPlanePointAndNormal
472     const vectorField& points,
473     const vector& origin,
474     const vector& axis
475 ) const
477     // LS in local CS
478     vector dir = (points[0] - origin);
479     dir -= axis*(axis&dir);
480     dir /= mag(dir);
481     coordinateSystem cs("cs", origin, axis, dir);
483     vectorField localPoints = cs.localPosition(points);
484     scalarField W = 1.0/(mag(points - origin) + SMALL);
486     scalarRectangularMatrix M
487     (
488         points.size(),
489         3,
490         0.0
491     );
493     for (label i=0; i<localPoints.size(); i++)
494     {
495         M[i][0] = localPoints[i].x();
496         M[i][1] = localPoints[i].y();
497         M[i][2] = 1.0;
498     }
500     scalarSquareMatrix MtM(3, 0.0);
501     for (label i = 0; i < MtM.n(); i++)
502     {
503         for (label j = 0; j < MtM.m(); j++)
504         {
505             for (label k = 0; k < M.n(); k++)
506             {
507                 MtM[i][j] += M[k][i]*M[k][j]*W[k];
508             }
509         }
510     }
512     scalarField MtR(3, 0);
513     for (label i = 0; i < MtR.size(); i++)
514     {
515         for (label j = 0; j < M.n(); j++)
516         {
517             MtR[i] += M[j][i]*localPoints[j].z()*W[j];
518         }
519     }
521     scalarSquareMatrix::LUsolve(MtM, MtR);
523     vector n0 = vector(-MtR[0], -MtR[1], 1);
524     n0 = cs.globalVector(n0);
525     n0 /= mag(n0);
527     vector p0 = vector(0, 0, MtR[2]);
528     p0 = cs.globalPosition(p0);
530     tmp<vectorField> pointAndNormal
531     (
532         new vectorField(2, vector::zero)
533     );
535     pointAndNormal()[0] = p0;
536     pointAndNormal()[1] = n0;
538     return pointAndNormal;
542 // tmp<vectorField> freeSurface::pointDisplacementSM()
543 // {
544 //     const pointField& points = aMesh().patch().localPoints();
545 //     const labelListList& pointFaces = aMesh().patch().pointFaces();
548 //     tmp<vectorField> tdisplacement
549 //     (
550 //         new vectorField
551 //         (
552 //             points.size(),
553 //             vector::zero
554 //         )
555 //     );
557 //     vectorField& displacement = tdisplacement();
560 //     forAll (pointFaces, pointI)
561 //     {
562 //         scalar weightsSum = 0.0;
564 //         const labelList& curPointFaces = pointFaces[pointI];
566 //         forAll (curPointFaces, faceI)
567 //         {
568 //             label curFace = curPointFaces[faceI];
570 //             scalar weight = 1.0/mag
571 //             (
572 //                 points[pointI]
573 //               - controlPoints()[curFace]
574 //             );
576 //             displacement[pointI] += weight*controlPoints()[curFace];
578 //             weightsSum += weight;
579 //         }
581 //         displacement[pointI] /= weightsSum;
583 //         displacement[pointI] -= points[pointI];
584 //     }
587 //     displacement = motionPointsMask()*
588 //         (pointsDisplacementDir()&displacement)*
589 //         pointsDisplacementDir();
592 //     return tdisplacement;
593 // }
596 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
598 } // End namespace Foam
600 // ************************************************************************* //