PageRenderTime 53ms CodeModel.GetById 26ms RepoModel.GetById 0ms app.codeStats 0ms

/applications/test/extendedStencil/testExtendedStencil.C

https://github.com/xyuan/OpenFOAM-1.7.x
C | 498 lines | 96 code | 54 blank | 348 comment | 1 complexity | 7d86af17f47de2670b2cd6ad156ca712 MD5 | raw file
  1. /*---------------------------------------------------------------------------*\
  2. ========= |
  3. \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
  4. \\ / O peration |
  5. \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
  6. \\/ M anipulation |
  7. -------------------------------------------------------------------------------
  8. License
  9. This file is part of OpenFOAM.
  10. OpenFOAM is free software: you can redistribute it and/or modify it
  11. under the terms of the GNU General Public License as published by
  12. the Free Software Foundation, either version 3 of the License, or
  13. (at your option) any later version.
  14. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  15. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16. FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  17. for more details.
  18. You should have received a copy of the GNU General Public License
  19. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
  20. Application
  21. testExtendedStencil
  22. Description
  23. Test app for determining extended stencil.
  24. \*---------------------------------------------------------------------------*/
  25. #include "argList.H"
  26. #include "fvMesh.H"
  27. #include "volFields.H"
  28. #include "Time.H"
  29. #include "mapDistribute.H"
  30. #include "OFstream.H"
  31. #include "meshTools.H"
  32. //#include "FECCellToFaceStencil.H"
  33. //#include "CFCCellToFaceStencil.H"
  34. //#include "CPCCellToFaceStencil.H"
  35. //#include "CECCellToFaceStencil.H"
  36. //#include "extendedCentredCellToFaceStencil.H"
  37. //#include "extendedUpwindCellToFaceStencil.H"
  38. //#include "centredCFCCellToFaceStencilObject.H"
  39. //#include "centredFECCellToFaceStencilObject.H"
  40. //#include "centredCPCCellToFaceStencilObject.H"
  41. //#include "centredCECCellToFaceStencilObject.H"
  42. //#include "upwindFECCellToFaceStencilObject.H"
  43. //#include "upwindCPCCellToFaceStencilObject.H"
  44. //#include "upwindCECCellToFaceStencilObject.H"
  45. //#include "upwindCFCCellToFaceStencilObject.H"
  46. #include "centredCFCFaceToCellStencilObject.H"
  47. using namespace Foam;
  48. // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49. void writeStencilOBJ
  50. (
  51. const fileName& fName,
  52. const point& fc,
  53. const List<point>& stencilCc
  54. )
  55. {
  56. OFstream str(fName);
  57. label vertI = 0;
  58. meshTools::writeOBJ(str, fc);
  59. vertI++;
  60. forAll(stencilCc, i)
  61. {
  62. meshTools::writeOBJ(str, stencilCc[i]);
  63. vertI++;
  64. str << "l 1 " << vertI << nl;
  65. }
  66. }
  67. // Stats
  68. void writeStencilStats(const labelListList& stencil)
  69. {
  70. label sumSize = 0;
  71. label nSum = 0;
  72. label minSize = labelMax;
  73. label maxSize = labelMin;
  74. forAll(stencil, i)
  75. {
  76. const labelList& sCells = stencil[i];
  77. if (sCells.size() > 0)
  78. {
  79. sumSize += sCells.size();
  80. nSum++;
  81. minSize = min(minSize, sCells.size());
  82. maxSize = max(maxSize, sCells.size());
  83. }
  84. }
  85. reduce(sumSize, sumOp<label>());
  86. reduce(nSum, sumOp<label>());
  87. sumSize /= nSum;
  88. reduce(minSize, minOp<label>());
  89. reduce(maxSize, maxOp<label>());
  90. Info<< "Stencil size :" << nl
  91. << " average : " << sumSize << nl
  92. << " min : " << minSize << nl
  93. << " max : " << maxSize << nl
  94. << endl;
  95. }
  96. // Main program:
  97. int main(int argc, char *argv[])
  98. {
  99. # include "addTimeOptions.H"
  100. # include "setRootCase.H"
  101. # include "createTime.H"
  102. // Get times list
  103. instantList Times = runTime.times();
  104. # include "checkTimeOptions.H"
  105. runTime.setTime(Times[startTime], startTime);
  106. # include "createMesh.H"
  107. // Force calculation of extended edge addressing
  108. const labelListList& edgeFaces = mesh.edgeFaces();
  109. const labelListList& edgeCells = mesh.edgeCells();
  110. const labelListList& pointCells = mesh.pointCells();
  111. Info<< "dummy:" << edgeFaces.size() + edgeCells.size() + pointCells.size()
  112. << endl;
  113. // Centred, semi-extended stencil (edge cells only)
  114. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  115. // {
  116. // //const FECCellToFaceStencil cfcStencil(mesh);
  117. // //const extendedCentredCellToFaceStencil addressing
  118. // //(
  119. // // cfcStencil
  120. // //);
  121. // const extendedCentredStencil& addressing =
  122. // centredFECCellToFaceStencilObject::New
  123. // (
  124. // mesh
  125. // );
  126. //
  127. // Info<< "faceEdgeCell:" << endl;
  128. // writeStencilStats(addressing.stencil());
  129. //
  130. // // Collect stencil cell centres
  131. // List<List<point> > stencilPoints(mesh.nFaces());
  132. // addressing.collectData
  133. // (
  134. // mesh.C(),
  135. // stencilPoints
  136. // );
  137. //
  138. // forAll(stencilPoints, faceI)
  139. // {
  140. // writeStencilOBJ
  141. // (
  142. // runTime.path()/"faceEdgeCell" + Foam::name(faceI) + ".obj",
  143. // mesh.faceCentres()[faceI],
  144. // stencilPoints[faceI]
  145. // );
  146. // }
  147. // }
  148. // // Centred, face stencil
  149. // // ~~~~~~~~~~~~~~~~~~~~~
  150. //
  151. // {
  152. // const extendedCentredCellToFaceStencil& addressing =
  153. // centredCFCCellToFaceStencilObject::New
  154. // (
  155. // mesh
  156. // );
  157. //
  158. // Info<< "cellFaceCell:" << endl;
  159. // writeStencilStats(addressing.stencil());
  160. //
  161. //
  162. // //// Do some interpolation.
  163. // //{
  164. // // const labelListList& stencil = addressing.stencil();
  165. // // List<List<scalar> > stencilWeights(stencil.size());
  166. // // forAll(stencil, faceI)
  167. // // {
  168. // // const labelList& fStencil = stencil[faceI];
  169. // //
  170. // // if (fStencil.size() > 0)
  171. // // {
  172. // // // Uniform weights
  173. // // stencilWeights[faceI] = scalarList
  174. // // (
  175. // // fStencil.size(),
  176. // // 1.0/fStencil.size()
  177. // // );
  178. // // }
  179. // // }
  180. // //
  181. // // tmp<surfaceVectorField> tfc
  182. // // (
  183. // // addressing.weightedSum(mesh.C(), stencilWeights)
  184. // // );
  185. // //}
  186. //
  187. //
  188. // // Collect stencil cell centres
  189. // List<List<point> > stencilPoints(mesh.nFaces());
  190. // addressing.collectData
  191. // (
  192. // mesh.C(),
  193. // stencilPoints
  194. // );
  195. //
  196. // forAll(stencilPoints, faceI)
  197. // {
  198. // if (stencilPoints[faceI].size() >= 15)
  199. // {
  200. // writeStencilOBJ
  201. // (
  202. // runTime.path()/"centredFace" + Foam::name(faceI) + ".obj",
  203. // mesh.faceCentres()[faceI],
  204. // stencilPoints[faceI]
  205. // );
  206. // }
  207. // }
  208. // }
  209. // // Centred, point stencil
  210. // // ~~~~~~~~~~~~~~~~~~~~~~
  211. //
  212. // {
  213. // //const extendedCentredCellToFaceStencil& addressing =
  214. // //centredCPCStencilObject::New
  215. // //(
  216. // // mesh
  217. // //);
  218. // //
  219. // //Info<< "cellPointCell:" << endl;
  220. // //writeStencilStats(addressing.stencil());
  221. // //
  222. // //
  223. // //// Collect stencil cell centres
  224. // //List<List<point> > stencilPoints(mesh.nFaces());
  225. // //addressing.collectData
  226. // //(
  227. // // mesh.C(),
  228. // // stencilPoints
  229. // //);
  230. // //
  231. // //forAll(stencilPoints, faceI)
  232. // //{
  233. // // writeStencilOBJ
  234. // // (
  235. // // runTime.path()/"centredPoint" + Foam::name(faceI) + ".obj",
  236. // // mesh.faceCentres()[faceI],
  237. // // stencilPoints[faceI]
  238. // // );
  239. // //}
  240. // }
  241. // // Centred, edge stencil
  242. // // ~~~~~~~~~~~~~~~~~~~~~~
  243. //
  244. // {
  245. // //const extendedCentredCellToFaceStencil& addressing =
  246. // //centredCECStencilObject::New
  247. // //(
  248. // // mesh
  249. // //);
  250. // //
  251. // //Info<< "cellEdgeCell:" << endl;
  252. // //writeStencilStats(addressing.stencil());
  253. // //
  254. // //
  255. // //// Collect stencil cell centres
  256. // //List<List<point> > stencilPoints(mesh.nFaces());
  257. // //addressing.collectData
  258. // //(
  259. // // mesh.C(),
  260. // // stencilPoints
  261. // //);
  262. // //
  263. // //forAll(stencilPoints, faceI)
  264. // //{
  265. // // writeStencilOBJ
  266. // // (
  267. // // runTime.path()/"centredEdge" + Foam::name(faceI) + ".obj",
  268. // // mesh.faceCentres()[faceI],
  269. // // stencilPoints[faceI]
  270. // // );
  271. // //}
  272. // }
  273. // Upwind, semi-extended stencil
  274. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  275. //{
  276. // const extendedUpwindCellToFaceStencil& addressing =
  277. // upwindFECCellToFaceStencilObject::New
  278. // (
  279. // mesh,
  280. // 0.5
  281. // );
  282. //
  283. // Info<< "upwind-faceEdgeCell:" << endl;
  284. // writeStencilStats(addressing.ownStencil());
  285. //
  286. // {
  287. // // Collect stencil cell centres
  288. // List<List<point> > ownPoints(mesh.nFaces());
  289. // addressing.collectData
  290. // (
  291. // addressing.ownMap(),
  292. // addressing.ownStencil(),
  293. // mesh.C(),
  294. // ownPoints
  295. // );
  296. //
  297. // forAll(ownPoints, faceI)
  298. // {
  299. // writeStencilOBJ
  300. // (
  301. // runTime.path()/"ownFEC" + Foam::name(faceI) + ".obj",
  302. // mesh.faceCentres()[faceI],
  303. // ownPoints[faceI]
  304. // );
  305. // }
  306. // }
  307. // {
  308. // // Collect stencil cell centres
  309. // List<List<point> > neiPoints(mesh.nFaces());
  310. // addressing.collectData
  311. // (
  312. // addressing.neiMap(),
  313. // addressing.neiStencil(),
  314. // mesh.C(),
  315. // neiPoints
  316. // );
  317. //
  318. // forAll(neiPoints, faceI)
  319. // {
  320. // writeStencilOBJ
  321. // (
  322. // runTime.path()/"neiFEC" + Foam::name(faceI) + ".obj",
  323. // mesh.faceCentres()[faceI],
  324. // neiPoints[faceI]
  325. // );
  326. // }
  327. // }
  328. //}
  329. // Upwind, extended stencil
  330. // ~~~~~~~~~~~~~~~~~~~~~~~~
  331. //{
  332. // const extendedUpwindCellToFaceStencil& addressing =
  333. // upwindCFCCellToFaceStencilObject::New
  334. // (
  335. // mesh,
  336. // 0.5
  337. // );
  338. //
  339. // Info<< "upwind-cellFaceCell:" << endl;
  340. // writeStencilStats(addressing.ownStencil());
  341. //
  342. // {
  343. // // Collect stencil cell centres
  344. // List<List<point> > ownPoints(mesh.nFaces());
  345. // addressing.collectData
  346. // (
  347. // addressing.ownMap(),
  348. // addressing.ownStencil(),
  349. // mesh.C(),
  350. // ownPoints
  351. // );
  352. //
  353. // forAll(ownPoints, faceI)
  354. // {
  355. // writeStencilOBJ
  356. // (
  357. // runTime.path()/"ownCFC" + Foam::name(faceI) + ".obj",
  358. // mesh.faceCentres()[faceI],
  359. // ownPoints[faceI]
  360. // );
  361. // }
  362. // }
  363. // {
  364. // // Collect stencil cell centres
  365. // List<List<point> > neiPoints(mesh.nFaces());
  366. // addressing.collectData
  367. // (
  368. // addressing.neiMap(),
  369. // addressing.neiStencil(),
  370. // mesh.C(),
  371. // neiPoints
  372. // );
  373. //
  374. // forAll(neiPoints, faceI)
  375. // {
  376. // writeStencilOBJ
  377. // (
  378. // runTime.path()/"neiCFC" + Foam::name(faceI) + ".obj",
  379. // mesh.faceCentres()[faceI],
  380. // neiPoints[faceI]
  381. // );
  382. // }
  383. // }
  384. //}
  385. //---- CELL CENTRED STENCIL -----
  386. // Centred, cell stencil
  387. // ~~~~~~~~~~~~~~~~~~~~~
  388. {
  389. const extendedCentredFaceToCellStencil& addressing =
  390. centredCFCFaceToCellStencilObject::New
  391. (
  392. mesh
  393. );
  394. Info<< "cellFaceCell:" << endl;
  395. writeStencilStats(addressing.stencil());
  396. // Collect stencil face centres
  397. List<List<point> > stencilPoints(mesh.nCells());
  398. addressing.collectData
  399. (
  400. mesh.Cf(),
  401. stencilPoints
  402. );
  403. forAll(stencilPoints, cellI)
  404. {
  405. writeStencilOBJ
  406. (
  407. runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
  408. mesh.cellCentres()[cellI],
  409. stencilPoints[cellI]
  410. );
  411. }
  412. }
  413. //XXXXXX
  414. // // Evaluate
  415. // List<List<scalar> > stencilData(faceStencils.size());
  416. // collectStencilData
  417. // (
  418. // distMap,
  419. // faceStencils,
  420. // vf,
  421. // stencilData
  422. // );
  423. // for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
  424. // {
  425. // const scalarList& stData = stencilData[faceI];
  426. // const scalarList& stWeight = fit[faceI];
  427. //
  428. // forAll(stData, i)
  429. // {
  430. // sf[faceI] += stWeight[i]*stData[i];
  431. // }
  432. // }
  433. // See finiteVolume/lnInclude/leastSquaresGrad.C
  434. //XXXXXX
  435. Pout<< "End\n" << endl;
  436. return 0;
  437. }
  438. // ************************************************************************* //