Loading...
Searching...
No Matches
meshRefinementTemplates.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2015 OpenFOAM Foundation
9 Copyright (C) 2019-2023 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "meshRefinement.H"
30#include "fvMesh.H"
31#include "globalIndex.H"
32#include "syncTools.H"
33
34// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
35
36template<class T>
38(
39 const labelList& newToOld,
40 const T& nullValue,
41 List<T>& elems
42)
43{
44 List<T> newElems(newToOld.size(), nullValue);
45
46 forAll(newElems, i)
47 {
48 const label oldI = newToOld[i];
49
50 if (oldI >= 0)
51 {
52 newElems[i] = elems[oldI];
53 }
54 }
55
56 elems.transfer(newElems);
57}
58
59
60template<class T>
62(
63 const bitSet& isMasterElem,
64 const UList<T>& values
65)
66{
67 if (values.size() != isMasterElem.size())
68 {
70 << "Number of elements in list " << values.size()
71 << " does not correspond to number of elements in isMasterElem "
72 << isMasterElem.size()
73 << exit(FatalError);
74 }
75
76 T sum = T(Zero);
77 label n = 0;
78
79 forAll(values, i)
80 {
81 if (isMasterElem.test(i))
82 {
83 sum += values[i];
84 n++;
85 }
86 }
87
90
91 if (n > 0)
92 {
93 return sum/n;
94 }
95 else
96 {
98 }
99}
100
101
102// Compare two lists over all boundary faces
103template<class T>
105(
106 const scalar tol,
107 const string& msg,
108 const UList<T>& faceData,
109 const UList<T>& syncedFaceData
110) const
111{
112 const label nBFaces = mesh_.nBoundaryFaces();
113
114 if (faceData.size() != nBFaces || syncedFaceData.size() != nBFaces)
115 {
117 << "Boundary faces:" << nBFaces
118 << " faceData:" << faceData.size()
119 << " syncedFaceData:" << syncedFaceData.size()
120 << abort(FatalError);
121 }
122
123 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
124
125 forAll(patches, patchi)
126 {
127 const polyPatch& pp = patches[patchi];
128
129 label bFacei = pp.start() - mesh_.nInternalFaces();
130
131 forAll(pp, i)
132 {
133 const T& data = faceData[bFacei];
134 const T& syncData = syncedFaceData[bFacei];
135
136 if (mag(data - syncData) > tol)
137 {
138 label facei = pp.start()+i;
139
141 << msg
142 << "patchFace:" << i
143 << " face:" << facei
144 << " fc:" << mesh_.faceCentres()[facei]
145 << " patch:" << pp.name()
146 << " faceData:" << data
147 << " syncedFaceData:" << syncData
148 << " diff:" << mag(data - syncData)
149 << abort(FatalError);
150 }
151
152 bFacei++;
153 }
155}
156
157
158// Print list sorted by coordinates. Used for comparing non-parallel v.s.
159// parallel operation
160template<class T>
162(
163 const UList<point>& points,
164 const UList<T>& data
165)
166{
167 const globalIndex globalPoints(points.size());
168
169 pointField allPoints(globalPoints.gather(points));
170 List<T> allData(globalPoints.gather(data));
171
172
173 scalarField magAllPoints(mag(allPoints-point(-0.317, 0.117, 0.501)));
174
175 labelList visitOrder(sortedOrder(magAllPoints));
176 for (const label allPointi : visitOrder)
177 {
178 Info<< allPoints[allPointi] << " : " << allData[allPointi]
179 << endl;
180 }
181}
182
183
184template<class GeoField>
185void Foam::meshRefinement::addPatchFields
186(
187 fvMesh& mesh,
188 const word& patchFieldType
189)
190{
191 // Field order is irrelevant...
192 for (GeoField& fld : mesh.objectRegistry::objects<GeoField>())
193 {
194 auto& bfld = fld.boundaryFieldRef();
195
196 const label newPatchi = bfld.size();
197 bfld.resize(newPatchi+1);
198
199 bfld.set
200 (
201 newPatchi,
202 GeoField::Patch::New
203 (
204 patchFieldType,
205 mesh.boundary()[newPatchi],
206 fld.internalField()
207 )
208 );
209 }
210}
211
212
213template<class GeoField>
214void Foam::meshRefinement::reorderPatchFields
215(
216 fvMesh& mesh,
217 const labelList& oldToNew
218)
219{
220 // Field order is irrelevant...
221 for (GeoField& fld : mesh.objectRegistry::objects<GeoField>())
223 fld.boundaryFieldRef().reorder(oldToNew);
224 }
225}
226
227
228template<class EnumContainer>
230(
231 const EnumContainer& namedEnum,
232 const wordList& words
233)
234{
235 int flags = 0;
236
237 for (const word& w : words)
238 {
239 flags |= namedEnum[w];
241
242 return flags;
243}
244
245
246template<class Type>
248(
249 const polyMesh& mesh,
250 const bitSet& isMasterEdge,
251 const labelList& meshPoints,
252 const edgeList& edges,
253 const scalarField& edgeWeights,
254 const Field<Type>& pointData,
255 Field<Type>& sum
256)
257{
258 if
259 (
260 edges.size() != isMasterEdge.size()
261 || edges.size() != edgeWeights.size()
262 || meshPoints.size() != pointData.size()
263 )
264 {
266 << "Inconsistent sizes for edge or point data:"
267 << " isMasterEdge:" << isMasterEdge.size()
268 << " edgeWeights:" << edgeWeights.size()
269 << " edges:" << edges.size()
270 << " pointData:" << pointData.size()
271 << " meshPoints:" << meshPoints.size()
272 << abort(FatalError);
273 }
274
275 sum.setSize(meshPoints.size());
276 sum = Type(Zero);
277
278 forAll(edges, edgeI)
279 {
280 if (isMasterEdge.test(edgeI))
281 {
282 const edge& e = edges[edgeI];
283
284 scalar eWeight = edgeWeights[edgeI];
285
286 label v0 = e[0];
287 label v1 = e[1];
288
289 sum[v0] += eWeight*pointData[v1];
290 sum[v1] += eWeight*pointData[v0];
291 }
292 }
293
295 (
296 mesh,
297 meshPoints,
298 sum,
300 Type(Zero) // null value
301 );
302}
303
304
305template<class Type>
307(
308 const dictionary& dict,
309 const word& keyword,
310 const bool noExit,
311 enum keyType::option matchOpt,
312 const Type& deflt
313)
314{
315 // If noExit, emit FatalIOError message without exiting
317 (
318 noExit
320 );
321
322 Type val(deflt);
323
324 if (!dict.readEntry(keyword, val, matchOpt, readOpt))
325 {
327 << "Entry '" << keyword << "' not found in dictionary "
328 << dict.name() << nl;
329 }
330
331 return val;
332}
333
334
335// ************************************************************************* //
label n
Info<< nl;Info<< "Write faMesh in vtk format:"<< nl;{ vtk::uindirectPatchWriter writer(aMesh.patch(), fileName(aMesh.time().globalPath()/vtkBaseFileName));writer.writeGeometry();globalIndex procAddr(aMesh.nFaces());labelList cellIDs;if(UPstream::master()) { cellIDs.resize(procAddr.totalSize());for(const labelRange &range :procAddr.ranges()) { auto slice=cellIDs.slice(range);slice=identity(range);} } writer.beginCellData(4);writer.writeProcIDs();writer.write("cellID", cellIDs);writer.write("area", aMesh.S().field());writer.write("normal", aMesh.faceAreaNormals());writer.beginPointData(1);writer.write("normal", aMesh.pointAreaNormals());Info<< " "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.time().globalPath()/(vtkBaseFileName+"-edges")));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
Generic templated field type that is much like a Foam::List except that it is expected to hold numeri...
Definition Field.H:172
readOption
Enumeration defining read preferences.
@ READ_IF_PRESENT
Reading is optional [identical to LAZY_READ].
@ MUST_READ
Reading required.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition List.H:72
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
Definition List.C:347
label size() const noexcept
Number of entries.
Definition PackedList.H:392
const Field< point_type > & faceCentres() const
Return face centres for patch.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition UList.H:89
void size(const label n)
Older name for setAddressableSize.
Definition UList.H:118
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition bitSet.H:61
bool test(const label pos) const
Test for True value at specified position, never auto-vivify entries.
Definition bitSet.H:334
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition dictionary.H:133
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
Definition edge.H:62
Mesh data needed to do the Finite Volume discretisation.
Definition fvMesh.H:85
const fvBoundaryMesh & boundary() const noexcept
Return reference to boundary mesh.
Definition fvMesh.H:395
Calculates a non-overlapping list of offsets based on an input size (eg, number of cells) from differ...
Definition globalIndex.H:77
Calculates points shared by more than two processor patches or cyclic patches.
option
Enumeration for the data type and search/match modes (bitmask).
Definition keyType.H:80
static T gAverage(const bitSet &isMasterElem, const UList< T > &values)
Helper: calculate average.
static void collectAndPrint(const UList< point > &points, const UList< T > &data)
Print list according to (collected and) sorted coordinate.
static Type get(const dictionary &dict, const word &keyword, const bool noExit, enum keyType::option matchOpt=keyType::REGEX, const Type &deflt=Zero)
Wrapper around dictionary::get which does not exit.
static int readFlags(const EnumContainer &namedEnum, const wordList &words)
Helper: convert wordList into bit pattern using provided Enum.
static void weightedSum(const polyMesh &mesh, const bitSet &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.
const fvMesh & mesh() const
Reference to mesh.
void testSyncBoundaryFaceList(const scalar mergeDistance, const string &, const UList< T > &, const UList< T > &) const
Compare two lists over all boundary faces.
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
A traits class, which is primarily used for primitives and vector-space.
Definition pTraits.H:64
A polyBoundaryMesh is a polyPatch list with registered IO, a reference to the associated polyMesh,...
Mesh consisting of general polyhedral cells.
Definition polyMesh.H:79
A patch is a list of labels that address the faces in the global face list.
Definition polyPatch.H:73
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
A class for handling words, derived from Foam::string.
Definition word.H:66
const volScalarField & T
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition error.H:629
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
const pointField & points
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
Definition HashOps.H:164
List< edge > edgeList
List of edge.
Definition edgeList.H:32
List< word > wordList
List of word.
Definition fileName.H:60
List< label > labelList
A List of labels.
Definition List.H:62
messageStream Info
Information stream (stdout output on master, null elsewhere).
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition Ostream.H:519
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
void reduce(T &value, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce).
errorManip< error > abort(error &err)
Definition errorManip.H:139
vector point
Point is a vector.
Definition point.H:37
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1, const label comm)
static constexpr const zero Zero
Global zero (0).
Definition zero.H:127
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
vectorField pointField
pointField is a vectorField.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition errorManip.H:125
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
constexpr char nl
The newline '\n' character (0x0a).
Definition Ostream.H:50
dictionary dict
volScalarField & e
#define forAll(list, i)
Loop across all elements in list.
Definition stdFoam.H:299