elasticityMotionSolver.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) 2007-2019 PCOpt/NTUA
9 Copyright (C) 2013-2019 FOSS GP
10 Copyright (C) 2019 OpenCFD Ltd.
11-------------------------------------------------------------------------------
12License
13 This file is part of OpenFOAM.
14
15 OpenFOAM is free software: you can redistribute it and/or modify it
16 under the terms of the GNU General Public License as published by
17 the Free Software Foundation, either version 3 of the License, or
18 (at your option) any later version.
19
20 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27
28\*---------------------------------------------------------------------------*/
29
31#include "motionInterpolation.H"
32#include "wallDist.H"
34#include "fvMatrices.H"
35#include "fvcDiv.H"
36#include "fvmDiv.H"
37#include "fvmDiv.H"
38#include "fvmLaplacian.H"
40
41// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42
43namespace Foam
44{
46
48 (
52 );
53}
54
55
56// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
57
59{
60 // Adjust boundary conditions based on the steps to be executed
62 {
63 pointPatchVectorField& pointBCs =
65 if (isA<fixedValuePointPatchVectorField>(pointBCs))
66 {
67 auto& fixedValueBCs =
68 refCast<fixedValuePointPatchVectorField>(pointBCs);
69 fixedValueBCs == fixedValueBCs/scalar(nSteps_);
70 }
71 }
72
73 // Copy boundary conditions to internalField
74 // Needed for interpolation to faces
76
77 // Interpolate boundary conditions from points to faces
79 {
81 if (isA<fixedValueFvPatchVectorField>(bField))
82 {
84 const polyPatch& patch = mesh().boundaryMesh()[pI];
85 forAll(bField, fI)
86 {
87 bField[fI] = patch[fI].average(points, pointMotionU_);
88 }
89 }
90 }
91}
92
93
94// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
95
97(
98 const polyMesh& mesh,
99 const IOdictionary& dict
100)
101:
102 motionSolver(mesh, dict, typeName),
103 fvMesh_
104 (
105 const_cast<fvMesh&>
106 (
107 refCast<const fvMesh>(mesh)
108 )
109 ),
110 pointMotionU_
111 (
113 (
114 "pointMotionU",
115 mesh.time().timeName(),
116 mesh,
117 IOobject::READ_IF_PRESENT,
118 IOobject::AUTO_WRITE
119 ),
122 fixedValuePointPatchVectorField::typeName
123 ),
124 cellMotionU_
125 (
127 (
128 "cellMotionU",
129 mesh.time().timeName(),
130 mesh,
131 IOobject::READ_IF_PRESENT,
132 IOobject::AUTO_WRITE
133 ),
134 fvMesh_,
135 dimensionedVector(pointMotionU_.dimensions(), Zero),
136 pointMotionU_.boundaryField().types()
137 ),
138 interpolationPtr_
139 (
140 coeffDict().found("interpolation")
141 ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
142 : motionInterpolation::New(fvMesh_)
143 ),
144 E_
145 (
147 (
148 "mu",
149 mesh.time().timeName(),
150 mesh,
151 IOobject::NO_READ,
152 IOobject::NO_WRITE
153 ),
154 fvMesh_,
156 zeroGradientFvPatchScalarField::typeName
157 ),
158 exponent_(this->coeffDict().get<scalar>("exponent")),
159 nSteps_(this->coeffDict().get<label>("steps")),
160 nIters_(this->coeffDict().get<label>("iters")),
161 tolerance_(this->coeffDict().get<scalar>("tolerance"))
162{}
163
164
165// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
166
168{
169 tmp<pointField> tnewPoints(new pointField(mesh().points()));
170
171 return tnewPoints;
172}
173
174
176{
177 // Re-init to zero
178 cellMotionU_.primitiveFieldRef() = vector::zero;
179
180 // Adjust boundary conditions based on the number of steps to be executed
181 // and interpolate to faces
182 setBoundaryConditions();
183
184 // Solve the elasticity equations in a stepped manner
185 for (label istep = 0; istep < nSteps_; ++istep)
186 {
187 Info<< "Step " << istep << endl;
188
189 // Update diffusivity
190 const scalarField& vols = mesh().cellVolumes();
191 E_.primitiveFieldRef() = 1./pow(vols, exponent_);
192 E_.correctBoundaryConditions();
193
194 for (label iter = 0; iter < nIters_; ++iter)
195 {
196 Info<< "Iteration " << iter << endl;
197 cellMotionU_.storePrevIter();
198 fvVectorMatrix dEqn
199 (
200 fvm::laplacian(2*E_, cellMotionU_)
201 + fvc::div(2*E_*T(fvc::grad(cellMotionU_)))
202 - fvc::div(E_*fvc::div(cellMotionU_)*tensor::I)
203 );
204
205 scalar residual = mag(dEqn.solve().initialResidual());
206 cellMotionU_.relax();
207
208 // Print execution time
209 fvMesh_.time().printExecutionTime(Info);
210
211 // Check convergence
212 if (residual < tolerance_)
213 {
214 Info<< "\n***Reached mesh movement convergence limit for step "
215 << istep
216 << " iteration " << iter << "***\n\n";
217 break;
218 }
219 }
220
221 // Interpolate from cells to points
222 interpolationPtr_->interpolate(cellMotionU_, pointMotionU_);
223 vectorField newPoints
224 (
225 mesh().points() + pointMotionU_.primitiveFieldRef()
226 );
227
228 // Move points and check mesh
229 fvMesh_.movePoints(newPoints);
230 fvMesh_.checkMesh(true);
231
232 if (debug)
233 {
234 Info<< " Writing new mesh points " << endl;
236 (
238 (
239 "points",
240 mesh().pointsInstance(),
241 mesh().meshSubDir,
242 mesh(),
245 ),
246 mesh().points()
247 );
248 points.write();
249 }
250 }
251}
252
253
255{
256 // Do nothing
257}
258
259
261{
262 // Do nothing
263}
264
265
266// ************************************************************************* //
bool found
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
void updateCoeffs()
Update the boundary condition coefficients.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
static const Tensor I
Definition: Tensor.H:81
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Mesh deformation based on the linear elasticity equations. The boundary displacement is set as a boun...
virtual tmp< pointField > curPoints() const
Return point location. Mesh is actually moved in solve()
label nSteps_
Intermediate steps to solve the PDEs.
void setBoundaryConditions()
Set boundary conditions of cellMotionU based on pointMotionU.
virtual void solve()
Solve for motion.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
Base class for interpolation of cell displacement fields, generated by fvMotionSolvers,...
void movePoints()
Update for new mesh geometry.
void updateMesh()
Update for new mesh topology.
Virtual base class for mesh motion solver.
Definition: motionSolver.H:61
const polyMesh & mesh() const
Return reference to mesh.
Definition: motionSolver.H:144
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:55
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
const scalarField & cellVolumes() const
Lookup type of boundary radiation properties.
Definition: lookup.H:66
A class for managing temporary objects.
Definition: tmp.H:65
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & T
dynamicFvMesh & mesh
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Calculate the divergence of the given field.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
const pointField & points
word timeName
Definition: getTimeIndex.H:3
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
Namespace for OpenFOAM.
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:131
const dimensionSet dimless
Dimensionless.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
A non-counting (dummy) refCount.
Definition: refCount.H:59