Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 #include <OpenFOAM/argList.H>
00071 #include <OpenFOAM/Time.H>
00072 #include <OpenFOAM/dimensionedTypes.H>
00073 #include <OpenFOAM/IFstream.H>
00074 #include "faceMesh.H"
00075 #include <dynamicMesh/polyTopoChange.H>
00076 #include <dynamicMesh/polyTopoChanger.H>
00077 #include <dynamicMesh/edgeCollapser.H>
00078 #include <OpenFOAM/mathematicalConstants.H>
00079 #include <OpenFOAM/globalMeshData.H>
00080 #include <dynamicMesh/perfectInterface.H>
00081
00082 #include "extrudedMesh/extrudedMesh.H"
00083 #include <extrudeModel/extrudeModel.H>
00084
00085 using namespace Foam;
00086
00087
00088
00089
00090 int main(int argc, char *argv[])
00091 {
00092 #include <OpenFOAM/setRootCase.H>
00093 #include "createTimeExtruded.H"
00094
00095 autoPtr<extrudedMesh> meshPtr(NULL);
00096
00097 IOdictionary dict
00098 (
00099 IOobject
00100 (
00101 "extrudeProperties",
00102 runTimeExtruded.constant(),
00103 runTimeExtruded,
00104 IOobject::MUST_READ
00105 )
00106 );
00107
00108 autoPtr<extrudeModel> model(extrudeModel::New(dict));
00109
00110 const word sourceType(dict.lookup("constructFrom"));
00111
00112 autoPtr<faceMesh> fMesh;
00113
00114 if (sourceType == "patch")
00115 {
00116 fileName sourceCasePath(dict.lookup("sourceCase"));
00117 sourceCasePath.expand();
00118 fileName sourceRootDir = sourceCasePath.path();
00119 fileName sourceCaseDir = sourceCasePath.name();
00120 word patchName(dict.lookup("sourcePatch"));
00121
00122 Info<< "Extruding patch " << patchName
00123 << " on mesh " << sourceCasePath << nl
00124 << endl;
00125
00126 Time runTime
00127 (
00128 Time::controlDictName,
00129 sourceRootDir,
00130 sourceCaseDir
00131 );
00132 #include <OpenFOAM/createPolyMesh.H>
00133
00134 label patchID = mesh.boundaryMesh().findPatchID(patchName);
00135
00136 if (patchID == -1)
00137 {
00138 FatalErrorIn(args.executable())
00139 << "Cannot find patch " << patchName
00140 << " in the source mesh.\n"
00141 << "Valid patch names are " << mesh.boundaryMesh().names()
00142 << exit(FatalError);
00143 }
00144
00145 const polyPatch& pp = mesh.boundaryMesh()[patchID];
00146 fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints()));
00147
00148 {
00149 fileName surfName(runTime.path()/patchName + ".sMesh");
00150 Info<< "Writing patch as surfaceMesh to "
00151 << surfName << nl << endl;
00152 OFstream os(surfName);
00153 os << fMesh() << nl;
00154 }
00155 }
00156 else if (sourceType == "surface")
00157 {
00158
00159 fileName surfName(dict.lookup("surface"));
00160
00161 Info<< "Extruding surfaceMesh read from file " << surfName << nl
00162 << endl;
00163
00164 IFstream is(surfName);
00165
00166 fMesh.reset(new faceMesh(is));
00167
00168 Info<< "Read patch from file " << surfName << nl
00169 << endl;
00170 }
00171 else
00172 {
00173 FatalErrorIn(args.executable())
00174 << "Illegal 'constructFrom' specification. Should either be "
00175 << "patch or surface." << exit(FatalError);
00176 }
00177
00178 Switch flipNormals(dict.lookup("flipNormals"));
00179
00180 if (flipNormals)
00181 {
00182 Info<< "Flipping faces." << nl << endl;
00183
00184 faceList faces(fMesh().size());
00185 forAll(faces, i)
00186 {
00187 faces[i] = fMesh()[i].reverseFace();
00188 }
00189 fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
00190 }
00191
00192
00193 Info<< "Extruding patch with :" << nl
00194 << " points : " << fMesh().points().size() << nl
00195 << " faces : " << fMesh().size() << nl
00196 << " normals[0] : " << fMesh().faceNormals()[0]
00197 << nl
00198 << endl;
00199
00200 extrudedMesh mesh
00201 (
00202 IOobject
00203 (
00204 extrudedMesh::defaultRegion,
00205 runTimeExtruded.constant(),
00206 runTimeExtruded
00207 ),
00208 fMesh(),
00209 model()
00210 );
00211
00212
00213 const boundBox& bb = mesh.globalData().bb();
00214 const vector span = bb.span();
00215 const scalar mergeDim = 1E-4 * bb.minDim();
00216
00217 Info<< "Mesh bounding box : " << bb << nl
00218 << " with span : " << span << nl
00219 << "Merge distance : " << mergeDim << nl
00220 << endl;
00221
00222 const polyBoundaryMesh& patches = mesh.boundaryMesh();
00223
00224 const label origPatchID = patches.findPatchID("originalPatch");
00225 const label otherPatchID = patches.findPatchID("otherSide");
00226
00227 if (origPatchID == -1 || otherPatchID == -1)
00228 {
00229 FatalErrorIn(args.executable())
00230 << "Cannot find patch originalPatch or otherSide." << nl
00231 << "Valid patches are " << patches.names() << exit(FatalError);
00232 }
00233
00234
00235
00236
00237 {
00238 Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
00239
00240
00241 edgeCollapser collapser(mesh);
00242
00243 const edgeList& edges = mesh.edges();
00244 const pointField& points = mesh.points();
00245
00246 forAll(edges, edgeI)
00247 {
00248 const edge& e = edges[edgeI];
00249
00250 scalar d = e.mag(points);
00251
00252 if (d < mergeDim)
00253 {
00254 Info<< "Merging edge " << e << " since length " << d
00255 << " << " << mergeDim << nl;
00256
00257
00258 collapser.collapseEdge(edgeI, e[0]);
00259 }
00260 }
00261
00262
00263 polyTopoChange meshMod(mesh);
00264
00265 bool anyChange = collapser.setRefinement(meshMod);
00266
00267 if (anyChange)
00268 {
00269
00270 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
00271
00272
00273 mesh.updateMesh(map);
00274
00275
00276 if (map().hasMotionPoints())
00277 {
00278 mesh.movePoints(map().preMotionPoints());
00279 }
00280 }
00281 }
00282
00283
00284
00285
00286
00287 Switch mergeFaces(dict.lookup("mergeFaces"));
00288 if (mergeFaces)
00289 {
00290 Info<< "Assuming full 360 degree axisymmetric case;"
00291 << " stitching faces on patches "
00292 << patches[origPatchID].name() << " and "
00293 << patches[otherPatchID].name() << " together ..." << nl << endl;
00294
00295 polyTopoChanger stitcher(mesh);
00296 stitcher.setSize(1);
00297
00298
00299 labelList isf(patches[origPatchID].size());
00300
00301 forAll (isf, i)
00302 {
00303 isf[i] = patches[origPatchID].start() + i;
00304 }
00305
00306 const word cutZoneName("originalCutFaceZone");
00307
00308 List<faceZone*> fz
00309 (
00310 1,
00311 new faceZone
00312 (
00313 cutZoneName,
00314 isf,
00315 boolList(isf.size(), false),
00316 0,
00317 mesh.faceZones()
00318 )
00319 );
00320
00321 mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
00322
00323
00324 stitcher.set
00325 (
00326 0,
00327 new perfectInterface
00328 (
00329 "couple",
00330 0,
00331 stitcher,
00332 cutZoneName,
00333 patches[origPatchID].name(),
00334 patches[otherPatchID].name()
00335 )
00336 );
00337
00338
00339 autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
00340
00341 mesh.movePoints(morphMap->preMotionPoints());
00342 }
00343
00344 if (!mesh.write())
00345 {
00346 FatalErrorIn(args.executable()) << "Failed writing mesh"
00347 << exit(FatalError);
00348 }
00349
00350 Info << "End\n" << endl;
00351
00352 return 0;
00353 }
00354
00355
00356