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
27 \*---------------------------------------------------------------------------*/
29 #include "enrichedPatch.H"
30 #include "DynamicList.H"
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 const Foam::label Foam::enrichedPatch::enrichedFaceRatio_ = 3;
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 void Foam::enrichedPatch::calcEnrichedFaces
41 const labelListList& pointsIntoMasterEdges,
42 const labelListList& pointsIntoSlaveEdges,
43 const pointField& projectedSlavePoints
46 if (enrichedFacesPtr_)
50 "void enrichedPatch::calcEnrichedFaces\n"
52 " const labelListList& pointsIntoMasterEdges,\n"
53 " const labelListList& pointsIntoSlaveEdges,\n"
54 " const pointField& projectedSlavePoints\n"
56 ) << "Enriched faces already calculated."
60 // Create a list of enriched faces
62 // 1) Grab the original face and start from point zero.
63 // 2) If the point has been merged away, grab the merge label;
64 // otherwise, keep the original label.
65 // 3) Go to the next edge. Collect all the points to be added along
66 // the edge; order them in the necessary direction and insert onto the
68 // 4) Grab the next point and return on step 2.
69 enrichedFacesPtr_ = new faceList(masterPatch_.size() + slavePatch_.size());
70 faceList& enrichedFaces = *enrichedFacesPtr_;
72 label nEnrichedFaces = 0;
74 const pointField& masterLocalPoints = masterPatch_.localPoints();
75 const faceList& masterLocalFaces = masterPatch_.localFaces();
76 const labelListList& masterFaceEdges = masterPatch_.faceEdges();
78 const faceList& slaveLocalFaces = slavePatch_.localFaces();
79 const labelListList& slaveFaceEdges = slavePatch_.faceEdges();
81 // For correct functioning of the enrichedPatch class, the slave
82 // faces need to be inserted first. See comments in
85 // Get reference to the point merge map
86 const Map<label>& pmm = pointMergeMap();
88 // Add slave faces into the enriched faces list
90 forAll (slavePatch_, faceI)
92 const face oldFace = slavePatch_[faceI];
93 const face oldLocalFace = slaveLocalFaces[faceI];
94 // Info << "old slave face " << faceI << ": " << oldFace << endl;
95 const labelList& curEdges = slaveFaceEdges[faceI];
97 DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
99 // Note: The number of points and edges in a face is always identical
100 // so both can be done is the same loop
104 Map<label>::const_iterator mpIter =
105 pmm.find(oldFace[i]);
107 if (mpIter == pmm.end())
110 newFace.append(oldFace[i]);
112 // Add the projected point into the patch support
115 oldFace[i], // Global label of point
116 projectedSlavePoints[oldLocalFace[i]] // Projected position
122 newFace.append(mpIter());
124 // Add the projected point into the patch support
127 mpIter(), // Merged global label of point
128 projectedSlavePoints[oldLocalFace[i]] // Projected position
132 // Grab the edge points
134 const labelList& slavePointsOnEdge =
135 pointsIntoSlaveEdges[curEdges[i]];
136 // Info << "slavePointsOnEdge for " << curEdges[i] << ": " << slavePointsOnEdge << endl;
137 // If there are no points on the edge, skip everything
138 // If there is only one point, no need for sorting
139 if (slavePointsOnEdge.size() > 0)
141 // Sort edge points in order
142 scalarField edgePointWeights(slavePointsOnEdge.size());
143 const point& startPoint = projectedSlavePoints[oldLocalFace[i]];
146 projectedSlavePoints[oldLocalFace.nextLabel(i)]
149 scalar magSqrE = magSqr(e);
159 "void enrichedPatch::calcEnrichedFaces\n"
161 " const labelListList& pointsIntoMasterEdges,\n"
162 " const labelListList& pointsIntoSlaveEdges,\n"
163 " const pointField& projectedSlavePoints\n"
165 ) << "Zero length edge in slave patch for face " << i
166 << ". This is not allowed."
167 << abort(FatalError);
170 pointField slavePosOnEdge(slavePointsOnEdge.size());
172 forAll (slavePointsOnEdge, edgePointI)
174 slavePosOnEdge[edgePointI] =
175 pointMap().find(slavePointsOnEdge[edgePointI])();
177 edgePointWeights[edgePointI] =
178 (e & (slavePosOnEdge[edgePointI] - startPoint));
183 // Check weights: all new points should be on the edge
184 if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
188 "void enrichedPatch::calcEnrichedFaces\n"
190 " const labelListList& pointsIntoMasterEdges,\n"
191 " const labelListList& pointsIntoSlaveEdges,\n"
192 " const pointField& projectedSlavePoints\n"
194 ) << "Invalid point edge weights. Some of points are"
195 << " not on the edge for edge " << curEdges[i]
196 << " of face " << faceI << " in slave patch." << nl
197 << "Min weight: " << min(edgePointWeights)
198 << " Max weight: " << max(edgePointWeights)
199 << abort(FatalError);
203 // Go through the points and collect them based on
204 // weights from lower to higher. This gives the
205 // correct order of points along the edge.
206 for (label passI = 0; passI < edgePointWeights.size(); passI++)
208 // Max weight can only be one, so the sorting is
209 // done by elimination.
210 label nextPoint = -1;
213 forAll (edgePointWeights, wI)
215 if (edgePointWeights[wI] < dist)
217 dist = edgePointWeights[wI];
222 // Insert the next point and reset its weight to exclude it
224 newFace.append(slavePointsOnEdge[nextPoint]);
225 edgePointWeights[nextPoint] = GREAT;
227 // Add the point into patch support
230 slavePointsOnEdge[nextPoint],
231 slavePosOnEdge[nextPoint]
236 // Info << "New slave face " << faceI << ": " << newFace << endl;
238 // Add the new face to the list
239 enrichedFaces[nEnrichedFaces].transfer(newFace.shrink());
243 // Add master faces into the enriched faces list
245 forAll (masterPatch_, faceI)
247 const face& oldFace = masterPatch_[faceI];
248 const face& oldLocalFace = masterLocalFaces[faceI];
249 // Info << "old master face: " << oldFace << endl;
250 const labelList& curEdges = masterFaceEdges[faceI];
252 DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
254 // Note: The number of points and edges in a face is always identical
255 // so both can be done is the same loop
259 Map<label>::const_iterator mpIter =
260 pmm.find(oldFace[i]);
262 if (mpIter == pmm.end())
265 newFace.append(oldFace[i]);
267 // Add the point into patch support
271 masterLocalPoints[oldLocalFace[i]]
277 newFace.append(mpIter());
279 // Add the point into support
280 pointMap().insert(mpIter(), masterLocalPoints[oldLocalFace[i]]);
283 // Grab the edge points
285 const labelList& masterPointsOnEdge =
286 pointsIntoMasterEdges[curEdges[i]];
288 // If there are no points on the edge, skip everything
289 // If there is only one point, no need for sorting
290 if (masterPointsOnEdge.size() > 0)
292 // Sort edge points in order
293 scalarField edgePointWeights(masterPointsOnEdge.size());
294 const point& startPoint = masterLocalPoints[oldLocalFace[i]];
297 masterLocalPoints[oldLocalFace.nextLabel(i)]
300 scalar magSqrE = magSqr(e);
310 "void enrichedPatch::calcEnrichedFaces\n"
312 " const labelListList& pointsIntoMasterEdges,\n"
313 " const labelListList& pointsIntoSlaveEdges,\n"
314 " const pointField& projectedSlavePoints\n"
316 ) << "Zero length edge in master patch for face " << i
317 << ". This is not allowed."
318 << abort(FatalError);
321 pointField masterPosOnEdge(masterPointsOnEdge.size());
323 forAll (masterPointsOnEdge, edgePointI)
325 masterPosOnEdge[edgePointI] =
326 pointMap().find(masterPointsOnEdge[edgePointI])();
328 edgePointWeights[edgePointI] =
329 (e & (masterPosOnEdge[edgePointI] - startPoint));
334 // Check weights: all new points should be on the edge
335 if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
339 "void enrichedPatch::calcEnrichedFaces\n"
341 " const labelListList& pointsIntoMasterEdges,\n"
342 " const labelListList& pointsIntoSlaveEdges,\n"
343 " const pointField& projectedSlavePoints\n"
345 ) << "Invalid point edge weights. Some of points are"
346 << " not on the edge for edge " << curEdges[i]
347 << " of face " << faceI << " in master patch." << nl
348 << "Min weight: " << min(edgePointWeights)
349 << " Max weight: " << max(edgePointWeights)
350 << abort(FatalError);
354 // Go through the points and collect them based on
355 // weights from lower to higher. This gives the
356 // correct order of points along the edge.
357 for (label pass = 0; pass < edgePointWeights.size(); pass++)
359 // Max weight can only be one, so the sorting is
360 // done by elimination.
361 label nextPoint = -1;
364 forAll (edgePointWeights, wI)
366 if (edgePointWeights[wI] < dist)
368 dist = edgePointWeights[wI];
373 // Insert the next point and reset its weight to exclude it
375 newFace.append(masterPointsOnEdge[nextPoint]);
376 edgePointWeights[nextPoint] = GREAT;
378 // Add the point into patch support
381 masterPointsOnEdge[nextPoint],
382 masterPosOnEdge[nextPoint]
387 // Info << "New master face: " << newFace << endl;
389 // Add the new face to the list
390 enrichedFaces[nEnrichedFaces].transfer(newFace.shrink());
394 // Check the support for the enriched patch
399 Info<< "Enriched patch support OK. Slave faces: "
400 << slavePatch_.size() << " Master faces: "
401 << masterPatch_.size() << endl;
407 "void enrichedPatch::calcEnrichedFaces\n"
409 " const labelListList& pointsIntoMasterEdges,\n"
410 " const labelListList& pointsIntoSlaveEdges,\n"
411 " const pointField& projectedSlavePoints\n"
413 ) << "Error in enriched patch support"
414 << abort(FatalError);
420 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
422 const Foam::faceList& Foam::enrichedPatch::enrichedFaces() const
424 if (!enrichedFacesPtr_)
426 FatalErrorIn("const faceList& enrichedPatch::enrichedFaces() const")
427 << "Enriched faces not available yet. Please use "
428 << "void enrichedPatch::calcEnrichedFaces\n"
430 << " const labelListList& pointsIntoMasterEdges,\n"
431 << " const labelListList& pointsIntoSlaveEdges,\n"
432 << " const pointField& projectedSlavePoints\n"
434 << " before trying to access faces."
435 << abort(FatalError);
438 return *enrichedFacesPtr_;
442 // ************************************************************************* //