v0.14.0
Searching...
No Matches
ArcLengthTools.hpp
Go to the documentation of this file.
1/** \file ArcLengthTools.cpp
2 * \ingroup arc_length_control
3 *
4 * Implementation of arc-length control method
5 *
6 * FIXME: Some variables not comply with naming convention, need to be fixed.
7 *
8 */
9
10
11
12#ifndef __ARC_LENGTH_TOOLS_HPP__
13#define __ARC_LENGTH_TOOLS_HPP__
14
15/**
16 * \brief Store variables for ArcLength analysis
17 *
18 * \ingroup arc_length_control
19
20 The constrain function if given by
21 \f[
22 r_\lambda = f_\lambda(\mathbf{x},\lambda) - s^2
23 \f]
24 where \f$f_\lambda(\mathbf{x},\lambda)\f$ is some constrain function, which has
25 general form given by
26 \f[
27 f_\lambda(\mathbf{x},\lambda) = \alpha f(\Delta\mathbf{x}) +
28 \beta^2 \Delta\lambda^2 \| \mathbf{F}_{\lambda} \|^2
29 \f]
30 where for example \f$f(\mathbf{x})=\|\Delta\mathbf{x}\|^2\f$ is some user
31 defined
32 function evaluating
33 increments vector of degrees of freedom \f$\Delta\mathbf{x}\f$. The increment
34 vector is
35 \f[
36 \Delta \mathbf{x} = \mathbf{x}-\mathbf{x}_0
37 \f]
38 and
39 \f[
40 \Delta \lambda = \lambda-\lambda_0.
41 \f]
42
43 For convenience we assume that
44 \f[
45 \frac{\partial f}{\partial \mathbf{x}}\Delta \mathbf{x}
46 =
47 \textrm{d}\mathbf{b} \Delta \mathbf{x},
48 \f]
49 as result linearised constrain equation takes form
50 \f[
51 \textrm{d}\mathbf{b} \delta \Delta x +
52 D \delta \Delta\lambda - r_{\lambda} = 0
53 \f]
54 where
55 \f[
56 D = 2\beta^2 \Delta\lambda \| \mathbf{F}_{\lambda} \|^2.
57 \f]
58
59 User need to implement functions calculating \f$f(\mathbf{x},\lambda)\f$, i.e.
60 function
61 \f$f(\|\Delta\mathbf{x}\|^2)\f$ and its derivative,
62 \f$\textrm{d}\mathbf{b}\f$.
63
64 */
66
67 virtual ~ArcLengthCtx() = default;
68
70
71 double s; ///< arc length radius
72 double beta; ///< force scaling factor
73 double alpha; ///< displacement scaling factor
74
75 SmartPetscObj<Vec> ghosTdLambda;
76 double dLambda; ///< increment of load factor
77 SmartPetscObj<Vec> ghostDiag;
78 double dIag; ///< diagonal value
79
80 double dx2; ///< inner_prod(dX,dX)
81 double F_lambda2; ///< inner_prod(F_lambda,F_lambda);
82 double res_lambda; ///< f_lambda - s
83 SmartPetscObj<Vec> F_lambda; ///< F_lambda reference load vector
84 SmartPetscObj<Vec>
85 db; ///< db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
86 SmartPetscObj<Vec> xLambda; ///< solution of eq. K*xLambda = F_lambda
87 SmartPetscObj<Vec> x0; ///< displacement vector at beginning of step
88 SmartPetscObj<Vec> dx; ///< dx = x-x0
89
90 /**
91 * \brief set arc radius
92 */
93 MoFEMErrorCode setS(double s);
94
95 /**
96 * \brief set parameters controlling arc-length equations
97 * alpha controls off diagonal therms
98 * beta controls diagonal therm
99 */
100 MoFEMErrorCode setAlphaBeta(double alpha, double beta);
101
102 ArcLengthCtx(MoFEM::Interface &m_field, const std::string &problem_name,
103 const std::string &field_name = "LAMBDA");
104
105 /** \brief Get global index of load factor
106 */
108 return arcDofRawPtr->getPetscGlobalDofIdx();
109 };
110
111 /** \brief Get local index of load factor
112 */
113 DofIdx getPetscLocalDofIdx() { return arcDofRawPtr->getPetscLocalDofIdx(); };
114
115 /** \brief Get value of load factor
116 */
117 FieldData &getFieldData() { return arcDofRawPtr->getFieldData(); }
118
119 /** \brief Get proc owning lambda dof
120 */
121 int getPart() { return arcDofRawPtr->getPart(); };
122
123private:
124 NumeredDofEntity *arcDofRawPtr;
125};
126
127#ifdef __SNESCTX_HPP__
128
129/**
130 * \brief It is ctx structure passed to SNES solver
131 * \ingroup arc_length_control
132 */
133struct ArcLengthSnesCtx : public SnesCtx {
134 ArcLengthCtx *arcPtrRaw;
135
136 ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
137 ArcLengthCtx *arc_ptr_raw)
138 : SnesCtx(m_field, problem_name), arcPtrRaw(arc_ptr_raw) {}
139
140 ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
141 boost::shared_ptr<ArcLengthCtx> arc_ptr)
142 : SnesCtx(m_field, problem_name), arcPtrRaw(arc_ptr.get()),
143 arcPtr(arc_ptr) {}
144
145private:
146 boost::shared_ptr<ArcLengthCtx> arcPtr;
147};
148
149#endif //__SNESCTX_HPP__
150
151#ifdef __TSCTX_HPP__
152
153/**
154 * \brief It is ctx structure passed to SNES solver
155 * \ingroup arc_length_control
156 */
157struct ArcLengthTsCtx : public TsCtx {
158 ArcLengthCtx *arcPtrRaw;
159
160 /// \deprecated use constructor with shared ptr
161 DEPRECATED ArcLengthTsCtx(MoFEM::Interface &m_field,
162 const std::string &problem_name,
163 ArcLengthCtx *arc_ptr_raw)
164 : TsCtx(m_field, problem_name), arcPtrRaw(arc_ptr_raw) {}
165
166 ArcLengthTsCtx(MoFEM::Interface &m_field, const std::string &problem_name,
167 boost::shared_ptr<ArcLengthCtx> arc_ptr)
168 : TsCtx(m_field, problem_name), arcPtrRaw(arc_ptr.get()),
169 arcPtr(arc_ptr) {}
170
171private:
172 boost::shared_ptr<ArcLengthCtx> arcPtr;
173};
174
175#endif // __TSCTX_HPP__
176
177/** \brief shell matrix for arc-length method
178 *
179 * \ingroup arc_length_control
180
181 Shell matrix which has structure:
182 \f[
183 \left[
184 \begin{array}{cc}
185 \mathbf{K} & -\mathbf{F}_\lambda \\
186 \textrm{d}\mathbf{b} & D
187 \end{array}
188 \right]
189 \left\{
190 \begin{array}{c}
191 \delta \Delta \mathbf{x} \\
192 \delta \Delta \lambda
193 \end{array}
194 \right\}
195 =
196 \left[
197 \begin{array}{c}
198 -\mathbf{f}_\textrm{int} \\
199 -r_\lambda
200 \end{array}
201 \right]
202 \f]
203
204 */
206
207 SmartPetscObj<Mat> Aij;
209 ArcLengthCtx *arcPtrRaw; // this is for back compatibility
210
211 ArcLengthMatShell(Mat aij, boost::shared_ptr<ArcLengthCtx> arc_ptr,
212 string problem_name);
213 virtual ~ArcLengthMatShell() = default;
214
215 /// \deprecated use constructor with shared_ptr
216 DEPRECATED ArcLengthMatShell(Mat aij, ArcLengthCtx *arc_ptr_raw,
217 string problem_name);
218
219 MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode);
220 friend MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f);
221
222private:
223 boost::shared_ptr<ArcLengthCtx> arcPtr;
224};
225
226/**
227 * mult operator for Arc Length Shell Mat
228 */
229MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f);
230
231/**
232 * \brief structure for Arc Length pre-conditioner
233 * \ingroup arc_length_control
234 */
236
237 SmartPetscObj<KSP> kSP;
238 SmartPetscObj<PC> pC;
239 SmartPetscObj<Mat> shellAij;
240 SmartPetscObj<Mat> Aij;
241
242 PCArcLengthCtx(Mat shell_Aij, Mat aij,
243 boost::shared_ptr<ArcLengthCtx> arc_ptr);
244
245 PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
246 boost::shared_ptr<ArcLengthCtx> arc_ptr);
247
248 ArcLengthCtx *arcPtrRaw; // this is for back compatibility
249
250 /// \deprecated use with shared_ptr
251 DEPRECATED PCArcLengthCtx(Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr_raw);
252 /// \deprecated use with shared_ptr
253 DEPRECATED PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
254 ArcLengthCtx *arc_ptr_raw);
255
256 friend MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x);
257 friend MoFEMErrorCode PCSetupArcLength(PC pc);
258
259private:
260 boost::shared_ptr<ArcLengthCtx> arcPtr;
261};
262
263/**
264 * apply operator for Arc Length pre-conditioner
265 * solves K*pc_x = pc_f
266 * solves K*xLambda = -dF_lambda
267 * solves ddlambda = ( res_lambda - db*xLambda )/( diag + db*pc_x )
268 * calculate pc_x = pc_x + ddlambda*xLambda
269 */
270MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x);
271
272/**
273 * set up structure for Arc Length pre-conditioner
274
275 * it sets pre-conditioner for matrix K
276 */
277MoFEMErrorCode PCSetupArcLength(PC pc);
278
279/**
280 * \brief Zero F_lambda
281 *
282 */
283struct ZeroFLmabda : public FEMethod {
284
285 boost::shared_ptr<ArcLengthCtx> arcPtr;
286
287 ZeroFLmabda(boost::shared_ptr<ArcLengthCtx> arc_ptr);
288
289 MoFEMErrorCode preProcess();
290};
291
292#ifdef __DIRICHLET_HPP__
293
294/**
295 * \brief Assemble F_lambda into the right hand side
296 *
297 * postProcess - assembly F_lambda
298 *
299 */
300struct AssembleFlambda : public FEMethod {
301
302 boost::shared_ptr<ArcLengthCtx> arcPtr;
303
304 AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
305 boost::shared_ptr<DirichletDisplacementBc> bc =
306 boost::shared_ptr<DirichletDisplacementBc>());
307
308 MoFEMErrorCode preProcess();
309 MoFEMErrorCode operator()();
310 MoFEMErrorCode postProcess();
311
312 inline void pushDirichletBC(boost::shared_ptr<DirichletDisplacementBc> bc) {
313 bCs.push_back(bc);
314 }
315
316private:
317 std::vector<boost::shared_ptr<DirichletDisplacementBc>> bCs;
318};
319
320#endif
321
322/**
323 * |brief Simple arc-length control of force
324 *
325 * This is added for testing, it simply control force, i.e.
326 *
327 * \f[
328 * \lambda = s
329 * \f]
330 *
331 * Constructor takes one argument,
332 * @param arc_ptr Pointer to arc-length CTX.
333 */
335
336 boost::shared_ptr<ArcLengthCtx> arcPtr;
337 const bool aSsemble;
338
339 SimpleArcLengthControl(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
340 const bool assemble = false);
342
343 MoFEMErrorCode preProcess();
344 MoFEMErrorCode operator()();
345 MoFEMErrorCode postProcess();
346
347 /** \brief Calculate internal lambda
348 */
349 double calculateLambdaInt();
350
351 /** \brief Calculate db
352 */
353 MoFEMErrorCode calculateDb();
354
355 MoFEMErrorCode calculateDxAndDlambda(Vec x);
356};
357
358/** \brief Implementation of spherical arc-length method
359 * \ingroup arc_length_control
360
361 \f[
362 \alpha \| \Delta\mathbf{x} \|^2
363 + \Delta\lambda^2 \beta^2 \| \mathbf{F}_\lambda \|^2
364 = s^2
365 \f]
366
367 This is particular implementation of ArcLength control, i.e. spherical arc
368 length control. If beta is set to 0 and alpha is non-zero it is cylindrical
369 arc-length control. Works well with general problem with non-linear
370 geometry. It not guarantee dissipative loading path in case of physical
371 nonlinearities.
372
373 */
375
376 ArcLengthCtx *arcPtrRaw; // this is for back compatibility
377
378 /// \deprecated use constructor with shared_ptr
380
381 SphericalArcLengthControl(boost::shared_ptr<ArcLengthCtx> &arc_ptr);
383
384 MoFEMErrorCode preProcess();
385 MoFEMErrorCode operator()();
386 MoFEMErrorCode postProcess();
387
388 /** \brief Calculate f_lambda(dx,lambda)
389
390 \f[
391 f_\lambda(\Delta\mathbf{x},\lambda) =
392 \alpha \| \Delta\mathbf{x} \|^2
393 + \Delta\lambda^2 \beta^2 \| \mathbf{F}_\lambda \|^2
394 \f]
395
396 */
397 virtual double calculateLambdaInt();
398
399 /** \brief Calculate db
400
401 \f[
402 \textrm{d}\mathbf{b} = 2 \alpha \Delta\mathbf{x}
403 \f]
404
405 */
406 virtual MoFEMErrorCode calculateDb();
407 virtual MoFEMErrorCode calculateDxAndDlambda(Vec x);
408 virtual MoFEMErrorCode calculateInitDlambda(double *dlambda);
409 virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda);
410
411private:
412 boost::shared_ptr<ArcLengthCtx> arcPtr;
413};
414
415#endif // __ARC_LENGTH_TOOLS_HPP__
416
417/**
418 \defgroup arc_length_control Arc-Length control
419 \ingroup user_modules
420**/
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
MoFEMErrorCode PCSetupArcLength(PC pc)
#define DEPRECATED
Definition: definitions.h:17
static double lambda
constexpr AssemblyType A
constexpr auto field_name
Store variables for ArcLength analysis.
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
double s
double res_lambda
f_lambda - s
DofIdx getPetscLocalDofIdx()
Get local index of load factor.
SmartPetscObj< Vec > dx
dx = x-x0
MoFEMErrorCode setAlphaBeta(double alpha, double beta)
set parameters controlling arc-length equations alpha controls off diagonal therms beta controls diag...
double alpha
displacement scaling factor
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
MoFEMErrorCode setS(double s)
MoFEM::Interface & mField
double F_lambda2
inner_prod(F_lambda,F_lambda);
double dIag
diagonal value
double dx2
inner_prod(dX,dX)
virtual ~ArcLengthCtx()=default
FieldData & getFieldData()
Get value of load factor.
SmartPetscObj< Vec > ghostDiag
SmartPetscObj< Vec > x0
displacement vector at beginning of step
double dLambda
increment of load factor
double beta
force scaling factor
SmartPetscObj< Vec > ghosTdLambda
int getPart()
Get proc owning lambda dof.
NumeredDofEntity * arcDofRawPtr
DofIdx getPetscGlobalDofIdx()
Get global index of load factor.
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
shell matrix for arc-length method
ArcLengthCtx * arcPtrRaw
friend MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
boost::shared_ptr< ArcLengthCtx > arcPtr
SmartPetscObj< Mat > Aij
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
virtual ~ArcLengthMatShell()=default
Deprecated interface functions.
structure for Arc Length pre-conditioner
ArcLengthCtx * arcPtrRaw
friend MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
boost::shared_ptr< ArcLengthCtx > arcPtr
friend MoFEMErrorCode PCSetupArcLength(PC pc)
SmartPetscObj< Mat > Aij
SmartPetscObj< Mat > shellAij
SmartPetscObj< KSP > kSP
SmartPetscObj< PC > pC
MoFEMErrorCode preProcess()
MoFEMErrorCode calculateDb()
Calculate db.
MoFEMErrorCode operator()()
double calculateLambdaInt()
Calculate internal lambda.
MoFEMErrorCode calculateDxAndDlambda(Vec x)
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode postProcess()
Implementation of spherical arc-length method.
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode postProcess()
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
virtual MoFEMErrorCode calculateDb()
Calculate db.
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Zero F_lambda.
boost::shared_ptr< ArcLengthCtx > arcPtr
MoFEMErrorCode preProcess()