v0.10.0
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  /** \brief Get global index of load factor
116  */
118  return arcDofRawPtr->getPetscGlobalDofIdx();
119  };
120 
121  /** \brief Get local index of load factor
122  */
123  DofIdx getPetscLocalDofIdx() { return arcDofRawPtr->getPetscLocalDofIdx(); };
124 
125  /** \brief Get value of load factor
126  */
127  FieldData &getFieldData() { return arcDofRawPtr->getFieldData(); }
128 
129  /** \brief Get proc owning lambda dof
130  */
131  int getPart() { return arcDofRawPtr->getPart(); };
132 
133 private:
134  NumeredDofEntity *arcDofRawPtr;
135 };
136 
137 #ifdef __SNESCTX_HPP__
138 
139 /**
140  * \brief It is ctx structure passed to SNES solver
141  * \ingroup arc_length_control
142  */
143 struct ArcLengthSnesCtx : public SnesCtx {
144  ArcLengthCtx *arcPtrRaw;
145 
146  ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
147  ArcLengthCtx *arc_ptr_raw)
148  : SnesCtx(m_field, problem_name), arcPtrRaw(arc_ptr_raw) {}
149 
150  ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
151  boost::shared_ptr<ArcLengthCtx> arc_ptr)
152  : SnesCtx(m_field, problem_name), arcPtrRaw(arc_ptr.get()),
153  arcPtr(arc_ptr) {}
154 
155 private:
156  boost::shared_ptr<ArcLengthCtx> arcPtr;
157 };
158 
159 #endif //__SNESCTX_HPP__
160 
161 #ifdef __TSCTX_HPP__
162 
163 /**
164  * \brief It is ctx structure passed to SNES solver
165  * \ingroup arc_length_control
166  */
167 struct ArcLengthTsCtx : public TsCtx {
168  ArcLengthCtx *arcPtrRaw;
169 
170  /// \deprecated use constructor with shared ptr
171  DEPRECATED ArcLengthTsCtx(MoFEM::Interface &m_field,
172  const std::string &problem_name,
173  ArcLengthCtx *arc_ptr_raw)
174  : TsCtx(m_field, problem_name), arcPtrRaw(arc_ptr_raw) {}
175 
176  ArcLengthTsCtx(MoFEM::Interface &m_field, const std::string &problem_name,
177  boost::shared_ptr<ArcLengthCtx> arc_ptr)
178  : TsCtx(m_field, problem_name), arcPtrRaw(arc_ptr.get()),
179  arcPtr(arc_ptr) {}
180 
181 private:
182  boost::shared_ptr<ArcLengthCtx> arcPtr;
183 };
184 
185 #endif // __TSCTX_HPP__
186 
187 /** \brief shell matrix for arc-length method
188  *
189  * \ingroup arc_length_control
190 
191  Shell matrix which has structure:
192  \f[
193  \left[
194  \begin{array}{cc}
195  \mathbf{K} & -\mathbf{F}_\lambda \\
196  \textrm{d}\mathbf{b} & D
197  \end{array}
198  \right]
199  \left\{
200  \begin{array}{c}
201  \delta \Delta \mathbf{x} \\
202  \delta \Delta \lambda
203  \end{array}
204  \right\}
205  =
206  \left[
207  \begin{array}{c}
208  -\mathbf{f}_\textrm{int} \\
209  -r_\lambda
210  \end{array}
211  \right]
212  \f]
213 
214  */
216 
217  SmartPetscObj<Mat> Aij;
218  string problemName;
219  ArcLengthCtx *arcPtrRaw; // this is for back compatibility
220 
221  ArcLengthMatShell(Mat aij, boost::shared_ptr<ArcLengthCtx> arc_ptr,
222  string problem_name);
223 
224  /// \deprecated use constructor with shared_ptr
225  DEPRECATED ArcLengthMatShell(Mat aij, ArcLengthCtx *arc_ptr_raw,
226  string problem_name);
227 
228  MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode);
229  friend MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f);
230 
231 private:
232  boost::shared_ptr<ArcLengthCtx> arcPtr;
233 };
234 
235 /**
236  * mult operator for Arc Length Shell Mat
237  */
238 MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f);
239 
240 /**
241  * \brief structure for Arc Length pre-conditioner
242  * \ingroup arc_length_control
243  */
245 
246  SmartPetscObj<KSP> kSP;
247  SmartPetscObj<PC> pC;
248  SmartPetscObj<Mat> shellAij;
249  SmartPetscObj<Mat> Aij;
250 
251  PCArcLengthCtx(Mat shell_Aij, Mat aij,
252  boost::shared_ptr<ArcLengthCtx> &arc_ptr);
253 
254  PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
255  boost::shared_ptr<ArcLengthCtx> &arc_ptr);
256 
257  ArcLengthCtx *arcPtrRaw; // this is for back compatibility
258  /// \deprecated use with shared_ptr
259  DEPRECATED PCArcLengthCtx(Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr_raw);
260 
261  friend MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x);
262  friend MoFEMErrorCode PCSetupArcLength(PC pc);
263 
264 private:
265  boost::shared_ptr<ArcLengthCtx> arcPtr;
266 };
267 
268 /**
269  * apply operator for Arc Length pre-conditioner
270  * solves K*pc_x = pc_f
271  * solves K*xLambda = -dF_lambda
272  * solves ddlambda = ( res_lambda - db*xLambda )/( diag + db*pc_x )
273  * calculate pc_x = pc_x + ddlambda*xLambda
274  */
275 MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x);
276 
277 /**
278  * set up structure for Arc Length pre-conditioner
279 
280  * it sets pre-conditioner for matrix K
281  */
283 
284 /**
285  * \brief Zero F_lambda
286  *
287  */
288 struct ZeroFLmabda : public FEMethod {
289 
290  boost::shared_ptr<ArcLengthCtx> arcPtr;
291 
292  ZeroFLmabda(boost::shared_ptr<ArcLengthCtx> arc_ptr);
293 
295 };
296 
297 #ifdef __DIRICHLET_HPP__
298 
299 /**
300  * \brief Assemble F_lambda into the right hand side
301  *
302  * postProcess - assembly F_lambda
303  *
304  */
305 struct AssembleFlambda : public FEMethod {
306 
307  boost::shared_ptr<ArcLengthCtx> arcPtr;
308 
309  AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
310  boost::shared_ptr<DirichletDisplacementBc> bc =
311  boost::shared_ptr<DirichletDisplacementBc>());
312 
313  MoFEMErrorCode preProcess();
314  MoFEMErrorCode operator()();
315  MoFEMErrorCode postProcess();
316 
317  inline void pushDirichletBC(boost::shared_ptr<DirichletDisplacementBc> bc) {
318  bCs.push_back(bc);
319  }
320 
321 private:
322  std::vector<boost::shared_ptr<DirichletDisplacementBc>> bCs;
323 };
324 
325 #endif
326 
327 /**
328  * |brief Simple arc-length control of force
329  *
330  * This is added for testing, it simply control force, i.e.
331  *
332  * \f[
333  * \lambda = s
334  * \f]
335  *
336  * Constructor takes one argument,
337  * @param arc_ptr Pointer to arc-length CTX.
338  */
340 
341  boost::shared_ptr<ArcLengthCtx> arcPtr;
342  const bool aSsemble;
343 
344  SimpleArcLengthControl(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
345  const bool assemble = false);
347 
351 
352  /** \brief Calculate internal lambda
353  */
354  double calculateLambdaInt();
355 
356  /** \brief Calculate db
357  */
359 
361 };
362 
363 /** \brief Implementation of spherical arc-length method
364  * \ingroup arc_length_control
365 
366  \f[
367  \alpha \| \Delta\mathbf{x} \|^2
368  + \Delta\lambda^2 \beta^2 \| \mathbf{F}_\lambda \|^2
369  = s^2
370  \f]
371 
372  This is particular implementation of ArcLength control, i.e. spherical arc
373  length control. If beta is set to 0 and alpha is non-zero it is cylindrical
374  arc-length control. Works well with general problem with non-linear
375  geometry. It not guarantee dissipative loading path in case of physical
376  nonlinearities.
377 
378  */
380 
381  ArcLengthCtx *arcPtrRaw; // this is for back compatibility
382 
383  /// \deprecated use constructor with shared_ptr
385 
386  SphericalArcLengthControl(boost::shared_ptr<ArcLengthCtx> &arc_ptr);
387  virtual ~SphericalArcLengthControl();
388 
392 
393  /** \brief Calculate f_lambda(dx,lambda)
394 
395  \f[
396  f_\lambda(\Delta\mathbf{x},\lambda) =
397  \alpha \| \Delta\mathbf{x} \|^2
398  + \Delta\lambda^2 \beta^2 \| \mathbf{F}_\lambda \|^2
399  \f]
400 
401  */
402  virtual double calculateLambdaInt();
403 
404  /** \brief Calculate db
405 
406  \f[
407  \textrm{d}\mathbf{b} = 2 \alpha \Delta\mathbf{x}
408  \f]
409 
410  */
411  virtual MoFEMErrorCode calculateDb();
412  virtual MoFEMErrorCode calculateDxAndDlambda(Vec x);
413  virtual MoFEMErrorCode calculateInitDlambda(double *dlambda);
414  virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda);
415 
416 private:
417  boost::shared_ptr<ArcLengthCtx> arcPtr;
418 };
419 
420 #endif // __ARC_LENGTH_TOOLS_HPP__
421 
422 /**
423  \defgroup arc_length_control Arc-Length control
424  \ingroup user_modules
425 **/
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
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:67
NumeredDofEntity * arcDofRawPtr
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)