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