patchToPatchStabilisation.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) 2023 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 \*---------------------------------------------------------------------------*/
25 
27 #include "PatchEdgeFacePointData.H"
28 #include "PatchEdgeFaceWave.H"
29 #include "SubField.H"
30 #include "globalIndex.H"
31 #include "OBJstream.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
38 }
39 
40 
41 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
42 
44 :
45  stabilisation_(false),
46  localStabilisationCells_(),
47  stabilisationMapPtr_(nullptr)
48 {}
49 
50 
51 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
52 
54 {}
55 
56 
57 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
58 
60 (
61  const polyPatch& patch,
62  const PackedBoolList& faceCoupleds
63 )
64 {
65  // Determine whether or not stabilisation is necessary
66  stabilisation_ = false;
67  forAll(faceCoupleds, facei)
68  {
69  if (!faceCoupleds[facei])
70  {
71  stabilisation_ = true;
72  break;
73  }
74  }
75  reduce(stabilisation_, orOp<bool>());
76 
77  // Quick return if nothing is to be done
78  if (!stabilisation_) return;
79 
80  // Construct initial edges. All edges that border a coupled face are added
81  // here. The wave will propagate everywhere for just the first iteration.
82  // Then most paths will end and subsequent iterations will propagate only
83  // through the uncoupled faces. This is a bit odd, but it is easier than
84  // doing the necessary synchronisation to determine which edges lie
85  // in-between coupled and non-coupled faces.
86  typedef PatchEdgeFacePointData<remote> info;
87  DynamicList<label> initialEdges(patch.nEdges());
88  DynamicList<info> initialEdgeInfos(patch.nEdges());
89  forAll(patch.edgeFaces(), edgei)
90  {
91  forAll(patch.edgeFaces()[edgei], edgeFacei)
92  {
93  const label facei = patch.edgeFaces()[edgei][edgeFacei];
94 
95  if (faceCoupleds[facei])
96  {
97  initialEdges.append(edgei);
98  initialEdgeInfos.append
99  (
100  info
101  (
102  remote(Pstream::myProcNo(), facei),
103  patch.edges()[edgei].centre(patch.localPoints()),
104  0
105  )
106  );
107  break;
108  }
109  }
110  }
111 
112  // Wave the information about the nearby coupled faces into the un-coupled
113  // faces. Base this wave on distance to the cut face. Initialise coupled
114  // faces to have a distance of zero, so that we do not waste time waving
115  // into coupled regions of the patch.
116  List<info> edgeInfos(patch.nEdges()), faceInfos(patch.size());
117  forAll(faceCoupleds, facei)
118  {
119  if (faceCoupleds[facei])
120  {
121  faceInfos[facei] =
122  info
123  (
124  remote(Pstream::myProcNo(), facei),
125  patch.faceCentres()[facei],
126  0
127  );
128  }
129  }
131  (
132  patch.boundaryMesh().mesh(),
133  patch,
134  initialEdges,
135  initialEdgeInfos,
136  edgeInfos,
137  faceInfos,
138  returnReduce(patch.nEdges(), sumOp<label>())
139  );
140 
141  // Check that the wave connected to all un-mapped faces
142  forAll(faceCoupleds, facei)
143  {
144  if (!faceCoupleds[facei] && !faceInfos[facei].valid(wave.data()))
145  {
147  << "Un-mapped face " << facei << " of patch " << patch.name()
148  << " on processor " << Pstream::myProcNo() << " with centre "
149  << "at " << patch.faceCentres()[facei] << " was not connected "
150  << "to a mapped cell by the stabilisation wave. This "
151  << "indicates that an entire non-contiguous region of patch "
152  << "lies outside of the other patch being mapped to. This is "
153  << "not recoverable." << exit(FatalError);
154  }
155  }
156 
157  // Construct the cell to local stabilisation cell map
158  const globalIndex cellGlobalIndex(patch.size());
159  localStabilisationCells_.resize(patch.size());
160  forAll(faceCoupleds, facei)
161  {
162  const remote& r = faceInfos[facei].data();
163  localStabilisationCells_[facei] =
164  faceCoupleds[facei]
165  ? cellGlobalIndex.toGlobal(facei)
166  : cellGlobalIndex.toGlobal(r.proci, r.elementi);
167  }
168 
169  // Construct the distribution map, if necessary
170  if (Pstream::parRun())
171  {
172  List<Map<label>> compactMap;
173  stabilisationMapPtr_.reset
174  (
175  new distributionMap
176  (
177  cellGlobalIndex,
178  localStabilisationCells_,
179  compactMap
180  )
181  );
182  }
183 
184  // Write out stabilisation connections
185  if (debug)
186  {
187  OBJstream obj
188  (
189  typeName + "_" + patch.name()
190  + (Pstream::parRun() ? "_proc" + name(Pstream::myProcNo()) : "")
191  + "_connections.obj"
192  );
193 
194  const pointField fcs(patch.faceCentres());
195  pointField sfcs(fcs);
196  stabilise(sfcs);
197 
198  forAll(fcs, celli)
199  {
200  const point& c = fcs[celli];
201  if (magSqr(c - fcs[celli]) == 0) continue;
202  obj.write(linePointRef(fcs[celli], c));
203  }
204  }
205 }
206 
207 
208 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
OFstream which keeps track of vertices.
Definition: OBJstream.H:56
virtual Ostream & write(const char)
Write character.
Definition: OBJstream.C:81
A bit-packed bool list.
Transport of nearest point location, plus data, for use in PatchEdgeFaceWave.
Wave propagation of information along patch. Every iteration information goes through one layer of fa...
label nEdges() const
Return number of edges in patch.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelListList & edgeFaces() const
Return edge-face addressing.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
T * data()
Return a pointer to the first data element,.
Definition: UListI.H:149
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Class containing processor-to-processor mapping information.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:64
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
const word & name() const
Return name.
Stabilisation data and routines for patch-to-patch interpolations.
void update(const polyPatch &patch, const PackedBoolList &faceCoupleds)
Compute the stabilisation addressing if necessary.
const polyMesh & mesh() const
Return the mesh reference.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:270
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:276
Struct for keeping processor, element (cell, face, point) index.
Definition: remote.H:57
label elementi
Element index.
Definition: remote.H:70
label proci
Processor index.
Definition: remote.H:67
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
bool valid(const PtrList< ModelType > &l)
const dimensionedScalar c
Speed of light in a vacuum.
label wave(const fvMesh &mesh, const List< labelPair > &changedPatchAndFaces, const label nCorrections, GeometricField< scalar, PatchField, GeoMesh > &distance, TrackingData &td, GeometricField< DataType, PatchField, GeoMesh > &... data)
Wave distance (and maybe additional) data from faces. If nCorrections is.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
defineTypeNameAndDebug(combustionModel, 0)
error FatalError
dimensioned< scalar > magSqr(const dimensioned< Type > &)