1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
26 Calculating cut faces of the enriched patch, together with the addressing
27 into master and slave patch.
29 \*---------------------------------------------------------------------------*/
31 #include "enrichedPatch.H"
33 #include "DynamicList.H"
34 #include "labelPair.H"
35 #include "primitiveMesh.H"
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 // If the cut face gets more than this number of points, it will be checked
41 const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 // Index of debug signs:
51 void Foam::enrichedPatch::calcCutFaces() const
53 if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
55 FatalErrorIn("void enrichedPatch::calcCutFaces() const")
56 << "Cut faces addressing already calculated."
60 const faceList& enFaces = enrichedFaces();
61 const labelList& mp = meshPoints();
62 const faceList& lf = localFaces();
63 const pointField& lp = localPoints();
64 const labelListList& pp = pointPoints();
65 // Pout << "enFaces: " << enFaces << endl;
66 // Pout << "lf: " << lf << endl;
67 // Pout << "lp: " << lp << endl;
68 // Pout << "pp: " << pp << endl;
69 const Map<labelList>& masterPointFaceAddr = masterPointFaces();
71 // Prepare the storage
72 DynamicList<face> cf(2*lf.size());
73 DynamicList<label> cfMaster(2*lf.size());
74 DynamicList<label> cfSlave(2*lf.size());
77 // Go through all the faces
78 // 1) For each face, start at point zero and grab the edge zero.
79 // Both points are added into the cut face. Mark the first edge
80 // as used and position the current point as the end of the first
81 // edge and previous point as point zero.
82 // 2) Grab all possible points from the current point. Excluding
83 // the previous point find the point giving the biggest right
84 // turn. Add the point to the face and mark the edges as used.
85 // Continue doing this until the face is complete, i.e. the next point
86 // to add is the first point of the face.
87 // 3) Once the facet is completed, register its owner from the face
88 // it has been created from (remember that the first lot of faces
89 // inserted into the enriched faces list are the slaves, to allow
90 // matching of the other side).
91 // 4) If the facet is created from an original slave face, go
92 // through the master patch and try to identify the master face
93 // this facet belongs to. This is recognised by the fact that the
94 // faces consists exclusively of the points of the master face and
95 // the points projected onto the face.
97 // Create a set of edge usage parameters
98 HashSet<edge, Hash<edge> > edgesUsedOnce(pp.size());
99 HashSet<edge, Hash<edge> > edgesUsedTwice
100 (pp.size()*primitiveMesh::edgesPerPoint_);
105 const face& curLocalFace = lf[faceI];
106 const face& curGlobalFace = enFaces[faceI];
108 // Pout<< "Doing face " << faceI << " local: " << curLocalFace << " or " << curGlobalFace << endl;
109 // if (faceI < slavePatch_.size())
111 // Pout<< "original slave: " << slavePatch_[faceI]
112 // << " local: " << slavePatch_.localFaces()[faceI] << endl;
116 // Pout<< "original master: "
117 // << masterPatch_[faceI - slavePatch_.size()] << " "
118 // << masterPatch_.localFaces()[faceI - slavePatch_.size()]
122 // pointField facePoints = curLocalFace.points(lp);
123 // forAll (curLocalFace, pointI)
125 // Pout << "v " << facePoints[pointI].x() << " "
126 // << facePoints[pointI].y() << " "
127 // << facePoints[pointI].z() << endl;
131 // Track the usage of face edges. When all edges are used, the
132 // face decomposition is complete.
133 // The seed edges include all the edges of the original face + all edges
134 // of other faces that have been used in the construction of the
135 // facet. Edges from other faces can be considered as
136 // internal to the current face if used only once.
138 // Track the edge usage to avoid duplicate faces and reset it to unused
139 boolList usedFaceEdges(curLocalFace.size(), false);
141 SLList<edge> edgeSeeds;
143 // Insert the edges of current face into the seed list.
144 edgeList cfe = curLocalFace.edges();
145 forAll (curLocalFace, edgeI)
147 edgeSeeds.append(cfe[edgeI]);
151 vector normal = curLocalFace.normal(lp);
152 normal /= mag(normal);
154 while (edgeSeeds.size() > 0)
156 // Pout << "edgeSeeds.size(): " << edgeSeeds.size() << endl;
157 const edge curEdge = edgeSeeds.removeHead();
159 // Locate the edge in current face
160 const label curEdgeWhich = curLocalFace.which(curEdge.start());
162 // Check if the edge is in current face and if it has been
163 // used already. If so, skip it
167 && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
170 // Edge is in the starting face.
171 // If unused, mark it as used; if used, skip it
172 if (usedFaceEdges[curEdgeWhich]) continue;
174 usedFaceEdges[curEdgeWhich] = true;
177 // If the edge has already been used twice, skip it
178 if (edgesUsedTwice.found(curEdge)) continue;
179 // Pout << "Trying new edge (" << mp[curEdge.start()] << ", " << mp[curEdge.end()] << ") seed: " << curEdge << " used: " << edgesUsedTwice.found(curEdge) << endl;
181 // Estimate the size of cut face as twice the size of original face
182 DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
183 DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
185 // Found unused edge.
186 label prevPointLabel = curEdge.start();
187 cutFaceGlobalPoints.append(mp[prevPointLabel]);
188 cutFaceLocalPoints.append(prevPointLabel);
189 // Pout << "prevPointLabel: " << mp[prevPointLabel] << endl;
190 // Grab current point and append it to the list
191 label curPointLabel = curEdge.end();
192 point curPoint = lp[curPointLabel];
193 cutFaceGlobalPoints.append(mp[curPointLabel]);
194 cutFaceLocalPoints.append(curPointLabel);
196 bool completedCutFace = false;
198 label faceSizeDebug = maxFaceSizeDebug_;
202 // Grab the next point options
203 // Pout << "curPointLabel: " << mp[curPointLabel] << endl;
204 const labelList& nextPoints = pp[curPointLabel];
205 // Pout << "nextPoints: " << IndirectList<label>(mp, nextPoints) << endl;
206 // Get the vector along the edge and the right vector
207 vector ahead = curPoint - lp[prevPointLabel];
208 ahead -= normal*(normal & ahead);
211 vector right = normal ^ ahead;
213 // Pout<< "normal: " << normal << " ahead: " << ahead << " right: " << right << endl;
215 scalar atanTurn = -GREAT;
216 label bestAtanPoint = -1;
218 forAll (nextPoints, nextI)
220 // Exclude the point we are coming from; there will always
221 // be more than one edge, so this is safe
222 if (nextPoints[nextI] != prevPointLabel)
224 // Pout << "cur point: " << curPoint << " trying for point: " << mp[nextPoints[nextI]] << " " << lp[nextPoints[nextI]];
225 vector newDir = lp[nextPoints[nextI]] - curPoint;
226 // Pout << " newDir: " << newDir << " mag: " << mag(newDir) << flush;
227 newDir -= normal*(normal & newDir);
228 scalar magNewDir = mag(newDir);
229 // Pout << " corrected: " << newDir << " mag: " << mag(newDir) << flush;
231 if (magNewDir < SMALL)
235 "void enrichedPatch::"
236 "calcCutFaces() const"
237 ) << "Zero length edge detected. Probable "
238 << "projection error: slave patch probably "
239 << "does not project onto master. "
240 << "Please switch on "
241 << "enriched patch debug for more info"
242 << abort(FatalError);
248 atan2(newDir & right, newDir & ahead);
250 // Pout << " atan: " << curAtanTurn << endl;
252 if (curAtanTurn > atanTurn)
254 bestAtanPoint = nextPoints[nextI];
255 atanTurn = curAtanTurn;
257 } // end of prev point skip
258 } // end of next point selection
259 // Pout<< " bestAtanPoint: " << bestAtanPoint << " or "
260 // << mp[bestAtanPoint] << endl;
261 // Selected next best point.
262 // Pout << "cutFaceGlobalPoints: " << cutFaceGlobalPoints << endl;
263 // Check if the edge about to be added has been used
264 // in the current face or twice in other faces. If
265 // so, the face is bad.
266 const label whichNextPoint = curLocalFace.which(curPointLabel);
271 && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
272 && usedFaceEdges[whichNextPoint]
275 // This edge is already used in current face
276 // face cannot be good; start on a new one
277 // Pout << "Double usage in current face, cannot be good" << endl;
278 completedCutFace = true;
282 if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
284 // This edge is already used -
285 // face cannot be good; start on a new one
286 completedCutFace = true;
287 // Pout << "Double usage elsewhere, cannot be good" << endl;
290 if (completedCutFace) continue;
292 // Insert the next best point into the list
293 if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
295 // Face is completed. Save it and move on to the next
297 completedCutFace = true;
301 Pout<< " local: " << cutFaceLocalPoints
302 << " one side: " << faceI;
307 cutFaceGlobal.transfer(cutFaceGlobalPoints.shrink());
309 cf.append(cutFaceGlobal);
312 cutFaceLocal.transfer(cutFaceLocalPoints.shrink());
313 // Pout << "\ncutFaceLocal: " << cutFaceLocal.points(lp) << endl;
314 // Go through all edges of the cut faces.
315 // If the edge corresponds to a starting face edge,
316 // mark the starting face edge as true
318 forAll (cutFaceLocal, cutI)
320 const edge curCutFaceEdge
323 cutFaceLocal.nextLabel(cutI)
326 // Increment the usage count using two hash sets
327 HashSet<edge, Hash<edge> >::iterator euoIter =
328 edgesUsedOnce.find(curCutFaceEdge);
330 if (euoIter == edgesUsedOnce.end())
332 // Pout << "Found edge not used before: "<< curCutFaceEdge << endl;
333 edgesUsedOnce.insert(curCutFaceEdge);
337 // Pout << "Found edge used once: " << curCutFaceEdge << endl;
338 edgesUsedOnce.erase(euoIter);
339 edgesUsedTwice.insert(curCutFaceEdge);
342 const label curCutFaceEdgeWhich = curLocalFace.which(curCutFaceEdge.start());
346 curCutFaceEdgeWhich > -1
347 && curLocalFace.nextLabel(curCutFaceEdgeWhich) == curCutFaceEdge.end()
350 // Found edge in original face
351 // Pout << "Found edge in orig face: " << curCutFaceEdge << ": " << curCutFaceEdgeWhich << endl;
352 usedFaceEdges[curCutFaceEdgeWhich] = true;
356 // Edge not in original face. Add it to seeds
357 // Pout << "Found new edge seed: " << curCutFaceEdge << endl;
358 edgeSeeds.append(curCutFaceEdge.reverseEdge());
363 // Find out what the other side is
366 // If the face is in the slave half of the
367 // enrichedFaces lists, it may be matched against
368 // the master face. It can be recognised by the
369 // fact that all its points belong to one master
370 // face. For this purpose:
371 // 1) Grab the master faces around the point zero
372 // of the cut face and collect all master faces
374 // 2) Loop through the rest of cut face points and
375 // if the point refers to one of the faces marked
376 // by point zero, increment its count.
377 // 3) When completed, the face whose count is
378 // equal to the number of points in the cut face
379 // is the other side. If this is not the case, there is no
380 // face on the other side.
382 if (faceI < slavePatch_.size())
384 Map<labelList>::const_iterator mpfAddrIter =
385 masterPointFaceAddr.find(cutFaceGlobal[0]);
387 bool otherSideFound = false;
389 if (mpfAddrIter != masterPointFaceAddr.end())
393 // Create the label count using point zero
394 const labelList& masterFacesOfPZero = mpfAddrIter();
396 labelList hits(masterFacesOfPZero.size(), 1);
401 pointI < cutFaceGlobal.size();
405 Map<labelList>::const_iterator
407 masterPointFaceAddr.find
409 cutFaceGlobal[pointI]
415 == masterPointFaceAddr.end()
418 // Point is off the master patch. Skip
423 const labelList& curMasterFaces =
426 // For every current face, try to find it in the
428 forAll (curMasterFaces, i)
430 forAll (masterFacesOfPZero, j)
435 == masterFacesOfPZero[j]
445 // If all point are found attempt matching
448 forAll (hits, pointI)
450 if (hits[pointI] == cutFaceGlobal.size())
453 otherSideFound = true;
457 masterFacesOfPZero[pointI]
460 cfSlave.append(faceI);
462 // Reverse the face such that it
463 // points out of the master patch
465 cf[cf.size() - 1].reverseFace();
469 Pout<< " other side: "
470 << masterFacesOfPZero[pointI]
474 } // end of for all hits
478 if (!otherSideFound || miss)
482 Pout << " solo slave A" << endl;
486 cfSlave.append(faceI);
491 // First point not in master patch
494 Pout << " solo slave B" << endl;
498 cfSlave.append(faceI);
505 Pout << " master side" << endl;
508 cfMaster.append(faceI - slavePatch_.size());
514 // Face not completed. Prepare for the next point search
515 prevPointLabel = curPointLabel;
516 curPointLabel = bestAtanPoint;
517 curPoint = lp[curPointLabel];
518 cutFaceGlobalPoints.append(mp[curPointLabel]);
519 cutFaceLocalPoints.append(curPointLabel);
521 if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
525 // Check for duplicate points in the face
526 forAll (cutFaceGlobalPoints, checkI)
530 label checkJ = checkI + 1;
531 checkJ < cutFaceGlobalPoints.size();
537 cutFaceGlobalPoints[checkI]
538 == cutFaceGlobalPoints[checkJ]
541 // Shrink local points for debugging
542 cutFaceLocalPoints.shrink();
546 if (faceI < slavePatch_.size())
548 origFace = slavePatch_[faceI];
550 slavePatch_.localFaces()[faceI];
556 [faceI - slavePatch_.size()];
559 masterPatch_.localFaces()
560 [faceI - slavePatch_.size()];
565 "void enrichedPatch::"
566 "calcCutFaces() const"
567 ) << "Duplicate point found in cut face. "
568 << "Error in the face cutting "
569 << "algorithm for global face "
570 << origFace << " local face "
571 << origFaceLocal << nl
572 << "Slave size: " << slavePatch_.size()
574 << masterPatch_.size()
575 << " index: " << faceI << ".\n"
576 << "Face: " << curGlobalFace << nl
577 << "Cut face: " << cutFaceGlobalPoints
578 << " local: " << cutFaceLocalPoints
580 << face(cutFaceLocalPoints).points(lp)
581 << abort(FatalError);
587 } while (!completedCutFace);
588 } // end of used edges
592 Pout << " Finished face " << faceI << endl;
595 } // end of local faces
597 // Re-pack the list into compact storage
598 cutFacesPtr_ = new faceList();
599 cutFacesPtr_->transfer(cf.shrink());
601 cutFaceMasterPtr_ = new labelList();
602 cutFaceMasterPtr_->transfer(cfMaster.shrink());
604 cutFaceSlavePtr_ = new labelList();
605 cutFaceSlavePtr_->transfer(cfSlave.shrink());
609 void Foam::enrichedPatch::clearCutFaces()
611 deleteDemandDrivenData(cutFacesPtr_);
613 deleteDemandDrivenData(cutFaceMasterPtr_);
614 deleteDemandDrivenData(cutFaceSlavePtr_);
618 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
620 const Foam::faceList& Foam::enrichedPatch::cutFaces() const
627 return *cutFacesPtr_;
631 const Foam::labelList& Foam::enrichedPatch::cutFaceMaster() const
633 if (!cutFaceMasterPtr_)
638 return *cutFaceMasterPtr_;
642 const Foam::labelList& Foam::enrichedPatch::cutFaceSlave() const
644 if (!cutFaceSlavePtr_)
649 return *cutFaceSlavePtr_;
653 // ************************************************************************* //