applications/utilities/mesh/manipulation/splitMesh/splitMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Splits mesh by making internal faces external. Uses attachDetach.
  27
  28     Generates a meshModifier of the form:
  29
  30     Splitter
  31     {
  32         type                       attachDetach;
  33         faceZoneName               membraneFaces;
  34         masterPatchName            masterPatch;
  35         slavePatchName             slavePatch;
  36         triggerTimes               runTime.value();
  37     }
  38
  39     so will detach at the current time and split all faces in membraneFaces
  40     into masterPatch and slavePatch (which have to be present but of 0 size)
  41
  42 \*---------------------------------------------------------------------------*/
  43
  44 #include "argList.H"
  45 #include "polyMesh.H"
  46 #include "Time.H"
  47 #include "polyTopoChange.H"
  48 #include "mapPolyMesh.H"
  49 #include "faceSet.H"
  50 #include "attachDetach.H"
  51 #include "polyTopoChanger.H"
  52 #include "regionSide.H"
  53 #include "primitiveFacePatch.H"
  54
  55 using namespace Foam;
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59 // Find edge between points v0 and v1.
  60 label findEdge(const primitiveMesh& mesh, const label v0, const label v1)
  61 {
  62     const labelList& pEdges = mesh.pointEdges()[v0];
  63
  64     forAll(pEdges, pEdgeI)
  65     {
  66         label edgeI = pEdges[pEdgeI];
  67
  68         const edge& e = mesh.edges()[edgeI];
  69
  70         if (e.otherVertex(v0) == v1)
  71         {
  72             return edgeI;
  73         }
  74     }
  75
  76     FatalErrorIn
  77     (
  78         "findEdge(const primitiveMesh&, const label, const label)"
  79     )   << "Cannot find edge between mesh points " << v0 << " and " << v1
  80         << abort(FatalError);
  81
  82     return -1;
  83 }
  84
  85
  86 // Checks whether patch present
  87 void checkPatch(const polyBoundaryMesh& bMesh, const word& name)
  88 {
  89     label patchI = bMesh.findPatchID(name);
  90
  91     if (patchI == -1)
  92     {
  93         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
  94             << "Cannot find patch " << name << endl
  95             << "It should be present but of zero size" << endl
  96             << "Valid patches are " << bMesh.names()
  97             << exit(FatalError);
  98     }
  99
 100     if (bMesh[patchI].size())
 101     {
 102         FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
 103             << "Patch " << name << " is present but non-zero size"
 104             << exit(FatalError);
 105     }
 106 }
 107
 108
 109 // Main program:
 110
 111 int main(int argc, char *argv[])
 112 {
 113     Foam::argList::noParallel();
 114
 115     Foam::argList::validArgs.append("faceSet");
 116     Foam::argList::validArgs.append("masterPatch");
 117     Foam::argList::validArgs.append("slavePatch");
 118     Foam::argList::validOptions.insert("overwrite", "");
 119
 120 #   include "setRootCase.H"
 121 #   include "createTime.H"
 122     runTime.functionObjects().off();
 123 #   include "createPolyMesh.H"
 124     const word oldInstance = mesh.pointsInstance();
 125
 126     word setName(args.additionalArgs()[0]);
 127     word masterPatch(args.additionalArgs()[1]);
 128     word slavePatch(args.additionalArgs()[2]);
 129     bool overwrite = args.optionFound("overwrite");
 130
 131     // List of faces to split
 132     faceSet facesSet(mesh, setName);
 133
 134     Info<< "Read " << facesSet.size() << " faces to split" << endl << endl;
 135
 136
 137     // Convert into labelList and check
 138
 139     labelList faces(facesSet.toc());
 140
 141     forAll(faces, i)
 142     {
 143         if (!mesh.isInternalFace(faces[i]))
 144         {
 145             FatalErrorIn(args.executable())
 146             << "Face " << faces[i] << " in faceSet " << setName
 147             << " is not an internal face."
 148             << exit(FatalError);
 149         }
 150     }
 151
 152
 153     // Check for empty master and slave patches
 154     checkPatch(mesh.boundaryMesh(), masterPatch);
 155     checkPatch(mesh.boundaryMesh(), slavePatch);
 156
 157
 158     //
 159     // Find 'side' of all faces on splitregion. Uses regionSide which needs
 160     // set of edges on side of this region. Use PrimitivePatch to find these.
 161     //
 162
 163     IndirectList<face> zoneFaces(mesh.faces(), faces);
 164
 165     // Calculation engine for set of faces in a mesh
 166     typedef PrimitivePatch<face, List, const pointField&> facePatch;
 167
 168
 169     // Addressing on faces only in mesh vertices.
 170     facePatch fPatch(zoneFaces(), mesh.points());
 171
 172     const labelList& meshPoints = fPatch.meshPoints();
 173
 174     // Mark all fence edges : edges on boundary of fPatch but not on boundary
 175     // of polyMesh
 176     labelHashSet fenceEdges(fPatch.size());
 177
 178     const labelListList& allEdgeFaces = fPatch.edgeFaces();
 179
 180     forAll(allEdgeFaces, patchEdgeI)
 181     {
 182         if (allEdgeFaces[patchEdgeI].size() == 1)
 183         {
 184             const edge& e = fPatch.edges()[patchEdgeI];
 185
 186             label edgeI =
 187                 findEdge
 188                 (
 189                     mesh,
 190                     meshPoints[e.start()],
 191                     meshPoints[e.end()]
 192                 );
 193
 194             fenceEdges.insert(edgeI);
 195         }
 196     }
 197
 198     // Find sides reachable from 0th face of faceSet
 199     label startFaceI = faces[0];
 200
 201     regionSide regionInfo
 202     (
 203         mesh,
 204         facesSet,
 205         fenceEdges,
 206         mesh.faceOwner()[startFaceI],
 207         startFaceI
 208     );
 209
 210     // Determine flip state for all faces in faceSet
 211     boolList zoneFlip(faces.size());
 212
 213     forAll(faces, i)
 214     {
 215         zoneFlip[i] = !regionInfo.sideOwner().found(faces[i]);
 216     }
 217
 218
 219     // Create and add face zones and mesh modifiers
 220     List<pointZone*> pz(0);
 221     List<faceZone*> fz(1);
 222     List<cellZone*> cz(0);
 223
 224     fz[0] =
 225         new faceZone
 226         (
 227             "membraneFaces",
 228             faces,
 229             zoneFlip,
 230             0,
 231             mesh.faceZones()
 232         );
 233
 234     Info << "Adding point and face zones" << endl;
 235     mesh.addZones(pz, fz, cz);
 236
 237     polyTopoChanger splitter(mesh);
 238     splitter.setSize(1);
 239
 240     // Add the sliding interface mesh modifier to start working at current
 241     // time
 242     splitter.set
 243     (
 244         0,
 245         new attachDetach
 246         (
 247             "Splitter",
 248             0,
 249             splitter,
 250             "membraneFaces",
 251             masterPatch,
 252             slavePatch,
 253             scalarField(1, runTime.value())
 254         )
 255     );
 256
 257     Info<< nl << "Constructed topologyModifier:" << endl;
 258     splitter[0].writeDict(Info);
 259
 260     if (!overwrite)
 261     {
 262         runTime++;
 263     }
 264
 265     splitter.changeMesh();
 266
 267     Info<< "Writing mesh to " << runTime.timeName() << endl;
 268     if (!mesh.write())
 269     {
 270         FatalErrorIn(args.executable())
 271             << "Failed writing polyMesh."
 272             << exit(FatalError);
 273     }
 274
 275     Info<< nl << "end" << endl;
 276     return 0;
 277 }
 278
 279
 280 // ************************************************************************* //