perfectInterface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  Best thing is probably to look at attachDetach which does almost exactly
26  the same but for the geometric matching of points and face centres.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "perfectInterface.H"
31 #include "polyTopoChanger.H"
32 #include "polyMesh.H"
33 #include "polyTopoChange.H"
35 #include "polyTopoChangeMap.H"
36 #include "matchPoints.H"
37 #include "polyModifyFace.H"
38 #include "polyRemovePoint.H"
39 #include "polyRemoveFace.H"
40 #include "indirectPrimitivePatch.H"
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 namespace Foam
45 {
46  defineTypeNameAndDebug(perfectInterface, 0);
48  (
49  polyMeshModifier,
50  perfectInterface,
51  dictionary
52  );
53 }
54 
55 
56 // Tolerance used as fraction of minimum edge length.
57 const Foam::scalar Foam::perfectInterface::tol_ = 1e-3;
58 
59 
60 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 
62 Foam::pointField Foam::perfectInterface::calcFaceCentres
63 (
64  const indirectPrimitivePatch& pp
65 )
66 {
67  const pointField& points = pp.points();
68 
69  pointField ctrs(pp.size());
70 
71  forAll(ctrs, patchFacei)
72  {
73  ctrs[patchFacei] = pp[patchFacei].centre(points);
74  }
75 
76  return ctrs;
77 }
78 
79 
80 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
81 
82 // Construct from components
84 (
85  const word& name,
86  const label index,
87  const polyTopoChanger& mme,
88  const word& faceZoneName,
89  const word& masterPatchName,
90  const word& slavePatchName
91 )
92 :
93  polyMeshModifier(name, index, mme, true),
94  faceZoneID_(faceZoneName, mme.mesh().faceZones()),
95  masterPatchID_(masterPatchName, mme.mesh().boundaryMesh()),
96  slavePatchID_(slavePatchName, mme.mesh().boundaryMesh())
97 {}
98 
99 
100 // Construct from dictionary
102 (
103  const word& name,
104  const dictionary& dict,
105  const label index,
106  const polyTopoChanger& mme
107 )
108 :
109  polyMeshModifier(name, index, mme, readBool(dict.lookup("active"))),
110  faceZoneID_
111  (
112  dict.lookup("faceZoneName"),
113  mme.mesh().faceZones()
114  ),
115  masterPatchID_
116  (
117  dict.lookup("masterPatchName"),
118  mme.mesh().boundaryMesh()
119  ),
120  slavePatchID_
121  (
122  dict.lookup("slavePatchName"),
123  mme.mesh().boundaryMesh()
124  )
125 {}
126 
127 
128 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
129 
131 {}
132 
133 
134 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
135 
137 {
138  // If modifier is inactive, skip change
139  if (!active())
140  {
141  if (debug)
142  {
143  Pout<< "bool perfectInterface::changeTopology() const "
144  << "for object " << name() << " : "
145  << "Inactive" << endl;
146  }
147 
148  return false;
149  }
150  else
151  {
152  // I want topo change every time step.
153  return true;
154  }
155 }
156 
157 
159 (
160  const indirectPrimitivePatch& pp0,
161  const indirectPrimitivePatch& pp1,
162  polyTopoChange& ref
163 ) const
164 {
165  const polyMesh& mesh = topoChanger().mesh();
166 
167  const polyBoundaryMesh& patches = mesh.boundaryMesh();
168 
169  // Some aliases
170  const edgeList& edges0 = pp0.edges();
171  const pointField& pts0 = pp0.localPoints();
172  const pointField& pts1 = pp1.localPoints();
173  const labelList& meshPts0 = pp0.meshPoints();
174  const labelList& meshPts1 = pp1.meshPoints();
175 
176 
177  // Get local dimension as fraction of minimum edge length
178 
179  scalar minLen = great;
180 
181  forAll(edges0, edgeI)
182  {
183  minLen = min(minLen, edges0[edgeI].mag(pts0));
184  }
185  scalar typDim = tol_*minLen;
186 
187  if (debug)
188  {
189  Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
190  << " pts0:" << pts0.size() << " pts1:" << pts1.size()
191  << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
192  }
193 
194 
195  // Determine pointMapping in mesh point labels. Uses geometric
196  // comparison to find correspondence between patch points.
197 
198  labelList renumberPoints(mesh.points().size());
199  forAll(renumberPoints, i)
200  {
201  renumberPoints[i] = i;
202  }
203  {
204  labelList from1To0Points(pts1.size());
205 
206  bool matchOk = matchPoints
207  (
208  pts1,
209  pts0,
210  scalarField(pts1.size(), typDim), // tolerance
211  true, // verbose
212  from1To0Points
213  );
214 
215  if (!matchOk)
216  {
218  << "Points on patch sides do not match to within tolerance "
219  << typDim << exit(FatalError);
220  }
221 
222  forAll(pts1, i)
223  {
224  renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
225  }
226  }
227 
228 
229 
230  // Calculate correspondence between patch faces
231 
232  labelList from0To1Faces(pp1.size());
233 
234  bool matchOk = matchPoints
235  (
236  calcFaceCentres(pp0),
237  calcFaceCentres(pp1),
238  scalarField(pp0.size(), typDim), // tolerance
239  true, // verbose
240  from0To1Faces
241  );
242 
243  if (!matchOk)
244  {
246  << "Face centres of patch sides do not match to within tolerance "
247  << typDim << exit(FatalError);
248  }
249 
250 
251 
252  // Now
253  // - renumber faces using pts1 (except patch1 faces)
254  // - remove patch1 faces. Remember cell label on owner side.
255  // - modify patch0 faces to be internal.
256 
257  // 1. Get faces to be renumbered
258  labelHashSet affectedFaces(2*pp1.size());
259  forAll(meshPts1, i)
260  {
261  label meshPointi = meshPts1[i];
262 
263  if (meshPointi != renumberPoints[meshPointi])
264  {
265  const labelList& pFaces = mesh.pointFaces()[meshPointi];
266 
267  forAll(pFaces, pFacei)
268  {
269  affectedFaces.insert(pFaces[pFacei]);
270  }
271  }
272  }
273  forAll(pp1, i)
274  {
275  affectedFaces.erase(pp1.addressing()[i]);
276  }
277  // Remove patch0 from renumbered faces. Should not be necessary since
278  // patch0 and 1 should not share any point (if created by mergeMeshing)
279  // so affectedFaces should not contain any patch0 faces but you can
280  // never be sure what the user is doing.
281  forAll(pp0, i)
282  {
283  label facei = pp0.addressing()[i];
284 
285  if (affectedFaces.erase(facei))
286  {
288  << "Found face " << facei << " vertices "
289  << mesh.faces()[facei] << " whose points are"
290  << " used both by master patch and slave patch" << endl;
291  }
292  }
293 
294 
295  // 2. Renumber (non patch0/1) faces.
296  forAllConstIter(labelHashSet, affectedFaces, iter)
297  {
298  const label facei = iter.key();
299  const face& f = mesh.faces()[facei];
300 
301  face newFace(f.size());
302 
303  forAll(newFace, fp)
304  {
305  newFace[fp] = renumberPoints[f[fp]];
306  }
307 
308  label nbr = -1;
309 
310  label patchi = -1;
311 
312  if (mesh.isInternalFace(facei))
313  {
314  nbr = mesh.faceNeighbour()[facei];
315  }
316  else
317  {
318  patchi = patches.whichPatch(facei);
319  }
320 
321  label zoneID = mesh.faceZones().whichZone(facei);
322 
323  bool zoneFlip = false;
324 
325  if (zoneID >= 0)
326  {
327  const faceZone& fZone = mesh.faceZones()[zoneID];
328 
329  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
330  }
331 
332  ref.setAction
333  (
335  (
336  newFace, // modified face
337  facei, // label of face being modified
338  mesh.faceOwner()[facei], // owner
339  nbr, // neighbour
340  false, // face flip
341  patchi, // patch for face
342  false, // remove from zone
343  zoneID, // zone for face
344  zoneFlip // face flip in zone
345  )
346  );
347  }
348 
349 
350  // 3. Remove patch1 points
351  forAll(meshPts1, i)
352  {
353  label meshPointi = meshPts1[i];
354 
355  if (meshPointi != renumberPoints[meshPointi])
356  {
357  ref.setAction(polyRemovePoint(meshPointi));
358  }
359  }
360 
361 
362  // 4. Remove patch1 faces
363  forAll(pp1, i)
364  {
365  label facei = pp1.addressing()[i];
366  ref.setAction(polyRemoveFace(facei));
367  }
368 
369 
370  // 5. Modify patch0 faces for new points (not really necessary; see
371  // comment above about patch1 and patch0 never sharing points) and
372  // becoming internal.
373  const boolList& mfFlip =
374  mesh.faceZones()[faceZoneID_.index()].flipMap();
375 
376  forAll(pp0, i)
377  {
378  label facei = pp0.addressing()[i];
379 
380  const face& f = mesh.faces()[facei];
381 
382  face newFace(f.size());
383 
384  forAll(newFace, fp)
385  {
386  newFace[fp] = renumberPoints[f[fp]];
387  }
388 
389  label own = mesh.faceOwner()[facei];
390 
391  label pp1Facei = pp1.addressing()[from0To1Faces[i]];
392 
393  label nbr = mesh.faceOwner()[pp1Facei];
394 
395  if (own < nbr)
396  {
397  ref.setAction
398  (
400  (
401  newFace, // modified face
402  facei, // label of face being modified
403  own, // owner
404  nbr, // neighbour
405  false, // face flip
406  -1, // patch for face
407  false, // remove from zone
408  faceZoneID_.index(), // zone for face
409  mfFlip[i] // face flip in zone
410  )
411  );
412  }
413  else
414  {
415  ref.setAction
416  (
418  (
419  newFace.reverseFace(), // modified face
420  facei, // label of face being modified
421  nbr, // owner
422  own, // neighbour
423  true, // face flip
424  -1, // patch for face
425  false, // remove from zone
426  faceZoneID_.index(), // zone for face
427  !mfFlip[i] // face flip in zone
428  )
429  );
430  }
431  }
432 }
433 
434 
436 {
437  if (debug)
438  {
439  Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
440  << "for object " << name() << " : "
441  << "masterPatchID_:" << masterPatchID_
442  << " slavePatchID_:" << slavePatchID_
443  << " faceZoneID_:" << faceZoneID_ << endl;
444  }
445 
446  if
447  (
448  masterPatchID_.active()
449  && slavePatchID_.active()
450  && faceZoneID_.active()
451  )
452  {
453  const polyMesh& mesh = topoChanger().mesh();
454 
455  const polyBoundaryMesh& patches = mesh.boundaryMesh();
456  const polyPatch& patch0 = patches[masterPatchID_.index()];
457  const polyPatch& patch1 = patches[slavePatchID_.index()];
458 
459 
460  labelList pp0Labels(identity(patch0.size())+patch0.start());
462  (
463  IndirectList<face>(mesh.faces(), pp0Labels),
464  mesh.points()
465  );
466 
467  labelList pp1Labels(identity(patch1.size())+patch1.start());
469  (
470  IndirectList<face>(mesh.faces(), pp1Labels),
471  mesh.points()
472  );
473 
474  setRefinement(pp0, pp1, ref);
475  }
476 }
477 
478 
480 {
481  // Update only my points. Nothing to be done here as points already
482  // shared by now.
483 }
484 
485 
487 {
488  // Mesh has changed topologically. Update local topological data
489  const polyMesh& mesh = topoChanger().mesh();
490 
491  faceZoneID_.update(mesh.faceZones());
492  masterPatchID_.update(mesh.boundaryMesh());
493  slavePatchID_.update(mesh.boundaryMesh());
494 }
495 
496 
498 {
499  os << nl << type() << nl
500  << name()<< nl
501  << faceZoneID_.name() << nl
502  << masterPatchID_.name() << nl
503  << slavePatchID_.name() << endl;
504 }
505 
506 
508 {
509  os << nl << name() << nl << token::BEGIN_BLOCK << nl
510 
511  << " type " << type()
513 
514  << " active " << active()
516 
517  << " faceZoneName " << faceZoneID_.name()
519 
520  << " masterPatchName " << masterPatchID_.name()
522 
523  << " slavePatchName " << slavePatchID_.name()
525 
526  << token::END_BLOCK << endl;
527 }
528 
529 
530 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
531 
532 
533 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
534 
535 
536 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
537 
538 
539 // ************************************************************************* //
const fvPatchList & patches
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
virtual void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
Class containing data for face removal.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool update()
Update the mesh for topology change, mesh to mesh mapping.
Definition: fvMesh.C:584
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Class describing modification of a face.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
virtual void modifyMotionPoints(pointField &motionPoints) const
Modify motion points to comply with the topological change.
virtual void write(Ostream &) const
Write.
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:307
bool readBool(Istream &)
Definition: boolIO.C:60
const Field< PointType > & localPoints() const
Return pointField of points in patch.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
fvMesh & mesh
Determine correspondence between points. See below.
Macros for easy insertion into run-time selection tables.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
A list of faces which address into the list of points.
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:33
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List of mesh modifiers defining the mesh dynamics.
A class for handling words, derived from string.
Definition: word.H:59
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
virtual void setRefinement(polyTopoChange &) const
Insert the layer addition/removal instructions.
Virtual base class for mesh modifiers.
Foam::polyBoundaryMesh.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
virtual ~perfectInterface()
Destructor.
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: MeshZones.C:221
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
bool update(const point &, const FvWallInfoDataBase< WallInfo, Type, Derived > &w2, const scalar tol, TrackingData &td)
Evaluate distance to point. Update distSqr, origin from whomever.
virtual void writeDict(Ostream &) const
Write dictionary.
const labelListList & pointFaces() const
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
const polyMesh & mesh() const
Return the mesh reference.
virtual bool changeTopology() const
Check for topology change.
perfectInterface(const word &name, const label index, const polyTopoChanger &mme, const word &faceZoneName, const word &masterPatchName, const word &slavePatchName)
Construct from components.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Class containing data for point removal.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.