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