regionModel1D.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-2020 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 "regionModel1D.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 namespace regionModels
33 {
34  defineTypeNameAndDebug(regionModel1D, 0);
35 }
36 }
37 
38 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
39 
40 void Foam::regionModels::regionModel1D::constructMeshObjects()
41 {
42  nMagSfPtr_.reset
43  (
45  (
46  IOobject
47  (
48  "nMagSf",
49  time().timeName(),
50  regionMesh(),
53  ),
54  regionMesh(),
56  )
57  );
58 }
59 
60 
61 void Foam::regionModels::regionModel1D::initialise()
62 {
64 
65  // Calculate boundaryFaceFaces and boundaryFaceCells
66 
67  DynamicList<label> faceIDs;
68  DynamicList<label> cellIDs;
69 
70  label localPyrolysisFacei = 0;
71 
72  const polyBoundaryMesh& rbm = regionMesh().boundaryMesh();
73 
74  forAll(intCoupledPatchIDs_, i)
75  {
76  const label patchi = intCoupledPatchIDs_[i];
77  const polyPatch& ppCoupled = rbm[patchi];
78  forAll(ppCoupled, localFacei)
79  {
80  label facei = ppCoupled.start() + localFacei;
81  label celli = -1;
82  label nFaces = 0;
83  label nCells = 0;
84  do
85  {
86  label ownCelli = regionMesh().faceOwner()[facei];
87  if (ownCelli != celli)
88  {
89  celli = ownCelli;
90  }
91  else
92  {
93  celli = regionMesh().faceNeighbour()[facei];
94  }
95  nCells++;
96  cellIDs.append(celli);
97  const cell& cFaces = regionMesh().cells()[celli];
98  facei = cFaces.opposingFaceLabel(facei, regionMesh().faces());
99  faceIDs.append(facei);
100  nFaces++;
101  } while (regionMesh().isInternalFace(facei));
102 
103  boundaryFaceOppositeFace_[localPyrolysisFacei] = facei;
104  faceIDs.remove(); // remove boundary face.
105  nFaces--;
106 
107  boundaryFaceFaces_[localPyrolysisFacei].transfer(faceIDs);
108  boundaryFaceCells_[localPyrolysisFacei].transfer(cellIDs);
109 
110  localPyrolysisFacei++;
111  nLayers_ = nCells;
112  }
113  }
114 
115  boundaryFaceOppositeFace_.setSize(localPyrolysisFacei);
116  boundaryFaceFaces_.setSize(localPyrolysisFacei);
117  boundaryFaceCells_.setSize(localPyrolysisFacei);
118 
119  surfaceScalarField& nMagSf = nMagSfPtr_();
120  surfaceScalarField::Boundary& nMagSfBf = nMagSf.boundaryFieldRef();
121 
122  localPyrolysisFacei = 0;
123 
124  forAll(intCoupledPatchIDs_, i)
125  {
126  const label patchi = intCoupledPatchIDs_[i];
127  const polyPatch& ppCoupled = rbm[patchi];
128  const vectorField& pNormals = ppCoupled.faceNormals();
129  nMagSfBf[patchi] = regionMesh().Sf().boundaryField()[patchi] & pNormals;
130  forAll(pNormals, localFacei)
131  {
132  const vector& n = pNormals[localFacei];
133  const labelList& faces = boundaryFaceFaces_[localPyrolysisFacei++];
134  forAll(faces, facei)
135  {
136  const label faceID = faces[facei];
137  nMagSf[faceID] = regionMesh().Sf()[faceID] & n;
138  }
139  }
140  }
141 }
142 
143 
144 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
145 
147 {
148  if (regionModel::read())
149  {
150  moveMesh_.readIfPresent("moveMesh", coeffs_);
151 
152  return true;
153  }
154  else
155  {
156  return false;
157  }
158 }
159 
160 
162 {
163  if (regionModel::read(dict))
164  {
165  moveMesh_.readIfPresent("moveMesh", coeffs_);
166 
167  return true;
168  }
169  else
170  {
171  return false;
172  }
173 }
174 
175 
177 (
178  const scalarList& deltaV,
179  const scalar minDelta
180 )
181 {
182  tmp<labelField> tcellMoveMap(new labelField(regionMesh().nCells(), 0));
183  labelField& cellMoveMap = tcellMoveMap.ref();
184 
185  if (!moveMesh_)
186  {
187  return cellMoveMap;
188  }
189 
190  pointField oldPoints = regionMesh().points();
191  pointField newPoints = oldPoints;
192 
193  const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
194 
195  label totalFaceId = 0;
196  forAll(intCoupledPatchIDs_, localPatchi)
197  {
198  label patchi = intCoupledPatchIDs_[localPatchi];
199  const polyPatch pp = bm[patchi];
200  const vectorField& cf = regionMesh().Cf().boundaryField()[patchi];
201 
202  forAll(pp, patchFacei)
203  {
204  const labelList& faces = boundaryFaceFaces_[totalFaceId];
205  const labelList& cells = boundaryFaceCells_[totalFaceId];
206 
207  const vector n = pp.faceNormals()[patchFacei];
208  const vector sf = pp.faceAreas()[patchFacei];
209 
210  List<point> oldCf(faces.size() + 1);
211  oldCf[0] = cf[patchFacei];
212  forAll(faces, i)
213  {
214  oldCf[i + 1] = regionMesh().faceCentres()[faces[i]];
215  }
216 
217  vector newDelta = Zero;
218  point nbrCf = oldCf[0];
219 
220  forAll(faces, i)
221  {
222  const label facei = faces[i];
223  const label celli = cells[i];
224 
225  const face f = regionMesh().faces()[facei];
226 
227  newDelta += (deltaV[celli]/mag(sf))*n;
228 
229  vector localDelta = Zero;
230  forAll(f, pti)
231  {
232  const label pointi = f[pti];
233 
234  if
235  (
236  mag((nbrCf - (oldPoints[pointi] + newDelta)) & n)
237  > minDelta
238  )
239  {
240  newPoints[pointi] = oldPoints[pointi] + newDelta;
241  localDelta = newDelta;
242  cellMoveMap[celli] = 1;
243  }
244  }
245  nbrCf = oldCf[i + 1] + localDelta;
246  }
247  // Modify boundary
248  const label bFacei = boundaryFaceOppositeFace_[totalFaceId];
249  const face f = regionMesh().faces()[bFacei];
250  const label celli = cells[cells.size() - 1];
251  newDelta += (deltaV[celli]/mag(sf))*n;
252  forAll(f, pti)
253  {
254  const label pointi = f[pti];
255  if
256  (
257  mag((nbrCf - (oldPoints[pointi] + newDelta)) & n)
258  > minDelta
259  )
260  {
261  newPoints[pointi] = oldPoints[pointi] + newDelta;
262  cellMoveMap[celli] = 1;
263  }
264  }
265  totalFaceId ++;
266  }
267  }
268  // Move points
269  regionMesh().movePoints(newPoints);
270 
271  return tcellMoveMap;
272 }
273 
274 
275 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
276 
277 Foam::regionModels::regionModel1D::regionModel1D
278 (
279  const fvMesh& mesh,
280  const word& regionType
281 )
282 :
283  regionModel(mesh, regionType),
284  boundaryFaceFaces_(),
285  boundaryFaceCells_(),
286  boundaryFaceOppositeFace_(),
287  nLayers_(0),
288  nMagSfPtr_(nullptr),
289  moveMesh_(false)
290 {}
291 
292 
293 Foam::regionModels::regionModel1D::regionModel1D
294 (
295  const fvMesh& mesh,
296  const word& regionType,
297  const word& modelName,
298  bool readFields
299 )
300 :
301  regionModel(mesh, regionType, modelName, false),
302  boundaryFaceFaces_(regionMesh().nCells()),
303  boundaryFaceCells_(regionMesh().nCells()),
304  boundaryFaceOppositeFace_(regionMesh().nCells()),
305  nLayers_(0),
306  nMagSfPtr_(nullptr),
307  moveMesh_(true)
308 {
309  constructMeshObjects();
310  initialise();
311 
312  if (readFields)
313  {
314  read();
315  }
316 }
317 
318 
319 Foam::regionModels::regionModel1D::regionModel1D
320 (
321  const fvMesh& mesh,
322  const word& regionType,
323  const word& modelName,
324  const dictionary& dict,
325  bool readFields
326 )
327 :
328  regionModel(mesh, regionType, modelName, dict, readFields),
329  boundaryFaceFaces_(regionMesh().nCells()),
330  boundaryFaceCells_(regionMesh().nCells()),
331  boundaryFaceOppositeFace_(regionMesh().nCells()),
332  nLayers_(0),
333  nMagSfPtr_(nullptr),
334  moveMesh_(false)
335 {
336  constructMeshObjects();
337  initialise();
338 
339  if (readFields)
340  {
341  read(dict);
342  }
343 }
344 
345 
346 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
347 
349 {}
350 
351 
352 // ************************************************************************* //
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
#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 dimensionSet dimArea
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
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
tmp< labelField > moveMesh(const scalarList &deltaV, const scalar minDelta=0.0)
Move mesh points according to change in cell volumes.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
defineTypeNameAndDebug(regionModel, 0)
bool read(const char *, int32_t &)
Definition: int32IO.C:85
const cellShapeList & cells
virtual ~regionModel1D()
Destructor.
A class for handling words, derived from string.
Definition: word.H:59
virtual bool read()
Read control parameters from dictionary.
Definition: regionModel.C:141
#define DebugInFunction
Report an information message using Foam::Info.
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
const Field< PointType > & faceNormals() const
Return face normals for patch.
Foam::polyBoundaryMesh.
const vectorField::subField faceAreas() const
Return face areas.
Definition: polyPatch.C:303
labelList f(nPoints)
volScalarField sf(fieldObject, mesh)
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual bool read()
Read control parameters from dictionary.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Base class for region models.
Definition: regionModel.H:55
Field< vector > vectorField
Specialisation of Field<T> for vector.
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Namespace for OpenFOAM.