Test-ExtendedStencil.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) 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 Application
25  testExtendedStencil
26 
27 Description
28  Test app for determining extended stencil.
29 
30 \*---------------------------------------------------------------------------*/
31 
32 #include "argList.H"
33 #include "fvMesh.H"
34 #include "volFields.H"
35 #include "Time.H"
36 //#include "mapDistribute.H"
37 #include "OFstream.H"
38 #include "meshTools.H"
39 //#include "FECCellToFaceStencil.H"
40 //#include "CFCCellToFaceStencil.H"
41 //#include "CPCCellToFaceStencil.H"
42 //#include "CECCellToFaceStencil.H"
43 //#include "extendedCentredCellToFaceStencil.H"
44 //#include "extendedUpwindCellToFaceStencil.H"
45 
46 //#include "centredCFCCellToFaceStencilObject.H"
47 //#include "centredFECCellToFaceStencilObject.H"
48 //#include "centredCPCCellToFaceStencilObject.H"
49 //#include "centredCECCellToFaceStencilObject.H"
50 
51 //#include "upwindFECCellToFaceStencilObject.H"
52 //#include "upwindCPCCellToFaceStencilObject.H"
53 //#include "upwindCECCellToFaceStencilObject.H"
54 
55 //#include "upwindCFCCellToFaceStencilObject.H"
56 //#include "centredCFCFaceToCellStencilObject.H"
57 
61 
62 using namespace Foam;
63 
64 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
65 
66 void writeStencilOBJ
67 (
68  const fileName& fName,
69  const point& fc,
70  const List<point>& stencilCc
71 )
72 {
73  OFstream str(fName);
74  label vertI = 0;
75 
76  meshTools::writeOBJ(str, fc);
77  vertI++;
78 
79  forAll(stencilCc, i)
80  {
81  meshTools::writeOBJ(str, stencilCc[i]);
82  vertI++;
83  str << "l 1 " << vertI << nl;
84  }
85 }
86 
87 
88 // Stats
89 void writeStencilStats(const labelListList& stencil)
90 {
91  label sumSize = 0;
92  label nSum = 0;
93  label minSize = labelMax;
94  label maxSize = labelMin;
95 
96  forAll(stencil, i)
97  {
98  const labelList& sCells = stencil[i];
99 
100  if (sCells.size() > 0)
101  {
102  sumSize += sCells.size();
103  nSum++;
104  minSize = min(minSize, sCells.size());
105  maxSize = max(maxSize, sCells.size());
106  }
107  }
108  reduce(sumSize, sumOp<label>());
109  reduce(nSum, sumOp<label>());
110  sumSize /= nSum;
111 
112  reduce(minSize, minOp<label>());
113  reduce(maxSize, maxOp<label>());
114 
115  Info<< "Stencil size :" << nl
116  << " average : " << sumSize << nl
117  << " min : " << minSize << nl
118  << " max : " << maxSize << nl
119  << endl;
120 }
121 
122 
123 // Main program:
124 
125 int main(int argc, char *argv[])
126 {
127  #include "addTimeOptions.H"
128  #include "setRootCase.H"
129  #include "createTime.H"
130 
131  // Get times list
132  instantList Times = runTime.times();
133  #include "checkTimeOptions.H"
134  runTime.setTime(Times[startTime], startTime);
135  #include "createMesh.H"
136 
137  // Force calculation of extended edge addressing
138  const labelListList& edgeFaces = mesh.edgeFaces();
139  const labelListList& edgeCells = mesh.edgeCells();
140  const labelListList& pointCells = mesh.pointCells();
141  Info<< "dummy:" << edgeFaces.size() + edgeCells.size() + pointCells.size()
142  << endl;
143 
144 
145  // Centred, semi-extended stencil (edge cells only)
146  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
147 
148 // {
149 // //const FECCellToFaceStencil cfcStencil(mesh);
150 // //const extendedCentredCellToFaceStencil addressing
151 // //(
152 // // cfcStencil
153 // //);
154 // const extendedCentredStencil& addressing =
155 // centredFECCellToFaceStencilObject::New
156 // (
157 // mesh
158 // );
159 //
160 // Info<< "faceEdgeCell:" << endl;
161 // writeStencilStats(addressing.stencil());
162 //
163 // // Collect stencil cell centres
164 // List<List<point> > stencilPoints(mesh.nFaces());
165 // addressing.collectData
166 // (
167 // mesh.C(),
168 // stencilPoints
169 // );
170 //
171 // forAll(stencilPoints, faceI)
172 // {
173 // writeStencilOBJ
174 // (
175 // runTime.path()/"faceEdgeCell" + Foam::name(faceI) + ".obj",
176 // mesh.faceCentres()[faceI],
177 // stencilPoints[faceI]
178 // );
179 // }
180 // }
181 
182 
183 
184 
185 // // Centred, face stencil
186 // // ~~~~~~~~~~~~~~~~~~~~~
187 //
188 // {
189 // const extendedCentredCellToFaceStencil& addressing =
190 // centredCFCCellToFaceStencilObject::New
191 // (
192 // mesh
193 // );
194 //
195 // Info<< "cellFaceCell:" << endl;
196 // writeStencilStats(addressing.stencil());
197 //
198 //
199 // //// Do some interpolation.
200 // //{
201 // // const labelListList& stencil = addressing.stencil();
202 // // List<List<scalar> > stencilWeights(stencil.size());
203 // // forAll(stencil, faceI)
204 // // {
205 // // const labelList& fStencil = stencil[faceI];
206 // //
207 // // if (fStencil.size() > 0)
208 // // {
209 // // // Uniform weights
210 // // stencilWeights[faceI] = scalarList
211 // // (
212 // // fStencil.size(),
213 // // 1.0/fStencil.size()
214 // // );
215 // // }
216 // // }
217 // //
218 // // tmp<surfaceVectorField> tfc
219 // // (
220 // // addressing.weightedSum(mesh.C(), stencilWeights)
221 // // );
222 // //}
223 //
224 //
225 // // Collect stencil cell centres
226 // List<List<point> > stencilPoints(mesh.nFaces());
227 // addressing.collectData
228 // (
229 // mesh.C(),
230 // stencilPoints
231 // );
232 //
233 // forAll(stencilPoints, faceI)
234 // {
235 // if (stencilPoints[faceI].size() >= 15)
236 // {
237 // writeStencilOBJ
238 // (
239 // runTime.path()/"centredFace" + Foam::name(faceI) + ".obj",
240 // mesh.faceCentres()[faceI],
241 // stencilPoints[faceI]
242 // );
243 // }
244 // }
245 // }
246 
247 
248 // // Centred, point stencil
249 // // ~~~~~~~~~~~~~~~~~~~~~~
250 //
251 // {
252 // //const extendedCentredCellToFaceStencil& addressing =
253 // //centredCPCStencilObject::New
254 // //(
255 // // mesh
256 // //);
257 // //
258 // //Info<< "cellPointCell:" << endl;
259 // //writeStencilStats(addressing.stencil());
260 // //
261 // //
262 // //// Collect stencil cell centres
263 // //List<List<point> > stencilPoints(mesh.nFaces());
264 // //addressing.collectData
265 // //(
266 // // mesh.C(),
267 // // stencilPoints
268 // //);
269 // //
270 // //forAll(stencilPoints, faceI)
271 // //{
272 // // writeStencilOBJ
273 // // (
274 // // runTime.path()/"centredPoint" + Foam::name(faceI) + ".obj",
275 // // mesh.faceCentres()[faceI],
276 // // stencilPoints[faceI]
277 // // );
278 // //}
279 // }
280 
281 
282 
283 // // Centred, edge stencil
284 // // ~~~~~~~~~~~~~~~~~~~~~~
285 //
286 // {
287 // //const extendedCentredCellToFaceStencil& addressing =
288 // //centredCECStencilObject::New
289 // //(
290 // // mesh
291 // //);
292 // //
293 // //Info<< "cellEdgeCell:" << endl;
294 // //writeStencilStats(addressing.stencil());
295 // //
296 // //
297 // //// Collect stencil cell centres
298 // //List<List<point> > stencilPoints(mesh.nFaces());
299 // //addressing.collectData
300 // //(
301 // // mesh.C(),
302 // // stencilPoints
303 // //);
304 // //
305 // //forAll(stencilPoints, faceI)
306 // //{
307 // // writeStencilOBJ
308 // // (
309 // // runTime.path()/"centredEdge" + Foam::name(faceI) + ".obj",
310 // // mesh.faceCentres()[faceI],
311 // // stencilPoints[faceI]
312 // // );
313 // //}
314 // }
315 
316 
317 
318  // Upwind, semi-extended stencil
319  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
320 
321  //{
322  // const extendedUpwindCellToFaceStencil& addressing =
323  // upwindFECCellToFaceStencilObject::New
324  // (
325  // mesh,
326  // 0.5
327  // );
328  //
329  // Info<< "upwind-faceEdgeCell:" << endl;
330  // writeStencilStats(addressing.ownStencil());
331  //
332  // {
333  // // Collect stencil cell centres
334  // List<List<point> > ownPoints(mesh.nFaces());
335  // addressing.collectData
336  // (
337  // addressing.ownMap(),
338  // addressing.ownStencil(),
339  // mesh.C(),
340  // ownPoints
341  // );
342  //
343  // forAll(ownPoints, faceI)
344  // {
345  // writeStencilOBJ
346  // (
347  // runTime.path()/"ownFEC" + Foam::name(faceI) + ".obj",
348  // mesh.faceCentres()[faceI],
349  // ownPoints[faceI]
350  // );
351  // }
352  // }
353  // {
354  // // Collect stencil cell centres
355  // List<List<point> > neiPoints(mesh.nFaces());
356  // addressing.collectData
357  // (
358  // addressing.neiMap(),
359  // addressing.neiStencil(),
360  // mesh.C(),
361  // neiPoints
362  // );
363  //
364  // forAll(neiPoints, faceI)
365  // {
366  // writeStencilOBJ
367  // (
368  // runTime.path()/"neiFEC" + Foam::name(faceI) + ".obj",
369  // mesh.faceCentres()[faceI],
370  // neiPoints[faceI]
371  // );
372  // }
373  // }
374  //}
375 
376 
377 
378  // Upwind, extended stencil
379  // ~~~~~~~~~~~~~~~~~~~~~~~~
380 
381  //{
382  // const extendedUpwindCellToFaceStencil& addressing =
383  // upwindCFCCellToFaceStencilObject::New
384  // (
385  // mesh,
386  // 0.5
387  // );
388  //
389  // Info<< "upwind-cellFaceCell:" << endl;
390  // writeStencilStats(addressing.ownStencil());
391  //
392  // {
393  // // Collect stencil cell centres
394  // List<List<point> > ownPoints(mesh.nFaces());
395  // addressing.collectData
396  // (
397  // addressing.ownMap(),
398  // addressing.ownStencil(),
399  // mesh.C(),
400  // ownPoints
401  // );
402  //
403  // forAll(ownPoints, faceI)
404  // {
405  // writeStencilOBJ
406  // (
407  // runTime.path()/"ownCFC" + Foam::name(faceI) + ".obj",
408  // mesh.faceCentres()[faceI],
409  // ownPoints[faceI]
410  // );
411  // }
412  // }
413  // {
414  // // Collect stencil cell centres
415  // List<List<point> > neiPoints(mesh.nFaces());
416  // addressing.collectData
417  // (
418  // addressing.neiMap(),
419  // addressing.neiStencil(),
420  // mesh.C(),
421  // neiPoints
422  // );
423  //
424  // forAll(neiPoints, faceI)
425  // {
426  // writeStencilOBJ
427  // (
428  // runTime.path()/"neiCFC" + Foam::name(faceI) + ".obj",
429  // mesh.faceCentres()[faceI],
430  // neiPoints[faceI]
431  // );
432  // }
433  // }
434  //}
435 
436 
437 
438  //---- CELL CENTRED STENCIL -----
439 
440  // Centred, cell stencil
441  // ~~~~~~~~~~~~~~~~~~~~~
442 
443  {
444  const extendedCentredCellToCellStencil& addressing =
446  (
447  mesh
448  );
449 
450  Info<< "cellCellCell:" << endl;
451  writeStencilStats(addressing.stencil());
452 
453  // Collect stencil cell centres
454  List<List<point> > stencilPoints(mesh.nCells());
455  addressing.collectData
456  (
457  mesh.C(),
458  stencilPoints
459  );
460 
461  forAll(stencilPoints, cellI)
462  {
464  (
465  runTime.path()/"centredCECCell" + Foam::name(cellI) + ".obj",
466  mesh.cellCentres()[cellI],
467  stencilPoints[cellI]
468  );
469  }
470  }
471  {
472  const extendedCentredCellToCellStencil& addressing =
474  (
475  mesh
476  );
477 
478  Info<< "cellCellCell:" << endl;
479  writeStencilStats(addressing.stencil());
480 
481  // Collect stencil cell centres
482  List<List<point> > stencilPoints(mesh.nCells());
483  addressing.collectData
484  (
485  mesh.C(),
486  stencilPoints
487  );
488 
489  forAll(stencilPoints, cellI)
490  {
492  (
493  runTime.path()/"centredCPCCell" + Foam::name(cellI) + ".obj",
494  mesh.cellCentres()[cellI],
495  stencilPoints[cellI]
496  );
497  }
498  }
499  {
500  const extendedCentredCellToCellStencil& addressing =
502  (
503  mesh
504  );
505 
506  Info<< "cellCellCell:" << endl;
507  writeStencilStats(addressing.stencil());
508 
509  // Collect stencil cell centres
510  List<List<point> > stencilPoints(mesh.nCells());
511  addressing.collectData
512  (
513  mesh.C(),
514  stencilPoints
515  );
516 
517  forAll(stencilPoints, cellI)
518  {
520  (
521  runTime.path()/"centredCFCCell" + Foam::name(cellI) + ".obj",
522  mesh.cellCentres()[cellI],
523  stencilPoints[cellI]
524  );
525  }
526  }
527 
528 
529 //XXXXXX
530 // // Evaluate
531 // List<List<scalar> > stencilData(faceStencils.size());
532 // collectStencilData
533 // (
534 // distMap,
535 // faceStencils,
536 // vf,
537 // stencilData
538 // );
539 // for (label faci = 0; faci < mesh.nInternalFaces(); faci++)
540 // {
541 // const scalarList& stData = stencilData[faceI];
542 // const scalarList& stWeight = fit[faceI];
543 //
544 // forAll(stData, i)
545 // {
546 // sf[faceI] += stWeight[i]*stData[i];
547 // }
548 // }
549 // See finiteVolume/lnInclude/leastSquaresGrad.C
550 //XXXXXX
551 
552  Pout<< "End\n" << endl;
553 
554  return 0;
555 }
556 
557 
558 // ************************************************************************* //
volFields.H
Foam::maxOp
Definition: ops.H:172
meshTools.H
centredCFCCellToCellStencilObject.H
Foam::fileName
A class for handling file names.
Definition: fileName.H:69
writeStencilStats
void writeStencilStats(const labelListList &stencil)
Definition: Test-ExtendedStencil.C:86
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
Foam::meshTools::writeOBJ
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Foam::primitiveMesh::edgeFaces
const labelListList & edgeFaces() const
Definition: primitiveMeshEdgeFaces.C:31
Foam::labelMax
static const label labelMax
Definition: label.H:62
Foam::minOp
Definition: ops.H:173
Foam::extendedCentredCellToCellStencil
Definition: extendedCentredCellToCellStencil.H:50
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
addTimeOptions.H
Foam::MeshObject< fvMesh, TopologicalMeshObject, centredCECCellToCellStencilObject >::New
static const centredCECCellToCellStencilObject & New(const fvMesh &mesh)
Definition: MeshObject.C:44
OFstream.H
Foam::primitiveMesh::nCells
label nCells() const
Definition: primitiveMeshI.H:64
main
int main(int argc, char *argv[])
Definition: Test-ExtendedStencil.C:122
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:43
Foam::labelMin
static const label labelMin
Definition: label.H:61
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
writeStencilOBJ
void writeStencilOBJ(const fileName &fName, const point &fc, const List< point > &stencilCc)
Definition: Test-ExtendedStencil.C:64
Foam::nl
static const char nl
Definition: Ostream.H:260
Foam::Info
messageStream Info
centredCPCCellToCellStencilObject.H
argList.H
Foam::primitiveMesh::edgeCells
const labelListList & edgeCells() const
Definition: primitiveMeshEdgeCells.C:32
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
fvMesh.H
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Foam::extendedCentredCellToCellStencil::stencil
const labelListList & stencil() const
Return reference to the stencil.
Definition: extendedCentredCellToCellStencil.H:92
checkTimeOptions.H
Foam::OFstream
Output to file stream.
Definition: OFstream.H:81
Foam::max
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
setRootCase.H
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:211
centredCECCellToCellStencilObject.H
Foam::fvMesh::C
const volVectorField & C() const
Return cell centres as volVectorField.
Definition: fvMeshGeometry.C:369
Foam::sumOp
Definition: ops.H:162
Foam::Pout
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Foam::Vector< scalar >
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
createMesh.H
Foam::primitiveMesh::pointCells
const labelListList & pointCells() const
Definition: primitiveMeshPointCells.C:108
createTime.H
startTime
Foam::label startTime
Definition: checkTimeOptions.H:5
Foam::List::size
void size(const label)
Override size to be inconsistent with allocated storage.
Foam::extendedCentredCellToCellStencil::collectData
void collectData(const GeometricField< Type, fvPatchField, volMesh > &fld, List< List< Type > > &stencilFld) const
Use map to get the data into stencil order.
Definition: extendedCentredCellToCellStencil.H:103
Foam::min
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47