meshRefinementTemplates.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-2021 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 
26 #include "meshRefinement.H"
27 #include "fvMesh.H"
28 #include "globalIndex.H"
29 #include "syncTools.H"
30 
31 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
32 
33 template<class T> void Foam::meshRefinement::updateList
34 (
35  const labelList& newToOld,
36  const T& nullValue,
37  List<T>& elems
38 )
39 {
40  List<T> newElems(newToOld.size(), nullValue);
41 
42  forAll(newElems, i)
43  {
44  const label oldi = newToOld[i];
45 
46  if (oldi >= 0)
47  {
48  newElems[i] = elems[oldi];
49  }
50  }
51 
52  elems.transfer(newElems);
53 }
54 
55 
56 template<class T>
58 (
59  const PackedBoolList& isMasterElem,
60  const UList<T>& values
61 )
62 {
63  if (values.size() != isMasterElem.size())
64  {
66  << "Number of elements in list " << values.size()
67  << " does not correspond to number of elements in isMasterElem "
68  << isMasterElem.size()
69  << exit(FatalError);
70  }
71 
72  T sum = Zero;
73  label n = 0;
74 
75  forAll(values, i)
76  {
77  if (isMasterElem[i])
78  {
79  sum += values[i];
80  n++;
81  }
82  }
83 
84  reduce(sum, sumOp<T>());
85  reduce(n, sumOp<label>());
86 
87  if (n > 0)
88  {
89  return sum/n;
90  }
91  else
92  {
93  return pTraits<T>::max;
94  }
95 }
96 
97 
98 template<class T>
100 (
101  const scalar tol,
102  const string& msg,
103  const UList<T>& faceData,
104  const UList<T>& syncedFaceData
105 ) const
106 {
107  const label nBFaces = mesh_.nFaces() - mesh_.nInternalFaces();
108 
109  if (faceData.size() != nBFaces || syncedFaceData.size() != nBFaces)
110  {
112  << "Boundary faces:" << nBFaces
113  << " faceData:" << faceData.size()
114  << " syncedFaceData:" << syncedFaceData.size()
115  << abort(FatalError);
116  }
117 
118  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
119 
120  forAll(patches, patchi)
121  {
122  const polyPatch& pp = patches[patchi];
123 
124  label bFacei = pp.start() - mesh_.nInternalFaces();
125 
126  forAll(pp, i)
127  {
128  const T& data = faceData[bFacei];
129  const T& syncData = syncedFaceData[bFacei];
130 
131  if (mag(data - syncData) > tol)
132  {
133  const label facei = pp.start() + i;
134 
136  << msg
137  << "patchFace:" << i
138  << " face:" << facei
139  << " fc:" << mesh_.faceCentres()[facei]
140  << " patch:" << pp.name()
141  << " faceData:" << data
142  << " syncedFaceData:" << syncData
143  << " diff:" << mag(data - syncData)
144  << abort(FatalError);
145  }
146 
147  bFacei++;
148  }
149  }
150 }
151 
152 
153 // Print list sorted by coordinates. Used for comparing non-parallel v.s.
154 // parallel operation
155 template<class T>
157 (
158  const UList<point>& points,
159  const UList<T>& data
160 )
161 {
162  globalIndex globalPoints(points.size());
163 
165  globalPoints.gather
166  (
167  Pstream::worldComm,
168  identity(Pstream::nProcs()),
169  points,
170  allPoints,
171  UPstream::msgType(),
172  Pstream::commsTypes::blocking
173  );
174 
175  List<T> allData;
176  globalPoints.gather
177  (
178  Pstream::worldComm,
179  identity(Pstream::nProcs()),
180  data,
181  allData,
182  UPstream::msgType(),
183  Pstream::commsTypes::blocking
184  );
185 
186 
187  scalarField magAllPoints(mag(allPoints-point(-0.317, 0.117, 0.501)));
188 
189  labelList visitOrder;
190  sortedOrder(magAllPoints, visitOrder);
191  forAll(visitOrder, i)
192  {
193  label allPointi = visitOrder[i];
194  Info<< allPoints[allPointi] << " : " << allData[allPointi]
195  << endl;
196  }
197 }
198 
199 
200 template<class Enum>
202 (
203  const Enum& namedEnum,
204  const wordList& words
205 )
206 {
207  int flags = 0;
208 
209  forAll(words, i)
210  {
211  const int index = namedEnum[words[i]];
212  const int val = 1<<index;
213  flags |= val;
214  }
215  return flags;
216 }
217 
218 
219 template<class Type>
221 (
222  const polyMesh& mesh,
223  const PackedBoolList& isMasterEdge,
224  const labelList& meshPoints,
225  const edgeList& edges,
226  const scalarField& edgeWeights,
227  const Field<Type>& pointData,
228  Field<Type>& sum
229 )
230 {
231  if
232  (
233  edges.size() != isMasterEdge.size()
234  || edges.size() != edgeWeights.size()
235  || meshPoints.size() != pointData.size()
236  )
237  {
239  << "Inconsistent sizes for edge or point data:"
240  << " isMasterEdge:" << isMasterEdge.size()
241  << " edgeWeights:" << edgeWeights.size()
242  << " edges:" << edges.size()
243  << " pointData:" << pointData.size()
244  << " meshPoints:" << meshPoints.size()
245  << abort(FatalError);
246  }
247 
248  sum.setSize(meshPoints.size());
249  sum = Zero;
250 
251  forAll(edges, edgei)
252  {
253  if (isMasterEdge[edgei])
254  {
255  const edge& e = edges[edgei];
256 
257  const scalar eWeight = edgeWeights[edgei];
258 
259  const label v0 = e[0];
260  const label v1 = e[1];
261 
262  sum[v0] += eWeight*pointData[v1];
263  sum[v1] += eWeight*pointData[v0];
264  }
265  }
266 
267  syncTools::syncPointList
268  (
269  mesh,
270  meshPoints,
271  sum,
272  plusEqOp<Type>(),
273  Type(Zero) // null value
274  );
275 }
276 
277 
278 // ************************************************************************* //
const fvPatchList & patches
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const word & name() const
Return name.
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Traits class for primitives.
Definition: pTraits.H:50
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
static void weightedSum(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, const scalarField &edgeWeights, const Field< Type > &data, Field< Type > &sum)
Helper: weighted sum (over all subset of mesh points) by.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
void testSyncBoundaryFaceList(const scalar mergeDistance, const string &, const UList< T > &, const UList< T > &) const
Compare two lists over all boundary faces.
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Foam::polyBoundaryMesh.
static int readFlags(const Enum &namedEnum, const wordList &)
Helper: convert wordList into bit pattern using provided.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:100
void setSize(const label)
Reset size of List.
Definition: List.C:281
static void collectAndPrint(const UList< point > &points, const UList< T > &data)
Print list according to (collected and) sorted coordinate.
const volScalarField & T
A bit-packed bool list.
label patchi
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
vector point
Point is a vector.
Definition: point.H:41
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
label size() const
Number of entries.
Definition: PackedListI.H:711
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105