FreeFOAM The Cross-Platform CFD Toolkit
surfaceToPatch.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 Application
25  surfaceToPatch
26 
27 Description
28  Reads surface and applies surface regioning to a mesh.
29 
30  Uses boundaryMesh to do the hard work.
31 
32 Usage
33 
34  - surfaceToPatch [OPTIONS] <Foam surface file>
35 
36  @param <Foam surface file> \n
37  @todo Detailed description of argument.
38 
39  @param -tol <fraction of mesh size>\n
40  Search toleracne.
41 
42  @param -faceSet <faceSet name>\n
43  Only apply to named face set.
44 
45  @param -case <dir>\n
46  Case directory.
47 
48  @param -help \n
49  Display help message.
50 
51  @param -doc \n
52  Display Doxygen API documentation page for this application.
53 
54  @param -srcDoc \n
55  Display Doxygen source documentation page for this application.
56 
57 \*---------------------------------------------------------------------------*/
58 
59 #include <OpenFOAM/argList.H>
60 #include <OpenFOAM/Time.H>
61 #include "boundaryMesh.H"
62 #include <OpenFOAM/polyMesh.H>
63 #include <meshTools/faceSet.H>
67 
68 using namespace Foam;
69 
70 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
71 
72 // Adds empty patch if not yet there. Returns patchID.
73 label addPatch(polyMesh& mesh, const word& patchName)
74 {
75  label patchI = mesh.boundaryMesh().findPatchID(patchName);
76 
77  if (patchI == -1)
78  {
79  const polyBoundaryMesh& patches = mesh.boundaryMesh();
80 
81  List<polyPatch*> newPatches(patches.size() + 1);
82 
83  patchI = 0;
84 
85  // Copy all old patches
86  forAll(patches, i)
87  {
88  const polyPatch& pp = patches[i];
89 
90  newPatches[patchI] =
91  pp.clone
92  (
93  patches,
94  patchI,
95  pp.size(),
96  pp.start()
97  ).ptr();
98 
99  patchI++;
100  }
101 
102  // Add zero-sized patch
103  newPatches[patchI] =
104  new polyPatch
105  (
106  patchName,
107  0,
108  mesh.nFaces(),
109  patchI,
110  patches
111  );
112 
113  mesh.removeBoundary();
114  mesh.addPatches(newPatches);
115 
116  Pout<< "Created patch " << patchName << " at " << patchI << endl;
117  }
118  else
119  {
120  Pout<< "Reusing patch " << patchName << " at " << patchI << endl;
121  }
122 
123  return patchI;
124 }
125 
126 
127 // Repatch single face. Return true if patch changed.
128 bool repatchFace
129 (
130  const polyMesh& mesh,
131  const boundaryMesh& bMesh,
132  const labelList& nearest,
133  const labelList& surfToMeshPatch,
134  const label faceI,
135  polyTopoChange& meshMod
136 )
137 {
138  bool changed = false;
139 
140  label bFaceI = faceI - mesh.nInternalFaces();
141 
142  if (nearest[bFaceI] != -1)
143  {
144  // Use boundary mesh one.
145  label bMeshPatchID = bMesh.whichPatch(nearest[bFaceI]);
146 
147  label patchID = surfToMeshPatch[bMeshPatchID];
148 
149  if (patchID != mesh.boundaryMesh().whichPatch(faceI))
150  {
151  label own = mesh.faceOwner()[faceI];
152 
153  label zoneID = mesh.faceZones().whichZone(faceI);
154 
155  bool zoneFlip = false;
156 
157  if (zoneID >= 0)
158  {
159  const faceZone& fZone = mesh.faceZones()[zoneID];
160 
161  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
162  }
163 
164  meshMod.setAction
165  (
167  (
168  mesh.faces()[faceI],// modified face
169  faceI, // label of face being modified
170  own, // owner
171  -1, // neighbour
172  false, // face flip
173  patchID, // patch for face
174  false, // remove from zone
175  zoneID, // zone for face
176  zoneFlip // face flip in zone
177  )
178  );
179 
180  changed = true;
181  }
182  }
183  else
184  {
185  changed = false;
186  }
187  return changed;
188 }
189 
190 
191 // Main program:
192 
193 int main(int argc, char *argv[])
194 {
196  argList::validArgs.append("surface file");
197  argList::validOptions.insert("faceSet", "faceSet name");
198  argList::validOptions.insert("tol", "fraction of mesh size");
199 
200 # include <OpenFOAM/setRootCase.H>
201 # include <OpenFOAM/createTime.H>
202 # include <OpenFOAM/createPolyMesh.H>
203 
204  fileName surfName(args.additionalArgs()[0]);
205 
206  Info<< "Reading surface from " << surfName << " ..." << endl;
207 
208  bool readSet = args.optionFound("faceSet");
209  word setName;
210 
211  if (readSet)
212  {
213  setName = args.option("faceSet");
214 
215  Info<< "Repatching only the faces in faceSet " << setName
216  << " according to nearest surface triangle ..." << endl;
217  }
218  else
219  {
220  Info<< "Patching all boundary faces according to nearest surface"
221  << " triangle ..." << endl;
222  }
223 
224  scalar searchTol = 1E-3;
225  args.optionReadIfPresent("tol", searchTol);
226 
227  // Get search box. Anything not within this box will not be considered.
228  const boundBox& meshBb = mesh.globalData().bb();
229  const vector searchSpan = searchTol * meshBb.span();
230 
231  Info<< "All boundary faces further away than " << searchTol
232  << " of mesh bounding box " << meshBb
233  << " will keep their patch label ..." << endl;
234 
235 
236  Info<< "Before patching:" << nl
237  << " patch\tsize" << endl;
238 
239  forAll(mesh.boundaryMesh(), patchI)
240  {
241  Info<< " " << mesh.boundaryMesh()[patchI].name() << '\t'
242  << mesh.boundaryMesh()[patchI].size() << endl;
243  }
244  Info<< endl;
245 
246 
247 
249 
250  // Load in the surface.
251  bMesh.readTriSurface(surfName);
252 
253  // Add all the boundaryMesh patches to the mesh.
254  const PtrList<boundaryPatch>& bPatches = bMesh.patches();
255 
256  // Map from surface patch ( = boundaryMesh patch) to polyMesh patch
257  labelList patchMap(bPatches.size());
258 
259  forAll(bPatches, i)
260  {
261  patchMap[i] = addPatch(mesh, bPatches[i].name());
262  }
263 
264  // Obtain nearest face in bMesh for each boundary face in mesh that
265  // is within search span.
266  // Note: should only determine for faceSet if working with that.
267  labelList nearest(bMesh.getNearest(mesh, searchSpan));
268 
269  {
270  // Dump unmatched faces to faceSet for debugging.
271  faceSet unmatchedFaces(mesh, "unmatchedFaces", nearest.size()/100);
272 
273  forAll(nearest, bFaceI)
274  {
275  if (nearest[bFaceI] == -1)
276  {
277  unmatchedFaces.insert(mesh.nInternalFaces() + bFaceI);
278  }
279  }
280 
281  Pout<< "Writing all " << unmatchedFaces.size()
282  << " unmatched faces to faceSet "
283  << unmatchedFaces.name()
284  << endl;
285 
286  unmatchedFaces.write();
287  }
288 
289 
290  polyTopoChange meshMod(mesh);
291 
292  label nChanged = 0;
293 
294  if (readSet)
295  {
296  faceSet faceLabels(mesh, setName);
297  Info<< "Read " << faceLabels.size() << " faces to repatch ..." << endl;
298 
299  forAllConstIter(faceSet, faceLabels, iter)
300  {
301  label faceI = iter.key();
302 
303  if (repatchFace(mesh, bMesh, nearest, patchMap, faceI, meshMod))
304  {
305  nChanged++;
306  }
307  }
308  }
309  else
310  {
311  forAll(nearest, bFaceI)
312  {
313  label faceI = mesh.nInternalFaces() + bFaceI;
314 
315  if (repatchFace(mesh, bMesh, nearest, patchMap, faceI, meshMod))
316  {
317  nChanged++;
318  }
319  }
320  }
321 
322  Pout<< "Changed " << nChanged << " boundary faces." << nl << endl;
323 
324  if (nChanged > 0)
325  {
326  meshMod.changeMesh(mesh, false);
327 
328  Info<< "After patching:" << nl
329  << " patch\tsize" << endl;
330 
331  forAll(mesh.boundaryMesh(), patchI)
332  {
333  Info<< " " << mesh.boundaryMesh()[patchI].name() << '\t'
334  << mesh.boundaryMesh()[patchI].size() << endl;
335  }
336  Info<< endl;
337 
338 
339  runTime++;
340 
341  // Write resulting mesh
342  Info << "Writing modified mesh to time " << runTime.value() << endl;
343  mesh.write();
344  }
345 
346 
347  Info<< "End\n" << endl;
348 
349  return 0;
350 }
351 
352 
353 // ************************ vim: set sw=4 sts=4 et: ************************ //