FreeFOAM The Cross-Platform CFD Toolkit
slidingInterfaceProjectPoints.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "slidingInterface.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/line.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 const Foam::scalar Foam::slidingInterface::pointMergeTolDefault_ = 0.05;
34 const Foam::scalar Foam::slidingInterface::edgeMergeTolDefault_ = 0.01;
35 const Foam::label Foam::slidingInterface::nFacesPerSlaveEdgeDefault_ = 5;
36 const Foam::label Foam::slidingInterface::edgeFaceEscapeLimitDefault_ = 10;
37 
38 const Foam::scalar Foam::slidingInterface::integralAdjTolDefault_ = 0.05;
39 const Foam::scalar Foam::slidingInterface::edgeMasterCatchFractionDefault_ = 0.4;
40 const Foam::scalar Foam::slidingInterface::edgeEndCutoffTolDefault_ = 0.0001;
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 // Index of debug signs:
46 // a - integral match adjustment: point adjusted
47 // n - integral match adjustment: point not adjusted
48 // . - loop of the edge-to-face interaction detection
49 // x - reversal of direction in edge-to-face interaction detection
50 // + - complete edge-to-face interaction detection
51 // z - incomplete edge-to-face interaction detection. This may be OK for edges
52 // crossing from one to the other side of multiply connected master patch
53 // * - colinear triangle: adjusting projection with slave face normal
54 // m - master point inserted into the edge
55 
56 bool Foam::slidingInterface::projectPoints() const
57 {
58  if (debug)
59  {
60  Pout<< "bool slidingInterface::projectPoints() : "
61  << " for object " << name() << " : "
62  << "Projecting slave points onto master surface." << endl;
63  }
64 
65  // Algorithm:
66  // 1) Go through all the points of the master and slave patch and calculate
67  // minimum edge length coming from the point. Calculate the point
68  // merge tolerance as the fraction of mimimum edge length.
69  // 2) Project all the slave points onto the master patch
70  // in the normal direction.
71  // 3) If some points have missed and the match is integral, the
72  // points in question will be adjusted. Find the nearest face for
73  // those points and continue with the procedure.
74  // 4) For every point, check all the points of the face it has hit.
75  // For every pair of points find if their distance is smaller than
76  // both the master and slave merge tolerance. If so, the slave point
77  // is moved to the location of the master point. Remember the master
78  // point index.
79  // 5) For every unmerged slave point, check its distance to all the
80  // edges of the face it has hit. If the distance is smaller than the
81  // edge merge tolerance, the point will be moved onto the edge.
82  // Remember the master edge index.
83  // 6) The remaning slave points will be projected into faces. Remember the
84  // master face index.
85  // 7) For the points that miss the master patch, grab the nearest face
86  // on the master and leave the slave point where it started
87  // from and the miss is recorded.
88 
89  const polyMesh& mesh = topoChanger().mesh();
90 
91  const primitiveFacePatch& masterPatch =
92  mesh.faceZones()[masterFaceZoneID_.index()]();
93 
94  const primitiveFacePatch& slavePatch =
95  mesh.faceZones()[slaveFaceZoneID_.index()]();
96 
97  // Get references to local points, local edges and local faces
98  // for master and slave patch
99 // const labelList& masterMeshPoints = masterPatch.meshPoints();
100  const pointField& masterLocalPoints = masterPatch.localPoints();
101  const faceList& masterLocalFaces = masterPatch.localFaces();
102  const edgeList& masterEdges = masterPatch.edges();
103  const labelListList& masterEdgeFaces = masterPatch.edgeFaces();
104  const labelListList& masterFaceEdges = masterPatch.faceEdges();
105  const labelListList& masterFaceFaces = masterPatch.faceFaces();
106 // Pout<< "Master patch. Local points: " << masterLocalPoints << nl
107 // << "Master patch. Mesh points: " << masterPatch.meshPoints() << nl
108 // << "Local faces: " << masterLocalFaces << nl
109 // << "local edges: " << masterEdges << endl;
110 
111 // const labelList& slaveMeshPoints = slavePatch.meshPoints();
112  const pointField& slaveLocalPoints = slavePatch.localPoints();
113  const edgeList& slaveEdges = slavePatch.edges();
114  const labelListList& slaveEdgeFaces = slavePatch.edgeFaces();
115  const vectorField& slavePointNormals = slavePatch.pointNormals();
116 // const vectorField& slaveFaceNormals = slavePatch.faceNormals();
117 // Pout<< "Slave patch. Local points: " << slaveLocalPoints << nl
118 // << "Slave patch. Mesh points: " << slavePatch.meshPoints() << nl
119 // << "Local faces: " << slavePatch.localFaces() << nl
120 // << "local edges: " << slaveEdges << endl;
121 
122  // Calculate min edge distance for points and faces
123 
124  // Calculate min edge length for the points and faces of master patch
125  scalarField minMasterPointLength(masterLocalPoints.size(), GREAT);
126  scalarField minMasterFaceLength(masterPatch.size(), GREAT);
127 
128  forAll (masterEdges, edgeI)
129  {
130  const edge& curEdge = masterEdges[edgeI];
131 
132  const scalar curLength =
133  masterEdges[edgeI].mag(masterLocalPoints);
134 
135  // Do points
136  minMasterPointLength[curEdge.start()] =
137  min
138  (
139  minMasterPointLength[curEdge.start()],
140  curLength
141  );
142 
143  minMasterPointLength[curEdge.end()] =
144  min
145  (
146  minMasterPointLength[curEdge.end()],
147  curLength
148  );
149 
150  // Do faces
151  const labelList& curFaces = masterEdgeFaces[edgeI];
152 
153  forAll (curFaces, faceI)
154  {
155  minMasterFaceLength[curFaces[faceI]] =
156  min
157  (
158  minMasterFaceLength[curFaces[faceI]],
159  curLength
160  );
161  }
162  }
163 
164 // Pout << "min length for master points: " << minMasterPointLength << endl
165 // << "min length for master faces: " << minMasterFaceLength << endl;
166 
167  // Calculate min edge length for the points and faces of slave patch
168  scalarField minSlavePointLength(slaveLocalPoints.size(), GREAT);
169  scalarField minSlaveFaceLength(slavePatch.size(), GREAT);
170 
171  forAll (slaveEdges, edgeI)
172  {
173  const edge& curEdge = slaveEdges[edgeI];
174 
175  const scalar curLength =
176  slaveEdges[edgeI].mag(slaveLocalPoints);
177 
178  // Do points
179  minSlavePointLength[curEdge.start()] =
180  min
181  (
182  minSlavePointLength[curEdge.start()],
183  curLength
184  );
185 
186  minSlavePointLength[curEdge.end()] =
187  min
188  (
189  minSlavePointLength[curEdge.end()],
190  curLength
191  );
192 
193  // Do faces
194  const labelList& curFaces = slaveEdgeFaces[edgeI];
195 
196  forAll (curFaces, faceI)
197  {
198  minSlaveFaceLength[curFaces[faceI]] =
199  min
200  (
201  minSlaveFaceLength[curFaces[faceI]],
202  curLength
203  );
204  }
205  }
206 
207 // Pout << "min length for slave points: " << minSlavePointLength << endl
208 // << "min length for slave faces: " << minSlaveFaceLength << endl;
209 
210  // Project slave points onto the master patch
211 
212  // Face hit by the slave point
213  List<objectHit> slavePointFaceHits =
214  slavePatch.projectPoints
215  (
216  masterPatch,
217  slavePointNormals,
218  projectionAlgo_
219  );
220 
221 // Pout << "USING N-SQAURED!!!" << endl;
222 // List<objectHit> slavePointFaceHits =
223 // projectPointsNSquared<face, List, const pointField&>
224 // (
225 // slavePatch,
226 // masterPatch,
227 // slavePointNormals,
228 // projectionAlgo_
229 // );
230 
231  if (debug)
232  {
233  label nHits = 0;
234 
235  forAll (slavePointFaceHits, pointI)
236  {
237  if (slavePointFaceHits[pointI].hit())
238  {
239  nHits++;
240  }
241  }
242 
243  Pout<< "Number of hits in point projection: " << nHits
244  << " out of " << slavePointFaceHits.size() << " points."
245  << endl;
246  }
247 
248  // Projected slave points are stored for mesh motion correction
249  if (projectedSlavePointsPtr_) delete projectedSlavePointsPtr_;
250 
251  projectedSlavePointsPtr_ =
252  new pointField(slavePointFaceHits.size(), vector::zero);
253  pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
254 
255  // Adjust projection to type of match
256 
257  label nAdjustedPoints = 0;
258 
259  // If the match is integral and some points have missed,
260  // find nearest master face
261  if (matchType_ == INTEGRAL)
262  {
263  if (debug)
264  {
265  Pout<< "bool slidingInterface::projectPoints() for object "
266  << name() << " : "
267  << "Adjusting point projection for integral match: ";
268  }
269 
270  forAll (slavePointFaceHits, pointI)
271  {
272  if (slavePointFaceHits[pointI].hit())
273  {
274  // Grab the hit point
275  projectedSlavePoints[pointI] =
276  masterLocalFaces
277  [slavePointFaceHits[pointI].hitObject()].ray
278  (
279  slaveLocalPoints[pointI],
280  slavePointNormals[pointI],
281  masterLocalPoints,
282  projectionAlgo_
283  ).hitPoint();
284  }
285  else
286  {
287  // Grab the nearest point on the face (edge)
288  pointHit missAdjust =
289  masterLocalFaces[slavePointFaceHits[pointI].hitObject()].ray
290  (
291  slaveLocalPoints[pointI],
292  slavePointNormals[pointI],
293  masterLocalPoints,
294  projectionAlgo_
295  );
296 
297  const point nearPoint = missAdjust.missPoint();
298  const point missPoint =
299  slaveLocalPoints[pointI]
300  + slavePointNormals[pointI]*missAdjust.distance();
301 
302  // Calculate the tolerance
303  const scalar mergeTol =
304  integralAdjTol_*minSlavePointLength[pointI];
305 
306  // Adjust the hit
307  if (mag(nearPoint - missPoint) < mergeTol)
308  {
309  if (debug)
310  {
311  Pout << "a";
312  }
313 
314 // Pout<< "Moving slave point in integral adjustment "
315 // << pointI << " miss point: " << missPoint
316 // << " near point: " << nearPoint
317 // << " mergeTol: " << mergeTol
318 // << " dist: " << mag(nearPoint - missPoint) << endl;
319 
320  projectedSlavePoints[pointI] = nearPoint;
321 
322  slavePointFaceHits[pointI] =
323  objectHit(true, slavePointFaceHits[pointI].hitObject());
324 
325  nAdjustedPoints++;
326  }
327  else
328  {
329  projectedSlavePoints[pointI] = slaveLocalPoints[pointI];
330 
331  if (debug)
332  {
333  Pout << "n";
334  }
335  }
336  }
337  }
338 
339  if (debug)
340  {
341  Pout << " done." << endl;
342  }
343  }
344  else if (matchType_ == PARTIAL)
345  {
346  forAll (slavePointFaceHits, pointI)
347  {
348  if (slavePointFaceHits[pointI].hit())
349  {
350  // Grab the hit point
351  projectedSlavePoints[pointI] =
352  masterLocalFaces
353  [slavePointFaceHits[pointI].hitObject()].ray
354  (
355  slaveLocalPoints[pointI],
356  slavePointNormals[pointI],
357  masterLocalPoints,
358  projectionAlgo_
359  ).hitPoint();
360  }
361  else
362  {
363  // The point remains where it started from
364  projectedSlavePoints[pointI] = slaveLocalPoints[pointI];
365  }
366  }
367  }
368  else
369  {
371  (
372  "bool slidingInterface::projectPoints() const"
373  ) << "Problem in point projection. Unknown sliding match type "
374  << " for object " << name()
375  << abort(FatalError);
376  }
377 
378  if (debug)
379  {
380  Pout<< "Number of adjusted points in projection: "
381  << nAdjustedPoints << endl;
382 
383  // Check for zero-length edges in slave projection
384  scalar minEdgeLength = GREAT;
385  scalar el = 0;
386  label nShortEdges = 0;
387 
388  forAll (slaveEdges, edgeI)
389  {
390  el = slaveEdges[edgeI].mag(projectedSlavePoints);
391 
392  if (el < SMALL)
393  {
394  Pout<< "Point projection problems for edge: "
395  << slaveEdges[edgeI] << ". Length = " << el
396  << endl;
397 
398  nShortEdges++;
399  }
400 
401  minEdgeLength = min(minEdgeLength, el);
402  }
403 
404  if (nShortEdges > 0)
405  {
407  (
408  "bool slidingInterface::projectPoints() const"
409  ) << "Problem in point projection. " << nShortEdges
410  << " short projected edges "
411  << "after adjustment for object " << name()
412  << abort(FatalError);
413  }
414  else
415  {
416  Pout << " ... projection OK." << endl;
417  }
418  }
419 // scalarField magDiffs(mag(slaveLocalPoints - projectedSlavePoints));
420 // Pout<< "Slave point face hits: " << slavePointFaceHits << nl
421 // << "slave points: " << slaveLocalPoints << nl
422 // << "Projected slave points: " << projectedSlavePoints
423 // << "diffs: " << magDiffs << endl;
424 
425  // Merge projected points against master points
426 
427  labelList slavePointPointHits(slaveLocalPoints.size(), -1);
428  labelList masterPointPointHits(masterLocalPoints.size(), -1);
429 
430  // Go through all the slave points and compare them against all the points
431  // in the master face they hit. If the distance between the projected point
432  // and any of the master face points is smaller than both tolerances,
433  // merge the projected point by:
434  // 1) adjusting the projected point to coincide with the
435  // master point it merges with
436  // 2) remembering the hit master point index in slavePointPointHits
437  // 3) resetting the hit face to -1
438  // 4) If a master point has been hit directly, it cannot be merged
439  // into the edge. Mark it as used in the reverse list
440 
441  label nMergedPoints = 0;
442 
443  forAll (projectedSlavePoints, pointI)
444  {
445  if (slavePointFaceHits[pointI].hit())
446  {
447  // Taking a non-const reference so the point can be adjusted
448  point& curPoint = projectedSlavePoints[pointI];
449 
450  // Get the hit face
451  const face& hitFace =
452  masterLocalFaces[slavePointFaceHits[pointI].hitObject()];
453 
454  label mergePoint = -1;
455  scalar mergeDist = GREAT;
456 
457  // Try all point before deciding on best fit.
458  forAll (hitFace, hitPointI)
459  {
460  scalar dist =
461  mag(masterLocalPoints[hitFace[hitPointI]] - curPoint);
462 
463  // Calculate the tolerance
464  const scalar mergeTol =
465  pointMergeTol_*
466  min
467  (
468  minSlavePointLength[pointI],
469  minMasterPointLength[hitFace[hitPointI]]
470  );
471 
472  if (dist < mergeTol && dist < mergeDist)
473  {
474  mergePoint = hitFace[hitPointI];
475  mergeDist = dist;
476 
477 // Pout<< "Merging slave point "
478 // << slavePatch.meshPoints()[pointI] << " at "
479 // << slaveLocalPoints[pointI] << " with master "
480 // << masterPatch.meshPoints()[mergePoint] << " at "
481 // << masterLocalPoints[mergePoint]
482 // << ". dist: " << mergeDist
483 // << " mergeTol: " << mergeTol << endl;
484  }
485  }
486 
487  if (mergePoint > -1)
488  {
489  // Point is to be merged with master point
490  nMergedPoints++;
491 
492  slavePointPointHits[pointI] = mergePoint;
493  curPoint = masterLocalPoints[mergePoint];
494  masterPointPointHits[mergePoint] = pointI;
495  }
496  }
497  }
498 
499 // Pout<< "slavePointPointHits: " << slavePointPointHits << nl
500 // << "masterPointPointHits: " << masterPointPointHits << endl;
501 
502  if (debug)
503  {
504  // Check for zero-length edges in slave projection
505  scalar minEdgeLength = GREAT;
506  scalar el = 0;
507 
508  forAll (slaveEdges, edgeI)
509  {
510  el = slaveEdges[edgeI].mag(projectedSlavePoints);
511 
512  if (el < SMALL)
513  {
514  Pout<< "Point projection problems for edge: "
515  << slaveEdges[edgeI] << ". Length = " << el
516  << endl;
517  }
518 
519  minEdgeLength = min(minEdgeLength, el);
520  }
521 
522  if (minEdgeLength < SMALL)
523  {
525  (
526  "bool slidingInterface::projectPoints() const"
527  ) << "Problem in point projection. Short projected edge "
528  << " after point merge for object " << name()
529  << abort(FatalError);
530  }
531  else
532  {
533  Pout << " ... point merge OK." << endl;
534  }
535  }
536 
537  // Merge projected points against master edges
538 
539  labelList slavePointEdgeHits(slaveLocalPoints.size(), -1);
540 
541  label nMovedPoints = 0;
542 
543  forAll (projectedSlavePoints, pointI)
544  {
545  // Eliminate the points merged into points
546  if (slavePointPointHits[pointI] < 0)
547  {
548  // Get current point position
549  point& curPoint = projectedSlavePoints[pointI];
550 
551  // Get the hit face
552  const labelList& hitFaceEdges =
553  masterFaceEdges[slavePointFaceHits[pointI].hitObject()];
554 
555  // Calculate the tolerance
556  const scalar mergeLength =
557  min
558  (
559  minSlavePointLength[pointI],
560  minMasterFaceLength[slavePointFaceHits[pointI].hitObject()]
561  );
562 
563  const scalar mergeTol = pointMergeTol_*mergeLength;
564 
565  scalar minDistance = GREAT;
566 
567  forAll (hitFaceEdges, edgeI)
568  {
569  const edge& curEdge = masterEdges[hitFaceEdges[edgeI]];
570 
571  pointHit edgeHit =
572  curEdge.line(masterLocalPoints).nearestDist(curPoint);
573 
574  if (edgeHit.hit())
575  {
576  scalar dist =
577  mag(edgeHit.hitPoint() - projectedSlavePoints[pointI]);
578 
579  if (dist < mergeTol && dist < minDistance)
580  {
581  // Point is to be moved onto master edge
582  nMovedPoints++;
583 
584  slavePointEdgeHits[pointI] = hitFaceEdges[edgeI];
585  projectedSlavePoints[pointI] = edgeHit.hitPoint();
586 
587  minDistance = dist;
588 
589 // Pout<< "Moving slave point "
590 // << slavePatch.meshPoints()[pointI]
591 // << " (" << pointI
592 // << ") at " << slaveLocalPoints[pointI]
593 // << " onto master edge " << hitFaceEdges[edgeI]
594 // << " or ("
595 // << masterLocalPoints[curEdge.start()]
596 // << masterLocalPoints[curEdge.end()]
597 // << ") hit: " << edgeHit.hitPoint()
598 // << ". dist: " << dist
599 // << " mergeTol: " << mergeTol << endl;
600  }
601  }
602  } // end of hit face edges
603 
604  if (slavePointEdgeHits[pointI] > -1)
605  {
606  // Projected slave point has moved. Re-attempt merge with
607  // master points of the edge
608  point& curPoint = projectedSlavePoints[pointI];
609 
610  const edge& hitMasterEdge =
611  masterEdges[slavePointEdgeHits[pointI]];
612 
613  label mergePoint = -1;
614  scalar mergeDist = GREAT;
615 
616  forAll (hitMasterEdge, hmeI)
617  {
618  scalar hmeDist =
619  mag(masterLocalPoints[hitMasterEdge[hmeI]] - curPoint);
620 
621  // Calculate the tolerance
622  const scalar mergeTol =
623  pointMergeTol_*
624  min
625  (
626  minSlavePointLength[pointI],
627  minMasterPointLength[hitMasterEdge[hmeI]]
628  );
629 
630  if (hmeDist < mergeTol && hmeDist < mergeDist)
631  {
632  mergePoint = hitMasterEdge[hmeI];
633  mergeDist = hmeDist;
634 
635 // Pout<< "Merging slave point; SECOND TRY "
636 // << slavePatch.meshPoints()[pointI] << " local "
637 // << pointI << " at "
638 // << slaveLocalPoints[pointI] << " with master "
639 // << masterPatch.meshPoints()[mergePoint] << " at "
640 // << masterLocalPoints[mergePoint]
641 // << ". hmeDist: " << mergeDist
642 // << " mergeTol: " << mergeTol << endl;
643  }
644  }
645 
646  if (mergePoint > -1)
647  {
648  // Point is to be merged with master point
649  slavePointPointHits[pointI] = mergePoint;
650  curPoint = masterLocalPoints[mergePoint];
651  masterPointPointHits[mergePoint] = pointI;
652 
653  slavePointFaceHits[pointI] =
654  objectHit(true, slavePointFaceHits[pointI].hitObject());
655 
656 
657  // Disable edge merge
658  slavePointEdgeHits[pointI] = -1;
659  }
660  }
661  }
662  }
663 
664  if (debug)
665  {
666  Pout<< "Number of merged master points: " << nMergedPoints << nl
667  << "Number of adjusted slave points: " << nMovedPoints << endl;
668 
669  // Check for zero-length edges in slave projection
670  scalar minEdgeLength = GREAT;
671  scalar el = 0;
672 
673  forAll (slaveEdges, edgeI)
674  {
675  el = slaveEdges[edgeI].mag(projectedSlavePoints);
676 
677  if (el < SMALL)
678  {
679  Pout<< "Point projection problems for edge: "
680  << slaveEdges[edgeI] << ". Length = " << el
681  << endl;
682  }
683 
684  minEdgeLength = min(minEdgeLength, el);
685  }
686 
687  if (minEdgeLength < SMALL)
688  {
690  (
691  "bool slidingInterface::projectPoints() const"
692  ) << "Problem in point projection. Short projected edge "
693  << " after correction for object " << name()
694  << abort(FatalError);
695  }
696  }
697 
698 // Pout << "slavePointEdgeHits: " << slavePointEdgeHits << endl;
699 
700  // Insert the master points into closest slave edge if appropriate
701 
702  // Algorithm:
703  // The face potentially interacts with all the points of the
704  // faces covering the path between its two ends. Since we are
705  // examining an arbitrary shadow of the edge on a non-Euclidian
706  // surface, it is typically quite hard to do a geometric point
707  // search (under a shadow of a straight line). Therefore, the
708  // search will be done topologically.
709  //
710  // I) Point collection
711  // -------------------
712  // 1) Grab the master faces to which the end points of the edge
713  // have been projected.
714  // 2) Starting from the face containing the edge start, grab all
715  // its points and put them into the point lookup map. Put the
716  // face onto the face lookup map.
717  // 3) If the face of the end point is on the face lookup, complete
718  // the point collection step (this is checked during insertion.
719  // 4) Start a new round of insertion. Visit all faces in the face
720  // lookup and add their neighbours which are not already on the
721  // map. Every time the new neighbour is found, add its points to
722  // the point lookup. If the face of the end point is inserted,
723  // continue with the current roundof insertion and stop the
724  // algorithm.
725  //
726  // II) Point check
727  // ---------------
728  // Grab all the points from the point collection and check them
729  // against the current edge. If the edge-to-point distance is
730  // smaller than the smallest tolerance in the game (min of
731  // master point tolerance and left and right slave face around
732  // the edge tolerance) and the nearest point hits the edge, the
733  // master point will break the slave edge. Check the actual
734  // distance of the original master position from the edge. If
735  // the distance is smaller than a fraction of slave edge
736  // length, the hit is considered valid. Store the slave edge
737  // index for every master point.
738 
739  labelList masterPointEdgeHits(masterLocalPoints.size(), -1);
740  scalarField masterPointEdgeDist(masterLocalPoints.size(), GREAT);
741 
742  // Note. "Processing slave edges" code is repeated twice in the
743  // sliding intergace coupling in order to allow the point
744  // projection to be done separately from the actual cutting.
745  // Please change consistently with coupleSlidingInterface.C
746  //
747 
748  if (debug)
749  {
750  Pout << "Processing slave edges " << endl;
751  }
752 
753  // Create a map of faces the edge can interact with
754  labelHashSet curFaceMap
755  (
756  nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
757  );
758 
760 
761  forAll (slaveEdges, edgeI)
762  {
763  const edge& curEdge = slaveEdges[edgeI];
764 
765  //HJ: check for all edges even if both ends have missed
766  // Experimental.
767 // if
768 // (
769 // slavePointFaceHits[curEdge.start()].hit()
770 // || slavePointFaceHits[curEdge.end()].hit()
771 // )
772  {
773  // Clear the maps
774  curFaceMap.clear();
775  addedFaces.clear();
776 
777  // Grab the faces for start and end points
778  const label startFace =
779  slavePointFaceHits[curEdge.start()].hitObject();
780  const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
781 
782  // Insert the start face into the list
783  curFaceMap.insert(startFace);
784  addedFaces.insert(startFace);
785 
786 // Pout << "Doing edge " << edgeI << " or " << curEdge << " start: " << slavePointFaceHits[curEdge.start()].hitObject() << " end " << slavePointFaceHits[curEdge.end()].hitObject() << endl;
787  // If the end face is on the list, the face collection is finished
788  label nSweeps = 0;
789  bool completed = false;
790 
791  while (nSweeps < edgeFaceEscapeLimit_)
792  {
793  nSweeps++;
794 
795  if (addedFaces.found(endFace))
796  {
797  completed = true;
798  }
799 
800  // Add all face neighbours of face in the map
801  const labelList cf = addedFaces.toc();
802  addedFaces.clear();
803 
804  forAll (cf, cfI)
805  {
806  const labelList& curNbrs = masterFaceFaces[cf[cfI]];
807 
808  forAll (curNbrs, nbrI)
809  {
810  if (!curFaceMap.found(curNbrs[nbrI]))
811  {
812  curFaceMap.insert(curNbrs[nbrI]);
813  addedFaces.insert(curNbrs[nbrI]);
814  }
815  }
816  }
817 
818  if (completed) break;
819 
820  if (debug)
821  {
822  Pout << ".";
823  }
824  }
825 
826  if (!completed)
827  {
828  if (debug)
829  {
830  Pout << "x";
831  }
832 
833  // It is impossible to reach the end from the start, probably
834  // due to disconnected domain. Do search in opposite direction
835 
836  label nReverseSweeps = 0;
837 
838  addedFaces.clear();
839  curFaceMap.insert(endFace);
840  addedFaces.insert(endFace);
841 
842  while (nReverseSweeps < edgeFaceEscapeLimit_)
843  {
844  nReverseSweeps++;
845 
846  if (addedFaces.found(startFace))
847  {
848  completed = true;
849  }
850 
851  // Add all face neighbours of face in the map
852  const labelList cf = addedFaces.toc();
853  addedFaces.clear();
854 
855  forAll (cf, cfI)
856  {
857  const labelList& curNbrs = masterFaceFaces[cf[cfI]];
858 
859  forAll (curNbrs, nbrI)
860  {
861  if (!curFaceMap.found(curNbrs[nbrI]))
862  {
863  curFaceMap.insert(curNbrs[nbrI]);
864  addedFaces.insert(curNbrs[nbrI]);
865  }
866  }
867  }
868 
869  if (completed) break;
870 
871  if (debug)
872  {
873  Pout << ".";
874  }
875  }
876  }
877 
878  if (completed)
879  {
880  if (debug)
881  {
882  Pout << "+ ";
883  }
884  }
885  else
886  {
887  if (debug)
888  {
889  Pout << "z ";
890  }
891  }
892 
893  // Collect the points
894 
895  // Create a map of points the edge can interact with
896  labelHashSet curPointMap
897  (
898  nFacesPerSlaveEdge_*primitiveMesh::pointsPerFace_
899  );
900 
901  const labelList curFaces = curFaceMap.toc();
902 // Pout << "curFaces: " << curFaces << endl;
903  forAll (curFaces, faceI)
904  {
905  const face& f = masterLocalFaces[curFaces[faceI]];
906 
907  forAll (f, pointI)
908  {
909  curPointMap.insert(f[pointI]);
910  }
911  }
912 
913  const labelList curMasterPoints = curPointMap.toc();
914 
915  // Check all the points against the edge.
916 
917  linePointRef edgeLine = curEdge.line(projectedSlavePoints);
918 
919  const vector edgeVec = edgeLine.vec();
920  const scalar edgeMag = edgeLine.mag();
921 
922  // Calculate actual distance involved in projection. This
923  // is used to reject master points out of reach.
924  // Calculated as a combination of travel distance in projection and
925  // edge length
926  scalar slaveCatchDist =
927  edgeMasterCatchFraction_*edgeMag
928  + 0.5*
929  (
930  mag
931  (
932  projectedSlavePoints[curEdge.start()]
933  - slaveLocalPoints[curEdge.start()]
934  )
935  + mag
936  (
937  projectedSlavePoints[curEdge.end()]
938  - slaveLocalPoints[curEdge.end()]
939  )
940  );
941 
942  // The point merge distance needs to be measured in the
943  // plane of the slave edge. The unit vector is calculated
944  // as a cross product of the edge vector and the edge
945  // projection direction. When checking for the distance
946  // in plane, a minimum of the master-to-edge and
947  // projected-master-to-edge distance is used, to avoid
948  // problems with badly defined master planes. HJ,
949  // 17/Oct/2004
950  vector edgeNormalInPlane =
951  edgeVec
952  ^ (
953  slavePointNormals[curEdge.start()]
954  + slavePointNormals[curEdge.end()]
955  );
956 
957  edgeNormalInPlane /= mag(edgeNormalInPlane);
958 
959  forAll (curMasterPoints, pointI)
960  {
961  const label cmp = curMasterPoints[pointI];
962 
963  // Skip the current point if the edge start or end has
964  // been adjusted onto in
965  if
966  (
967  slavePointPointHits[curEdge.start()] == cmp
968  || slavePointPointHits[curEdge.end()] == cmp
969  || masterPointPointHits[cmp] > -1
970  )
971  {
972 // Pout << "Edge already snapped to point. Skipping." << endl;
973  continue;
974  }
975 
976  // Check if the point actually hits the edge within bounds
977  pointHit edgeLineHit =
978  edgeLine.nearestDist(masterLocalPoints[cmp]);
979 
980  if (edgeLineHit.hit())
981  {
982  // If the distance to the line is smaller than
983  // the tolerance the master point needs to be
984  // inserted into the edge
985 
986  // Strict checking of slave cut to avoid capturing
987  // end points.
988  scalar cutOnSlave =
989  ((edgeLineHit.hitPoint() - edgeLine.start()) & edgeVec)
990  /sqr(edgeMag);
991 
992  scalar distInEdgePlane =
993  min
994  (
995  edgeLineHit.distance(),
996  mag
997  (
998  (
999  masterLocalPoints[cmp]
1000  - edgeLineHit.hitPoint()
1001  )
1002  & edgeNormalInPlane
1003  )
1004  );
1005 // Pout << "master point: " << cmp
1006 // << " cutOnSlave " << cutOnSlave
1007 // << " distInEdgePlane: " << distInEdgePlane
1008 // << " tol1: " << pointMergeTol_*edgeMag
1009 // << " hitDist: " << edgeLineHit.distance()
1010 // << " tol2: " <<
1011 // min
1012 // (
1013 // slaveCatchDist,
1014 // masterPointEdgeDist[cmp]
1015 // ) << endl;
1016 
1017  // Not a point hit, check for edge
1018  if
1019  (
1020  cutOnSlave > edgeEndCutoffTol_
1021  && cutOnSlave < 1.0 - edgeEndCutoffTol_ // check edge cut
1022  && distInEdgePlane < edgeMergeTol_*edgeMag // merge plane
1023  && edgeLineHit.distance()
1024  < min
1025  (
1026  slaveCatchDist,
1027  masterPointEdgeDist[cmp]
1028  )
1029  )
1030  {
1031  if (debug)
1032  {
1033  if (masterPointEdgeHits[cmp] == -1)
1034  {
1035  // First hit
1036  Pout << "m";
1037  }
1038  else
1039  {
1040  // Repeat hit
1041  Pout << "M";
1042  }
1043  }
1044 
1045  // Snap to point onto edge
1046  masterPointEdgeHits[cmp] = edgeI;
1047  masterPointEdgeDist[cmp] = edgeLineHit.distance();
1048 
1049 // Pout<< "Inserting master point "
1050 // << masterPatch.meshPoints()[cmp]
1051 // << " (" << cmp
1052 // << ") at " << masterLocalPoints[cmp]
1053 // << " into slave edge " << edgeI
1054 // << " " << curEdge
1055 // << " cutOnSlave: " << cutOnSlave
1056 // << " distInEdgePlane: " << distInEdgePlane
1057 // << ". dist: " << edgeLineHit.distance()
1058 // << " mergeTol: "
1059 // << edgeMergeTol_*edgeMag
1060 // << " other tol: " <<
1061 // min
1062 // (
1063 // slaveCatchDist,
1064 // masterPointEdgeDist[cmp]
1065 // )
1066 // << endl;
1067  }
1068  }
1069  }
1070  } // End if both ends missing
1071  } // End all slave edges
1072 
1073  if (debug)
1074  {
1075  Pout << endl;
1076  }
1077 // Pout << "masterPointEdgeHits: " << masterPointEdgeHits << endl;
1078 
1079  if (debug)
1080  {
1081  Pout<< "bool slidingInterface::projectPoints() for object "
1082  << name() << " : "
1083  << "Finished projecting points. Topology = ";
1084  }
1085 
1086 // Pout<< "slavePointPointHits: " << slavePointPointHits << nl
1087 // << "slavePointEdgeHits: " << slavePointEdgeHits << nl
1088 // << "slavePointFaceHits: " << slavePointFaceHits << nl
1089 // << "masterPointEdgeHits: " << masterPointEdgeHits << endl;
1090 
1091  // The four lists:
1092  // - slavePointPointHits
1093  // - slavePointEdgeHits
1094  // - slavePointFaceHits
1095  // - masterPointEdgeHits
1096  // define how the two patches will be merged topologically.
1097  // If the lists have not changed since the last merge, the
1098  // sliding interface changes only geometrically and simple mesh
1099  // motion will suffice. Otherwise, a topological change is
1100  // required.
1101 
1102  // Compare the result with the current state
1103  if (!attached_)
1104  {
1105  // The mesh needs to change topologically
1106  trigger_ = true;
1107 
1108  // Store the addressing arrays and projected points
1109  deleteDemandDrivenData(slavePointPointHitsPtr_);
1110  slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1111 
1112  deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1113  slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1114 
1115  deleteDemandDrivenData(slavePointFaceHitsPtr_);
1116  slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1117 
1118  deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1119  masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1120 
1121  if (debug)
1122  {
1123  Pout << "(Detached interface) changing." << endl;
1124  }
1125  }
1126  else
1127  {
1128  // Compare the lists against the stored lists. If any of them
1129  // has changed, topological change will be executed.
1130  trigger_ = false;
1131 
1132  if
1133  (
1134  !slavePointPointHitsPtr_
1135  || !slavePointEdgeHitsPtr_
1136  || !slavePointFaceHitsPtr_
1137  || !masterPointEdgeHitsPtr_
1138  )
1139  {
1140  // Must be restart. Force topological change.
1141  deleteDemandDrivenData(slavePointPointHitsPtr_);
1142  slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1143 
1144  deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1145  slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1146 
1147  deleteDemandDrivenData(slavePointFaceHitsPtr_);
1148  slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1149 
1150  deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1151  masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1152 
1153  if (debug)
1154  {
1155  Pout << "(Attached interface restart) changing." << endl;
1156  }
1157 
1158  trigger_ = true;
1159  return trigger_;
1160  }
1161 
1162  if (slavePointPointHits != (*slavePointPointHitsPtr_))
1163  {
1164  if (debug)
1165  {
1166  Pout << "(Point projection) ";
1167  }
1168 
1169  trigger_ = true;
1170  }
1171 
1172  if (slavePointEdgeHits != (*slavePointEdgeHitsPtr_))
1173  {
1174  if (debug)
1175  {
1176  Pout << "(Edge projection) ";
1177  }
1178 
1179  trigger_ = true;
1180  }
1181 
1182  // Comparison for face hits needs to exclude points that have hit
1183  // another point or edge
1184  bool faceHitsDifferent = false;
1185 
1186  const List<objectHit>& oldPointFaceHits = *slavePointFaceHitsPtr_;
1187 
1188  forAll (slavePointFaceHits, pointI)
1189  {
1190  if
1191  (
1192  slavePointPointHits[pointI] < 0
1193  && slavePointEdgeHits[pointI] < 0
1194  )
1195  {
1196  // This is a straight face hit
1197  if (slavePointFaceHits[pointI] != oldPointFaceHits[pointI])
1198  {
1199  // Two lists are different
1200  faceHitsDifferent = true;
1201  break;
1202  }
1203  }
1204  }
1205 
1206  if (faceHitsDifferent)
1207  {
1208  if (debug)
1209  {
1210  Pout << "(Face projection) ";
1211  }
1212 
1213  trigger_ = true;
1214 
1215  }
1216 
1217  if (masterPointEdgeHits != (*masterPointEdgeHitsPtr_))
1218  {
1219  if (debug)
1220  {
1221  Pout << "(Master point projection) ";
1222  }
1223 
1224  trigger_ = true;
1225  }
1226 
1227 
1228  if (trigger_)
1229  {
1230  // Grab new data
1231  deleteDemandDrivenData(slavePointPointHitsPtr_);
1232  slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1233 
1234  deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1235  slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1236 
1237  deleteDemandDrivenData(slavePointFaceHitsPtr_);
1238  slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1239 
1240  deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1241  masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1242 
1243  if (debug)
1244  {
1245  Pout << "changing." << endl;
1246  }
1247  }
1248  else
1249  {
1250  if (debug)
1251  {
1252  Pout << "preserved." << endl;
1253  }
1254  }
1255  }
1256 
1257  return trigger_;
1258 }
1259 
1260 
1261 void Foam::slidingInterface::clearPointProjection() const
1262 {
1263  deleteDemandDrivenData(slavePointPointHitsPtr_);
1264  deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1265  deleteDemandDrivenData(slavePointFaceHitsPtr_);
1266  deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1267 
1268  deleteDemandDrivenData(projectedSlavePointsPtr_);
1269 }
1270 
1271 
1272 // ************************ vim: set sw=4 sts=4 et: ************************ //