ccmReaderMesh.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2016-2021 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*----------------------------------------------------------------------------*/
27 
28 #include "ccmReader.H"
29 
30 #include "emptyPolyPatch.H"
31 #include "symmetryPolyPatch.H"
32 #include "wallPolyPatch.H"
33 #include "Fstream.H"
34 #include "IOdictionary.H"
35 
36 #include "ccmBoundaryInfo.H"
38 #include "SortableList.H"
39 #include "mergePoints.H"
40 #include "bitSet.H"
41 #include "ListOps.H"
42 
43 #include "ccmInternal.H" // include last to avoid any strange interactions
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 Foam::labelList Foam::ccm::reader::patchStartList(label initial) const
48 {
49  labelList startLst(patchSizes_.size(), Zero);
50 
51  label patchFaceI = initial;
52  forAll(patchSizes_, patchI)
53  {
54  startLst[patchI] = patchFaceI;
55  patchFaceI += patchSizes_[patchI];
56  }
57 
58  return startLst;
59 }
60 
61 
62 void Foam::ccm::reader::printSizes() const
63 {
64  Info<<"nPoints:" << nPoints_
65  << " nCells:" << nCells_
66  << " nFaces:" << nFaces_
67  << " nInternalFaces:" << nInternalFaces_
68  << endl;
69 }
70 
71 
72 // Determine if the geometry looks good.
73 // We'll insist that everything be on the first ("default") state,
74 // and is on the same file, and has a problem description
75 //
76 // The detection reads the problem description as well
77 //
78 bool Foam::ccm::reader::detectGeometry()
79 {
80  // Call once
81  if (geometryStatus_ != UNKNOWN)
82  {
83  return (geometryStatus_ == OKAY || geometryStatus_ == READ);
84  }
85 
86  // Explicitly restricted to 'default' state node
87  ccmID stateNode;
88 
89  // Get associated problem description
90  ccmID probNode;
91 
92  // Only check the first processor
93  ccmID processorNode;
94  int procI = 0;
95 
96  // Geometry needs vertices and topology
97  ccmID verticesNode, topoNode;
98 
99  if
100  (
101  CCMIOGetState
102  (
103  nullptr,
104  (globalState_->root),
105  "default",
106  &probNode,
107  &stateNode
108  )
109  == kCCMIONoErr
110 
111  && CCMIONextEntity
112  (
113  nullptr,
114  stateNode,
115  kCCMIOProcessor,
116  &procI,
117  &processorNode
118  )
119  == kCCMIONoErr
120 
121  && CCMIOReadProcessor
122  (
123  nullptr,
124  processorNode,
125  &verticesNode,
126  &topoNode,
127  nullptr, // Ignore initialField
128  nullptr // Ignore solutionNode
129  )
130  == kCCMIONoErr
131 
132  && CCMIOIsFromSameFile((globalState_->root), verticesNode)
133  && CCMIOIsFromSameFile((globalState_->root), topoNode)
134  && CCMIOIsValidEntity(probNode)
135  )
136  {
137  readProblemDescription(probNode);
138  geometryStatus_ = OKAY;
139  }
140  else
141  {
142  // Missing/incomplete geometry node and/or problem node
143  geometryStatus_ = BAD;
144  }
145 
146  return (geometryStatus_ == OKAY || geometryStatus_ == READ);
147 }
148 
149 
150 // read the geometry without any sorting
151 void Foam::ccm::reader::readMeshTopology
152 (
153  const scalar scaleFactor
154 )
155 {
156  ccmID verticesNode, topoNode;
157 
158  // Use first ("default") state node
159  ccmID stateNode;
160  int stateI = 0;
161 
162  // Use the first processor to find the mesh nodes
163  ccmID processorNode;
164  int procI = 0;
165 
166  if
167  (
168  CCMIONextEntity
169  (
170  nullptr,
171  (globalState_->root),
172  kCCMIOState,
173  &stateI,
174  &stateNode
175  )
176  == kCCMIONoErr
177 
178  && CCMIONextEntity
179  (
180  nullptr,
181  stateNode,
182  kCCMIOProcessor,
183  &procI,
184  &processorNode
185  )
186  == kCCMIONoErr
187 
188  && CCMIOReadProcessor
189  (
190  nullptr,
191  processorNode,
192  &verticesNode,
193  &topoNode,
194  nullptr, // Ignore initialField
195  nullptr // Ignore solutionNode
196  )
197  == kCCMIONoErr
198  )
199  {
200  labelList origPointId = readVertices(verticesNode, scaleFactor);
201  readCells(topoNode);
202  readMonitoring(topoNode);
203 
204  // Renumber vertex labels (ccm -> Foam)
205  {
206  label maxId = max(origPointId);
207  labelList mapToFoam(invert(maxId+1, origPointId));
208 
209  for (face& f : faces_)
210  {
211  inplaceRenumber(mapToFoam, f);
212  }
213  }
214  origPointId.clear();
215 
216  // Renumber owners/neighbours cell labels (ccm -> Foam)
217  {
218  label maxId = max(origCellId_);
219  labelList mapToFoam(invert(maxId+1, origCellId_));
220 
221  inplaceRenumber(mapToFoam, faceOwner_);
222  inplaceRenumber(mapToFoam, faceNeighbour_);
223  }
224 
225  // Juggle solids into fluid as required
226  juggleSolids();
227 
228  // Report sizes
229  printSizes();
230 
231  // Remove unwanted fluid/porous/solid types
232  removeUnwanted();
233 
234  // Collapse interfaces between domains (eg, fluid|porosity)
235  cleanupInterfaces();
236 
237  // Use point merge to join conformal interfaces (STARCCM)
238  mergeInplaceInterfaces();
239  }
240 }
241 
242 
243 // readVertices:
244 // 1) read the vertex data
245 // 2) read the map (which maps the index into the data array with the Id number)
246 //
247 // returns the original point Id
248 Foam::labelList Foam::ccm::reader::readVertices
249 (
250  const ccmID& verticesNode,
251  const scalar scaleFactor
252 )
253 {
254  int dims = 1;
255  float scale = 1.0;
256 
257  CCMIOSize size = 0;
258  ccmID mapId;
259 
260 #ifdef DEBUG_CCMIOREAD
261  Info<< "readVertices()" << endl;
262 #endif
263 
264  // Determine dimensions and mapId
265  CCMIOEntitySize
266  (
267  &(globalState_->error),
268  verticesNode,
269  &size,
270  nullptr
271  );
272 
273  CCMIOReadVerticesd
274  (
275  &(globalState_->error),
276  verticesNode,
277  &dims,
278  &scale,
279  &mapId,
280  nullptr,
281  kCCMIOStart,
282  kCCMIOEnd
283  );
284  assertNoError("problem finding 'Vertices' node");
285 
286  if (dims != 3)
287  {
289  << "can only handle 3-dimensional vertices"
290  << exit(FatalError);
291  }
292 
293  nPoints_ = size;
294  labelList origPointId(nPoints_);
295 
296  readMap
297  (
298  mapId,
299  origPointId
300  );
301 
302  // Temporary storage for the points - reading piecemeal is much slower
303  List<scalar> vrts(3*nPoints_);
304 
305  // The ccm data is double precision too
306  CCMIOReadVerticesd
307  (
308  &(globalState_->error),
309  verticesNode,
310  nullptr,
311  nullptr,
312  nullptr,
313  vrts.data(),
314  kCCMIOStart,
315  kCCMIOEnd
316  );
317  assertNoError("problem reading 'Vertices' node");
318 
319  // Convert to foam Points
320  points_.setSize(nPoints_);
321 
322  scalar effectiveScale = scale * scaleFactor;
323  forAll(points_, i)
324  {
325  points_[i].x() = effectiveScale * vrts[i*3];
326  points_[i].y() = effectiveScale * vrts[i*3+1];
327  points_[i].z() = effectiveScale * vrts[i*3+2];
328  }
329 
330  vrts.clear();
331 
332 #ifdef DEBUG_CCMIOREAD
333  Info<< "readVertices: " << nPoints_ << endl;
334 #endif
335 
336  return origPointId;
337 }
338 
339 
340 // readCells:
341 // - read faces, faceOwner and faceNeighbour
342 // - finally read interfaces
343 void Foam::ccm::reader::readCells
344 (
345  const ccmID& topoNode
346 )
347 {
348  CCMIOSize size = 0;
349  ccmID cellsNode, mapId;
350  ccmID nodeId;
351 
352 #ifdef DEBUG_CCMIOREAD
353  Info<< "readCells()" << endl;
354 #endif
355 
356  // Determine dimensions and mapId information for 'Cells'
357  CCMIOGetEntity
358  (
359  &(globalState_->error),
360  topoNode,
361  kCCMIOCells,
362  0,
363  &cellsNode
364  );
365 
366  CCMIOEntitySize
367  (
368  &(globalState_->error),
369  cellsNode,
370  &size,
371  nullptr
372  );
373  assertNoError("cannot get 'Cells' node");
374 
375  nCells_ = size;
376 
377 #ifdef DEBUG_CCMIOREAD
378  Info<< "readCells: " << nCells_ << endl;
379 #endif
380 
381  // Store cell ids so that we know which cells are which
382  origCellId_.setSize(nCells_);
383  cellTableId_.setSize(nCells_);
384 
385  CCMIOReadCells
386  (
387  &(globalState_->error),
388  cellsNode,
389  &mapId,
390  cellTableId_.data(),
391  kCCMIOStart,
392  kCCMIOEnd
393  );
394  assertNoError("Error reading 'Cells' node");
395 
396  readMap
397  (
398  mapId,
399  origCellId_
400  );
401 
402  // Determine dimensions and mapId information for 'InternalFaces'
403  CCMIOGetEntity
404  (
405  &(globalState_->error),
406  topoNode,
407  kCCMIOInternalFaces,
408  0,
409  &nodeId
410  );
411  CCMIOEntitySize
412  (
413  &(globalState_->error),
414  nodeId,
415  &size,
416  nullptr
417  );
418  nInternalFaces_ = size;
419 
420  nFaces_ = nInternalFaces_;
421 
422  // First pass:
423  //
424  // Determine patch sizes before reading internal faces
425  // also determine the original boundary regions
426 
427  label nPatches = 0;
428  DynamicList<ccmBoundaryInfo> bndInfo;
429 
430  // Number of children in the parent node is more than number of patches,
431  // but is a good start for allocation
432  if
433  (
434  CCMIOGetNumberOfChildren
435  (
436  nullptr,
437  topoNode.node,
438  &nPatches
439  ) == kCCMIONoErr
440  )
441  {
442  bndInfo.setCapacity(nPatches);
443  }
444 
445  for
446  (
447  int index = 0;
448  CCMIONextEntity
449  (
450  nullptr,
451  topoNode,
452  kCCMIOBoundaryFaces,
453  &index,
454  &nodeId
455  ) == kCCMIONoErr
456  && CCMIOEntitySize
457  (
458  &(globalState_->error),
459  nodeId,
460  &size,
461  nullptr
462  ) == kCCMIONoErr;
463  /* nop */
464  )
465  {
466  ccmBoundaryInfo info;
467  info.size = size;
468  nFaces_ += size;
469 
470  CCMIOGetEntityIndex
471  (
472  &(globalState_->error),
473  nodeId,
474  &(info.ccmIndex)
475  );
476 
477  // Name directly from the node (eg, STARCCM)
478  info.setPatchName(ccmReadOptstr("Label", nodeId));
479 
480  // Lookup the name, type from boundary region info:
481  auto dictIter = boundaryRegion_.find(info.ccmIndex);
482  if (dictIter.found())
483  {
484  dictionary& dict = dictIter.val();
485 
486  const word patchName(dict.get<word>("Label"));
487  const word patchType(dict.get<word>("BoundaryType"));
488 
489  if (!patchName.empty())
490  {
491  info.patchName = patchName;
492  }
493 
494  if (!patchType.empty())
495  {
496  info.patchType = patchType;
497  }
498 
499  // Optional, but potentially useful information:
500  dict.add("BoundaryIndex", info.ccmIndex);
501  dict.add("size", info.size);
502  }
503 
504  bndInfo.append(info);
505  }
506 
507  // Redimension lists according to the overall sizes
508  faces_.setSize(nFaces_);
509  faceOwner_.setSize(nFaces_);
510  faceNeighbour_.setSize(nFaces_);
511  origFaceId_.setSize(nFaces_);
512 
513 
514  // May be too large, but is a good place start size
515  patchSizes_.setSize(bndInfo.size());
516  origBndId_.setSize(bndInfo.size());
517  patchSizes_ = 0;
518  origBndId_ = -1;
519  nPatches = 0;
520 
521  for (ccmBoundaryInfo& info : bndInfo)
522  {
523  if (info.patchId != -1)
524  {
525  // Already inserted - eg, as interface
526  origBndId_[info.patchId] = info.ccmIndex;
527  }
528  else
529  {
530  info.patchId = nPatches++;
531  origBndId_[info.patchId] = info.ccmIndex;
532  }
533 
534  patchSizes_[info.patchId] += info.size;
535  }
536 
537  // Shrink to sizes actually used
538  patchSizes_.setSize(nPatches);
539  origBndId_.setSize(nPatches);
540 
541  // Boundary info indices flattened and sorted by patchId
542 
543  IndirectList<ccmBoundaryInfo> ccmLookupOrder(bndInfo, labelList());
544  {
545  DynamicList<label> addr(bndInfo.size());
546  for (int patchI = 0; patchI < nPatches; ++patchI)
547  {
548  forAll(bndInfo, infoI)
549  {
550  if (bndInfo[infoI].patchId == patchI)
551  {
552  addr.append(infoI);
553  }
554  }
555  }
556 
557  ccmLookupOrder.addressing() = std::move(addr);
558  }
559 
560 
561  //
562  // Now we are ready to do the reading
563  //
564  CCMIOGetEntity
565  (
566  &(globalState_->error),
567  topoNode,
568  kCCMIOInternalFaces,
569  0,
570  &nodeId
571  );
572 
573  // get allocation sizes
574  CCMIOReadFaces
575  (
576  &(globalState_->error),
577  nodeId,
578  kCCMIOInternalFaces,
579  &mapId,
580  &size,
581  nullptr,
582  kCCMIOStart,
583  kCCMIOEnd
584  );
585 
586  List<label> mapData(nInternalFaces_);
587  List<label> faceCells(2*nInternalFaces_);
588  List<label> ccmFaces(size);
589 
590  CCMIOReadFaces
591  (
592  &(globalState_->error),
593  nodeId,
594  kCCMIOInternalFaces,
595  nullptr,
596  nullptr,
597  ccmFaces.data(),
598  kCCMIOStart,
599  kCCMIOEnd
600  );
601 
602  CCMIOReadFaceCells
603  (
604  &(globalState_->error),
605  nodeId,
606  kCCMIOInternalFaces,
607  faceCells.data(),
608  kCCMIOStart,
609  kCCMIOEnd
610  );
611  assertNoError("Error reading internal faces");
612 
613  readMap
614  (
615  mapId,
616  mapData
617  );
618 
619  // Copy into Foam list
620  // ccmFaces are organized as [nVert vrt1 .. vrtN]
621  unsigned int pos = 0;
622  for (label faceI = 0; faceI < nInternalFaces_; ++faceI)
623  {
624  origFaceId_[faceI] = mapData[faceI];
625  faceOwner_[faceI] = faceCells[2*faceI];
626  faceNeighbour_[faceI] = faceCells[2*faceI+1];
627  face& f = faces_[faceI];
628 
629  f.setSize(ccmFaces[pos++]);
630  forAll(f, fp)
631  {
632  f[fp] = ccmFaces[pos++];
633  }
634  }
635 
636  // Read the boundary faces
637  // ~~~~~~~~~~~~~~~~~~~~~~~
638  label patchFaceI = nInternalFaces_;
639 
640  for (const ccmBoundaryInfo& info : ccmLookupOrder)
641  {
642  const unsigned int patchSize = info.size;
643 
644  if
645  (
646  CCMIOGetEntity
647  (
648  &(globalState_->error),
649  topoNode,
650  kCCMIOBoundaryFaces,
651  info.ccmIndex,
652  &nodeId
653  ) == kCCMIONoErr
654  )
655  {
656  CCMIOReadFaces
657  (
658  &(globalState_->error),
659  nodeId,
660  kCCMIOBoundaryFaces,
661  &mapId,
662  &size, // size needed for the faces
663  nullptr,
664  kCCMIOStart,
665  kCCMIOEnd
666  );
667 
668  mapData.setSize(patchSize);
669  faceCells.setSize(patchSize);
670  ccmFaces.setSize(size);
671 
672  readMap
673  (
674  mapId,
675  mapData
676  );
677 
678  CCMIOReadFaces
679  (
680  &(globalState_->error),
681  nodeId,
682  kCCMIOBoundaryFaces,
683  nullptr,
684  nullptr,
685  ccmFaces.data(),
686  kCCMIOStart,
687  kCCMIOEnd
688  );
689  CCMIOReadFaceCells
690  (
691  &(globalState_->error),
692  nodeId,
693  kCCMIOBoundaryFaces,
694  faceCells.data(),
695  kCCMIOStart,
696  kCCMIOEnd
697  );
698  assertNoError
699  (
700  "Error reading boundary face cells - index "
701  + ::Foam::name(info.ccmIndex)
702  );
703 
704  // Copy into Foam list
705  // ccmFaces are organized as [nVert vrt1 .. vrtN]
706  unsigned int pos = 0;
707  for (unsigned int i = 0; i < patchSize; ++i, ++patchFaceI)
708  {
709  origFaceId_[patchFaceI] = mapData[i];
710  faceOwner_[patchFaceI] = faceCells[i];
711  faceNeighbour_[patchFaceI] = -1;
712 
713  face& f = faces_[patchFaceI];
714 
715  f.setSize(ccmFaces[pos++]);
716  forAll(f, fp)
717  {
718  f[fp] = ccmFaces[pos++];
719  }
720  }
721  }
722  else
723  {
724  assertNoError
725  (
726  "Error reading boundary faces - index "
727  + ::Foam::name(info.ccmIndex)
728  );
729  }
730  }
731 
732  readInterfaces(cellsNode);
733 }
734 
735 
736 // Read any interfaces
737 // 1) interfaces between domains (fluid/solid, fluid/porosity)
738 // 2) PROSTAR baffles
739 //
740 void Foam::ccm::reader::readInterfaces
741 (
742  const ccmID& cellsNode
743 )
744 {
745 #ifdef DEBUG_CCMIOREAD
746  Info<< "readInterfaces()" << endl;
747 #endif
748 
749  label nBaffleInterface = 0, nInterfaceTotal = 0;
750  CCMIOIndex size = 0, dims = 0;
751 
752  bafInterfaces_.clear();
753  domInterfaces_.clear();
754 
755  ccmID interfaceNode;
756 
757  if
758  (
759  CCMIOGetEntity
760  (
761  nullptr,
762  cellsNode,
763  kCCMIOInterfaces,
764  0,
765  &interfaceNode
766  )
767  == kCCMIONoErr
768 
769  && CCMIOGetOptInfo
770  (
771  nullptr,
772  interfaceNode,
773  "FaceIds",
774  nullptr,
775  &size,
776  &dims,
777  nullptr
778  )
779  == kCCMIONoErr
780 
781  && size > 0 && dims == 2
782  )
783  {
784  nInterfaceTotal = size;
785  }
786  else
787  {
788  return;
789  }
790 
791  // Get ProstarBaffles
792  // formatted as [ prostarId faceId1 faceId2 celltableId ]
793  if
794  (
795  CCMIOGetOptInfo
796  (
797  nullptr,
798  interfaceNode,
799  "ProstarBaffles",
800  nullptr,
801  &size,
802  nullptr,
803  nullptr
804  )
805  == kCCMIONoErr
806 
807  && size > 0
808  )
809  {
810  // Be paranoid - force an integral value of 4
811  nBaffleInterface = (size - size % 4) / 4;
812  }
813 
814  // Determine sizes
815  label nDomainInterface = nInterfaceTotal - nBaffleInterface;
816 
817  // Maximum dimension
818  List<int> mapData(max(2 * nInterfaceTotal, 4 * nBaffleInterface), -1);
819 
820  bafInterfaces_.setSize(nBaffleInterface);
821  domInterfaces_.setSize(nDomainInterface);
822 
823  // Face number mapping
824  label maxId = max(origFaceId_);
825  labelList toFoamFaces(invert(maxId+1, origFaceId_));
826 
827  if (nBaffleInterface > 0)
828  {
829  mapData = -1;
830 
831  CCMIOReadOpt1i
832  (
833  &(globalState_->error),
834  interfaceNode,
835  "ProstarBaffles",
836  mapData.data(),
837  kCCMIOStart,
838  kCCMIOEnd
839  );
840  assertNoError("problem reading interface 'ProstarBaffles'");
841 
842 
843  // Copy/compress the desired entries (cannot use a SubList)
844  // Transform
845  // from [ prostarId faceId1 faceId2 celltableId ]
846  // to [ faceId1 faceId2 ]
847  for (label i=0; i < nBaffleInterface; ++i)
848  {
849  mapData[i*2] = mapData[i*4+1];
850  mapData[i*2+1] = mapData[i*4+2];
851  }
852 
853  for (label i=2*nBaffleInterface; i < 4*nBaffleInterface; ++i)
854  {
855  mapData[i] = -1;
856  }
857 
858  // Translate ccm faceId -> foam faceId
859  inplaceRenumber(toFoamFaces, mapData);
860 
861  forAll(bafInterfaces_, i)
862  {
863  bafInterfaces_[i][0] = mapData[i*2];
864  bafInterfaces_[i][1] = mapData[i*2+1];
865  }
866  }
867 
868  // Baffles are not domInterfaces, use hash to skip them
869  SubList<label> subMap(mapData, 2*nBaffleInterface);
870  labelHashSet hashedFace(subMap);
871 
872  mapData = -1;
873 
874  CCMIOReadOpt2i
875  (
876  &(globalState_->error),
877  interfaceNode,
878  "FaceIds",
879  mapData.data(),
880  kCCMIOStart,
881  kCCMIOEnd
882  );
883  assertNoError("problem reading interface 'FaceIds'");
884 
885  // Translate ccm faceId -> foam faceId
886  // newer version handles negatives indices
887  inplaceRenumber(toFoamFaces, mapData);
888 
889  nDomainInterface = 0;
890  for (label i=0; i < nInterfaceTotal; ++i)
891  {
892  label face0 = mapData[i*2];
893  label face1 = mapData[i*2+1];
894 
895  if
896  (
897  !hashedFace.found(face0)
898  && !hashedFace.found(face1)
899  )
900  {
901  domInterfaces_[nDomainInterface][0] = face0;
902  domInterfaces_[nDomainInterface][1] = face1;
903  ++nDomainInterface;
904  if (nDomainInterface >= domInterfaces_.size())
905  {
906  break;
907  }
908  }
909  }
910 
911  // Truncate for extra safety
912  domInterfaces_.setSize(nDomainInterface);
913 }
914 
915 
916 // Read monitoring faces
917 //
918 void Foam::ccm::reader::readMonitoring
919 (
920  const ccmID& topoId
921 )
922 {
923 #ifdef DEBUG_CCMIOREAD
924  Info<< "readMonitoring()" << endl;
925 #endif
926 
927 #ifdef WITH_MONITORING
928  CCMIONode topoNode, monitorParent;
929  CCMIONode monitorNode;
930 
931  // CCMIOID -> CCMIONODE
932  if
933  (
934  CCMIOGetEntityNode
935  (
936  nullptr,
937  topoId,
938  &topoNode
939  )
940  == kCCMIONoErr
941 
942  // Get "/Meshes/FaceBasedTopology/MonitorBoundaryRegions"
943  && CCMIOGetNode
944  (
945  nullptr,
946  topoNode,
947  "MonitorBoundaryRegions",
948  &monitorParent
949  )
950  == kCCMIONoErr
951  )
952  {
953  labelList toFoamFaces(invert(nFaces_+1, origFaceId_));
954 
955  // Simulate CCMIONextEntity
956  for
957  (
958  int index = 0;
959  CCMIOGetNextChildWithLabel
960  (
961  nullptr,
962  monitorParent,
963  "boundaryFaces",
964  &index,
965  &monitorNode
966  ) == kCCMIONoErr;
967  /* nop */
968  )
969  {
970 #ifdef DEBUG_MONITORING
971  Info<< "index = " << index << endl;
972 #endif
973 
974  int nMonFaces = 0;
975  // Simulate EntitySize ?
976  CCMIOReadNodei
977  (
978  nullptr,
979  monitorNode,
980  "NumFaces",
981  &nMonFaces
982  );
983 
984  int ccmRegionId = ccmGetEntityIndex(monitorNode);
985 
986  ccmID mapId;
987 
988  int idVal;
989  CCMIOReadNodei
990  (
991  nullptr,
992  monitorNode,
993  "MapId",
994  &idVal
995  );
996 
997 #ifdef DEBUG_MONITORING
998  Info<< "monitoring mapId " << idVal
999  << " with nFaces = " << nMonFaces
1000  << endl;
1001 #endif
1002  // Is it risky changing parents here?
1003  CCMIOGetEntity
1004  (
1005  nullptr,
1006  (globalState_->root),
1007  kCCMIOMap,
1008  idVal,
1009  &mapId
1010  );
1011 
1012  List<label> mapData(nMonFaces);
1013  // mapData.setSize(nMonFaces);
1014  // faceCells.setSize(nMonFaces);
1015 
1016  readMap
1017  (
1018  mapId,
1019  mapData
1020  );
1021 
1022 #ifdef DEBUG_MONITORING
1023  Info<< "map: " << mapData << nl
1024  << "toFoam: " << toFoamFaces
1025  << endl;
1026 #endif
1027 
1028  // Translate ccm faceId -> foam faceId
1029  // newer version handles negatives indices
1030  inplaceRenumber(toFoamFaces, mapData);
1031 
1032 #ifdef DEBUG_MONITORING
1033  Info<< "map: " << mapData << nl
1034  << "ccmRegionId: " << ccmRegionId << endl;
1035 #endif
1036 
1037  auto iter = boundaryRegion_.cfind(ccmRegionId);
1038 
1039  word zoneName;
1040  if (iter.found())
1041  {
1042  iter().readEntry("Label", zoneName);
1043  }
1044  else
1045  {
1046  zoneName = "monitoring_" + Foam::name(ccmRegionId);
1047  }
1048 
1049  monitoringSets_.insert(zoneName, mapData);
1050 
1051  // CCMIONode subNode;
1052  //
1053  //- simulate ReadFaceCells with kCCMIOBoundaryFaces
1054  // CCMIOGetNode(nullptr, childNode, "Cells", &subNode);
1055  // CCMIORead1i
1056  // (
1057  // nullptr, subNode, faceCells.data(),
1058  // kCCMIOStart, kCCMIOEnd
1059  // );
1060  //
1061  // Info << "cells: " << faceCells << endl;
1062  }
1063  }
1064 #endif
1065 }
1066 
1067 
1068 // Move solid faces from Default_Boundary_Region -> Default_Boundary_Solid
1069 void Foam::ccm::reader::juggleSolids()
1070 {
1071  if (!option().keepSolid())
1072  {
1073  return;
1074  }
1075 
1076  // Find "Default_Boundary_Region"
1077  label defaultBoundaryRegion = boundaryRegion_.findIndex
1078  (
1079  defaultBoundaryName
1080  );
1081 
1082  // Find "Default_Boundary_Solid"
1083  label defaultBoundarySolid = boundaryRegion_.findIndex
1084  (
1085  defaultSolidBoundaryName
1086  );
1087 
1088  // Cannot do anything if Default_Boundary_Region does not exist
1089  // or if Default_Boundary_Solid already exists
1090  if
1091  (
1092  defaultBoundaryRegion < 0
1093  || defaultBoundarySolid >= 0
1094  )
1095  {
1096  return;
1097  }
1098 
1099  // Identify solid cells
1100  // ~~~~~~~~~~~~~~~~~~~~
1101  label nSolids = 0;
1102  bitSet solidCells(cellTableId_.size(), false);
1103  {
1104  Map<word> solidMap = cellTable_.solids();
1105 
1106  forAll(cellTableId_, cellI)
1107  {
1108  if (solidMap.found(cellTableId_[cellI]))
1109  {
1110  solidCells.set(cellI);
1111  ++nSolids;
1112  }
1113  }
1114  }
1115 
1116  if (!nSolids)
1117  {
1118  return;
1119  }
1120 
1121 
1122  // The corresponding Foam patch
1123  const label patchIndex = origBndId_.find(defaultBoundaryRegion);
1124  const label nPatchFaces = patchSizes_[patchIndex];
1125 
1126  labelList patchStarts(patchStartList(nInternalFaces_));
1127  label adjustPatch = 0;
1128  for (label i = 0; i < nPatchFaces; ++i)
1129  {
1130  label faceI = patchStarts[patchIndex] + i;
1131  label cellI = faceOwner_[faceI];
1132 
1133  if (solidCells.test(cellI))
1134  {
1135  ++adjustPatch;
1136  }
1137  }
1138 
1139  // No solid cells on the Default_Boundary_Region
1140  if (!adjustPatch)
1141  {
1142  return;
1143  }
1144 
1145 
1146  // Insert Default_Boundary_Solid immediately after Default_Boundary_Region
1147  // then we only need to adjust a single patch and can easily re-merge
1148  // later
1149  label nPatches = patchSizes_.size();
1150  patchStarts.setSize(nPatches+1, 0);
1151  patchSizes_.setSize(nPatches+1, 0);
1152  origBndId_.setSize(nPatches+1, 0);
1153 
1154  // make room for new entry
1155  for (label i = nPatches; i > patchIndex; --i)
1156  {
1157  patchStarts[i] = patchStarts[i-1];
1158  patchSizes_[i] = patchSizes_[i-1];
1159  origBndId_[i] = origBndId_[i-1];
1160  }
1161 
1162  // Adjust start and sizes
1163  patchSizes_[patchIndex] -= adjustPatch;
1164  patchSizes_[patchIndex+1] = adjustPatch;
1165  patchStarts[patchIndex+1] =
1166  patchStarts[patchIndex] + patchSizes_[patchIndex];
1167 
1168  origBndId_[patchIndex+1] = boundaryRegion_.append
1169  (
1170  dictionary
1171  (
1172  IStringStream
1173  (
1174  "BoundaryType wall;"
1175  "Label " + word(defaultSolidBoundaryName) + ";"
1176  )()
1177  )
1178  );
1179 
1180  label fluidFace = patchStarts[patchIndex];
1181  label solidFace = patchStarts[patchIndex+1];
1182 
1183  labelList oldToNew(identity(nFaces_));
1184  for (label i = 0; i < nPatchFaces; ++i)
1185  {
1186  label faceI = patchStarts[patchIndex] + i;
1187  label cellI = faceOwner_[faceI];
1188 
1189  if (solidCells.test(cellI))
1190  {
1191  oldToNew[faceI] = solidFace++;
1192  }
1193  else
1194  {
1195  oldToNew[faceI] = fluidFace++;
1196  }
1197  }
1198 
1199  // Re-order faces, owners/neighbours
1200  inplaceReorder(oldToNew, faces_);
1201  inplaceReorder(oldToNew, faceOwner_);
1202  inplaceReorder(oldToNew, faceNeighbour_);
1203  inplaceReorder(oldToNew, origFaceId_);
1204 
1205  renumberInterfaces(oldToNew);
1206 }
1207 
1208 
1209 // In CCM, separate mesh domains are used for fluid/porous/solid
1210 // Thus the fluid/porous/solid cells correspond uniquely to a face owner
1211 void Foam::ccm::reader::removeUnwanted()
1212 {
1213  // Identify fluid/porous/solid cells for removal
1214  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1215  label nRemove = 0;
1216  bitSet removeCells(cellTableId_.size(), false);
1217 
1218  {
1219  Map<word> fluidMap = cellTable_.fluids();
1220  Map<word> porousMap = selectPorous(cellTable_);
1221  Map<word> solidMap = cellTable_.solids();
1222  Map<word> removeMap;
1223 
1224  forAll(cellTableId_, cellI)
1225  {
1226  label tableId = cellTableId_[cellI];
1227 
1228  if
1229  (
1230  porousMap.found(tableId)
1231  ? !option().keepPorous()
1232  : fluidMap.found(tableId)
1233  ? !option().keepFluid()
1234  : solidMap.found(tableId)
1235  ? !option().keepSolid()
1236  : false
1237  )
1238  {
1239  removeCells.set(cellI);
1240  ++nRemove;
1241  removeMap.set(tableId, cellTable_.name(tableId));
1242  }
1243  }
1244 
1245  if (nRemove)
1246  {
1247  Map<word> keepMap;
1248 
1249  forAllConstIters(cellTable_, iter)
1250  {
1251  const label tableId = iter.key();
1252  if (!removeMap.found(tableId))
1253  {
1254  keepMap.set(tableId, cellTable_.name(tableId));
1255  }
1256  }
1257 
1258  Info<<"remove "<< nRemove << " cells in "
1259  << removeMap.size() << " unwanted cellZone(s)" << nl;
1260 
1261  forAllConstIters(removeMap, iter)
1262  {
1263  Info<< " zone "
1264  << iter.key() << " : " << iter.val() << nl;
1265  }
1266 
1267  Info<<"retain "<< (nCells_ - nRemove) << " cells in "
1268  << keepMap.size() << " cellZone(s)" << nl;
1269 
1270  forAllConstIters(keepMap, iter)
1271  {
1272  Info<< " zone "
1273  << iter.key() << " : " << iter.val() << nl;
1274  }
1275  }
1276  }
1277 
1278  if (!nRemove)
1279  {
1280  return;
1281  }
1282 
1283  // Remove all faces where the owner corresponds to a removed cell
1284  // Adjust the nInternalFaces and patch sizes accordingly
1285  label adjustInternal = 0;
1286  labelList adjustPatchSize(patchSizes_.size(), Zero);
1287 
1288  label newFaceI = 0;
1289  label oldFaceI = nFaces_ - 1;
1290  labelList oldToNew(nFaces_, -1);
1291  for (label faceI = 0; faceI < nFaces_; ++faceI)
1292  {
1293  label cellI = faceOwner_[faceI];
1294  if (removeCells.test(cellI))
1295  {
1296  if (faceI < nInternalFaces_)
1297  {
1298  ++adjustInternal;
1299  }
1300  else
1301  {
1302  // need to adjust the patch sizes
1303  label beg = nInternalFaces_;
1304  forAll(patchSizes_, patchI)
1305  {
1306  label end = beg + patchSizes_[patchI];
1307 
1308  if (faceI >= beg && faceI < end)
1309  {
1310  ++adjustPatchSize[patchI];
1311  break;
1312  }
1313 
1314  beg = end;
1315  }
1316  }
1317 
1318  // Put discarded faces at the end of the list
1319  oldToNew[faceI] = oldFaceI--;
1320  }
1321  else
1322  {
1323  if (newFaceI != faceI)
1324  {
1325  faces_[newFaceI] = faces_[faceI];
1326  faceOwner_[newFaceI] = faceOwner_[faceI];
1327  faceNeighbour_[newFaceI] = faceNeighbour_[faceI];
1328  origFaceId_[newFaceI] = origFaceId_[faceI];
1329  }
1330 
1331  // New position for the face
1332  oldToNew[faceI] = newFaceI++;
1333  }
1334  }
1335 
1336  nFaces_ = newFaceI;
1337 
1338  // Redimension
1339  faces_.setSize(nFaces_);
1340  faceOwner_.setSize(nFaces_);
1341  faceNeighbour_.setSize(nFaces_);
1342  origFaceId_.setSize(nFaces_);
1343 
1344  // Adjust internal faces and patch sizes
1345 
1346  nInternalFaces_ -= adjustInternal;
1347  forAll(patchSizes_, patchI)
1348  {
1349  patchSizes_[patchI] -= adjustPatchSize[patchI];
1350  }
1351 
1352  renumberInterfaces(oldToNew);
1353 
1354  // Renumber cells
1355  label nCell = 0;
1356  oldToNew.setSize(nCells_, -1);
1357  for (label cellI = 0; cellI < nCells_; ++cellI)
1358  {
1359  if (!removeCells.test(cellI))
1360  {
1361  if (nCell != cellI)
1362  {
1363  origCellId_[nCell] = origCellId_[cellI];
1364  cellTableId_[nCell] = cellTableId_[cellI];
1365  }
1366  oldToNew[cellI] = nCell;
1367  ++nCell;
1368  }
1369  }
1370 
1371  inplaceRenumber(oldToNew, faceOwner_);
1372  inplaceRenumber(oldToNew, faceNeighbour_);
1373 
1374  // Redimension
1375  nCells_ = nCell;
1376  origCellId_.setSize(nCells_);
1377  cellTableId_.setSize(nCells_);
1378 
1379 
1380  // Remove unused points - adjust points, faces accordingly
1381  oldToNew.setSize(nPoints_);
1382  oldToNew = -1;
1383 
1384  // Mark up all the used points
1385  for (const face& f : faces_)
1386  {
1387  for (const label pointi : f)
1388  {
1389  ++oldToNew[pointi];
1390  }
1391  }
1392 
1393  label nPointUsed = 0;
1394  forAll(oldToNew, ptI)
1395  {
1396  if (oldToNew[ptI] >= 0)
1397  {
1398  oldToNew[ptI] = nPointUsed;
1399  if (ptI != nPointUsed)
1400  {
1401  points_[nPointUsed] = points_[ptI];
1402  }
1403  ++nPointUsed;
1404  }
1405  }
1406 
1407  nPoints_ = nPointUsed;
1408  points_.setSize(nPoints_);
1409 
1410  for (face& f : faces_)
1411  {
1412  inplaceRenumber(oldToNew, f);
1413  }
1414 
1415  // Report sizes
1416  printSizes();
1417 }
1418 
1419 
1420 void Foam::ccm::reader::validateInterface
1421 (
1422  List<labelPair>& lst
1423 )
1424 {
1425  label nElem = 0;
1426  forAll(lst, elemI)
1427  {
1428  label face0 = lst[elemI][0];
1429  label face1 = lst[elemI][1];
1430 
1431  if (face0 < nFaces_ && face1 < nFaces_)
1432  {
1433  if (nElem != elemI)
1434  {
1435  lst[nElem][0] = face0;
1436  lst[nElem][1] = face1;
1437  }
1438  ++nElem;
1439  }
1440  }
1441  lst.setSize(nElem);
1442 }
1443 
1444 
1445 void Foam::ccm::reader::renumberInterfaces
1446 (
1447  const labelUList& oldToNew
1448 )
1449 {
1450  forAll(domInterfaces_, elemI)
1451  {
1452  domInterfaces_[elemI][0] = oldToNew[domInterfaces_[elemI][0]];
1453  domInterfaces_[elemI][1] = oldToNew[domInterfaces_[elemI][1]];
1454  }
1455 
1456  forAll(bafInterfaces_, elemI)
1457  {
1458  bafInterfaces_[elemI][0] = oldToNew[bafInterfaces_[elemI][0]];
1459  bafInterfaces_[elemI][1] = oldToNew[bafInterfaces_[elemI][1]];
1460  }
1461 
1462  validateInterface(domInterfaces_);
1463  validateInterface(bafInterfaces_);
1464 }
1465 
1466 
1467 //
1468 // 1) remove interfaces between domains (fluid/porosity; fluid/solid, etc)
1469 // 2) reorganize baffle interfaces into [0-N/2; N/2-N] lists at the beginning
1470 // of the corresponding patch
1471 //
1472 void Foam::ccm::reader::cleanupInterfaces()
1473 {
1474  validateInterface(bafInterfaces_);
1475  validateInterface(domInterfaces_);
1476 
1477  if (bafInterfaces_.size() <= 0 && domInterfaces_.size() <= 0)
1478  {
1479  Info<<"0 baffle interface pairs" << nl
1480  <<"0 domain interface pairs" << endl;
1481  return;
1482  }
1483 
1484 #ifdef DEBUG_BAFFLES
1485  Info<< "baffle Interfaces " << bafInterfaces_ << nl
1486  << "domain Interfaces " << domInterfaces_ << nl
1487  << "nCells:" << nCells_ << nl
1488  << "nFaces:" << nFaces_ << nl
1489  << "patchSizes:" << patchSizes_ << nl
1490  << "nInternalFaces:" << nInternalFaces_ << endl;
1491 
1492  forAll(domInterfaces_, elemI)
1493  {
1494  const label face0 = domInterfaces_[elemI][0];
1495  const label face1 = domInterfaces_[elemI][1];
1496 
1497  Info<< "interface [" << elemI << "] = "
1498  << face0 << " - " << face1 << " own/neigh = "
1499  << faceOwner_[face0] << "/" << faceNeighbour_[face0] << " "
1500  << faceOwner_[face1] << "/" << faceNeighbour_[face1] << endl;
1501  }
1502 #endif
1503 
1504  // Only reorder faces that need it
1505  labelList oldToNew(nFaces_, -1);
1506 
1507  // - move side0 face from domInterfaces to join the internal faces
1508  // - move the redundant side1 to the end of the list for later deletion
1509  label begOfList = nInternalFaces_;
1510  label endOfList = nFaces_ - 1;
1511 
1512  // The patch sizes (and the start) will definitely change
1513  const labelList origPatchStarts(patchStartList(nInternalFaces_));
1514  labelList adjustPatchSize(patchSizes_.size(), Zero);
1515  labelList bafflePatchCount(patchSizes_.size(), Zero);
1516 
1517  // The new dimensions after merging the domain interfaces:
1518  nInternalFaces_ += domInterfaces_.size();
1519  nFaces_ -= domInterfaces_.size();
1520 
1521  Info<< domInterfaces_.size() << " domain interface pairs";
1522  if (domInterfaces_.size())
1523  {
1524  Info<<" to merge" << endl;
1525  printSizes();
1526  }
1527  else
1528  {
1529  Info<< endl;
1530  }
1531 
1532  forAll(domInterfaces_, elemI)
1533  {
1534  label face0 = domInterfaces_[elemI][0];
1535  label face1 = domInterfaces_[elemI][1];
1536 
1537  oldToNew[face0] = begOfList++;
1538  oldToNew[face1] = endOfList--; // End of list for truncationx
1539 
1540  // face0 gets a new neighbour, face1 loses its owner
1541  faceNeighbour_[face0] = faceOwner_[face1];
1542  faceOwner_[face1] = -1;
1543 
1544  // Need to adjust the patch sizes
1545  forAll(patchSizes_, patchI)
1546  {
1547  label beg = origPatchStarts[patchI];
1548  label end = beg + patchSizes_[patchI];
1549 
1550  if (face0 >= beg && face0 < end)
1551  {
1552  ++adjustPatchSize[patchI];
1553  }
1554  if (face1 >= beg && face1 < end)
1555  {
1556  ++adjustPatchSize[patchI];
1557  }
1558  }
1559  }
1560 
1561  // Count the number of baffles per patch
1562  forAll(bafInterfaces_, elemI)
1563  {
1564  label face0 = bafInterfaces_[elemI][0];
1565  label face1 = bafInterfaces_[elemI][1];
1566 
1567  forAll(patchSizes_, patchI)
1568  {
1569  label beg = origPatchStarts[patchI];
1570  label end = beg + patchSizes_[patchI];
1571 
1572  if (face0 >= beg && face0 < end)
1573  {
1574  ++bafflePatchCount[patchI];
1575  }
1576  if (face1 >= beg && face1 < end)
1577  {
1578  ++bafflePatchCount[patchI];
1579  }
1580  }
1581  }
1582 
1583 
1584  if (option().removeBaffles())
1585  {
1586  // The new dimensions after merging the baffles:
1587  nInternalFaces_ += bafInterfaces_.size();
1588  nFaces_ -= bafInterfaces_.size();
1589  }
1590 
1591  Info<< bafInterfaces_.size() << " baffle interface pairs";
1592  if (bafInterfaces_.size())
1593  {
1594  if (option().removeBaffles())
1595  {
1596  Info<< " to merge" << endl;
1597  printSizes();
1598  }
1599  else
1600  {
1601  Info<< " to be sorted" << endl;
1602  }
1603  }
1604  else
1605  {
1606  Info<< endl;
1607  }
1608 
1609  if (option().removeBaffles())
1610  {
1611  forAll(bafInterfaces_, elemI)
1612  {
1613  label face0 = bafInterfaces_[elemI][0];
1614  label face1 = bafInterfaces_[elemI][1];
1615 
1616  oldToNew[face0] = begOfList++;
1617  oldToNew[face1] = endOfList--; // End of list for truncation
1618 
1619  // face0 gets a new neighbour, face1 loses its owner
1620  faceNeighbour_[face0] = faceOwner_[face1];
1621  faceOwner_[face1] = -1;
1622  }
1623 
1624 
1625  // Boundary faces continue from the new nInternalFaces
1626  label pos = nInternalFaces_;
1627  forAll(patchSizes_, patchI)
1628  {
1629  label beg = origPatchStarts[patchI];
1630  label end = beg + patchSizes_[patchI];
1631 
1632  // Fill in values for the remainder of the boundary patch
1633  for (label faceI = beg; faceI < end; ++faceI)
1634  {
1635  if (oldToNew[faceI] < 0)
1636  {
1637  oldToNew[faceI] = pos;
1638  ++pos;
1639  }
1640  }
1641  }
1642 
1643  // Baffles have been resolved - remove last traces
1644  bafInterfaces_.clear();
1645  }
1646  else
1647  {
1648  // This check is probably unnecessary
1649  forAll(bafflePatchCount, patchI)
1650  {
1651  if (bafflePatchCount[patchI] % 2)
1652  {
1653  Info<< "WARNING: patch " << patchI
1654  << " has an uneven number of baffles ("
1655  << bafflePatchCount[patchI] << ") expect strange results"
1656  << endl;
1657  }
1658  }
1659 
1660 
1661  // Reordered faces continue from the new nInternalFaces
1662  label pos = nInternalFaces_;
1663  forAll(patchSizes_, patchI)
1664  {
1665  const label beg = origPatchStarts[patchI];
1666  const label end = beg + patchSizes_[patchI];
1667 
1668  const label nsize = bafflePatchCount[patchI];
1669  if (nsize > 0)
1670  {
1671  // Reorganize baffle interfaces into [0-N/2; N/2-N] lists
1672  // at the beginning of the corresponding patch
1673  const label nsizeby2 = (nsize - nsize % 2) / 2;
1674  label nsorted = 0;
1675 
1676  // Renumber the normal (baffle) interfaces
1677  forAll(bafInterfaces_, elemI)
1678  {
1679  const label face0 = bafInterfaces_[elemI][0];
1680  const label face1 = bafInterfaces_[elemI][1];
1681 
1682  if
1683  (
1684  (face0 >= beg && face0 < end)
1685  || (face1 >= beg && face1 < end)
1686  )
1687  {
1688  oldToNew[face0] = pos + nsorted;
1689  oldToNew[face1] = pos + nsorted + nsizeby2;
1690 
1691  // Mark destination of the faces, but cannot renumber
1692  // yet. Use negative to potential overlap with other
1693  // patch regions
1694  bafInterfaces_[elemI][0] = -oldToNew[face0];
1695  bafInterfaces_[elemI][1] = -oldToNew[face1];
1696 
1697  ++nsorted;
1698  }
1699  }
1700  pos += 2*nsorted;
1701  }
1702 
1703  // Fill in values for the remainder of the boundary patch
1704  for (label faceI = beg; faceI < end; ++faceI)
1705  {
1706  if (oldToNew[faceI] < 0)
1707  {
1708  oldToNew[faceI] = pos;
1709  ++pos;
1710  }
1711  }
1712  }
1713 
1714  // Finalize new numbers for the normal (baffle) interfaces
1715  forAll(bafInterfaces_, elemI)
1716  {
1717  bafInterfaces_[elemI][0] = abs(bafInterfaces_[elemI][0]);
1718  bafInterfaces_[elemI][1] = abs(bafInterfaces_[elemI][1]);
1719  }
1720  }
1721 
1722 #ifdef DEBUG_BAFFLES
1723  Info<< "remap with " << oldToNew << nl
1724  << "owners:" << faceOwner_ << nl
1725  << "neighbours:" << faceNeighbour_ << nl
1726  << endl;
1727 #endif
1728 
1729  // Re-order faces, owners/neighbours
1730  inplaceReorder(oldToNew, faces_);
1731  inplaceReorder(oldToNew, faceOwner_);
1732  inplaceReorder(oldToNew, faceNeighbour_);
1733  inplaceReorder(oldToNew, origFaceId_);
1734 
1735  if (monitoringSets_.size())
1736  {
1737 #ifdef WITH_MONITORING
1738  // Modify oldToNew mapping to account for monitoring faces that
1739  // coincided with a domain interface
1740  //
1741  // TODO - should modify flip map as well
1742  forAll(domInterfaces_, elemI)
1743  {
1744  label face0 = domInterfaces_[elemI][0];
1745  label face1 = domInterfaces_[elemI][1];
1746  oldToNew[face1] = oldToNew[face0];
1747  }
1748 
1749  forAllIters(monitoringSets_, iter)
1750  {
1751  inplaceRenumber(oldToNew, iter.val());
1752  }
1753 #endif
1754  }
1755 
1756  // We can finally drop this information now
1757  domInterfaces_.clear();
1758 
1759  // Truncate lists
1760  faces_.setSize(nFaces_);
1761  faceOwner_.setSize(nFaces_);
1762  faceNeighbour_.setSize(nFaces_);
1763  origFaceId_.setSize(nFaces_);
1764 
1765  // Remove empty patches:
1766 
1767  // Fix patch sizes:
1768  oldToNew.setSize(patchSizes_.size());
1769  oldToNew = -1;
1770 
1771  label nPatches = 0;
1772  forAll(patchSizes_, patchI)
1773  {
1774  patchSizes_[patchI] -= adjustPatchSize[patchI];
1775  if (option().removeBaffles())
1776  {
1777  patchSizes_[patchI] -= bafflePatchCount[patchI];
1778  }
1779 
1780  if (patchSizes_[patchI])
1781  {
1782  oldToNew[patchI] = nPatches++;
1783  }
1784  }
1785 
1786  inplaceReorder(oldToNew, patchSizes_);
1787  inplaceReorder(oldToNew, origBndId_);
1788 
1789  patchSizes_.setSize(nPatches);
1790  origBndId_.setSize(nPatches);
1791 
1792 #ifdef DEBUG_BAFFLES
1793  Info<< "nCells:" << nCells_ << nl
1794  << "nFaces:" << nFaces_ << nl
1795  << "PatchSizes:" << patchSizes_ << nl
1796  << "nInternalFaces:" << nInternalFaces_ << nl
1797  << endl;
1798 #endif
1799 }
1800 
1801 
1802 //
1803 // Merge STARCCM in-place interfaces
1804 //
1805 void Foam::ccm::reader::mergeInplaceInterfaces()
1806 {
1807  if (interfaceDefinitions_.empty())
1808  {
1809  return;
1810  }
1811  if (!option().mergeInterfaces())
1812  {
1813  Info<< interfaceDefinitions_.size() << " interface definitions"
1814  << " - leaving unmerged" << endl;
1815  return;
1816  }
1817 
1818  // List of patch pairs that are interfaces
1819  DynamicList<labelPair> interfacePatches(interfaceDefinitions_.size());
1820 
1821  label nWarn = 0;
1822 
1823  forAllConstIters(interfaceDefinitions_, iter)
1824  {
1825  const interfaceEntry& ifentry = iter.val();
1826 
1827  labelPair patchPair
1828  (
1829  origBndId_.find(ifentry.bnd0),
1830  origBndId_.find(ifentry.bnd1)
1831  );
1832 
1833  if
1834  (
1835  patchPair[0] == patchPair[1]
1836  || patchPair[0] < 0
1837  || patchPair[1] < 0
1838  )
1839  {
1840  // This should not happen
1841  Info<<"Warning : bad interface " << ifentry.id << " " << ifentry
1842  <<" on patches " << patchPair << endl;
1843  }
1844  else if
1845  (
1846  patchSizes_[patchPair[0]] != patchSizes_[patchPair[1]]
1847  || patchSizes_[patchPair[0]] == 0
1848  || patchSizes_[patchPair[1]] == 0
1849  )
1850  {
1851  if (!nWarn++)
1852  {
1853  Info<<"Warning: skip interface with zero or different"
1854  << " number of faces" << nl;
1855  }
1856 
1857  Info<<" Interface:" << ifentry.id << " " << ifentry
1858  <<" patches " << patchPair
1859  <<" sizes ("
1860  << patchSizes_[patchPair[0]]
1861  << " " << patchSizes_[patchPair[1]] << ")"
1862  << nl;
1863  }
1864  else
1865  {
1866  interfacePatches.append(patchPair);
1867  }
1868  }
1869 
1870  if (interfacePatches.empty())
1871  {
1872  return;
1873  }
1874 
1875 
1876  // Local point mapping
1877  labelList mergedPointMap;
1878 
1879  // Global remapping
1880  labelList oldToNew(identity(points_.size()));
1881 
1882  const labelList origPatchStarts(patchStartList(nInternalFaces_));
1883 
1884  label nMergedTotal = 0;
1885 
1886  // Markup points to merge
1887  bitSet whichPoints(points_.size());
1888 
1889  Info<< "interface merge points (tol="
1890  << option().mergeTol() << "):" << endl;
1891 
1892  DynamicList<label> interfacesToMerge(interfacePatches.size());
1893  forAll(interfacePatches, interI)
1894  {
1895  const label patch0 = interfacePatches[interI][0];
1896  const label patch1 = interfacePatches[interI][1];
1897  const label nPatch0Faces = patchSizes_[patch0];
1898  const label nPatch1Faces = patchSizes_[patch1];
1899 
1900  // Markup points to merge
1901  whichPoints.reset();
1902  for (label local0FaceI = 0; local0FaceI < nPatch0Faces; ++local0FaceI)
1903  {
1904  const face& f = faces_[origPatchStarts[patch0] + local0FaceI];
1905 
1906  for (const label pointi : f)
1907  {
1908  // Simultaneously account for previous point merges
1909  whichPoints.set(oldToNew[pointi]);
1910  }
1911  }
1912  for (label local1FaceI = 0; local1FaceI < nPatch1Faces; ++local1FaceI)
1913  {
1914  const face& f = faces_[origPatchStarts[patch1] + local1FaceI];
1915 
1916  for (const label pointi : f)
1917  {
1918  // Simultaneously account for previous point merges
1919  whichPoints.set(oldToNew[pointi]);
1920  }
1921  }
1922 
1923  // The global addresses
1924  labelList addr(whichPoints.toc());
1925 
1926  const UIndirectList<point> pointsToMerge(points_, addr);
1927 
1928  Info<< " patch " << patch0 << ',' << patch1 << ": ("
1929  << nPatch0Faces << " and " << nPatch1Faces << " faces) " << flush;
1930 
1931  const label nMerged = mergePoints
1932  (
1933  pointsToMerge,
1934  option().mergeTol(),
1935  false,
1936  mergedPointMap
1937  );
1938 
1939  Info<< nMerged << " from " << pointsToMerge.size() << " points"
1940  << endl;
1941 
1942  if (nMerged)
1943  {
1944  // Two-steps:
1945  // * Identify duplicate points - without changing the order!
1946  // * Transcribe local to global addressing (oldToNew)
1947 
1948  forAll(mergedPointMap, pti)
1949  {
1950  const label mergedPti = mergedPointMap[pti];
1951 
1952  if (mergedPti < 0)
1953  {
1954  continue; // Already seen as duplicate
1955  }
1956 
1957  const label origPointi = oldToNew[addr[pti]];
1958 
1959  // Find further duplicate points
1960  for
1961  (
1962  label dupPti = pti+1;
1963  (dupPti = mergedPointMap.find(mergedPti, dupPti)) != -1;
1964  ++dupPti
1965  )
1966  {
1967  oldToNew[addr[dupPti]] = origPointi;
1968  mergedPointMap[dupPti] = -1; // Mark as already seen
1969  }
1970  }
1971 
1972  interfacesToMerge.append(interI);
1973  nMergedTotal += nMerged;
1974  }
1975  }
1976 
1977 
1978  //
1979  // Nothing to do
1980  //
1981  if (!nMergedTotal)
1982  {
1983  return;
1984  }
1985 
1986  // Update point references to account for point merge:
1987  for (face& f : faces_)
1988  {
1989  inplaceRenumber(oldToNew, f);
1990  }
1991 
1992  // Determine which points are actually in use:
1993  oldToNew.resize_nocopy(nPoints_);
1994  oldToNew = -1;
1995 
1996  // Mark up all the used points
1997  for (const face& f : faces_)
1998  {
1999  for (const label pointi : f)
2000  {
2001  ++oldToNew[pointi];
2002  }
2003  }
2004 
2005  label nPointUsed = 0;
2006  forAll(oldToNew, ptI)
2007  {
2008  if (oldToNew[ptI] >= 0)
2009  {
2010  oldToNew[ptI] = nPointUsed;
2011  if (ptI != nPointUsed)
2012  {
2013  points_[nPointUsed] = points_[ptI];
2014  }
2015  ++nPointUsed;
2016  }
2017  }
2018 
2019  // Info<< "merge " << nMergedTotal << " points from "
2020  // << nPoints_ << " to " << nPointUsed << endl;
2021 
2022  nPoints_ = nPointUsed;
2023  points_.resize(nPoints_);
2024 
2025  for (face& f : faces_)
2026  {
2027  inplaceRenumber(oldToNew, f);
2028  }
2029 
2030 
2031  //
2032  // Merge the faces as well
2033  //
2034  Info<< "interface merge faces:" << endl;
2035 
2036  nMergedTotal = 0;
2037  labelList adjustPatchSize(patchSizes_.size(), Zero);
2038  forAll(interfacesToMerge, mergeI)
2039  {
2040  const label patch0 = interfacePatches[interfacesToMerge[mergeI]][0];
2041  const label patch1 = interfacePatches[interfacesToMerge[mergeI]][1];
2042 
2043  labelList faceAddr0(patchSizes_[patch0]);
2044  labelList faceAddr1(patchSizes_[patch1]);
2045 
2046  forAll(faceAddr0, localFaceI)
2047  {
2048  faceAddr0[localFaceI] = origPatchStarts[patch0] + localFaceI;
2049  }
2050  forAll(faceAddr1, localFaceI)
2051  {
2052  faceAddr1[localFaceI] = origPatchStarts[patch1] + localFaceI;
2053  }
2054 
2055  if (faceAddr0.size() != faceAddr1.size())
2056  {
2057  // This should not occur, we avoided the same thing above
2058  continue;
2059  }
2060 
2061  // Improve comparison speed by sorting by distance
2062  SortableList<scalar> pts0MagSqr
2063  (
2064  magSqr
2065  (
2067  (
2068  UIndirectList<face>
2069  (
2070  faces_,
2071  faceAddr0
2072  ),
2073  points_
2074  ).faceCentres()
2075  )
2076  );
2077  SortableList<scalar> pts1MagSqr
2078  (
2079  magSqr
2080  (
2082  (
2083  UIndirectList<face>
2084  (
2085  faces_,
2086  faceAddr1
2087  ),
2088  points_
2089  ).faceCentres()
2090  )
2091  );
2092 
2093  label nMerged = 0;
2094 
2095  // Record which faces failed to merge - use slower ad hoc merging
2096  labelHashSet failed0, failed1;
2097  forAll(pts0MagSqr, sortI)
2098  {
2099  const label face0I = faceAddr0[pts0MagSqr.indices()[sortI]];
2100  const label face1I = faceAddr1[pts1MagSqr.indices()[sortI]];
2101 
2102  // This is what we expect
2103  if (face::compare(faces_[face0I], faces_[face1I]))
2104  {
2105  ++nMerged;
2106 
2107  // Moved from boundary patch to internal patch
2108  ++adjustPatchSize[patch0];
2109  ++adjustPatchSize[patch1];
2110 
2111  if (faceOwner_[face0I] < faceOwner_[face1I])
2112  {
2113  // keep 0, discard 1
2114  faceNeighbour_[face0I] = faceOwner_[face1I];
2115  faceNeighbour_[face1I] = faceOwner_[face1I] = -1;
2116  }
2117  else
2118  {
2119  // keep 1, discard 0
2120  faceNeighbour_[face1I] = faceOwner_[face0I];
2121  faceNeighbour_[face0I] = faceOwner_[face0I] = -1;
2122  }
2123  }
2124  else
2125  {
2126  failed0.set(face0I);
2127  failed1.set(face1I);
2128  }
2129  }
2130 
2131  // Perhaps some sorting issues, recheck
2132  if (failed0.size())
2133  {
2134  // Note which one were successful
2135  labelHashSet done(failed0.size());
2136 
2137  for (const label face0I : failed0)
2138  {
2139  for (const label face1I : failed1)
2140  {
2141  // This is what we expect
2142  if (face::compare(faces_[face0I], faces_[face1I]))
2143  {
2144  ++nMerged;
2145 
2146  // Moved from boundary patch to internal patch
2147  ++adjustPatchSize[patch0];
2148  ++adjustPatchSize[patch1];
2149 
2150  if (faceOwner_[face0I] < faceOwner_[face1I])
2151  {
2152  // keep 0, discard 1
2153  faceNeighbour_[face0I] = faceOwner_[face1I];
2154  faceNeighbour_[face1I] = faceOwner_[face1I] = -1;
2155  }
2156  else
2157  {
2158  // keep 1, discard 0
2159  faceNeighbour_[face1I] = faceOwner_[face0I];
2160  faceNeighbour_[face0I] = faceOwner_[face0I] = -1;
2161  }
2162 
2163  failed1.erase(face1I); // Never check again
2164  done.set(face0I); // Mark as done
2165  break; // Stop looking
2166  }
2167  }
2168  }
2169 
2170  // Transfer to note how many were successful
2171  failed0 = done;
2172  }
2173 
2174  Info<< " patch " << patch0 << ',' << patch1 << ": "
2175  << nMerged << " from " << faceAddr0.size() << " faces";
2176 
2177  if (failed0.size())
2178  {
2179  Info<< " (" << failed0.size() << " merged ad hoc)";
2180  }
2181  Info<< endl;
2182 
2183 
2184  nMergedTotal += nMerged;
2185  }
2186 
2187 
2188  // Nothing to do
2189  if (!nMergedTotal)
2190  {
2191  return;
2192  }
2193 
2194  // Info<< "merge " << nMergedTotal << " faces from "
2195  // << nFaces_ << " to " << (nFaces_ - nMergedTotal) << endl;
2196 
2197  oldToNew.setSize(nFaces_);
2198  oldToNew = -1;
2199 
2200  // Remaining external faces will be shifted here:
2201  label extFaceI = nInternalFaces_ + nMergedTotal;
2202 
2203  // Start over
2204  nInternalFaces_ = 0;
2205  label nFaceUsed = 0;
2206  for (label faceI = 0; faceI < nFaces_; ++faceI)
2207  {
2208  if (faceOwner_[faceI] != -1)
2209  {
2210  if (faceNeighbour_[faceI] != -1)
2211  {
2212  // Internal face
2213  oldToNew[faceI] = nInternalFaces_;
2214  ++nInternalFaces_;
2215  ++nFaceUsed;
2216  }
2217  else
2218  {
2219  // External face
2220  oldToNew[faceI] = extFaceI;
2221  ++extFaceI;
2222  ++nFaceUsed;
2223  }
2224  }
2225  }
2226 
2227  if (nFaceUsed != extFaceI)
2228  {
2230  << "coding error: used " << nFaceUsed
2231  << " faces, but expected to use " << extFaceI << " faces"
2232  << exit(FatalError);
2233  }
2234 
2235  // Re-order faces, owners/neighbours
2236  inplaceReorder(oldToNew, faces_, true);
2237  inplaceReorder(oldToNew, faceOwner_, true);
2238  inplaceReorder(oldToNew, faceNeighbour_, true);
2239  inplaceReorder(oldToNew, origFaceId_, true);
2240 
2241  nFaces_ = nFaceUsed;
2242 
2243  faces_.setSize(nFaces_);
2244  faceOwner_.setSize(nFaces_);
2245  faceNeighbour_.setSize(nFaces_);
2246  origFaceId_.setSize(nFaces_);
2247 
2248  // Fix patch sizes:
2249  oldToNew.setSize(patchSizes_.size());
2250  oldToNew = -1;
2251 
2252  label nPatches = 0;
2253  forAll(patchSizes_, patchI)
2254  {
2255  patchSizes_[patchI] -= adjustPatchSize[patchI];
2256  if (patchSizes_[patchI])
2257  {
2258  oldToNew[patchI] = nPatches++;
2259  }
2260  }
2261 
2262  inplaceReorder(oldToNew, patchSizes_);
2263  inplaceReorder(oldToNew, origBndId_);
2264 
2265  patchSizes_.setSize(nPatches);
2266  origBndId_.setSize(nPatches);
2267 
2268  // Report we are done
2269  Info<< ".." << endl;
2270 }
2271 
2272 
2273 //
2274 // Re-order mesh into Foam convention
2275 // - owner < neighbour
2276 // - face vertices such that normal points away from owner
2277 // - order faces: upper-triangular for internal faces;
2278 // boundary faces after internal faces
2279 //
2280 void Foam::ccm::reader::reorderMesh()
2281 {
2282  // Set owner/neighbour so owner < neighbour
2283  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2284 
2285  forAll(faceOwner_, faceI)
2286  {
2287  const label nbr = faceNeighbour_[faceI];
2288  const label own = faceOwner_[faceI];
2289 
2290  if (nbr >= cellTableId_.size() || own >= cellTableId_.size())
2291  {
2293  << "face:" << faceI
2294  << " nbr:" << nbr
2295  << " own:" << own
2296  << " nCells:" << cellTableId_.size()
2297  << exit(FatalError);
2298  }
2299 
2300  if (nbr >= 0 && nbr < own)
2301  {
2302  faceOwner_[faceI] = faceNeighbour_[faceI];
2303  faceNeighbour_[faceI] = own;
2304  faces_[faceI].flip();
2305  }
2306 
2307  // And check the face
2308  const face& f = faces_[faceI];
2309 
2310  for (const label pointi : f)
2311  {
2312  if (pointi < 0 || pointi >= points_.size())
2313  {
2315  << "face:" << faceI << " f:" << f
2316  << abort(FatalError);
2317  }
2318  }
2319  }
2320 
2321  // Do upper-triangular ordering
2322  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2323  labelList oldToNew(faceOwner_.size(), -1);
2324 
2325  // Create cells (inverse of face-to-cell addressing)
2326  cellList cellFaceAddr;
2327  primitiveMesh::calcCells
2328  (
2329  cellFaceAddr,
2330  faceOwner_,
2331  faceNeighbour_,
2332  nCells_
2333  );
2334 
2335  label newFaceI = 0;
2336  forAll(cellFaceAddr, cellI)
2337  {
2338  const labelList& cFaces = cellFaceAddr[cellI];
2339  SortableList<label> nbr(cFaces.size(), -1);
2340 
2341  forAll(cFaces, i)
2342  {
2343  label faceI = cFaces[i];
2344  label nbrCellI = faceNeighbour_[faceI];
2345 
2346  if (nbrCellI >= 0)
2347  {
2348  // Internal face. Get cell on other side
2349  if (nbrCellI == cellI)
2350  {
2351  nbrCellI = faceOwner_[faceI];
2352  }
2353 
2354  // cellI is master
2355  if (cellI < nbrCellI)
2356  {
2357  nbr[i] = nbrCellI;
2358  }
2359  }
2360  }
2361 
2362  nbr.sort();
2363 
2364  forAll(nbr, i)
2365  {
2366  if (nbr[i] >= 0)
2367  {
2368  oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
2369  }
2370  }
2371  }
2372 
2373  // Reorder faces accordingly
2374  inplaceReorder(oldToNew, faces_);
2375  inplaceReorder(oldToNew, faceOwner_);
2376  inplaceReorder(oldToNew, faceNeighbour_);
2377  inplaceReorder(oldToNew, origFaceId_);
2378 
2379  forAllIters(monitoringSets_, iter)
2380  {
2381  labelList& lst = iter.val();
2382  inplaceRenumber(oldToNew, lst);
2383 
2384  // disallow monitoring on boundaries
2385  label nElem = 0;
2386  forAll(lst, i)
2387  {
2388  if (lst[i] >= 0 && lst[i] < nInternalFaces_)
2389  {
2390  if (nElem != i)
2391  {
2392  lst[nElem] = lst[i];
2393  }
2394  ++nElem;
2395  }
2396  }
2397 
2398  if (nElem)
2399  {
2400  lst.setSize(nElem);
2401  }
2402  else
2403  {
2404  Info << "remove monitor " << iter.key() << endl;
2405  monitoringSets_.erase(iter);
2406  }
2407  }
2408 }
2409 
2410 
2411 // Attach patches
2412 //
2413 // - patchSizes_ : obvious
2414 // - origBndId_ : lookup for patch name/type
2415 //
2416 void Foam::ccm::reader::addPatches
2417 (
2418  polyMesh& mesh
2419 ) const
2420 {
2421  // Create patches
2422  // use patch types to determine what Foam types to generate
2423  List<polyPatch*> newPatches(origBndId_.size());
2424 
2425  label meshFaceI = nInternalFaces_;
2426  wordHashSet hashedNames(origBndId_.size());
2427 
2428  // lookup patch names/types from the problem description
2429  // provide some fallback values
2430  forAll(newPatches, patchI)
2431  {
2432  const word fallbackName(polyPatch::defaultName(patchI));
2433  word patchName;
2434  word patchType;
2435 
2436  auto citer = boundaryRegion_.cfind(origBndId_[patchI]);
2437 
2438  if (citer.found())
2439  {
2440  citer().readEntry("Label", patchName);
2441  citer().readEntry("BoundaryType", patchType);
2442  }
2443  else
2444  {
2445  patchName = fallbackName;
2446  patchType = "patch";
2447  }
2448 
2449  // Avoid duplicate names
2450  // - don't bother checking if the modified name is also a duplicate
2451  if (hashedNames.found(patchName))
2452  {
2453  Info<< "renamed patch " << patchName << " to ";
2454  patchName = fallbackName + "_" + patchName;
2455  Info<< patchName << endl;
2456  }
2457  hashedNames.insert(patchName);
2458 
2459  Info<< "patch " << patchI
2460  << " (start: " << meshFaceI << " size: " << patchSizes_[patchI]
2461  << ") name: " << patchName
2462  << endl;
2463 
2464  if (patchType == "wall")
2465  {
2466  newPatches[patchI] =
2467  new wallPolyPatch
2468  (
2469  patchName,
2470  patchSizes_[patchI],
2471  meshFaceI,
2472  patchI,
2473  mesh.boundaryMesh(),
2474  patchType
2475  );
2476  }
2477  else if (patchType == "symmetry")
2478  {
2479  newPatches[patchI] =
2480  new symmetryPolyPatch
2481  (
2482  patchName,
2483  patchSizes_[patchI],
2484  meshFaceI,
2485  patchI,
2486  mesh.boundaryMesh(),
2487  patchType
2488  );
2489  }
2490  else if (patchType == "empty")
2491  {
2492  // Note: not ccm name, may have been introduced by us
2493  newPatches[patchI] =
2494  new emptyPolyPatch
2495  (
2496  patchName,
2497  patchSizes_[patchI],
2498  meshFaceI,
2499  patchI,
2500  mesh.boundaryMesh(),
2501  patchType
2502  );
2503  }
2504  else
2505  {
2506  // All other ccm types become straight polyPatch:
2507  // 'inlet', 'outlet', 'pressure'.
2508  newPatches[patchI] =
2509  new polyPatch
2510  (
2511  patchName,
2512  patchSizes_[patchI],
2513  meshFaceI,
2514  patchI,
2515  mesh.boundaryMesh(),
2516  patchType
2517  );
2518  }
2519 
2520  meshFaceI += patchSizes_[patchI];
2521  }
2522 
2523  if (meshFaceI != mesh.nFaces())
2524  {
2526  << "meshFaceI:" << meshFaceI << " nFaces:" << mesh.nFaces()
2527  << abort(FatalError);
2528  }
2529 
2530  mesh.addPatches(newPatches);
2531 }
2532 
2533 
2534 
2535 // Attach faceZones based on the monitoring boundary conditions
2536 void Foam::ccm::reader::addFaceZones
2537 (
2538  polyMesh& mesh
2539 ) const
2540 {
2541  label nZone = monitoringSets_.size();
2542  mesh.faceZones().setSize(nZone);
2543 
2544  if (!nZone)
2545  {
2546  return;
2547  }
2548 
2549  nZone = 0;
2550  forAllConstIters(monitoringSets_, iter)
2551  {
2552  Info<< "faceZone " << nZone
2553  << " (size: " << iter().size() << ") name: "
2554  << iter.key() << endl;
2555 
2556  mesh.faceZones().set
2557  (
2558  nZone,
2559  new faceZone
2560  (
2561  iter.key(),
2562  iter(),
2563  false, // none are flipped
2564  nZone,
2565  mesh.faceZones()
2566  )
2567  );
2568 
2569  ++nZone;
2570  }
2571 
2572  mesh.faceZones().writeOpt(IOobject::AUTO_WRITE);
2573  warnDuplicates("faceZones", mesh.faceZones().names());
2574 }
2575 
2576 
2577 // Remove most of the ccm-specific information with the exception of information
2578 // auxiliary to the normal polyMesh:
2579 // (cellTableId_, origCellId_, origFaceId_, interfaces_, ...)
2581 {
2582  // allow re-reading the file
2583  if (geometryStatus_ == OKAY || geometryStatus_ == READ)
2584  {
2585  geometryStatus_ = UNKNOWN;
2586  }
2587 
2588  points_.clear();
2589  faces_.clear();
2590  faceOwner_.clear();
2591  faceNeighbour_.clear();
2592 
2593  origBndId_.clear();
2594  patchSizes_.clear();
2595 }
2596 
2597 
2598 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2599 
2602  const objectRegistry& registry,
2603  const fileName& remappingDictName
2604 )
2605 {
2606  if (!readGeometry())
2607  {
2608  return nullptr;
2609  }
2610 
2611  // merge cellTable and rename boundaryRegion
2612  remapMeshInfo(registry, remappingDictName);
2613 
2614  // Construct polyMesh
2615  // ~~~~~~~~~~~~~~~~~~
2617  (
2618  IOobject
2619  (
2621  "constant",
2622  registry,
2625  ),
2626  std::move(points_),
2627  std::move(faces_),
2628  std::move(faceOwner_),
2629  std::move(faceNeighbour_)
2630  );
2631  polyMesh& mesh = *meshPtr;
2632 
2633  addPatches(mesh);
2634 
2635  // Attach cellZones based on the cellTable Id
2636  // any other values can be extracted later from the cellTable dictionary
2637  cellTable_.addCellZones(mesh, cellTableId_);
2638  warnDuplicates("cellZones", mesh.cellZones().names());
2639 
2640  addFaceZones(mesh);
2641 
2642  warnDuplicates("boundaries", mesh.boundaryMesh().names());
2643  clearGeom();
2644 
2645  return meshPtr;
2646 }
2647 
2648 
2649 // ************************************************************************* //
Foam::autoPtr::New
static autoPtr< T > New(Args &&... args)
Construct autoPtr of T with forwarding arguments.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
meshPtr
Foam::autoPtr< Foam::fvMesh > meshPtr(nullptr)
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:194
Foam::fileName
A class for handling file names.
Definition: fileName.H:73
uindirectPrimitivePatch.H
Foam::polyMesh::defaultRegion
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:318
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
wallPolyPatch.H
symmetryPolyPatch.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
bitSet.H
Foam::ccm::reader::clearGeom
void clearGeom()
Clear out some information after obtaining a polyMesh.
Definition: ccmReaderMesh.C:2580
Foam::invert
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:61
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::face::compare
static int compare(const face &a, const face &b)
Compare faces.
Definition: face.C:281
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:54
Foam::objectRegistry
Registry of regIOobjects.
Definition: objectRegistry.H:60
Foam::flush
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:361
SortableList.H
nPatches
const label nPatches
Definition: printMeshSummary.H:30
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
ccmInternal.H
Internal bits for wrapping libccmio - do not use directly.
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
forAllIters
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:223
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
emptyPolyPatch.H
Foam::mergePoints
label mergePoints(const PointList &points, const scalar mergeTol, const bool verbose, labelList &pointMap, typename PointList::const_reference origin=PointList::value_type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
IOdictionary.H
Foam::autoPtr< Foam::polyMesh >
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::nl
constexpr char nl
Definition: Ostream.H:404
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
Fstream.H
f
labelList f(nPoints)
Foam::List< label >
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
mergePoints.H
Merge points. See below.
Foam::wordHashSet
HashSet< word, Hash< word > > wordHashSet
A HashSet of words, uses string hasher.
Definition: HashSet.H:77
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
ccmBoundaryInfo.H
Container for holding STARCCM boundary information.
Foam::inplaceReorder
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
Definition: ListOpsTemplates.C:124
ccmReader.H
patchId
label patchId(-1)
ListOps.H
Various functions to operate on Lists.
Foam::uindirectPrimitivePatch
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
Definition: uindirectPrimitivePatch.H:49
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Foam::patchIdentifier::defaultName
static word defaultName(const label n=-1)
Default patch name: "patch" or "patchN".
Definition: patchIdentifier.H:76
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
Foam::ccm::reader::mesh
autoPtr< polyMesh > mesh(const objectRegistry &registry, const fileName &remappingDictName=fileName::null)
Construct the polyMesh from the read geometry.
Definition: ccmReaderMesh.C:2601
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::HashSet::set
bool set(const Key &key)
Same as insert (no value to overwrite)
Definition: HashSet.H:197
Foam::pos
dimensionedScalar pos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:177