springRenumber.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 | Copyright (C) 2011-2012 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.subDict(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 // ************************************************************************* //
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
virtual labelList renumber(const pointField &) const
Return the order in which cells need to be visited, i.e.
dimensioned< scalar > mag(const dimensioned< Type > &)
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
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 size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
void shuffle(UList< T > &)
Definition: UList.C:135
label nCells() const
messageStream Info
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Namespace for OpenFOAM.
label readLabel(Istream &is)
Definition: label.H:64
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
Macros for easy insertion into run-time selection tables.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
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.
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Abstract base class for renumbering.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
scalar maxCo
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)