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-2019 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 // Add a T entry
34 template<class T> void Foam::meshRefinement::updateList
35 (
36  const labelList& newToOld,
37  const T& nullValue,
38  List<T>& elems
39 )
40 {
41  List<T> newElems(newToOld.size(), nullValue);
42 
43  forAll(newElems, i)
44  {
45  label oldI = newToOld[i];
46 
47  if (oldI >= 0)
48  {
49  newElems[i] = elems[oldI];
50  }
51  }
52 
53  elems.transfer(newElems);
54 }
55 
56 
57 template<class T>
59 (
60  const PackedBoolList& isMasterElem,
61  const UList<T>& values
62 )
63 {
64  if (values.size() != isMasterElem.size())
65  {
67  << "Number of elements in list " << values.size()
68  << " does not correspond to number of elements in isMasterElem "
69  << isMasterElem.size()
70  << exit(FatalError);
71  }
72 
73  T sum = Zero;
74  label n = 0;
75 
76  forAll(values, i)
77  {
78  if (isMasterElem[i])
79  {
80  sum += values[i];
81  n++;
82  }
83  }
84 
85  reduce(sum, sumOp<T>());
86  reduce(n, sumOp<label>());
87 
88  if (n > 0)
89  {
90  return sum/n;
91  }
92  else
93  {
94  return pTraits<T>::max;
95  }
96 }
97 
98 
99 // Compare two lists over all boundary faces
100 template<class T>
102 (
103  const scalar tol,
104  const string& msg,
105  const UList<T>& faceData,
106  const UList<T>& syncedFaceData
107 ) const
108 {
109  label nBFaces = mesh_.nFaces() - mesh_.nInternalFaces();
110 
111  if (faceData.size() != nBFaces || syncedFaceData.size() != nBFaces)
112  {
114  << "Boundary faces:" << nBFaces
115  << " faceData:" << faceData.size()
116  << " syncedFaceData:" << syncedFaceData.size()
117  << abort(FatalError);
118  }
119 
120  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
121 
122  forAll(patches, patchi)
123  {
124  const polyPatch& pp = patches[patchi];
125 
126  label bFacei = pp.start() - mesh_.nInternalFaces();
127 
128  forAll(pp, i)
129  {
130  const T& data = faceData[bFacei];
131  const T& syncData = syncedFaceData[bFacei];
132 
133  if (mag(data - syncData) > tol)
134  {
135  label facei = pp.start()+i;
136 
138  << msg
139  << "patchFace:" << i
140  << " face:" << facei
141  << " fc:" << mesh_.faceCentres()[facei]
142  << " patch:" << pp.name()
143  << " faceData:" << data
144  << " syncedFaceData:" << syncData
145  << " diff:" << mag(data - syncData)
146  << abort(FatalError);
147  }
148 
149  bFacei++;
150  }
151  }
152 }
153 
154 
155 // Print list sorted by coordinates. Used for comparing non-parallel v.s.
156 // parallel operation
157 template<class T>
159 (
160  const UList<point>& points,
161  const UList<T>& data
162 )
163 {
164  globalIndex globalPoints(points.size());
165 
167  globalPoints.gather
168  (
169  Pstream::worldComm,
170  identity(Pstream::nProcs()),
171  points,
172  allPoints,
173  UPstream::msgType(),
174  Pstream::commsTypes::blocking
175  );
176 
177  List<T> allData;
178  globalPoints.gather
179  (
180  Pstream::worldComm,
181  identity(Pstream::nProcs()),
182  data,
183  allData,
184  UPstream::msgType(),
185  Pstream::commsTypes::blocking
186  );
187 
188 
189  scalarField magAllPoints(mag(allPoints-point(-0.317, 0.117, 0.501)));
190 
191  labelList visitOrder;
192  sortedOrder(magAllPoints, visitOrder);
193  forAll(visitOrder, i)
194  {
195  label allPointi = visitOrder[i];
196  Info<< allPoints[allPointi] << " : " << allData[allPointi]
197  << endl;
198  }
199 }
200 
201 
202 template<class Enum>
204 (
205  const Enum& namedEnum,
206  const wordList& words
207 )
208 {
209  int flags = 0;
210 
211  forAll(words, i)
212  {
213  int index = namedEnum[words[i]];
214  int val = 1<<index;
215  flags |= val;
216  }
217  return flags;
218 }
219 
220 
221 template<class Type>
223 (
224  const polyMesh& mesh,
225  const PackedBoolList& isMasterEdge,
226  const labelList& meshPoints,
227  const edgeList& edges,
228  const scalarField& edgeWeights,
229  const Field<Type>& pointData,
230  Field<Type>& sum
231 )
232 {
233  if
234  (
235  edges.size() != isMasterEdge.size()
236  || edges.size() != edgeWeights.size()
237  || meshPoints.size() != pointData.size()
238  )
239  {
241  << "Inconsistent sizes for edge or point data:"
242  << " isMasterEdge:" << isMasterEdge.size()
243  << " edgeWeights:" << edgeWeights.size()
244  << " edges:" << edges.size()
245  << " pointData:" << pointData.size()
246  << " meshPoints:" << meshPoints.size()
247  << abort(FatalError);
248  }
249 
250  sum.setSize(meshPoints.size());
251  sum = Zero;
252 
253  forAll(edges, edgeI)
254  {
255  if (isMasterEdge[edgeI])
256  {
257  const edge& e = edges[edgeI];
258 
259  scalar eWeight = edgeWeights[edgeI];
260 
261  label v0 = e[0];
262  label v1 = e[1];
263 
264  sum[v0] += eWeight*pointData[v1];
265  sum[v1] += eWeight*pointData[v0];
266  }
267  }
268 
269  syncTools::syncPointList
270  (
271  mesh,
272  meshPoints,
273  sum,
274  plusEqOp<Type>(),
275  Type(Zero) // null value
276  );
277 }
278 
279 
280 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
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:319
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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
patches[0]
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
volScalarField & e
Elementary charge.
Definition: createFields.H:11
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.
const volScalarField & T
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.
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:303
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:74
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:299
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342