|
| 1 | +/*---------------------------------------------------------------------------*\ |
| 2 | + ========= | |
| 3 | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| 4 | + \\ / O peration | |
| 5 | + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation |
| 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 | +\*---------------------------------------------------------------------------*/ |
| 25 | + |
| 26 | +#include "fvCFD.H" |
| 27 | + |
| 28 | +#include "cellModeller.H" |
| 29 | + |
| 30 | +int main(int argc, char *argv[]) |
| 31 | +{ |
| 32 | + // set up the case |
| 33 | + #include "setRootCase.H" |
| 34 | + |
| 35 | + // create the run time object |
| 36 | + Info<< "Create time\n" << endl; |
| 37 | + Time runTime |
| 38 | + ( |
| 39 | + Time::controlDictName, |
| 40 | + args.rootPath(), |
| 41 | + args.caseName() |
| 42 | + ); |
| 43 | + |
| 44 | + // disable post-processing etc. |
| 45 | + runTime.functionObjects().off(); |
| 46 | + |
| 47 | + // --- |
| 48 | + // create cell types |
| 49 | + const cellModel& hex = *(cellModeller::lookup("hex")); |
| 50 | + const cellModel& prism = *(cellModeller::lookup("prism")); |
| 51 | + const cellModel& pyr = *(cellModeller::lookup("pyr")); |
| 52 | + const cellModel& tet = *(cellModeller::lookup("tet")); |
| 53 | + Info << "Created cell types" << endl; |
| 54 | + |
| 55 | + // create a point cloud between which the cells will stretch |
| 56 | + pointField points(12); |
| 57 | + points[ 0] = vector(0.0, 0, 0); |
| 58 | + points[ 1] = vector(0.5, 0, 0); |
| 59 | + points[ 2] = vector(1.0, 0, 0); |
| 60 | + points[ 3] = vector(0.0, 1, 0); |
| 61 | + points[ 4] = vector(0.5, 1, 0); |
| 62 | + points[ 5] = vector(1.0, 1, 0); |
| 63 | + points[ 6] = vector(0.0, 0, 0.1); |
| 64 | + points[ 7] = vector(0.5, 0, 0.1); |
| 65 | + points[ 8] = vector(1.0, 0, 0.1); |
| 66 | + points[ 9] = vector(0.0, 1, 0.1); |
| 67 | + points[10] = vector(0.5, 1, 0.1); |
| 68 | + points[11] = vector(1.0, 1, 0.1); |
| 69 | + Info << "Created points" << endl; |
| 70 | + |
| 71 | + // create the cells from a particular cell shape and list of points |
| 72 | + // TODO check if point ordering is important or not |
| 73 | + List<cellShape> cells; |
| 74 | + |
| 75 | + List<label> cellPoints(8); |
| 76 | + |
| 77 | + cellPoints[0] = 0; |
| 78 | + cellPoints[1] = 3; |
| 79 | + cellPoints[2] = 9; |
| 80 | + cellPoints[3] = 6; |
| 81 | + cellPoints[4] = 1; |
| 82 | + cellPoints[5] = 4; |
| 83 | + cellPoints[6] = 10; |
| 84 | + cellPoints[7] = 7; |
| 85 | + cells.append(cellShape(hex, cellPoints)); |
| 86 | + cellPoints[0] = 1; |
| 87 | + cellPoints[1] = 4; |
| 88 | + cellPoints[2] = 10; |
| 89 | + cellPoints[3] = 7; |
| 90 | + cellPoints[4] = 2; |
| 91 | + cellPoints[5] = 5; |
| 92 | + cellPoints[6] = 11; |
| 93 | + cellPoints[7] = 8; |
| 94 | + cells.append(cellShape(hex, cellPoints)); |
| 95 | + Info << "Created cells" << endl; |
| 96 | + |
| 97 | + // create patch definitions from lists of points defining individual faces |
| 98 | + // this is a list of List<List<label> > objects; each lowest-level list |
| 99 | + // contains indices of vertices making up the patch; these are grouped by |
| 100 | + // boundary index, i.e. patch number. The overall list holds all of those in |
| 101 | + // a single container. |
| 102 | + faceListList patchFaces; |
| 103 | + |
| 104 | + List<label> patchPoints(4); |
| 105 | + |
| 106 | + patchPoints[0] = 3; |
| 107 | + patchPoints[1] = 9; |
| 108 | + patchPoints[2] = 10; |
| 109 | + patchPoints[3] = 4; |
| 110 | + patchFaces.append(List<face>(1, face(patchPoints))); |
| 111 | + Info << "Created patches" << endl; |
| 112 | + |
| 113 | + // patch names |
| 114 | + List<word> boundaryPatchNames; |
| 115 | + boundaryPatchNames.append("side0"); |
| 116 | + |
| 117 | + // types and physical types |
| 118 | + // TODO figure out how these are supposed to look really, now just use defaults |
| 119 | + wordList boundaryPatchTypes(patchFaces.size(), polyPatch::typeName); |
| 120 | + wordList boundaryPatchPhysicalTypes |
| 121 | + ( |
| 122 | + patchFaces.size(), |
| 123 | + polyPatch::typeName |
| 124 | + ); |
| 125 | + |
| 126 | + // default values for other fields |
| 127 | + word regionName = polyMesh::defaultRegion; |
| 128 | + word defaultFacesName = "defaultFaces"; |
| 129 | + word defaultFacesType = emptyPolyPatch::typeName; |
| 130 | + Info << "Created physical and default values" << endl; |
| 131 | + |
| 132 | + // create the mesh |
| 133 | + polyMesh mesh |
| 134 | + ( |
| 135 | + IOobject |
| 136 | + ( |
| 137 | + regionName, |
| 138 | + runTime.constant(), |
| 139 | + runTime |
| 140 | + ), |
| 141 | + xferMove(points), |
| 142 | + cells, |
| 143 | + patchFaces, |
| 144 | + boundaryPatchNames, |
| 145 | + boundaryPatchTypes, |
| 146 | + defaultFacesName, |
| 147 | + defaultFacesType, |
| 148 | + boundaryPatchPhysicalTypes |
| 149 | + ); |
| 150 | + Info << "Created the mesh object" << endl; |
| 151 | + |
| 152 | + // --- |
| 153 | + // Write the grid |
| 154 | + Info << nl << "Writing extruded mesh to time = " << runTime.timeName() << nl << endl; |
| 155 | + mesh.write(); |
| 156 | + |
| 157 | +/* |
| 158 | + // These two create the time system (instance called runTime) and fvMesh (instance called mesh). |
| 159 | + #include "createTime.H" |
| 160 | + #include "createMesh.H" |
| 161 | +
|
| 162 | + // runTime and mesh are instances of objects (or classes). |
| 163 | + // If you are not familiar with what a class or object is, it is HIGHLY RECOMMENDED you visit this |
| 164 | + // website and only come back once you've read everything about classes, inheritance and polymorphism: |
| 165 | + // http://www.cplusplus.com/doc/tutorial/classes/ |
| 166 | + // Note how the next lines call functions .timeName(), .C() and .Cf() implemented in the objects. |
| 167 | + // It is also important to realise that mesh.C() and .Cf() return vector fields denoting centres of each |
| 168 | + // cell and internal face. |
| 169 | + // Calling the mesh.C().size() method therefore yields the total size of the mesh. |
| 170 | + Info << "Hello there, the most recent time folder found is " << runTime.timeName() << nl |
| 171 | + << "The mesh has " << mesh.C().size() << " cells and " << mesh.Cf().size() |
| 172 | + << " internal faces in it. Wubalubadubdub!" << nl << endl; |
| 173 | +
|
| 174 | + // It's possible to iterate over every cell in a standard C++ for loop |
| 175 | + for (label cellI = 0; cellI < mesh.C().size(); cellI++) |
| 176 | + if (cellI%20 == 0) // only show every twentieth cell not to spam the screen too much |
| 177 | + Info << "Cell " << cellI << " with centre at " << mesh.C()[cellI] << endl; |
| 178 | + Info << endl; // spacer |
| 179 | +
|
| 180 | + // Each cell is constructed of faces - these may either be internal or constitute a |
| 181 | + // boundary, or a patch in OpenFOAM terms; internal faces have an owner cell |
| 182 | + // and a neighbour. |
| 183 | + for (label faceI = 0; faceI < mesh.owner().size(); faceI++) |
| 184 | + if (faceI%40 == 0) |
| 185 | + Info << "Internal face " << faceI << " with centre at " << mesh.Cf()[faceI] |
| 186 | + << " with owner cell " << mesh.owner()[faceI] |
| 187 | + << " and neighbour " << mesh.neighbour()[faceI] << endl; |
| 188 | + Info << endl; |
| 189 | +
|
| 190 | + // Boundary conditions may be accessed through the boundaryMesh object. |
| 191 | + // In reality, each boundary face is also included in the constant/polyMesh/faces |
| 192 | + // description. But, in that file, the internal faces are defined first. |
| 193 | + // In addition, the constant/polyMesh/boundary file defines the starting faceI |
| 194 | + // indices from which boundary face definitions start. |
| 195 | + // OpenFOAM also provides a macro definition for for loops over all entries |
| 196 | + // in a field or a list, which saves up on the amount of typing. |
| 197 | + forAll(mesh.boundaryMesh(), patchI) |
| 198 | + Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with " |
| 199 | + << mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face " |
| 200 | + << mesh.boundary()[patchI].start() << endl; |
| 201 | + Info << endl; |
| 202 | +
|
| 203 | + // Faces adjacent to boundaries may be accessed as follows. |
| 204 | + // Also, a useful thing to know about a face is its normal vector and face area. |
| 205 | + label patchFaceI(0); |
| 206 | + forAll(mesh.boundaryMesh(), patchI) |
| 207 | + Info << "Patch " << patchI << " has its face " << patchFaceI << " adjacent to cell " |
| 208 | + << mesh.boundary()[patchI].patch().faceCells()[patchFaceI] |
| 209 | + << ". It has normal vector " << mesh.boundary()[patchI].Sf()[patchFaceI] |
| 210 | + << " and surface area " << mag(mesh.boundary()[patchI].Sf()[patchFaceI]) |
| 211 | + << endl; |
| 212 | + Info << endl; |
| 213 | +
|
| 214 | + // For internal faces, method .Sf() can be called directly on the mesh instance. |
| 215 | + // Moreover, there is a shorthand method .magSf() which returns the surface area |
| 216 | + // as a scalar. |
| 217 | + // For internal faces, the normal vector points from the owner to the neighbour |
| 218 | + // and the owner has a smaller cellI index than the neighbour. For boundary faces, |
| 219 | + // the normals always point outside of the domain (they have "imaginary" neighbours |
| 220 | + // which do not exist). |
| 221 | +
|
| 222 | + // It is possible to look at the points making up each face in more detail. |
| 223 | + // First, we define a few shorthands by getting references to the respective |
| 224 | + // objects in the mesh. These are defined as constants since we do not aim to |
| 225 | + // alter the mesh in any way. |
| 226 | + // NOTE: these lists refer to the physical definition of the mesh and thus |
| 227 | + // include boundary faces. Use can be made of the mesh.boundary()[patchI].Cf().size() |
| 228 | + // and mesh.boundary()[patchI].start() methods to check whether the face is internal |
| 229 | + // or lies on a boundary. |
| 230 | + const faceList& fcs = mesh.faces(); |
| 231 | + const List<point>& pts = mesh.points(); |
| 232 | + const List<point>& cents = mesh.faceCentres(); |
| 233 | +
|
| 234 | + forAll(fcs,faceI) |
| 235 | + if (faceI%80==0) |
| 236 | + { |
| 237 | + if (faceI<mesh.Cf().size()) |
| 238 | + Info << "Internal face "; |
| 239 | + else |
| 240 | + { |
| 241 | + forAll(mesh.boundary(),patchI) |
| 242 | + if ((mesh.boundary()[patchI].start()<= faceI) && |
| 243 | + (faceI < mesh.boundary()[patchI].start()+mesh.boundary()[patchI].Cf().size())) |
| 244 | + { |
| 245 | + Info << "Face on patch " << patchI << ", faceI "; |
| 246 | + break; // exit the forAll loop prematurely |
| 247 | + } |
| 248 | + } |
| 249 | +
|
| 250 | + Info << faceI << " with centre at " << cents[faceI] |
| 251 | + << " has " << fcs[faceI].size() << " vertices:"; |
| 252 | + forAll(fcs[faceI],vertexI) |
| 253 | + // Note how fcs[faceI] holds the indices of points whose coordinates |
| 254 | + // are stored in the pts list. |
| 255 | + Info << " " << pts[fcs[faceI][vertexI]]; |
| 256 | + Info << endl; |
| 257 | + } |
| 258 | + Info << endl; |
| 259 | +
|
| 260 | + // In the original cavity tutorial, on which the test case is based, |
| 261 | + // the frontAndBack boundary is defined as and "empty" type. This is a special |
| 262 | + // BC case which may cause unexpected behaviour as its .Cf() field has size of 0. |
| 263 | + // Type of a patch may be checked to avoid running into this problem if there |
| 264 | + // is a substantial risk that an empty patch type will appear |
| 265 | + label patchID(0); |
| 266 | + const polyPatch& pp = mesh.boundaryMesh()[patchID]; |
| 267 | + if (isA<emptyPolyPatch>(pp)) |
| 268 | + { |
| 269 | + // patch patchID is of type "empty". |
| 270 | + Info << "You will not see this." << endl; |
| 271 | + } |
| 272 | +
|
| 273 | + // Patches may also be retrieved from the mesh using their name. This could be |
| 274 | + // useful if the user were to refer to a particular patch from a dictionary |
| 275 | + // (like when you do when calculating forces on a particular patch). |
| 276 | + word patchName("movingWall"); |
| 277 | + patchID = mesh.boundaryMesh().findPatchID(patchName); |
| 278 | + Info << "Retrieved patch " << patchName << " at index " << patchID << " using its name only." << nl << endl; |
| 279 | +*/ |
| 280 | + Info<< "End\n" << endl; |
| 281 | + |
| 282 | + return 0; |
| 283 | +} |
| 284 | + |
| 285 | + |
| 286 | +// ************************************************************************* // |
0 commit comments