FreeFOAM The Cross-Platform CFD Toolkit
polyMeshAdder.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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 \*----------------------------------------------------------------------------*/
25 
26 #include "polyMeshAdder.H"
28 #include <OpenFOAM/IOobject.H>
29 #include "faceCoupleInfo.H"
31 #include <OpenFOAM/SortableList.H>
32 #include <OpenFOAM/Time.H>
34 #include <OpenFOAM/mergePoints.H>
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 //- Append all mapped elements of a list to a DynamicList
42 void Foam::polyMeshAdder::append
43 (
44  const labelList& map,
45  const labelList& lst,
46  DynamicList<label>& dynLst
47 )
48 {
49  dynLst.setCapacity(dynLst.size() + lst.size());
50 
51  forAll(lst, i)
52  {
53  label newElem = map[lst[i]];
54 
55  if (newElem != -1)
56  {
57  dynLst.append(newElem);
58  }
59  }
60 }
61 
62 
63 //- Append all mapped elements of a list to a DynamicList
65 (
66  const labelList& map,
67  const labelList& lst,
68  const SortableList<label>& sortedLst,
69  DynamicList<label>& dynLst
70 )
71 {
72  dynLst.setSize(dynLst.size() + lst.size());
73 
74  forAll(lst, i)
75  {
76  label newElem = map[lst[i]];
77 
78  if (newElem != -1 && findSortedIndex(sortedLst, newElem) == -1)
79  {
80  dynLst.append(newElem);
81  }
82  }
83 }
84 
85 
86 // Get index of patch in new set of patchnames/types
87 Foam::label Foam::polyMeshAdder::patchIndex
88 (
89  const polyPatch& p,
90  DynamicList<word>& allPatchNames,
91  DynamicList<word>& allPatchTypes
92 )
93 {
94  // Find the patch name on the list. If the patch is already there
95  // and patch types match, return index
96  const word& pType = p.type();
97  const word& pName = p.name();
98 
99  label patchI = findIndex(allPatchNames, pName);
100 
101  if (patchI == -1)
102  {
103  // Patch not found. Append to the list
104  allPatchNames.append(pName);
105  allPatchTypes.append(pType);
106 
107  return allPatchNames.size() - 1;
108  }
109  else if (allPatchTypes[patchI] == pType)
110  {
111  // Found name and types match
112  return patchI;
113  }
114  else
115  {
116  // Found the name, but type is different
117 
118  // Duplicate name is not allowed. Create a composite name from the
119  // patch name and case name
120  const word& caseName = p.boundaryMesh().mesh().time().caseName();
121 
122  allPatchNames.append(pName + "_" + caseName);
123  allPatchTypes.append(pType);
124 
125  Pout<< "label patchIndex(const polyPatch& p) : "
126  << "Patch " << p.index() << " named "
127  << pName << " in mesh " << caseName
128  << " already exists, but patch types"
129  << " do not match.\nCreating a composite name as "
130  << allPatchNames[allPatchNames.size() - 1] << endl;
131 
132  return allPatchNames.size() - 1;
133  }
134 }
135 
136 
137 // Get index of zone in new set of zone names
138 Foam::label Foam::polyMeshAdder::zoneIndex
139 (
140  const word& curName,
141  DynamicList<word>& names
142 )
143 {
144  label zoneI = findIndex(names, curName);
145 
146  if (zoneI != -1)
147  {
148  return zoneI;
149  }
150  else
151  {
152  // Not found. Add new name to the list
153  names.append(curName);
154 
155  return names.size() - 1;
156  }
157 }
158 
159 
160 void Foam::polyMeshAdder::mergePatchNames
161 (
162  const polyBoundaryMesh& patches0,
163  const polyBoundaryMesh& patches1,
164 
165  DynamicList<word>& allPatchNames,
166  DynamicList<word>& allPatchTypes,
167 
168  labelList& from1ToAllPatches,
169  labelList& fromAllTo1Patches
170 )
171 {
172  // Insert the mesh0 patches and zones
173  append(patches0.names(), allPatchNames);
174  append(patches0.types(), allPatchTypes);
175 
176 
177  // Patches
178  // ~~~~~~~
179  // Patches from 0 are taken over as is; those from 1 get either merged
180  // (if they share name and type) or appended.
181  // Empty patches are filtered out much much later on.
182 
183  // Add mesh1 patches and build map both ways.
184  from1ToAllPatches.setSize(patches1.size());
185 
186  forAll(patches1, patchI)
187  {
188  from1ToAllPatches[patchI] = patchIndex
189  (
190  patches1[patchI],
191  allPatchNames,
192  allPatchTypes
193  );
194  }
195  allPatchTypes.shrink();
196  allPatchNames.shrink();
197 
198  // Invert 1 to all patch map
199  fromAllTo1Patches.setSize(allPatchNames.size());
200  fromAllTo1Patches = -1;
201 
202  forAll(from1ToAllPatches, i)
203  {
204  fromAllTo1Patches[from1ToAllPatches[i]] = i;
205  }
206 }
207 
208 
209 Foam::labelList Foam::polyMeshAdder::getPatchStarts
210 (
211  const polyBoundaryMesh& patches
212 )
213 {
214  labelList patchStarts(patches.size());
215  forAll(patches, patchI)
216  {
217  patchStarts[patchI] = patches[patchI].start();
218  }
219  return patchStarts;
220 }
221 
222 
223 Foam::labelList Foam::polyMeshAdder::getPatchSizes
224 (
225  const polyBoundaryMesh& patches
226 )
227 {
228  labelList patchSizes(patches.size());
229  forAll(patches, patchI)
230  {
231  patchSizes[patchI] = patches[patchI].size();
232  }
233  return patchSizes;
234 }
235 
236 
237 Foam::List<Foam::polyPatch*> Foam::polyMeshAdder::combinePatches
238 (
239  const polyMesh& mesh0,
240  const polyMesh& mesh1,
241  const polyBoundaryMesh& allBoundaryMesh,
242  const label nAllPatches,
243  const labelList& fromAllTo1Patches,
244 
245  const label nInternalFaces,
246  const labelList& nFaces,
247 
248  labelList& from0ToAllPatches,
249  labelList& from1ToAllPatches
250 )
251 {
252  const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
253  const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
254 
255 
256  // Compacted new patch list.
257  DynamicList<polyPatch*> allPatches(nAllPatches);
258 
259 
260  // Map from 0 to all patches (since gets compacted)
261  from0ToAllPatches.setSize(patches0.size());
262  from0ToAllPatches = -1;
263 
264  label startFaceI = nInternalFaces;
265 
266  // Copy patches0 with new sizes. First patches always come from
267  // mesh0 and will always be present.
268  for (label patchI = 0; patchI < patches0.size(); patchI++)
269  {
270  // Originates from mesh0. Clone with new size & filter out empty
271  // patch.
272  label filteredPatchI;
273 
274  if (nFaces[patchI] == 0 && isA<processorPolyPatch>(patches0[patchI]))
275  {
276  //Pout<< "Removing zero sized mesh0 patch "
277  // << patches0[patchI].name() << endl;
278  filteredPatchI = -1;
279  }
280  else
281  {
282  filteredPatchI = allPatches.size();
283 
284  allPatches.append
285  (
286  patches0[patchI].clone
287  (
288  allBoundaryMesh,
289  filteredPatchI,
290  nFaces[patchI],
291  startFaceI
292  ).ptr()
293  );
294  startFaceI += nFaces[patchI];
295  }
296 
297  // Record new index in allPatches
298  from0ToAllPatches[patchI] = filteredPatchI;
299 
300  // Check if patch was also in mesh1 and update its addressing if so.
301  if (fromAllTo1Patches[patchI] != -1)
302  {
303  from1ToAllPatches[fromAllTo1Patches[patchI]] = filteredPatchI;
304  }
305  }
306 
307  // Copy unique patches of mesh1.
308  forAll(from1ToAllPatches, patchI)
309  {
310  label allPatchI = from1ToAllPatches[patchI];
311 
312  if (allPatchI >= patches0.size())
313  {
314  // Patch has not been merged with any mesh0 patch.
315 
316  label filteredPatchI;
317 
318  if
319  (
320  nFaces[allPatchI] == 0
321  && isA<processorPolyPatch>(patches1[patchI])
322  )
323  {
324  //Pout<< "Removing zero sized mesh1 patch "
325  // << patches1[patchI].name() << endl;
326  filteredPatchI = -1;
327  }
328  else
329  {
330  filteredPatchI = allPatches.size();
331 
332  allPatches.append
333  (
334  patches1[patchI].clone
335  (
336  allBoundaryMesh,
337  filteredPatchI,
338  nFaces[allPatchI],
339  startFaceI
340  ).ptr()
341  );
342  startFaceI += nFaces[allPatchI];
343  }
344 
345  from1ToAllPatches[patchI] = filteredPatchI;
346  }
347  }
348 
349  allPatches.shrink();
350 
351  return allPatches;
352 }
353 
354 
355 Foam::labelList Foam::polyMeshAdder::getFaceOrder
356 (
357  const cellList& cells,
358  const label nInternalFaces,
359  const labelList& owner,
360  const labelList& neighbour
361 )
362 {
363  labelList oldToNew(owner.size(), -1);
364 
365  // Leave boundary faces in order
366  for (label faceI = nInternalFaces; faceI < owner.size(); faceI++)
367  {
368  oldToNew[faceI] = faceI;
369  }
370 
371  // First unassigned face
372  label newFaceI = 0;
373 
374  forAll(cells, cellI)
375  {
376  const labelList& cFaces = cells[cellI];
377 
378  SortableList<label> nbr(cFaces.size());
379 
380  forAll(cFaces, i)
381  {
382  label faceI = cFaces[i];
383 
384  label nbrCellI = neighbour[faceI];
385 
386  if (nbrCellI != -1)
387  {
388  // Internal face. Get cell on other side.
389  if (nbrCellI == cellI)
390  {
391  nbrCellI = owner[faceI];
392  }
393 
394  if (cellI < nbrCellI)
395  {
396  // CellI is master
397  nbr[i] = nbrCellI;
398  }
399  else
400  {
401  // nbrCell is master. Let it handle this face.
402  nbr[i] = -1;
403  }
404  }
405  else
406  {
407  // External face. Do later.
408  nbr[i] = -1;
409  }
410  }
411 
412  nbr.sort();
413 
414  forAll(nbr, i)
415  {
416  if (nbr[i] != -1)
417  {
418  oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
419  }
420  }
421  }
422 
423 
424  // Check done all faces.
425  forAll(oldToNew, faceI)
426  {
427  if (oldToNew[faceI] == -1)
428  {
430  (
431  "polyMeshAdder::getFaceOrder"
432  "(const cellList&, const label, const labelList&"
433  ", const labelList&) const"
434  ) << "Did not determine new position"
435  << " for face " << faceI
436  << abort(FatalError);
437  }
438  }
439 
440  return oldToNew;
441 }
442 
443 
444 // Extends face f with split points. cutEdgeToPoints gives for every
445 // edge the points introduced inbetween the endpoints.
446 void Foam::polyMeshAdder::insertVertices
447 (
448  const edgeLookup& cutEdgeToPoints,
449  const Map<label>& meshToMaster,
450  const labelList& masterToCutPoints,
451  const face& masterF,
452 
453  DynamicList<label>& workFace,
454  face& allF
455 )
456 {
457  workFace.clear();
458 
459  // Check any edge for being cut (check on the cut so takes account
460  // for any point merging on the cut)
461 
462  forAll(masterF, fp)
463  {
464  label v0 = masterF[fp];
465  label v1 = masterF.nextLabel(fp);
466 
467  // Copy existing face point
468  workFace.append(allF[fp]);
469 
470  // See if any edge between v0,v1
471 
472  Map<label>::const_iterator v0Fnd = meshToMaster.find(v0);
473  if (v0Fnd != meshToMaster.end())
474  {
475  Map<label>::const_iterator v1Fnd = meshToMaster.find(v1);
476  if (v1Fnd != meshToMaster.end())
477  {
478  // Get edge in cutPoint numbering
479  edge cutEdge
480  (
481  masterToCutPoints[v0Fnd()],
482  masterToCutPoints[v1Fnd()]
483  );
484 
485  edgeLookup::const_iterator iter = cutEdgeToPoints.find(cutEdge);
486 
487  if (iter != cutEdgeToPoints.end())
488  {
489  const edge& e = iter.key();
490  const labelList& addedPoints = iter();
491 
492  // cutPoints first in allPoints so no need for renumbering
493  if (e[0] == cutEdge[0])
494  {
495  forAll(addedPoints, i)
496  {
497  workFace.append(addedPoints[i]);
498  }
499  }
500  else
501  {
502  forAllReverse(addedPoints, i)
503  {
504  workFace.append(addedPoints[i]);
505  }
506  }
507  }
508  }
509  }
510  }
511 
512  if (workFace.size() != allF.size())
513  {
514  allF.transfer(workFace);
515  }
516 }
517 
518 
519 // Adds primitives (cells, faces, points)
520 // Cells:
521 // - all of mesh0
522 // - all of mesh1
523 // Faces:
524 // - all uncoupled of mesh0
525 // - all coupled faces
526 // - all uncoupled of mesh1
527 // Points:
528 // - all coupled
529 // - all uncoupled of mesh0
530 // - all uncoupled of mesh1
531 void Foam::polyMeshAdder::mergePrimitives
532 (
533  const polyMesh& mesh0,
534  const polyMesh& mesh1,
535  const faceCoupleInfo& coupleInfo,
536 
537  const label nAllPatches, // number of patches in the new mesh
538  const labelList& fromAllTo1Patches,
539  const labelList& from1ToAllPatches,
540 
541  pointField& allPoints,
542  labelList& from0ToAllPoints,
543  labelList& from1ToAllPoints,
544 
545  faceList& allFaces,
546  labelList& allOwner,
547  labelList& allNeighbour,
548  label& nInternalFaces,
549  labelList& nFacesPerPatch,
550  label& nCells,
551 
552  labelList& from0ToAllFaces,
553  labelList& from1ToAllFaces,
554  labelList& from1ToAllCells
555 )
556 {
557  const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
558  const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
559 
560  const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
561  const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
562  const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();
563 
564 
565  // Points
566  // ~~~~~~
567 
568  // Storage for new points
569  allPoints.setSize(mesh0.nPoints() + mesh1.nPoints());
570  label allPointI = 0;
571 
572  from0ToAllPoints.setSize(mesh0.nPoints());
573  from0ToAllPoints = -1;
574  from1ToAllPoints.setSize(mesh1.nPoints());
575  from1ToAllPoints = -1;
576 
577  // Copy coupled points (on cut)
578  {
579  const pointField& cutPoints = coupleInfo.cutPoints();
580 
581  //const labelListList& cutToMasterPoints =
582  // coupleInfo.cutToMasterPoints();
583  labelListList cutToMasterPoints
584  (
586  (
587  cutPoints.size(),
588  coupleInfo.masterToCutPoints()
589  )
590  );
591 
592  //const labelListList& cutToSlavePoints =
593  // coupleInfo.cutToSlavePoints();
594  labelListList cutToSlavePoints
595  (
597  (
598  cutPoints.size(),
599  coupleInfo.slaveToCutPoints()
600  )
601  );
602 
603  forAll(cutPoints, i)
604  {
605  allPoints[allPointI] = cutPoints[i];
606 
607  // Mark all master and slave points referring to this point.
608 
609  const labelList& masterPoints = cutToMasterPoints[i];
610 
611  forAll(masterPoints, j)
612  {
613  label mesh0PointI = masterPatch.meshPoints()[masterPoints[j]];
614  from0ToAllPoints[mesh0PointI] = allPointI;
615  }
616 
617  const labelList& slavePoints = cutToSlavePoints[i];
618 
619  forAll(slavePoints, j)
620  {
621  label mesh1PointI = slavePatch.meshPoints()[slavePoints[j]];
622  from1ToAllPoints[mesh1PointI] = allPointI;
623  }
624  allPointI++;
625  }
626  }
627 
628  // Add uncoupled mesh0 points
629  forAll(mesh0.points(), pointI)
630  {
631  if (from0ToAllPoints[pointI] == -1)
632  {
633  allPoints[allPointI] = mesh0.points()[pointI];
634  from0ToAllPoints[pointI] = allPointI;
635  allPointI++;
636  }
637  }
638 
639  // Add uncoupled mesh1 points
640  forAll(mesh1.points(), pointI)
641  {
642  if (from1ToAllPoints[pointI] == -1)
643  {
644  allPoints[allPointI] = mesh1.points()[pointI];
645  from1ToAllPoints[pointI] = allPointI;
646  allPointI++;
647  }
648  }
649 
650  allPoints.setSize(allPointI);
651 
652 
653  // Faces
654  // ~~~~~
655 
656  // Sizes per patch
657  nFacesPerPatch.setSize(nAllPatches);
658  nFacesPerPatch = 0;
659 
660  // Storage for faces and owner/neighbour
661  allFaces.setSize(mesh0.nFaces() + mesh1.nFaces());
662  allOwner.setSize(allFaces.size());
663  allOwner = -1;
664  allNeighbour.setSize(allFaces.size());
665  allNeighbour = -1;
666  label allFaceI = 0;
667 
668  from0ToAllFaces.setSize(mesh0.nFaces());
669  from0ToAllFaces = -1;
670  from1ToAllFaces.setSize(mesh1.nFaces());
671  from1ToAllFaces = -1;
672 
673  // Copy mesh0 internal faces (always uncoupled)
674  for (label faceI = 0; faceI < mesh0.nInternalFaces(); faceI++)
675  {
676  allFaces[allFaceI] = renumber(from0ToAllPoints, mesh0.faces()[faceI]);
677  allOwner[allFaceI] = mesh0.faceOwner()[faceI];
678  allNeighbour[allFaceI] = mesh0.faceNeighbour()[faceI];
679  from0ToAllFaces[faceI] = allFaceI++;
680  }
681 
682  // Copy coupled faces. Every coupled face has an equivalent master and
683  // slave. Also uncount as boundary faces all the newly coupled faces.
684  const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
685  const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();
686 
687  forAll(cutFaces, i)
688  {
689  label masterFaceI = cutToMasterFaces[i];
690 
691  label mesh0FaceI = masterPatch.addressing()[masterFaceI];
692 
693  if (from0ToAllFaces[mesh0FaceI] == -1)
694  {
695  // First occurrence of face
696  from0ToAllFaces[mesh0FaceI] = allFaceI;
697 
698  // External face becomes internal so uncount
699  label patch0 = patches0.whichPatch(mesh0FaceI);
700  nFacesPerPatch[patch0]--;
701  }
702 
703  label slaveFaceI = cutToSlaveFaces[i];
704 
705  label mesh1FaceI = slavePatch.addressing()[slaveFaceI];
706 
707  if (from1ToAllFaces[mesh1FaceI] == -1)
708  {
709  from1ToAllFaces[mesh1FaceI] = allFaceI;
710 
711  label patch1 = patches1.whichPatch(mesh1FaceI);
712  nFacesPerPatch[from1ToAllPatches[patch1]]--;
713  }
714 
715  // Copy cut face (since cutPoints are copied first no renumbering
716  // nessecary)
717  allFaces[allFaceI] = cutFaces[i];
718  allOwner[allFaceI] = mesh0.faceOwner()[mesh0FaceI];
719  allNeighbour[allFaceI] = mesh1.faceOwner()[mesh1FaceI] + mesh0.nCells();
720 
721  allFaceI++;
722  }
723 
724  // Copy mesh1 internal faces (always uncoupled)
725  for (label faceI = 0; faceI < mesh1.nInternalFaces(); faceI++)
726  {
727  allFaces[allFaceI] = renumber(from1ToAllPoints, mesh1.faces()[faceI]);
728  allOwner[allFaceI] = mesh1.faceOwner()[faceI] + mesh0.nCells();
729  allNeighbour[allFaceI] = mesh1.faceNeighbour()[faceI] + mesh0.nCells();
730  from1ToAllFaces[faceI] = allFaceI++;
731  }
732 
733  nInternalFaces = allFaceI;
734 
735  // Copy (unmarked/uncoupled) external faces in new order.
736  for (label allPatchI = 0; allPatchI < nAllPatches; allPatchI++)
737  {
738  if (allPatchI < patches0.size())
739  {
740  // Patch is present in mesh0
741  const polyPatch& pp = patches0[allPatchI];
742 
743  nFacesPerPatch[allPatchI] += pp.size();
744 
745  label faceI = pp.start();
746 
747  forAll(pp, i)
748  {
749  if (from0ToAllFaces[faceI] == -1)
750  {
751  // Is uncoupled face since has not yet been dealt with
752  allFaces[allFaceI] = renumber
753  (
754  from0ToAllPoints,
755  mesh0.faces()[faceI]
756  );
757  allOwner[allFaceI] = mesh0.faceOwner()[faceI];
758  allNeighbour[allFaceI] = -1;
759 
760  from0ToAllFaces[faceI] = allFaceI++;
761  }
762  faceI++;
763  }
764  }
765  if (fromAllTo1Patches[allPatchI] != -1)
766  {
767  // Patch is present in mesh1
768  const polyPatch& pp = patches1[fromAllTo1Patches[allPatchI]];
769 
770  nFacesPerPatch[allPatchI] += pp.size();
771 
772  label faceI = pp.start();
773 
774  forAll(pp, i)
775  {
776  if (from1ToAllFaces[faceI] == -1)
777  {
778  // Is uncoupled face
779  allFaces[allFaceI] = renumber
780  (
781  from1ToAllPoints,
782  mesh1.faces()[faceI]
783  );
784  allOwner[allFaceI] =
785  mesh1.faceOwner()[faceI]
786  + mesh0.nCells();
787  allNeighbour[allFaceI] = -1;
788 
789  from1ToAllFaces[faceI] = allFaceI++;
790  }
791  faceI++;
792  }
793  }
794  }
795  allFaces.setSize(allFaceI);
796  allOwner.setSize(allFaceI);
797  allNeighbour.setSize(allFaceI);
798 
799 
800  // So now we have all ok for one-to-one mapping.
801  // For split slace faces:
802  // - mesh consistent with slave side
803  // - mesh not consistent with owner side. It is not zipped up, the
804  // original faces need edges split.
805 
806  // Use brute force to prevent having to calculate addressing:
807  // - get map from master edge to split edges.
808  // - check all faces to find any edge that is split.
809  {
810  // From two cut-points to labels of cut-points inbetween.
811  // (in order: from e[0] to e[1]
812  const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();
813 
814  // Get map of master face (in mesh labels) that are in cut. These faces
815  // do not need to be renumbered.
816  labelHashSet masterCutFaces(cutToMasterFaces.size());
817  forAll(cutToMasterFaces, i)
818  {
819  label meshFaceI = masterPatch.addressing()[cutToMasterFaces[i]];
820 
821  masterCutFaces.insert(meshFaceI);
822  }
823 
824  DynamicList<label> workFace(100);
825 
826  forAll(from0ToAllFaces, face0)
827  {
828  if (!masterCutFaces.found(face0))
829  {
830  label allFaceI = from0ToAllFaces[face0];
831 
832  insertVertices
833  (
834  cutEdgeToPoints,
835  masterPatch.meshPointMap(),
836  coupleInfo.masterToCutPoints(),
837  mesh0.faces()[face0],
838 
839  workFace,
840  allFaces[allFaceI]
841  );
842  }
843  }
844 
845  // Same for slave face
846 
847  labelHashSet slaveCutFaces(cutToSlaveFaces.size());
848  forAll(cutToSlaveFaces, i)
849  {
850  label meshFaceI = slavePatch.addressing()[cutToSlaveFaces[i]];
851 
852  slaveCutFaces.insert(meshFaceI);
853  }
854 
855  forAll(from1ToAllFaces, face1)
856  {
857  if (!slaveCutFaces.found(face1))
858  {
859  label allFaceI = from1ToAllFaces[face1];
860 
861  insertVertices
862  (
863  cutEdgeToPoints,
864  slavePatch.meshPointMap(),
865  coupleInfo.slaveToCutPoints(),
866  mesh1.faces()[face1],
867 
868  workFace,
869  allFaces[allFaceI]
870  );
871  }
872  }
873  }
874 
875  // Now we have a full facelist and owner/neighbour addressing.
876 
877 
878  // Cells
879  // ~~~~~
880 
881  from1ToAllCells.setSize(mesh1.nCells());
882  from1ToAllCells = -1;
883 
884  forAll(mesh1.cells(), i)
885  {
886  from1ToAllCells[i] = i + mesh0.nCells();
887  }
888 
889  // Make cells (= cell-face addressing)
890  nCells = mesh0.nCells() + mesh1.nCells();
891  cellList allCells(nCells);
892  primitiveMesh::calcCells(allCells, allOwner, allNeighbour, nCells);
893 
894  // Reorder faces for upper-triangular order.
895  labelList oldToNew
896  (
897  getFaceOrder
898  (
899  allCells,
900  nInternalFaces,
901  allOwner,
902  allNeighbour
903  )
904  );
905 
906  inplaceReorder(oldToNew, allFaces);
907  inplaceReorder(oldToNew, allOwner);
908  inplaceReorder(oldToNew, allNeighbour);
909  inplaceRenumber(oldToNew, from0ToAllFaces);
910  inplaceRenumber(oldToNew, from1ToAllFaces);
911 }
912 
913 
914 void Foam::polyMeshAdder::mergePointZones
915 (
916  const pointZoneMesh& pz0,
917  const pointZoneMesh& pz1,
918  const labelList& from0ToAllPoints,
919  const labelList& from1ToAllPoints,
920 
921  DynamicList<word>& zoneNames,
922  labelList& from1ToAll,
923  List<DynamicList<label> >& pzPoints
924 )
925 {
926  zoneNames.setCapacity(pz0.size() + pz1.size());
927 
928  // Names
929  append(pz0.names(), zoneNames);
930 
931  from1ToAll.setSize(pz1.size());
932 
933  forAll (pz1, zoneI)
934  {
935  from1ToAll[zoneI] = zoneIndex(pz1[zoneI].name(), zoneNames);
936  }
937  zoneNames.shrink();
938 
939  // Point labels per merged zone
940  pzPoints.setSize(zoneNames.size());
941 
942  forAll(pz0, zoneI)
943  {
944  append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
945  }
946 
947  // Get sorted zone contents for duplicate element recognition
948  PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
949  forAll(pzPoints, zoneI)
950  {
951  pzPointsSorted.set
952  (
953  zoneI,
954  new SortableList<label>(pzPoints[zoneI])
955  );
956  }
957 
958  // Now we have full addressing for points so do the pointZones of mesh1.
959  forAll(pz1, zoneI)
960  {
961  // Relabel all points of zone and add to correct pzPoints.
962  label allZoneI = from1ToAll[zoneI];
963 
964  append
965  (
966  from1ToAllPoints,
967  pz1[zoneI],
968  pzPointsSorted[allZoneI],
969  pzPoints[allZoneI]
970  );
971  }
972 
973  forAll(pzPoints, i)
974  {
975  pzPoints[i].shrink();
976  }
977 }
978 
979 
980 void Foam::polyMeshAdder::mergeFaceZones
981 (
982  const faceZoneMesh& fz0,
983  const faceZoneMesh& fz1,
984  const labelList& from0ToAllFaces,
985  const labelList& from1ToAllFaces,
986 
987  DynamicList<word>& zoneNames,
988  labelList& from1ToAll,
989  List<DynamicList<label> >& fzFaces,
990  List<DynamicList<bool> >& fzFlips
991 )
992 {
993  zoneNames.setCapacity(fz0.size() + fz1.size());
994 
995  append(fz0.names(), zoneNames);
996 
997  from1ToAll.setSize(fz1.size());
998 
999  forAll (fz1, zoneI)
1000  {
1001  from1ToAll[zoneI] = zoneIndex(fz1[zoneI].name(), zoneNames);
1002  }
1003  zoneNames.shrink();
1004 
1005 
1006  // Create temporary lists for faceZones.
1007  fzFaces.setSize(zoneNames.size());
1008  fzFlips.setSize(zoneNames.size());
1009  forAll(fz0, zoneI)
1010  {
1011  DynamicList<label>& newZone = fzFaces[zoneI];
1012  DynamicList<bool>& newFlip = fzFlips[zoneI];
1013 
1014  newZone.setCapacity(fz0[zoneI].size());
1015  newFlip.setCapacity(newZone.size());
1016 
1017  const labelList& addressing = fz0[zoneI];
1018  const boolList& flipMap = fz0[zoneI].flipMap();
1019 
1020  forAll(addressing, i)
1021  {
1022  label faceI = addressing[i];
1023 
1024  if (from0ToAllFaces[faceI] != -1)
1025  {
1026  newZone.append(from0ToAllFaces[faceI]);
1027  newFlip.append(flipMap[i]);
1028  }
1029  }
1030  }
1031 
1032  // Get sorted zone contents for duplicate element recognition
1033  PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
1034  forAll(fzFaces, zoneI)
1035  {
1036  fzFacesSorted.set
1037  (
1038  zoneI,
1039  new SortableList<label>(fzFaces[zoneI])
1040  );
1041  }
1042 
1043  // Now we have full addressing for faces so do the faceZones of mesh1.
1044  forAll(fz1, zoneI)
1045  {
1046  label allZoneI = from1ToAll[zoneI];
1047 
1048  DynamicList<label>& newZone = fzFaces[allZoneI];
1049  const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
1050  DynamicList<bool>& newFlip = fzFlips[allZoneI];
1051 
1052  newZone.setCapacity(newZone.size() + fz1[zoneI].size());
1053  newFlip.setCapacity(newZone.size());
1054 
1055  const labelList& addressing = fz1[zoneI];
1056  const boolList& flipMap = fz1[zoneI].flipMap();
1057 
1058  forAll(addressing, i)
1059  {
1060  label faceI = addressing[i];
1061  label allFaceI = from1ToAllFaces[faceI];
1062 
1063  if
1064  (
1065  allFaceI != -1
1066  && findSortedIndex(newZoneSorted, allFaceI) == -1
1067  )
1068  {
1069  newZone.append(allFaceI);
1070  newFlip.append(flipMap[i]);
1071  }
1072  }
1073  }
1074 
1075  forAll(fzFaces, i)
1076  {
1077  fzFaces[i].shrink();
1078  fzFlips[i].shrink();
1079  }
1080 }
1081 
1082 
1083 void Foam::polyMeshAdder::mergeCellZones
1084 (
1085  const cellZoneMesh& cz0,
1086  const cellZoneMesh& cz1,
1087  const labelList& from1ToAllCells,
1088 
1089  DynamicList<word>& zoneNames,
1090  labelList& from1ToAll,
1091  List<DynamicList<label> >& czCells
1092 )
1093 {
1094  zoneNames.setCapacity(cz0.size() + cz1.size());
1095 
1096  append(cz0.names(), zoneNames);
1097 
1098  from1ToAll.setSize(cz1.size());
1099  forAll (cz1, zoneI)
1100  {
1101  from1ToAll[zoneI] = zoneIndex(cz1[zoneI].name(), zoneNames);
1102  }
1103  zoneNames.shrink();
1104 
1105 
1106  // Create temporary lists for cellZones.
1107  czCells.setSize(zoneNames.size());
1108  forAll(cz0, zoneI)
1109  {
1110  // Insert mesh0 cells
1111  append(cz0[zoneI], czCells[zoneI]);
1112  }
1113 
1114 
1115  // Cell mapping is trivial.
1116  forAll(cz1, zoneI)
1117  {
1118  label allZoneI = from1ToAll[zoneI];
1119 
1120  append(from1ToAllCells, cz1[zoneI], czCells[allZoneI]);
1121  }
1122 
1123  forAll(czCells, i)
1124  {
1125  czCells[i].shrink();
1126  }
1127 }
1128 
1129 
1130 void Foam::polyMeshAdder::mergeZones
1131 (
1132  const polyMesh& mesh0,
1133  const polyMesh& mesh1,
1134  const labelList& from0ToAllPoints,
1135  const labelList& from0ToAllFaces,
1136 
1137  const labelList& from1ToAllPoints,
1138  const labelList& from1ToAllFaces,
1139  const labelList& from1ToAllCells,
1140 
1141  DynamicList<word>& pointZoneNames,
1142  List<DynamicList<label> >& pzPoints,
1143 
1144  DynamicList<word>& faceZoneNames,
1145  List<DynamicList<label> >& fzFaces,
1146  List<DynamicList<bool> >& fzFlips,
1147 
1148  DynamicList<word>& cellZoneNames,
1149  List<DynamicList<label> >& czCells
1150 )
1151 {
1152  labelList from1ToAllPZones;
1153  mergePointZones
1154  (
1155  mesh0.pointZones(),
1156  mesh1.pointZones(),
1157  from0ToAllPoints,
1158  from1ToAllPoints,
1159 
1160  pointZoneNames,
1161  from1ToAllPZones,
1162  pzPoints
1163  );
1164 
1165  labelList from1ToAllFZones;
1166  mergeFaceZones
1167  (
1168  mesh0.faceZones(),
1169  mesh1.faceZones(),
1170  from0ToAllFaces,
1171  from1ToAllFaces,
1172 
1173  faceZoneNames,
1174  from1ToAllFZones,
1175  fzFaces,
1176  fzFlips
1177  );
1178 
1179  labelList from1ToAllCZones;
1180  mergeCellZones
1181  (
1182  mesh0.cellZones(),
1183  mesh1.cellZones(),
1184  from1ToAllCells,
1185 
1186  cellZoneNames,
1187  from1ToAllCZones,
1188  czCells
1189  );
1190 }
1191 
1192 
1193 void Foam::polyMeshAdder::addZones
1194 (
1195  const DynamicList<word>& pointZoneNames,
1196  const List<DynamicList<label> >& pzPoints,
1197 
1198  const DynamicList<word>& faceZoneNames,
1199  const List<DynamicList<label> >& fzFaces,
1200  const List<DynamicList<bool> >& fzFlips,
1201 
1202  const DynamicList<word>& cellZoneNames,
1203  const List<DynamicList<label> >& czCells,
1204 
1205  polyMesh& mesh
1206 )
1207 {
1208  List<pointZone*> pZones(pzPoints.size());
1209  forAll(pZones, i)
1210  {
1211  pZones[i] = new pointZone
1212  (
1213  pointZoneNames[i],
1214  pzPoints[i],
1215  i,
1216  mesh.pointZones()
1217  );
1218  }
1219 
1220  List<faceZone*> fZones(fzFaces.size());
1221  forAll(fZones, i)
1222  {
1223  fZones[i] = new faceZone
1224  (
1225  faceZoneNames[i],
1226  fzFaces[i],
1227  fzFlips[i],
1228  i,
1229  mesh.faceZones()
1230  );
1231  }
1232 
1233  List<cellZone*> cZones(czCells.size());
1234  forAll(cZones, i)
1235  {
1236  cZones[i] = new cellZone
1237  (
1238  cellZoneNames[i],
1239  czCells[i],
1240  i,
1241  mesh.cellZones()
1242  );
1243  }
1244 
1245  mesh.addZones
1246  (
1247  pZones,
1248  fZones,
1249  cZones
1250  );
1251 }
1252 
1253 
1254 
1255 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1256 
1257 // Returns new mesh and sets
1258 // - map from new cell/face/point/patch to either mesh0 or mesh1
1259 //
1260 // mesh0Faces/mesh1Faces: corresponding faces on both meshes.
1263  const IOobject& io,
1264  const polyMesh& mesh0,
1265  const polyMesh& mesh1,
1266  const faceCoupleInfo& coupleInfo,
1268 )
1269 {
1270  const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1271  const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1272 
1273 
1274  DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1275  DynamicList<word> allPatchTypes(allPatchNames.size());
1276 
1277  // Patch maps
1278  labelList from1ToAllPatches(patches1.size());
1279  labelList fromAllTo1Patches(allPatchNames.size(), -1);
1280 
1281  mergePatchNames
1282  (
1283  patches0,
1284  patches1,
1285  allPatchNames,
1286  allPatchTypes,
1287  from1ToAllPatches,
1288  fromAllTo1Patches
1289  );
1290 
1291 
1292  // New points
1293  pointField allPoints;
1294 
1295  // Map from mesh0/1 points to allPoints.
1296  labelList from0ToAllPoints(mesh0.nPoints(), -1);
1297  labelList from1ToAllPoints(mesh1.nPoints(), -1);
1298 
1299  // New faces
1300  faceList allFaces;
1301  label nInternalFaces;
1302 
1303  // New cells
1304  labelList allOwner;
1305  labelList allNeighbour;
1306  label nCells;
1307 
1308  // Sizes per patch
1309  labelList nFaces(allPatchNames.size(), 0);
1310 
1311  // Maps
1312  labelList from0ToAllFaces(mesh0.nFaces(), -1);
1313  labelList from1ToAllFaces(mesh1.nFaces(), -1);
1314 
1315  // Map
1316  labelList from1ToAllCells(mesh1.nCells(), -1);
1317 
1318  mergePrimitives
1319  (
1320  mesh0,
1321  mesh1,
1322  coupleInfo,
1323 
1324  allPatchNames.size(),
1325  fromAllTo1Patches,
1326  from1ToAllPatches,
1327 
1328  allPoints,
1329  from0ToAllPoints,
1330  from1ToAllPoints,
1331 
1332  allFaces,
1333  allOwner,
1334  allNeighbour,
1335  nInternalFaces,
1336  nFaces,
1337  nCells,
1338 
1339  from0ToAllFaces,
1340  from1ToAllFaces,
1341  from1ToAllCells
1342  );
1343 
1344 
1345  // Zones
1346  // ~~~~~
1347 
1348  DynamicList<word> pointZoneNames;
1349  List<DynamicList<label> > pzPoints;
1350 
1351  DynamicList<word> faceZoneNames;
1352  List<DynamicList<label> > fzFaces;
1353  List<DynamicList<bool> > fzFlips;
1354 
1355  DynamicList<word> cellZoneNames;
1356  List<DynamicList<label> > czCells;
1357 
1358  mergeZones
1359  (
1360  mesh0,
1361  mesh1,
1362 
1363  from0ToAllPoints,
1364  from0ToAllFaces,
1365 
1366  from1ToAllPoints,
1367  from1ToAllFaces,
1368  from1ToAllCells,
1369 
1370  pointZoneNames,
1371  pzPoints,
1372 
1373  faceZoneNames,
1374  fzFaces,
1375  fzFlips,
1376 
1377  cellZoneNames,
1378  czCells
1379  );
1380 
1381 
1382  // Patches
1383  // ~~~~~~~
1384 
1385  // Map from 0 to all patches (since gets compacted)
1386  labelList from0ToAllPatches(patches0.size(), -1);
1387 
1388  List<polyPatch*> allPatches
1389  (
1390  combinePatches
1391  (
1392  mesh0,
1393  mesh1,
1394  patches0, // Should be boundaryMesh() on new mesh.
1395  allPatchNames.size(),
1396  fromAllTo1Patches,
1397  mesh0.nInternalFaces()
1398  + mesh1.nInternalFaces()
1399  + coupleInfo.cutFaces().size(),
1400  nFaces,
1401 
1402  from0ToAllPatches,
1403  from1ToAllPatches
1404  )
1405  );
1406 
1407 
1408  // Map information
1409  // ~~~~~~~~~~~~~~~
1410 
1411  mapPtr.reset
1412  (
1413  new mapAddedPolyMesh
1414  (
1415  mesh0.nPoints(),
1416  mesh0.nFaces(),
1417  mesh0.nCells(),
1418 
1419  mesh1.nPoints(),
1420  mesh1.nFaces(),
1421  mesh1.nCells(),
1422 
1423  from0ToAllPoints,
1424  from0ToAllFaces,
1425  identity(mesh0.nCells()),
1426 
1427  from1ToAllPoints,
1428  from1ToAllFaces,
1429  from1ToAllCells,
1430 
1431  from0ToAllPatches,
1432  from1ToAllPatches,
1433  getPatchSizes(patches0),
1434  getPatchStarts(patches0)
1435  )
1436  );
1437 
1438 
1439 
1440  // Now we have extracted all information from all meshes.
1441  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1442 
1443  // Construct mesh
1444  autoPtr<polyMesh> tmesh
1445  (
1446  new polyMesh
1447  (
1448  io,
1449  xferMove(allPoints),
1450  xferMove(allFaces),
1451  xferMove(allOwner),
1452  xferMove(allNeighbour)
1453  )
1454  );
1455  polyMesh& mesh = tmesh();
1456 
1457  // Add zones to new mesh.
1458  addZones
1459  (
1460  pointZoneNames,
1461  pzPoints,
1462 
1463  faceZoneNames,
1464  fzFaces,
1465  fzFlips,
1466 
1467  cellZoneNames,
1468  czCells,
1469  mesh
1470  );
1471 
1472  // Add patches to new mesh
1473  mesh.addPatches(allPatches);
1474 
1475  return tmesh;
1476 }
1477 
1478 
1479 // Inplace add mesh1 to mesh0
1482  polyMesh& mesh0,
1483  const polyMesh& mesh1,
1484  const faceCoupleInfo& coupleInfo,
1485  const bool validBoundary
1486 )
1487 {
1488  const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
1489  const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
1490 
1491  DynamicList<word> allPatchNames(patches0.size() + patches1.size());
1492  DynamicList<word> allPatchTypes(allPatchNames.size());
1493 
1494  // Patch maps
1495  labelList from1ToAllPatches(patches1.size());
1496  labelList fromAllTo1Patches(allPatchNames.size(), -1);
1497 
1498  mergePatchNames
1499  (
1500  patches0,
1501  patches1,
1502  allPatchNames,
1503  allPatchTypes,
1504  from1ToAllPatches,
1505  fromAllTo1Patches
1506  );
1507 
1508 
1509  // New points
1510  pointField allPoints;
1511 
1512  // Map from mesh0/1 points to allPoints.
1513  labelList from0ToAllPoints(mesh0.nPoints(), -1);
1514  labelList from1ToAllPoints(mesh1.nPoints(), -1);
1515 
1516  // New faces
1517  faceList allFaces;
1518  labelList allOwner;
1519  labelList allNeighbour;
1520  label nInternalFaces;
1521  // Sizes per patch
1522  labelList nFaces(allPatchNames.size(), 0);
1523  label nCells;
1524 
1525  // Maps
1526  labelList from0ToAllFaces(mesh0.nFaces(), -1);
1527  labelList from1ToAllFaces(mesh1.nFaces(), -1);
1528  // Map
1529  labelList from1ToAllCells(mesh1.nCells(), -1);
1530 
1531  mergePrimitives
1532  (
1533  mesh0,
1534  mesh1,
1535  coupleInfo,
1536 
1537  allPatchNames.size(),
1538  fromAllTo1Patches,
1539  from1ToAllPatches,
1540 
1541  allPoints,
1542  from0ToAllPoints,
1543  from1ToAllPoints,
1544 
1545  allFaces,
1546  allOwner,
1547  allNeighbour,
1548  nInternalFaces,
1549  nFaces,
1550  nCells,
1551 
1552  from0ToAllFaces,
1553  from1ToAllFaces,
1554  from1ToAllCells
1555  );
1556 
1557 
1558  // Zones
1559  // ~~~~~
1560 
1561  DynamicList<word> pointZoneNames;
1562  List<DynamicList<label> > pzPoints;
1563 
1564  DynamicList<word> faceZoneNames;
1565  List<DynamicList<label> > fzFaces;
1566  List<DynamicList<bool> > fzFlips;
1567 
1568  DynamicList<word> cellZoneNames;
1569  List<DynamicList<label> > czCells;
1570 
1571  mergeZones
1572  (
1573  mesh0,
1574  mesh1,
1575 
1576  from0ToAllPoints,
1577  from0ToAllFaces,
1578 
1579  from1ToAllPoints,
1580  from1ToAllFaces,
1581  from1ToAllCells,
1582 
1583  pointZoneNames,
1584  pzPoints,
1585 
1586  faceZoneNames,
1587  fzFaces,
1588  fzFlips,
1589 
1590  cellZoneNames,
1591  czCells
1592  );
1593 
1594 
1595  // Patches
1596  // ~~~~~~~
1597 
1598 
1599  // Store mesh0 patch info before modifying patches0.
1600  labelList mesh0PatchSizes(getPatchSizes(patches0));
1601  labelList mesh0PatchStarts(getPatchStarts(patches0));
1602 
1603  // Map from 0 to all patches (since gets compacted)
1604  labelList from0ToAllPatches(patches0.size(), -1);
1605 
1606  // Inplace extend mesh0 patches (note that patches0.size() now also
1607  // has changed)
1608  polyBoundaryMesh& allPatches =
1609  const_cast<polyBoundaryMesh&>(mesh0.boundaryMesh());
1610  allPatches.setSize(allPatchNames.size());
1611 
1612  label startFaceI = nInternalFaces;
1613 
1614  // Copy patches0 with new sizes. First patches always come from
1615  // mesh0 and will always be present.
1616  label allPatchI = 0;
1617 
1618  forAll(from0ToAllPatches, patch0)
1619  {
1620  // Originates from mesh0. Clone with new size & filter out empty
1621  // patch.
1622 
1623  if (nFaces[patch0] == 0 && isA<processorPolyPatch>(allPatches[patch0]))
1624  {
1625  //Pout<< "Removing zero sized mesh0 patch " << allPatchNames[patch0]
1626  // << endl;
1627  from0ToAllPatches[patch0] = -1;
1628  // Check if patch was also in mesh1 and update its addressing if so.
1629  if (fromAllTo1Patches[patch0] != -1)
1630  {
1631  from1ToAllPatches[fromAllTo1Patches[patch0]] = -1;
1632  }
1633  }
1634  else
1635  {
1636  // Clone.
1637  allPatches.set
1638  (
1639  allPatchI,
1640  allPatches[patch0].clone
1641  (
1642  allPatches,
1643  allPatchI,
1644  nFaces[patch0],
1645  startFaceI
1646  )
1647  );
1648 
1649  // Record new index in allPatches
1650  from0ToAllPatches[patch0] = allPatchI;
1651 
1652  // Check if patch was also in mesh1 and update its addressing if so.
1653  if (fromAllTo1Patches[patch0] != -1)
1654  {
1655  from1ToAllPatches[fromAllTo1Patches[patch0]] = allPatchI;
1656  }
1657 
1658  startFaceI += nFaces[patch0];
1659 
1660  allPatchI++;
1661  }
1662  }
1663 
1664  // Copy unique patches of mesh1.
1665  forAll(from1ToAllPatches, patch1)
1666  {
1667  label uncompactAllPatchI = from1ToAllPatches[patch1];
1668 
1669  if (uncompactAllPatchI >= from0ToAllPatches.size())
1670  {
1671  // Patch has not been merged with any mesh0 patch.
1672 
1673  if
1674  (
1675  nFaces[uncompactAllPatchI] == 0
1676  && isA<processorPolyPatch>(patches1[patch1])
1677  )
1678  {
1679  //Pout<< "Removing zero sized mesh1 patch "
1680  // << allPatchNames[uncompactAllPatchI] << endl;
1681  from1ToAllPatches[patch1] = -1;
1682  }
1683  else
1684  {
1685  // Clone.
1686  allPatches.set
1687  (
1688  allPatchI,
1689  patches1[patch1].clone
1690  (
1691  allPatches,
1692  allPatchI,
1693  nFaces[uncompactAllPatchI],
1694  startFaceI
1695  )
1696  );
1697 
1698  // Record new index in allPatches
1699  from1ToAllPatches[patch1] = allPatchI;
1700 
1701  startFaceI += nFaces[uncompactAllPatchI];
1702 
1703  allPatchI++;
1704  }
1705  }
1706  }
1707 
1708  allPatches.setSize(allPatchI);
1709 
1710 
1711  // Construct map information before changing mesh0 primitives
1712  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1713 
1715  (
1716  new mapAddedPolyMesh
1717  (
1718  mesh0.nPoints(),
1719  mesh0.nFaces(),
1720  mesh0.nCells(),
1721 
1722  mesh1.nPoints(),
1723  mesh1.nFaces(),
1724  mesh1.nCells(),
1725 
1726  from0ToAllPoints,
1727  from0ToAllFaces,
1728  identity(mesh0.nCells()),
1729 
1730  from1ToAllPoints,
1731  from1ToAllFaces,
1732  from1ToAllCells,
1733 
1734  from0ToAllPatches,
1735  from1ToAllPatches,
1736 
1737  mesh0PatchSizes,
1738  mesh0PatchStarts
1739  )
1740  );
1741 
1742 
1743 
1744  // Now we have extracted all information from all meshes.
1745  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1746 
1747  labelList patchSizes(getPatchSizes(allPatches));
1748  labelList patchStarts(getPatchStarts(allPatches));
1749 
1750  mesh0.resetMotion(); // delete any oldPoints.
1751  mesh0.resetPrimitives
1752  (
1753  xferMove(allPoints),
1754  xferMove(allFaces),
1755  xferMove(allOwner),
1756  xferMove(allNeighbour),
1757  patchSizes, // size
1758  patchStarts, // patchstarts
1759  validBoundary // boundary valid?
1760  );
1761 
1762  // Add zones to new mesh.
1763  mesh0.pointZones().clear();
1764  mesh0.faceZones().clear();
1765  mesh0.cellZones().clear();
1766  addZones
1767  (
1768  pointZoneNames,
1769  pzPoints,
1770 
1771  faceZoneNames,
1772  fzFaces,
1773  fzFlips,
1774 
1775  cellZoneNames,
1776  czCells,
1777  mesh0
1778  );
1779 
1780  return mapPtr;
1781 }
1782 
1783 
1786  const polyMesh& mesh,
1787  const scalar mergeDist
1788 )
1789 {
1790  const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
1791  const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
1792 
1793  // Because of adding the missing pieces e.g. when redistributing a mesh
1794  // it can be that there are multiple points on the same processor that
1795  // refer to the same shared point.
1796 
1797  // Invert point-to-shared addressing
1798  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1799 
1800  Map<labelList> sharedToMesh(sharedPointLabels.size());
1801 
1802  label nMultiple = 0;
1803 
1804  forAll(sharedPointLabels, i)
1805  {
1806  label pointI = sharedPointLabels[i];
1807 
1808  label sharedI = sharedPointAddr[i];
1809 
1810  Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
1811 
1812  if (iter != sharedToMesh.end())
1813  {
1814  // sharedI already used by other point. Add this one.
1815 
1816  nMultiple++;
1817 
1818  labelList& connectedPointLabels = iter();
1819 
1820  label sz = connectedPointLabels.size();
1821 
1822  // Check just to make sure.
1823  if (findIndex(connectedPointLabels, pointI) != -1)
1824  {
1825  FatalErrorIn("polyMeshAdder::findSharedPoints(..)")
1826  << "Duplicate point in sharedPoint addressing." << endl
1827  << "When trying to add point " << pointI << " on shared "
1828  << sharedI << " with connected points "
1829  << connectedPointLabels
1830  << abort(FatalError);
1831  }
1832 
1833  connectedPointLabels.setSize(sz+1);
1834  connectedPointLabels[sz] = pointI;
1835  }
1836  else
1837  {
1838  sharedToMesh.insert(sharedI, labelList(1, pointI));
1839  }
1840  }
1841 
1842 
1843  // Assign single master for every shared with multiple geometric points
1844  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1845 
1846  Map<label> pointToMaster(nMultiple);
1847 
1848  forAllConstIter(Map<labelList>, sharedToMesh, iter)
1849  {
1850  const labelList& connectedPointLabels = iter();
1851 
1852  //Pout<< "For shared:" << iter.key()
1853  // << " found points:" << connectedPointLabels
1854  // << " at coords:"
1855  // << pointField(mesh.points(), connectedPointLabels) << endl;
1856 
1857  if (connectedPointLabels.size() > 1)
1858  {
1859  const pointField connectedPoints
1860  (
1861  mesh.points(),
1862  connectedPointLabels
1863  );
1864 
1865  labelList toMergedPoints;
1866  pointField mergedPoints;
1867  bool hasMerged = Foam::mergePoints
1868  (
1869  connectedPoints,
1870  mergeDist,
1871  false,
1872  toMergedPoints,
1873  mergedPoints
1874  );
1875 
1876  if (hasMerged)
1877  {
1878  // Invert toMergedPoints
1879  const labelListList mergeSets
1880  (
1882  (
1883  mergedPoints.size(),
1884  toMergedPoints
1885  )
1886  );
1887 
1888  // Find master for valid merges
1889  forAll(mergeSets, setI)
1890  {
1891  const labelList& mergeSet = mergeSets[setI];
1892 
1893  if (mergeSet.size() > 1)
1894  {
1895  // Pick lowest numbered point
1896  label masterPointI = labelMax;
1897 
1898  forAll(mergeSet, i)
1899  {
1900  label pointI = connectedPointLabels[mergeSet[i]];
1901 
1902  masterPointI = min(masterPointI, pointI);
1903  }
1904 
1905  forAll(mergeSet, i)
1906  {
1907  label pointI = connectedPointLabels[mergeSet[i]];
1908 
1909  //Pout<< "Merging point " << pointI
1910  // << " at " << mesh.points()[pointI]
1911  // << " into master point "
1912  // << masterPointI
1913  // << " at " << mesh.points()[masterPointI]
1914  // << endl;
1915 
1916  pointToMaster.insert(pointI, masterPointI);
1917  }
1918  }
1919  }
1920  }
1921  }
1922  }
1923 
1924  //- Old: geometric merging. Causes problems for two close shared points.
1925  //labelList sharedToMerged;
1926  //pointField mergedPoints;
1927  //bool hasMerged = Foam::mergePoints
1928  //(
1929  // pointField
1930  // (
1931  // mesh.points(),
1932  // sharedPointLabels
1933  // ),
1934  // mergeDist,
1935  // false,
1936  // sharedToMerged,
1937  // mergedPoints
1938  //);
1939  //
1942  //
1943  //Map<label> pointToMaster(10*sharedToMerged.size());
1944  //
1945  //if (hasMerged)
1946  //{
1947  // labelListList mergeSets
1948  // (
1949  // invertOneToMany
1950  // (
1951  // sharedToMerged.size(),
1952  // sharedToMerged
1953  // )
1954  // );
1955  //
1956  // label nMergeSets = 0;
1957  //
1958  // forAll(mergeSets, setI)
1959  // {
1960  // const labelList& mergeSet = mergeSets[setI];
1961  //
1962  // if (mergeSet.size() > 1)
1963  // {
1964  // // Take as master the shared point with the lowest mesh
1965  // // point label. (rather arbitrarily - could use max or
1966  // // any other one of the points)
1967  //
1968  // nMergeSets++;
1969  //
1970  // label masterI = labelMax;
1971  //
1972  // forAll(mergeSet, i)
1973  // {
1974  // label sharedI = mergeSet[i];
1975  //
1976  // masterI = min(masterI, sharedPointLabels[sharedI]);
1977  // }
1978  //
1979  // forAll(mergeSet, i)
1980  // {
1981  // label sharedI = mergeSet[i];
1982  //
1983  // pointToMaster.insert(sharedPointLabels[sharedI], masterI);
1984  // }
1985  // }
1986  // }
1987  //
1988  // //if (debug)
1989  // //{
1990  // // Pout<< "polyMeshAdder : merging:"
1991  // // << pointToMaster.size() << " into " << nMergeSets
1992  // // << " sets." << endl;
1993  // //}
1994  //}
1995 
1996  return pointToMaster;
1997 }
1998 
1999 
2002  const polyMesh& mesh,
2003  const Map<label>& pointToMaster,
2004  polyTopoChange& meshMod
2005 )
2006 {
2007  // Remove all non-master points.
2008  forAll(mesh.points(), pointI)
2009  {
2010  Map<label>::const_iterator iter = pointToMaster.find(pointI);
2011 
2012  if (iter != pointToMaster.end())
2013  {
2014  if (iter() != pointI)
2015  {
2016  meshMod.removePoint(pointI, iter());
2017  }
2018  }
2019  }
2020 
2021  // Modify faces for points. Note: could use pointFaces here but want to
2022  // avoid addressing calculation.
2023  const faceList& faces = mesh.faces();
2024 
2025  forAll(faces, faceI)
2026  {
2027  const face& f = faces[faceI];
2028 
2029  bool hasMerged = false;
2030 
2031  forAll(f, fp)
2032  {
2033  label pointI = f[fp];
2034 
2035  Map<label>::const_iterator iter = pointToMaster.find(pointI);
2036 
2037  if (iter != pointToMaster.end())
2038  {
2039  if (iter() != pointI)
2040  {
2041  hasMerged = true;
2042  break;
2043  }
2044  }
2045  }
2046 
2047  if (hasMerged)
2048  {
2049  face newF(f);
2050 
2051  forAll(f, fp)
2052  {
2053  label pointI = f[fp];
2054 
2055  Map<label>::const_iterator iter = pointToMaster.find(pointI);
2056 
2057  if (iter != pointToMaster.end())
2058  {
2059  newF[fp] = iter();
2060  }
2061  }
2062 
2063  label patchID = mesh.boundaryMesh().whichPatch(faceI);
2064  label nei = (patchID == -1 ? mesh.faceNeighbour()[faceI] : -1);
2065  label zoneID = mesh.faceZones().whichZone(faceI);
2066  bool zoneFlip = false;
2067 
2068  if (zoneID >= 0)
2069  {
2070  const faceZone& fZone = mesh.faceZones()[zoneID];
2071  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
2072  }
2073 
2074  meshMod.setAction
2075  (
2077  (
2078  newF, // modified face
2079  faceI, // label of face
2080  mesh.faceOwner()[faceI], // owner
2081  nei, // neighbour
2082  false, // face flip
2083  patchID, // patch for face
2084  false, // remove from zone
2085  zoneID, // zone for face
2086  zoneFlip // face flip in zone
2087  )
2088  );
2089  }
2090  }
2091 }
2092 
2093 
2094 // ************************ vim: set sw=4 sts=4 et: ************************ //