FreeFOAM The Cross-Platform CFD Toolkit
splitMesh.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  splitMesh
26 
27 Description
28  Splits mesh by making internal faces external.
29 
30  Uses attachDetach. Generates a meshModifier of the form:
31 
32  Splitter
33  {
34  type attachDetach;
35  faceZoneName membraneFaces;
36  masterPatchName masterPatch;
37  slavePatchName slavePatch;
38  triggerTimes runTime.value();
39  }
40 
41  so will detach at the current time and split all faces in membraneFaces
42  into masterPatch and slavePatch (which have to be present but of 0 size)
43 
44 Usage
45 
46  - splitMesh [OPTIONS] <face set> <masterPatch> <slavePatch>
47 
48  @param <face set> \n
49  @todo Detailed description of argument.
50 
51  @param <masterPatch> \n
52  @todo Detailed description of argument.
53 
54  @param <slavePatch> \n
55  @todo Detailed description of argument.
56 
57  @param -overwrite \n
58  Overwrite existing data.
59 
60  @param -case <dir>\n
61  Case directory.
62 
63  @param -help \n
64  Display help message.
65 
66  @param -doc \n
67  Display Doxygen API documentation page for this application.
68 
69  @param -srcDoc \n
70  Display Doxygen source documentation page for this application.
71 
72 \*---------------------------------------------------------------------------*/
73 
74 #include <OpenFOAM/argList.H>
75 #include <OpenFOAM/polyMesh.H>
76 #include <OpenFOAM/Time.H>
78 #include <OpenFOAM/mapPolyMesh.H>
79 #include <meshTools/faceSet.H>
82 #include "regionSide.H"
84 
85 using namespace Foam;
86 
87 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
88 
89 // Find edge between points v0 and v1.
90 label findEdge(const primitiveMesh& mesh, const label v0, const label v1)
91 {
92  const labelList& pEdges = mesh.pointEdges()[v0];
93 
94  forAll(pEdges, pEdgeI)
95  {
96  label edgeI = pEdges[pEdgeI];
97 
98  const edge& e = mesh.edges()[edgeI];
99 
100  if (e.otherVertex(v0) == v1)
101  {
102  return edgeI;
103  }
104  }
105 
107  (
108  "findEdge(const primitiveMesh&, const label, const label)"
109  ) << "Cannot find edge between mesh points " << v0 << " and " << v1
110  << abort(FatalError);
111 
112  return -1;
113 }
114 
115 
116 // Checks whether patch present
117 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
118 {
119  label patchI = bMesh.findPatchID(name);
120 
121  if (patchI == -1)
122  {
123  FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
124  << "Cannot find patch " << name << endl
125  << "It should be present but of zero size" << endl
126  << "Valid patches are " << bMesh.names()
127  << exit(FatalError);
128  }
129 
130  if (bMesh[patchI].size())
131  {
132  FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
133  << "Patch " << name << " is present but non-zero size"
134  << exit(FatalError);
135  }
136 }
137 
138 
139 // Main program:
140 
141 int main(int argc, char *argv[])
142 {
144 
145  Foam::argList::validArgs.append("faceSet");
146  Foam::argList::validArgs.append("masterPatch");
147  Foam::argList::validArgs.append("slavePatch");
148  Foam::argList::validOptions.insert("overwrite", "");
149 
150 # include <OpenFOAM/setRootCase.H>
151 # include <OpenFOAM/createTime.H>
152  runTime.functionObjects().off();
153 # include <OpenFOAM/createPolyMesh.H>
154  const word oldInstance = mesh.pointsInstance();
155 
156  word setName(args.additionalArgs()[0]);
157  word masterPatch(args.additionalArgs()[1]);
158  word slavePatch(args.additionalArgs()[2]);
159  bool overwrite = args.optionFound("overwrite");
160 
161  // List of faces to split
162  faceSet facesSet(mesh, setName);
163 
164  Info<< "Read " << facesSet.size() << " faces to split" << endl << endl;
165 
166 
167  // Convert into labelList and check
168 
169  labelList faces(facesSet.toc());
170 
171  forAll(faces, i)
172  {
173  if (!mesh.isInternalFace(faces[i]))
174  {
176  << "Face " << faces[i] << " in faceSet " << setName
177  << " is not an internal face."
178  << exit(FatalError);
179  }
180  }
181 
182 
183  // Check for empty master and slave patches
184  checkPatch(mesh.boundaryMesh(), masterPatch);
185  checkPatch(mesh.boundaryMesh(), slavePatch);
186 
187 
188  //
189  // Find 'side' of all faces on splitregion. Uses regionSide which needs
190  // set of edges on side of this region. Use PrimitivePatch to find these.
191  //
192 
193  // Addressing on faces only in mesh vertices.
194  primitiveFacePatch fPatch
195  (
197  (
198  mesh.faces(),
199  faces
200  ),
201  mesh.points()
202  );
203 
204  const labelList& meshPoints = fPatch.meshPoints();
205 
206  // Mark all fence edges : edges on boundary of fPatch but not on boundary
207  // of polyMesh
208  labelHashSet fenceEdges(fPatch.size());
209 
210  const labelListList& allEdgeFaces = fPatch.edgeFaces();
211 
212  forAll(allEdgeFaces, patchEdgeI)
213  {
214  if (allEdgeFaces[patchEdgeI].size() == 1)
215  {
216  const edge& e = fPatch.edges()[patchEdgeI];
217 
218  label edgeI =
219  findEdge
220  (
221  mesh,
222  meshPoints[e.start()],
223  meshPoints[e.end()]
224  );
225 
226  fenceEdges.insert(edgeI);
227  }
228  }
229 
230  // Find sides reachable from 0th face of faceSet
231  label startFaceI = faces[0];
232 
233  regionSide regionInfo
234  (
235  mesh,
236  facesSet,
237  fenceEdges,
238  mesh.faceOwner()[startFaceI],
239  startFaceI
240  );
241 
242  // Determine flip state for all faces in faceSet
243  boolList zoneFlip(faces.size());
244 
245  forAll(faces, i)
246  {
247  zoneFlip[i] = !regionInfo.sideOwner().found(faces[i]);
248  }
249 
250 
251  // Create and add face zones and mesh modifiers
252  List<pointZone*> pz(0);
253  List<faceZone*> fz(1);
254  List<cellZone*> cz(0);
255 
256  fz[0] =
257  new faceZone
258  (
259  "membraneFaces",
260  faces,
261  zoneFlip,
262  0,
263  mesh.faceZones()
264  );
265 
266  Info << "Adding point and face zones" << endl;
267  mesh.addZones(pz, fz, cz);
268 
269  attachPolyTopoChanger splitter(mesh);
270  splitter.setSize(1);
271 
272  // Add the sliding interface mesh modifier to start working at current
273  // time
274  splitter.set
275  (
276  0,
277  new attachDetach
278  (
279  "Splitter",
280  0,
281  splitter,
282  "membraneFaces",
283  masterPatch,
284  slavePatch,
285  scalarField(1, runTime.value())
286  )
287  );
288 
289  Info<< nl << "Constructed topologyModifier:" << endl;
290  splitter[0].writeDict(Info);
291 
292  if (!overwrite)
293  {
294  runTime++;
295  }
296 
297  splitter.attach();
298 
299  if (overwrite)
300  {
301  mesh.setInstance(oldInstance);
302  }
303 
304  Info<< "Writing mesh to " << runTime.timeName() << endl;
305  if (!mesh.write())
306  {
308  << "Failed writing polyMesh."
309  << exit(FatalError);
310  }
311 
312  Info<< nl << "end" << endl;
313  return 0;
314 }
315 
316 
317 // ************************ vim: set sw=4 sts=4 et: ************************ //