viewFactorsGen.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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  viewFactorsGen
26 
27 Description
28  View factors are calculated based on a face agglomeration array
29  (finalAgglom generated by faceAgglomerate utility).
30 
31  Each view factor between the agglomerated faces i and j (Fij) is calculated
32  using a double integral of the sub-areas composing the agglomaration.
33 
34  The patches involved in the view factor calculation are taken from the qr
35  volScalarField (radiative flux) when is greyDiffusiveRadiationViewFactor
36  otherwise they are not included.
37 
38 \*---------------------------------------------------------------------------*/
39 
40 #include "argList.H"
41 #include "Time.H"
42 #include "fvMesh.H"
43 #include "singleCellFvMesh.H"
44 #include "volFields.H"
45 #include "surfaceFields.H"
48 #include "cyclicAMIPolyPatch.H"
49 #include "mapDistribute.H"
50 #include "meshTools.H"
52 #include "DynamicField.H"
53 #include "scalarMatrices.H"
54 #include "scalarListIOList.H"
55 
56 using namespace Foam;
57 
58 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
59 
60 triSurface triangulate
61 (
62  const polyBoundaryMesh& bMesh,
63  const labelHashSet& includePatches,
64  const labelListIOList& finalAgglom,
65  labelList& triSurfaceToAgglom,
66  const globalIndex& globalNumbering,
67  const polyBoundaryMesh& coarsePatches
68 )
69 {
70  const polyMesh& mesh = bMesh.mesh();
71 
72  // Storage for surfaceMesh. Size estimate.
73  DynamicList<labelledTri> triangles
74  (
75  mesh.nFaces() - mesh.nInternalFaces()
76  );
77 
78  label newPatchi = 0;
79  label localTriFacei = 0;
80 
81  forAllConstIter(labelHashSet, includePatches, iter)
82  {
83  const label patchi = iter.key();
84  const polyPatch& patch = bMesh[patchi];
85  const pointField& points = patch.points();
86 
87  label nTriTotal = 0;
88 
89  forAll(patch, patchFacei)
90  {
91  const face& f = patch[patchFacei];
92 
93  faceList triFaces(f.nTriangles(points));
94 
95  label nTri = 0;
96 
97  f.triangles(points, nTri, triFaces);
98 
99  forAll(triFaces, triFacei)
100  {
101  const face& f = triFaces[triFacei];
102 
103  triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
104 
105  nTriTotal++;
106 
107  triSurfaceToAgglom[localTriFacei++] = globalNumbering.toGlobal
108  (
110  finalAgglom[patchi][patchFacei]
111  + coarsePatches[patchi].start()
112  );
113  }
114  }
115 
116  newPatchi++;
117  }
118 
119  triSurfaceToAgglom.resize(localTriFacei);
120 
121  triangles.shrink();
122 
123  // Create globally numbered tri surface
124  triSurface rawSurface(triangles, mesh.points());
125 
126  // Create locally numbered tri surface
127  triSurface surface
128  (
129  rawSurface.localFaces(),
130  rawSurface.localPoints()
131  );
132 
133  // Add patch names to surface
134  surface.patches().setSize(newPatchi);
135 
136  newPatchi = 0;
137 
138  forAllConstIter(labelHashSet, includePatches, iter)
139  {
140  const label patchi = iter.key();
141  const polyPatch& patch = bMesh[patchi];
142 
143  surface.patches()[newPatchi].index() = patchi;
144  surface.patches()[newPatchi].name() = patch.name();
145  surface.patches()[newPatchi].geometricType() = patch.type();
146 
147  newPatchi++;
148  }
149 
150  return surface;
151 }
152 
153 
154 void writeRays
155 (
156  const fileName& fName,
157  const pointField& compactCf,
158  const pointField& myFc,
159  const labelListList& visibleFaceFaces
160 )
161 {
162  OFstream str(fName);
163  label vertI = 0;
164 
165  Pout<< "Dumping rays to " << str.name() << endl;
166 
167  forAll(myFc, facei)
168  {
169  const labelList visFaces = visibleFaceFaces[facei];
170  forAll(visFaces, faceRemote)
171  {
172  label compactI = visFaces[faceRemote];
173  const point& remoteFc = compactCf[compactI];
174 
175  meshTools::writeOBJ(str, myFc[facei]);
176  vertI++;
177  meshTools::writeOBJ(str, remoteFc);
178  vertI++;
179  str << "l " << vertI-1 << ' ' << vertI << nl;
180  }
181  }
182  string cmd("objToVTK " + fName + " " + fName.lessExt() + ".vtk");
183  Pout<< "cmd:" << cmd << endl;
184  system(cmd);
185 }
186 
187 
188 scalar calculateViewFactorFij
189 (
190  const vector& i,
191  const vector& j,
192  const vector& dAi,
193  const vector& dAj
194 )
195 {
196  vector r = i - j;
197  scalar rMag = mag(r);
198 
199  if (rMag > small)
200  {
201  scalar dAiMag = mag(dAi);
202  scalar dAjMag = mag(dAj);
203 
204  vector ni = dAi/dAiMag;
205  vector nj = dAj/dAjMag;
206  scalar cosThetaJ = mag(nj & r)/rMag;
207  scalar cosThetaI = mag(ni & r)/rMag;
208 
209  return
210  (
211  (cosThetaI*cosThetaJ*dAjMag*dAiMag)
213  );
214  }
215  else
216  {
217  return 0;
218  }
219 }
220 
221 
222 void insertMatrixElements
223 (
224  const globalIndex& globalNumbering,
225  const label fromProci,
226  const labelListList& globalFaceFaces,
227  const scalarListList& viewFactors,
228  scalarSquareMatrix& matrix
229 )
230 {
231  forAll(viewFactors, facei)
232  {
233  const scalarList& vf = viewFactors[facei];
234  const labelList& globalFaces = globalFaceFaces[facei];
235 
236  label globalI = globalNumbering.toGlobal(fromProci, facei);
237  forAll(globalFaces, i)
238  {
239  matrix[globalI][globalFaces[i]] = vf[i];
240  }
241  }
242 }
243 
244 
245 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
246 
247 int main(int argc, char *argv[])
248 {
249  #include "addRegionOption.H"
250  #include "setRootCase.H"
251  #include "createTime.H"
252  #include "createNamedMesh.H"
253 
254  // Read view factor dictionary
255  IOdictionary viewFactorDict
256  (
257  IOobject
258  (
259  "viewFactorsDict",
260  runTime.constant(),
261  mesh,
264  )
265  );
266 
267  const bool writeViewFactors =
268  viewFactorDict.lookupOrDefault<bool>("writeViewFactorMatrix", false);
269 
270  const bool dumpRays =
271  viewFactorDict.lookupOrDefault<bool>("dumpRays", false);
272 
273  const label debug = viewFactorDict.lookupOrDefault<label>("debug", 0);
274 
275  volScalarField qr
276  (
277  IOobject
278  (
279  "qr",
280  runTime.timeName(),
281  mesh,
284  ),
285  mesh
286  );
287 
288  // Read agglomeration map
289  labelListIOList finalAgglom
290  (
291  IOobject
292  (
293  "finalAgglom",
294  mesh.facesInstance(),
295  mesh,
298  false
299  )
300  );
301 
302  // Create the coarse mesh using agglomeration
303  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
304 
305  if (debug)
306  {
307  Pout << "\nCreating single cell mesh..." << endl;
308  }
309 
310  singleCellFvMesh coarseMesh
311  (
312  IOobject
313  (
314  "coarse:" + mesh.name(),
315  runTime.timeName(),
316  runTime,
319  ),
320  mesh,
321  finalAgglom
322  );
323 
324  if (debug)
325  {
326  Pout << "\nCreated single cell mesh..." << endl;
327  }
328 
329 
330  // Calculate total number of fine and coarse faces
331  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
332 
333  label nCoarseFaces = 0; // total number of coarse faces
334  label nFineFaces = 0; // total number of fine faces
335 
336  const polyBoundaryMesh& patches = mesh.boundaryMesh();
337  const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
338 
339  labelList viewFactorsPatches(patches.size());
340 
341  const volScalarField::Boundary& qrb = qr.boundaryField();
342 
343  label count = 0;
344  forAll(qrb, patchi)
345  {
346  const polyPatch& pp = patches[patchi];
347  const fvPatchScalarField& qrpI = qrb[patchi];
348 
349  if ((isA<fixedValueFvPatchScalarField>(qrpI)) && (pp.size() > 0))
350  {
351  viewFactorsPatches[count] = qrpI.patch().index();
352  nCoarseFaces += coarsePatches[patchi].size();
353  nFineFaces += patches[patchi].size();
354  count ++;
355  }
356  }
357 
358  viewFactorsPatches.resize(count);
359 
360  // total number of coarse faces
361  label totalNCoarseFaces = nCoarseFaces;
362 
363  reduce(totalNCoarseFaces, sumOp<label>());
364 
365  if (Pstream::master())
366  {
367  Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
368  }
369 
370  if (Pstream::master() && debug)
371  {
372  Pout << "\nView factor patches included in the calculation : "
373  << viewFactorsPatches << endl;
374  }
375 
376  // Collect local Cf and Sf on coarse mesh
377  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378 
379  DynamicList<point> localCoarseCf(nCoarseFaces);
380  DynamicList<point> localCoarseSf(nCoarseFaces);
381 
382  DynamicList<label> localAgg(nCoarseFaces);
383 
384  forAll(viewFactorsPatches, i)
385  {
386  const label patchID = viewFactorsPatches[i];
387 
388  const polyPatch& pp = patches[patchID];
389  const labelList& agglom = finalAgglom[patchID];
390  label nAgglom = max(agglom)+1;
391  labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
392  const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
393 
394  const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
395  const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
396 
398  includePatches.insert(patchID);
399 
400  forAll(coarseCf, facei)
401  {
402  point cf = coarseCf[facei];
403 
404  const label coarseFacei = coarsePatchFace[facei];
405  const labelList& fineFaces = coarseToFine[coarseFacei];
406  const label agglomI =
407  agglom[fineFaces[0]] + coarsePatches[patchID].start();
408 
409  // Construct single face
411  (
412  UIndirectList<face>(pp, fineFaces),
413  pp.points()
414  );
415 
416 
417  List<point> availablePoints
418  (
419  upp.faceCentres().size()
420  + upp.localPoints().size()
421  );
422 
424  (
425  availablePoints,
426  upp.faceCentres().size()
427  ) = upp.faceCentres();
428 
430  (
431  availablePoints,
432  upp.localPoints().size(),
433  upp.faceCentres().size()
434  ) = upp.localPoints();
435 
436  point cfo = cf;
437  scalar dist = great;
438  forAll(availablePoints, iPoint)
439  {
440  point cfFine = availablePoints[iPoint];
441  if (mag(cfFine-cfo) < dist)
442  {
443  dist = mag(cfFine-cfo);
444  cf = cfFine;
445  }
446  }
447 
448  point sf = coarseSf[facei];
449  localCoarseCf.append(cf);
450  localCoarseSf.append(sf);
451  localAgg.append(agglomI);
452  }
453  }
454 
455 
456  // Distribute local coarse Cf and Sf for shooting rays
457  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
458 
459  List<pointField> remoteCoarseCf(Pstream::nProcs());
460  List<pointField> remoteCoarseSf(Pstream::nProcs());
461  List<labelField> remoteCoarseAgg(Pstream::nProcs());
462 
463  remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
464  remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
465  remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
466 
467  Pstream::gatherList(remoteCoarseCf);
468  Pstream::scatterList(remoteCoarseCf);
469  Pstream::gatherList(remoteCoarseSf);
470  Pstream::scatterList(remoteCoarseSf);
471  Pstream::gatherList(remoteCoarseAgg);
472  Pstream::scatterList(remoteCoarseAgg);
473 
474 
475  globalIndex globalNumbering(nCoarseFaces);
476 
477  // Set up searching engine for obstacles
478  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
479  #include "searchingEngine.H"
480 
481 
482  // Determine rays between coarse face centres
483  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
484  DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
485 
486  DynamicList<label> rayEndFace(rayStartFace.size());
487 
488 
489  // Return rayStartFace in local index andrayEndFace in global index
490  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
491 
492  #include "shootRays.H"
493 
494  // Calculate number of visible faces from local index
495  labelList nVisibleFaceFaces(nCoarseFaces, 0);
496 
497  forAll(rayStartFace, i)
498  {
499  nVisibleFaceFaces[rayStartFace[i]]++;
500  }
501 
502  labelListList visibleFaceFaces(nCoarseFaces);
503 
504  label nViewFactors = 0;
505  forAll(nVisibleFaceFaces, facei)
506  {
507  visibleFaceFaces[facei].setSize(nVisibleFaceFaces[facei]);
508  nViewFactors += nVisibleFaceFaces[facei];
509  }
510 
511 
512  // - Construct compact numbering
513  // - return map from remote to compact indices
514  // (per processor (!= myProcNo) a map from remote index to compact index)
515  // - construct distribute map
516  // - renumber rayEndFace into compact addressing
517 
518  List<Map<label>> compactMap(Pstream::nProcs());
519 
520  mapDistribute map(globalNumbering, rayEndFace, compactMap);
521 
522  labelListIOList IOsubMap
523  (
524  IOobject
525  (
526  "subMap",
527  mesh.facesInstance(),
528  mesh,
531  false
532  ),
533  map.subMap()
534  );
535  IOsubMap.write();
536 
537 
538  labelListIOList IOconstructMap
539  (
540  IOobject
541  (
542  "constructMap",
543  mesh.facesInstance(),
544  mesh,
547  false
548  ),
549  map.constructMap()
550  );
551  IOconstructMap.write();
552 
553 
554  IOList<label> consMapDim
555  (
556  IOobject
557  (
558  "constructMapDim",
559  mesh.facesInstance(),
560  mesh,
563  false
564  ),
565  List<label>(1, map.constructSize())
566  );
567  consMapDim.write();
568 
569 
570  // visibleFaceFaces has:
571  // (local face, local viewed face) = compact viewed face
572  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
573 
574  nVisibleFaceFaces = 0;
575  forAll(rayStartFace, i)
576  {
577  label facei = rayStartFace[i];
578  label compactI = rayEndFace[i];
579  visibleFaceFaces[facei][nVisibleFaceFaces[facei]++] = compactI;
580  }
581 
582 
583  // Construct data in compact addressing
584  // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
585  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
586 
587  pointField compactCoarseCf(map.constructSize(), Zero);
588  pointField compactCoarseSf(map.constructSize(), Zero);
589  List<List<point>> compactFineSf(map.constructSize());
590  List<List<point>> compactFineCf(map.constructSize());
591 
592  DynamicList<label> compactPatchId(map.constructSize());
593 
594  // Insert my coarse local values
595  SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
596  SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
597 
598  // Insert my fine local values
599  label compactI = 0;
600  forAll(viewFactorsPatches, i)
601  {
602  label patchID = viewFactorsPatches[i];
603  const labelList& agglom = finalAgglom[patchID];
604  label nAgglom = max(agglom)+1;
605  labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
606  const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
607 
608  forAll(coarseToFine, coarseI)
609  {
610  compactPatchId.append(patchID);
611  List<point>& fineCf = compactFineCf[compactI];
612  List<point>& fineSf = compactFineSf[compactI++];
613 
614  const label coarseFacei = coarsePatchFace[coarseI];
615  const labelList& fineFaces = coarseToFine[coarseFacei];
616 
617  fineCf.setSize(fineFaces.size());
618  fineSf.setSize(fineFaces.size());
619 
620  fineCf = UIndirectList<point>
621  (
622  mesh.Cf().boundaryField()[patchID],
623  coarseToFine[coarseFacei]
624  );
625  fineSf = UIndirectList<point>
626  (
627  mesh.Sf().boundaryField()[patchID],
628  coarseToFine[coarseFacei]
629  );
630  }
631  }
632 
633  // Do all swapping
634  map.distribute(compactCoarseSf);
635  map.distribute(compactCoarseCf);
636  map.distribute(compactFineCf);
637  map.distribute(compactFineSf);
638 
639  map.distribute(compactPatchId);
640 
641 
642  // Plot all rays between visible faces.
643  if (dumpRays)
644  {
645  writeRays
646  (
647  runTime.path()/"allVisibleFaces.obj",
648  compactCoarseCf,
649  remoteCoarseCf[Pstream::myProcNo()],
650  visibleFaceFaces
651  );
652  }
653 
654 
655  // Fill local view factor matrix
656  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
657 
659  (
660  IOobject
661  (
662  "F",
663  mesh.facesInstance(),
664  mesh,
667  false
668  ),
669  nCoarseFaces
670  );
671 
672  label totalPatches = coarsePatches.size();
673  reduce(totalPatches, maxOp<label>());
674 
675  // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
676  scalarSquareMatrix sumViewFactorPatch
677  (
678  totalPatches,
679  0.0
680  );
681 
682  scalarList patchArea(totalPatches, 0.0);
683 
684  if (Pstream::master())
685  {
686  Info<< "\nCalculating view factors..." << endl;
687  }
688 
689  if (mesh.nSolutionD() == 3)
690  {
691  forAll(localCoarseSf, coarseFacei)
692  {
693  const List<point>& localFineSf = compactFineSf[coarseFacei];
694  const vector Ai = sum(localFineSf);
695  const List<point>& localFineCf = compactFineCf[coarseFacei];
696  const label fromPatchId = compactPatchId[coarseFacei];
697  patchArea[fromPatchId] += mag(Ai);
698 
699  const labelList& visCoarseFaces = visibleFaceFaces[coarseFacei];
700 
701  forAll(visCoarseFaces, visCoarseFacei)
702  {
703  F[coarseFacei].setSize(visCoarseFaces.size());
704  label compactJ = visCoarseFaces[visCoarseFacei];
705  const List<point>& remoteFineSj = compactFineSf[compactJ];
706  const List<point>& remoteFineCj = compactFineCf[compactJ];
707 
708  const label toPatchId = compactPatchId[compactJ];
709 
710  scalar Fij = 0;
711  forAll(localFineSf, i)
712  {
713  const vector& dAi = localFineSf[i];
714  const vector& dCi = localFineCf[i];
715 
716  forAll(remoteFineSj, j)
717  {
718  const vector& dAj = remoteFineSj[j];
719  const vector& dCj = remoteFineCj[j];
720 
721  scalar dIntFij = calculateViewFactorFij
722  (
723  dCi,
724  dCj,
725  dAi,
726  dAj
727  );
728 
729  Fij += dIntFij;
730  }
731  }
732  F[coarseFacei][visCoarseFacei] = Fij/mag(Ai);
733  sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
734  }
735  }
736  }
737  else if (mesh.nSolutionD() == 2)
738  {
739  const boundBox& box = mesh.bounds();
740  const Vector<label>& dirs = mesh.geometricD();
741  vector emptyDir = Zero;
742  forAll(dirs, i)
743  {
744  if (dirs[i] == -1)
745  {
746  emptyDir[i] = 1.0;
747  }
748  }
749 
750  scalar wideBy2 = (box.span() & emptyDir)*2.0;
751 
752  forAll(localCoarseSf, coarseFacei)
753  {
754  const vector& Ai = localCoarseSf[coarseFacei];
755  const vector& Ci = localCoarseCf[coarseFacei];
756  vector Ain = Ai/mag(Ai);
757  vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
758  vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
759 
760  const label fromPatchId = compactPatchId[coarseFacei];
761  patchArea[fromPatchId] += mag(Ai);
762 
763  const labelList& visCoarseFaces = visibleFaceFaces[coarseFacei];
764  forAll(visCoarseFaces, visCoarseFacei)
765  {
766  F[coarseFacei].setSize(visCoarseFaces.size());
767  label compactJ = visCoarseFaces[visCoarseFacei];
768  const vector& Aj = compactCoarseSf[compactJ];
769  const vector& Cj = compactCoarseCf[compactJ];
770 
771  const label toPatchId = compactPatchId[compactJ];
772 
773  vector Ajn = Aj/mag(Aj);
774  vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
775  vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
776 
777  scalar d1 = mag(R1i - R2j);
778  scalar d2 = mag(R2i - R1j);
779  scalar s1 = mag(R1i - R1j);
780  scalar s2 = mag(R2i - R2j);
781 
782  scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
783 
784  F[coarseFacei][visCoarseFacei] = Fij;
785  sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
786  }
787  }
788  }
789 
790  if (Pstream::master())
791  {
792  Info << "Writing view factor matrix..." << endl;
793  }
794 
795  // Write view factors matrix in listlist form
796  F.write();
797 
798  reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
799  reduce(patchArea, sumOp<scalarList>());
800 
801 
802  if (Pstream::master() && debug)
803  {
804  forAll(viewFactorsPatches, i)
805  {
806  label patchi = viewFactorsPatches[i];
807  forAll(viewFactorsPatches, i)
808  {
809  label patchJ = viewFactorsPatches[i];
810  Info << "F" << patchi << patchJ << ": "
811  << sumViewFactorPatch[patchi][patchJ]/patchArea[patchi]
812  << endl;
813  }
814  }
815  }
816 
817 
818  if (writeViewFactors)
819  {
820  volScalarField viewFactorField
821  (
822  IOobject
823  (
824  "viewFactorField",
825  mesh.time().timeName(),
826  mesh,
829  ),
830  mesh,
832  );
833 
834  volScalarField::Boundary& viewFactorFieldBf =
835  viewFactorField.boundaryFieldRef();
836 
837  label compactI = 0;
838  forAll(viewFactorsPatches, i)
839  {
840  label patchID = viewFactorsPatches[i];
841  const labelList& agglom = finalAgglom[patchID];
842  label nAgglom = max(agglom)+1;
843  labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
844  const labelList& coarsePatchFace =
845  coarseMesh.patchFaceMap()[patchID];
846 
847  forAll(coarseToFine, coarseI)
848  {
849  const scalar Fij = sum(F[compactI]);
850  const label coarseFaceID = coarsePatchFace[coarseI];
851  const labelList& fineFaces = coarseToFine[coarseFaceID];
852  forAll(fineFaces, fineId)
853  {
854  const label faceID = fineFaces[fineId];
855  viewFactorFieldBf[patchID][faceID] = Fij;
856  }
857  compactI++;
858  }
859  }
860  viewFactorField.write();
861  }
862 
863 
864  // Invert compactMap (from processor+localface to compact) to go
865  // from compact to processor+localface (expressed as a globalIndex)
866  // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
867  labelList compactToGlobal(map.constructSize());
868 
869  // Local indices first (note: are not in compactMap)
870  for (label i = 0; i < globalNumbering.localSize(); i++)
871  {
872  compactToGlobal[i] = globalNumbering.toGlobal(i);
873  }
874 
875 
876  forAll(compactMap, proci)
877  {
878  const Map<label>& localToCompactMap = compactMap[proci];
879 
880  forAllConstIter(Map<label>, localToCompactMap, iter)
881  {
882  compactToGlobal[iter()] = globalNumbering.toGlobal
883  (
884  proci,
885  iter.key()
886  );
887  }
888  }
889 
890 
891  if (Pstream::master())
892  {
893  scalarSquareMatrix Fmatrix(totalNCoarseFaces, 0.0);
894 
895  labelListList globalFaceFaces(visibleFaceFaces.size());
896 
897  // Create globalFaceFaces needed to insert view factors
898  // in F to the global matrix Fmatrix
899  forAll(globalFaceFaces, facei)
900  {
901  globalFaceFaces[facei] = renumber
902  (
903  compactToGlobal,
904  visibleFaceFaces[facei]
905  );
906  }
907 
908  labelListIOList IOglobalFaceFaces
909  (
910  IOobject
911  (
912  "globalFaceFaces",
913  mesh.facesInstance(),
914  mesh,
917  false
918  ),
919  globalFaceFaces
920  );
921  IOglobalFaceFaces.write();
922  }
923  else
924  {
925  labelListList globalFaceFaces(visibleFaceFaces.size());
926  forAll(globalFaceFaces, facei)
927  {
928  globalFaceFaces[facei] = renumber
929  (
930  compactToGlobal,
931  visibleFaceFaces[facei]
932  );
933  }
934 
935  labelListIOList IOglobalFaceFaces
936  (
937  IOobject
938  (
939  "globalFaceFaces",
940  mesh.facesInstance(),
941  mesh,
944  false
945  ),
946  globalFaceFaces
947  );
948 
949  IOglobalFaceFaces.write();
950  }
951 
952  Info<< "End\n" << endl;
953  return 0;
954 }
955 
956 
957 // ************************************************************************* //
Foam::surfaceFields.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:82
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
const dimensionedScalar & F
Faraday constant: default SI units: [C/mol].
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
fvMesh as subset of other mesh. Consists of one cell and all original bounday faces. Useful when manipulating boundary data. Single internal cell only needed to be able to manipulate in a standard way.
fileName path() const
Return path.
Definition: Time.H:265
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
const word & name() const
Return name.
const word & name() const
Return name.
Definition: IOobject.H:303
A class for handling file names.
Definition: fileName.H:79
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:808
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
label index() const
Return the index of this patch in the fvBoundaryMesh.
Definition: fvPatch.H:173
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
label nInternalFaces() const
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Output to file stream.
Definition: OFstream.H:82
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:622
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
label triangles(const pointField &points, label &triI, faceList &triFaces) const
Split into triangles using existing points.
Definition: face.C:724
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:814
A list of faces which address into the list of points.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
A List obtained as a section of another List.
Definition: SubList.H:53
dynamicFvMesh & mesh
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
const Field< PointType > & points() const
Return reference to global points.
const polyMesh & mesh() const
Return the mesh reference.
static const zero Zero
Definition: zero.H:97
Triangle with additional region number.
Definition: labelledTri.H:57
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:325
Foam::polyBoundaryMesh.
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
labelHashSet includePatches
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
Definition: polyMesh.C:842
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setSize(const label)
Reset size of List.
Definition: List.C:281
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:440
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
label patchi
Class containing processor-to-processor mapping information.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
fileName lessExt() const
Return file name without extension (part before last .)
Definition: fileName.C:272
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:84
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label nTriangles() const
Number of triangles after splitting.
Definition: faceI.H:130
virtual bool write(const bool write=true) const
Write using setting from DB.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
Triangulated surface description with patch information.
Definition: triSurface.H:66
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Namespace for OpenFOAM.
int system(const std::string &command)
Execute the specified command.
Definition: POSIX.C:1234
label localSize() const
My local size.
Definition: globalIndexI.H:60