springRenumber.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-2018 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 "springRenumber.H"
28 #include "decompositionMethod.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(springRenumber, 0);
35 
37  (
38  renumberMethod,
39  springRenumber,
40  dictionary
41  );
42 }
43 
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
47 Foam::springRenumber::springRenumber(const dictionary& renumberDict)
48 :
49  renumberMethod(renumberDict),
50  dict_(renumberDict.optionalSubDict(typeName+"Coeffs")),
51  maxCo_(readScalar(dict_.lookup("maxCo"))),
52  maxIter_(readLabel(dict_.lookup("maxIter"))),
53  freezeFraction_(readScalar(dict_.lookup("freezeFraction")))
54 {}
55 
56 
57 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
58 
60 (
61  const polyMesh& mesh,
62  const pointField& points
63 ) const
64 {
65  CompactListList<label> cellCells;
67  (
68  mesh,
69  identity(mesh.nCells()),
70  mesh.nCells(),
71  false, // local only
72  cellCells
73  );
74 
75  return renumber(cellCells(), points);
76 }
77 
78 
80 (
81  const labelListList& cellCells,
82  const pointField& points
83 ) const
84 {
85  // Look at cell index as a 1D position parameter.
86  // Move cells to the average 'position' of their neighbour.
87 
88  scalarField position(cellCells.size());
89  forAll(position, celli)
90  {
91  position[celli] = celli;
92  }
93 
94  labelList oldToNew(identity(cellCells.size()));
95 
96  scalar maxCo = maxCo_ * cellCells.size();
97 
98  for (label iter = 0; iter < maxIter_; iter++)
99  {
100  // Pout<< "Iteration : " << iter << nl
101  // << "------------"
102  // << endl;
103 
104  // Pout<< "Position :" << nl
105  // << " min : " << min(position) << nl
106  // << " max : " << max(position) << nl
107  // << " avg : " << average(position) << nl
108  // << endl;
109 
110  // Sum force per cell.
111  scalarField sumForce(cellCells.size(), 0.0);
112  forAll(cellCells, oldCelli)
113  {
114  const labelList& cCells = cellCells[oldCelli];
115  label celli = oldToNew[oldCelli];
116 
117  forAll(cCells, i)
118  {
119  label nbrCelli = oldToNew[cCells[i]];
120 
121  sumForce[celli] += (position[nbrCelli]-position[celli]);
122  }
123  }
124 
125  // Pout<< "Force :" << nl
126  // << " min : " << min(sumForce) << nl
127  // << " max : " << max(sumForce) << nl
128  // << " avgMag : " << average(mag(sumForce)) << nl
129  // << "DeltaT : " << deltaT << nl
130  // << endl;
131 
132  // Limit displacement
133  scalar deltaT = maxCo / max(mag(sumForce));
134 
135  Info<< "Iter:" << iter
136  << " maxCo:" << maxCo
137  << " deltaT:" << deltaT
138  << " average force:" << average(mag(sumForce)) << endl;
139 
140  // Determine displacement.
141  scalarField displacement(deltaT*sumForce);
142 
143  // Pout<< "Displacement :" << nl
144  // << " min : " << min(displacement) << nl
145  // << " max : " << max(displacement) << nl
146  // << " avgMag : " << average(mag(displacement)) << nl
147  // << endl;
148 
149  // Calculate new position and scale to be within original range
150  // (0..nCells-1) for ease of postprocessing.
151  position += displacement;
152  position -= min(position);
153  position *= (position.size()-1)/max(position);
154 
155  // Slowly freeze.
156  maxCo *= freezeFraction_;
157  }
158 
159  // writeOBJ("endPosition.obj", cellCells, position);
160 
161  // Move cells to new position
163  sortedOrder(position, shuffle);
164 
165  // Reorder oldToNew
166  inplaceReorder(shuffle, oldToNew);
167 
168  return invert(oldToNew.size(), oldToNew);
169 }
170 
171 
172 // ************************************************************************* //
void shuffle(UList< T > &)
Definition: UList.C:143
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual labelList renumber(const pointField &) const
Return the order in which cells need to be visited, i.e.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Macros for easy insertion into run-time selection tables.
static void calcCellCells(const polyMesh &mesh, const labelList &agglom, const label nLocalCoarse, const bool global, CompactListList< label > &cellCells)
Helper: determine (local or global) cellCells from mesh.
Abstract base class for renumbering.
stressControl lookup("compactNormalStress") >> compactNormalStress
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
label readLabel(Istream &is)
Definition: label.H:64
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
scalar maxCo
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Namespace for OpenFOAM.