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