v0.14.0
ArcLengthTools.cpp
Go to the documentation of this file.
1 /** \file ArcLengthTools.hpp
2 
3  Implementation of Arc Length element
4 
5  */
6 
7 
8 
9 #include <MoFEM.hpp>
10 using namespace MoFEM;
11 #include <ArcLengthTools.hpp>
12 
13 #define ArcFunctionBegin \
14  MoFEMFunctionBegin; \
15  MOFEM_LOG_CHANNEL("WORLD"); \
16  MOFEM_LOG_FUNCTION(); \
17  MOFEM_LOG_TAG("WORLD", "Arc");
18 
19 // ********************
20 // Arc-length ctx class
21 
22 MoFEMErrorCode ArcLengthCtx::setS(double s) {
23  ArcFunctionBegin;
24  this->s = s;
25  MOFEM_LOG_C("WORLD", Sev::inform, "\tSet s = %6.4e", this->s);
26  MoFEMFunctionReturn(0);
27 }
28 
29 MoFEMErrorCode ArcLengthCtx::setAlphaBeta(double alpha, double beta) {
30  ArcFunctionBegin;
31  this->alpha = alpha;
32  this->beta = beta;
33  MOFEM_LOG_C("WORLD", Sev::inform, "\tSet alpha = %6.4e beta = %6.4e",
34  this->alpha, this->beta);
35  MoFEMFunctionReturn(0);
36 }
37 
38 ArcLengthCtx::ArcLengthCtx(MoFEM::Interface &m_field,
39  const std::string &problem_name,
40  const std::string &field_name)
41  : mField(m_field), dx2(0), F_lambda2(0), res_lambda(0) {
42 
43  auto create_f_lambda = [&]() {
44  ArcFunctionBegin;
45  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(problem_name, ROW,
46  F_lambda);
47  CHKERR VecSetOption(F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
48  MoFEMFunctionReturn(0);
49  };
50 
51  auto vec_duplicate = [&]() {
52  MoFEMFunctionBegin;
53  db = vectorDuplicate(F_lambda);
54  xLambda = vectorDuplicate(F_lambda);
55  x0 = vectorDuplicate(F_lambda);
56  dx = vectorDuplicate(F_lambda);
57  MoFEMFunctionReturn(0);
58  };
59 
60  auto zero_vectors = [&]() {
61  MoFEMFunctionBegin;
62  CHKERR VecZeroEntries(F_lambda);
63  CHKERR VecZeroEntries(db);
64  CHKERR VecZeroEntries(xLambda);
65  CHKERR VecZeroEntries(x0);
66  CHKERR VecZeroEntries(dx);
67  MoFEMFunctionReturn(0);
68  };
69 
70  auto find_lambda_dof = [&]() {
71  MoFEMFunctionBegin;
72 
73  const Problem *problem_ptr;
74  CHKERR m_field.get_problem(problem_name, &problem_ptr);
75  boost::shared_ptr<NumeredDofEntity_multiIndex> dofs_ptr_no_const =
76  problem_ptr->getNumeredRowDofsPtr();
77  auto bit_number = m_field.get_field_bit_number(field_name);
78  auto dIt = dofs_ptr_no_const->get<Unique_mi_tag>().lower_bound(
79  FieldEntity::getLoBitNumberUId(bit_number));
80  auto hi_dit = dofs_ptr_no_const->get<Unique_mi_tag>().upper_bound(
81  FieldEntity::getHiBitNumberUId(bit_number));
82  if (std::distance(dIt, hi_dit) != 1)
83  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
84  ("can not find unique LAMBDA (load factor) but found " +
85  boost::lexical_cast<std::string>(std::distance(dIt, hi_dit)))
86  .c_str());
87  arcDofRawPtr = (*dIt).get();
88  MoFEMFunctionReturn(0);
89  };
90 
91  auto create_ghost_vecs = [&]() {
92  MoFEMFunctionBegin;
93  Vec ghost_d_lambda, ghost_diag;
94  if ((unsigned int)mField.get_comm_rank() == arcDofRawPtr->getPart()) {
95  CHKERR VecCreateGhostWithArray(mField.get_comm(), 1, 1, 0, PETSC_NULL,
96  &dLambda, &ghost_d_lambda);
97 
98  CHKERR VecCreateGhostWithArray(mField.get_comm(), 1, 1, 0, PETSC_NULL,
99  &dIag, &ghost_diag);
100  } else {
101  int one[] = {0};
102  CHKERR VecCreateGhostWithArray(mField.get_comm(), 0, 1, 1, one, &dLambda,
103  &ghost_d_lambda);
104  CHKERR VecCreateGhostWithArray(mField.get_comm(), 0, 1, 1, one, &dIag,
105  &ghost_diag);
106  }
107  dLambda = 0;
108  dIag = 0;
109  ghosTdLambda = SmartPetscObj<Vec>(ghost_d_lambda);
110  ghostDiag = SmartPetscObj<Vec>(ghost_diag);
111  MoFEMFunctionReturn(0);
112  };
113 
114  ierr = create_f_lambda();
115  CHKERRABORT(PETSC_COMM_SELF, ierr);
116 
117  ierr = vec_duplicate();
118  CHKERRABORT(PETSC_COMM_SELF, ierr);
119 
120  ierr = zero_vectors();
121  CHKERRABORT(PETSC_COMM_SELF, ierr);
122 
123  ierr = find_lambda_dof();
124  CHKERRABORT(PETSC_COMM_SELF, ierr);
125 
126  ierr = create_ghost_vecs();
127  CHKERRABORT(PETSC_COMM_SELF, ierr);
128 }
129 
130 // ***********************
131 // Arc-length shell matrix
132 
133 ArcLengthMatShell::ArcLengthMatShell(Mat aij, ArcLengthCtx *arc_ptr_raw,
134  string problem_name)
135  : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr_raw) {}
136 
137 ArcLengthMatShell::ArcLengthMatShell(Mat aij,
138  boost::shared_ptr<ArcLengthCtx> arc_ptr,
139  string problem_name)
140  : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr.get()),
141  arcPtr(arc_ptr) {}
142 
143 MoFEMErrorCode ArcLengthMatShell::setLambda(Vec ksp_x, double *lambda,
144  ScatterMode scattermode) {
145  MoFEMFunctionBegin;
146 
147  int part = arcPtrRaw->getPart();
148  int rank = arcPtrRaw->mField.get_comm_rank();
149 
150  switch (scattermode) {
151  case SCATTER_FORWARD: {
152  Vec lambda_ghost;
153  if (rank == part) {
154  CHKERR VecCreateGhostWithArray(arcPtrRaw->mField.get_comm(), 1, 1, 0,
155  PETSC_NULL, lambda, &lambda_ghost);
156  } else {
157  int one[] = {0};
158  CHKERR VecCreateGhostWithArray(arcPtrRaw->mField.get_comm(), 0, 1, 1, one,
159  lambda, &lambda_ghost);
160  }
161  int idx = arcPtrRaw->getPetscGlobalDofIdx();
162  if (part == rank) {
163  CHKERR VecGetValues(ksp_x, 1, &idx, lambda);
164  }
165  CHKERR VecGhostUpdateBegin(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
166  CHKERR VecGhostUpdateEnd(lambda_ghost, INSERT_VALUES, SCATTER_FORWARD);
167  CHKERR VecDestroy(&lambda_ghost);
168  } break;
169  case SCATTER_REVERSE: {
170  if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
171  PetscScalar *array;
172  CHKERR VecGetArray(ksp_x, &array);
173  array[arcPtrRaw->getPetscLocalDofIdx()] = *lambda;
174  CHKERR VecRestoreArray(ksp_x, &array);
175  }
176  } break;
177  default:
178  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
179  }
180 
181  MoFEMFunctionReturn(0);
182 }
183 
184 MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f) {
185  MoFEMFunctionBegin;
186  void *void_ctx;
187  CHKERR MatShellGetContext(A, &void_ctx);
188  ArcLengthMatShell *ctx = static_cast<ArcLengthMatShell *>(void_ctx);
189  CHKERR MatMult(ctx->Aij, x, f);
190  double lambda;
191  CHKERR ctx->setLambda(x, &lambda, SCATTER_FORWARD);
192  double db_dot_x;
193  CHKERR VecDot(ctx->arcPtrRaw->db, x, &db_dot_x);
194  double f_lambda;
195  f_lambda = ctx->arcPtrRaw->dIag * lambda + db_dot_x;
196  CHKERR ctx->setLambda(f, &f_lambda, SCATTER_REVERSE);
197  CHKERR VecAXPY(f, lambda, ctx->arcPtrRaw->F_lambda);
198  MoFEMFunctionReturn(0);
199 }
200 
201 // arc-length preconditioner
202 
203 PCArcLengthCtx::PCArcLengthCtx(Mat shell_Aij, Mat aij, ArcLengthCtx *arc_ptr)
204  : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr) {
205  auto comm = PetscObjectComm((PetscObject)aij);
206  pC = createPC(comm);
207  kSP = createKSP(comm);
208  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
209  CHKERRABORT(PETSC_COMM_WORLD, ierr);
210 }
211 
212 PCArcLengthCtx::PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
213  ArcLengthCtx *arc_ptr)
214  : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
215  arcPtrRaw(arc_ptr) {
216  auto comm = PetscObjectComm((PetscObject)aij);
217  kSP = createKSP(comm);
218  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
219  CHKERRABORT(PETSC_COMM_WORLD, ierr);
220 }
221 
222 PCArcLengthCtx::PCArcLengthCtx(Mat shell_Aij, Mat aij,
223  boost::shared_ptr<ArcLengthCtx> arc_ptr)
224  : shellAij(shell_Aij, true), Aij(aij, true), arcPtrRaw(arc_ptr.get()),
225  arcPtr(arc_ptr) {
226  auto comm = PetscObjectComm((PetscObject)aij);
227  pC = createPC(comm);
228  kSP = createKSP(comm);
229  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
230  CHKERRABORT(PETSC_COMM_WORLD, ierr);
231 }
232 
233 PCArcLengthCtx::PCArcLengthCtx(PC pc, Mat shell_Aij, Mat aij,
234  boost::shared_ptr<ArcLengthCtx> arc_ptr)
235  : pC(pc, true), shellAij(shell_Aij, true), Aij(aij, true),
236  arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {
237  auto comm = PetscObjectComm((PetscObject)aij);
238  kSP = createKSP(comm);
239  ierr = KSPAppendOptionsPrefix(kSP, "arc_length_");
240  CHKERRABORT(PETSC_COMM_WORLD, ierr);
241 }
242 
243 MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x) {
244  ArcFunctionBegin;
245  void *void_ctx;
246  CHKERR PCShellGetContext(pc, &void_ctx);
247  PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
248  void *void_MatCtx;
249  MatShellGetContext(ctx->shellAij, &void_MatCtx);
250  ArcLengthMatShell *mat_ctx = static_cast<ArcLengthMatShell *>(void_MatCtx);
251  PetscBool same;
252  PetscObjectTypeCompare((PetscObject)ctx->kSP, KSPPREONLY, &same);
253 
254  double res_lambda;
255  CHKERR mat_ctx->setLambda(pc_f, &res_lambda, SCATTER_FORWARD);
256 
257  // Solve residual
258  CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
259  CHKERR KSPSetInitialGuessKnoll(ctx->kSP, PETSC_FALSE);
260  CHKERR KSPSolve(ctx->kSP, pc_f, pc_x);
261  double db_dot_pc_x;
262  CHKERR VecDot(ctx->arcPtrRaw->db, pc_x, &db_dot_pc_x);
263 
264  // Solve for x_lambda
265  if (same != PETSC_TRUE) {
266  CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_TRUE);
267  } else {
268  CHKERR KSPSetInitialGuessNonzero(ctx->kSP, PETSC_FALSE);
269  }
270  CHKERR KSPSolve(ctx->kSP, ctx->arcPtrRaw->F_lambda, ctx->arcPtrRaw->xLambda);
271  double db_dot_x_lambda;
272  CHKERR VecDot(ctx->arcPtrRaw->db, ctx->arcPtrRaw->xLambda, &db_dot_x_lambda);
273 
274  // Calculate d_lambda
275  double denominator = ctx->arcPtrRaw->dIag + db_dot_x_lambda;
276  double ddlambda = -(res_lambda - db_dot_pc_x) / denominator;
277 
278  // Update solution vector
279  CHKERR VecAXPY(pc_x, -ddlambda, ctx->arcPtrRaw->xLambda);
280  CHKERR mat_ctx->setLambda(pc_x, &ddlambda, SCATTER_REVERSE);
281 
282  if (ddlambda != ddlambda || denominator == 0) {
283 
284  double nrm2_pc_f, nrm2_db, nrm2_pc_x, nrm2_xLambda;
285  CHKERR VecNorm(pc_f, NORM_2, &nrm2_pc_f);
286  CHKERR VecNorm(ctx->arcPtrRaw->db, NORM_2, &nrm2_db);
287  CHKERR VecNorm(pc_x, NORM_2, &nrm2_pc_x);
288  CHKERR VecNorm(ctx->arcPtrRaw->xLambda, NORM_2, &nrm2_xLambda);
289 
290  MOFEM_LOG("WORLD", Sev::error) << "problem with ddlambda=" << ddlambda;
291  MOFEM_LOG("WORLD", Sev::error) << "res_lambda=" << res_lambda;
292  MOFEM_LOG("WORLD", Sev::error) << "denominator=" << denominator;
293  MOFEM_LOG("WORLD", Sev::error) << "db_dot_pc_x=" << db_dot_pc_x;
294  MOFEM_LOG("WORLD", Sev::error) << "db_dot_x_lambda=" << db_dot_x_lambda;
295  MOFEM_LOG("WORLD", Sev::error) << "diag=" << ctx->arcPtrRaw->dIag;
296  MOFEM_LOG("WORLD", Sev::error) << "nrm2_db=" << nrm2_db;
297  MOFEM_LOG("WORLD", Sev::error) << "nrm2_pc_f=" << nrm2_pc_f;
298  MOFEM_LOG("WORLD", Sev::error) << "nrm2_pc_x=" << nrm2_pc_x;
299  MOFEM_LOG("WORLD", Sev::error) << "nrm2_xLambda=" << nrm2_xLambda;
300 
301  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
302  "Increment of lambda is not number");
303 
304  }
305 
306  // Debugging PC
307  if (0) {
308  Vec y;
309  CHKERR VecDuplicate(pc_x, &y);
310  CHKERR MatMult(ctx->shellAij, pc_x, y);
311  CHKERR VecAXPY(y, -1, pc_f);
312  double res_lambda_y;
313  CHKERR mat_ctx->setLambda(y, &res_lambda_y, SCATTER_FORWARD);
314  double zero;
315  CHKERR mat_ctx->setLambda(y, &zero, SCATTER_REVERSE);
316  double norm_y;
317  CHKERR VecNorm(y, NORM_2, &norm_y);
318  MOFEM_LOG_C("WORLD", Sev::noisy, "Debug res y = %3.4e res_lambda_y = %3.4e",
319  norm_y, res_lambda_y);
320  CHKERR VecDestroy(&y);
321  }
322 
323  MoFEMFunctionReturn(0);
324 }
325 
326 MoFEMErrorCode PCSetupArcLength(PC pc) {
327  MoFEMFunctionBegin;
328  void *void_ctx;
329  CHKERR PCShellGetContext(pc, &void_ctx);
330  PCArcLengthCtx *ctx = static_cast<PCArcLengthCtx *>(void_ctx);
331  auto get_pc_ops = [&](auto pc) {
332  MoFEMFunctionBegin;
333  Mat shell_aij_raw, aij_raw;
334  CHKERR PCGetOperators(pc, &shell_aij_raw, &aij_raw);
335  ctx->shellAij = SmartPetscObj<Mat>(shell_aij_raw, true);
336  ctx->Aij = SmartPetscObj<Mat>(aij_raw, true);
337  MoFEMFunctionReturn(0);
338  };
339  CHKERR get_pc_ops(pc);
340  CHKERR PCSetUseAmat(pc, PETSC_TRUE);
341  CHKERR PCSetOperators(ctx->pC, ctx->Aij, ctx->Aij);
342  CHKERR PCSetFromOptions(ctx->pC);
343  CHKERR PCSetUp(ctx->pC);
344 #if PETSC_VERSION_LT(3, 12, 0)
345  CHKERR KSPSetTabLevel(ctx->kSP, 3);
346 #else
347  CHKERR PetscObjectSetTabLevel((PetscObject)ctx->kSP, 3);
348 #endif
349  CHKERR KSPSetFromOptions(ctx->kSP);
350  CHKERR KSPSetOperators(ctx->kSP, ctx->Aij, ctx->Aij);
351  CHKERR KSPSetPC(ctx->kSP, ctx->pC);
352  CHKERR KSPSetUp(ctx->kSP);
353  MoFEMFunctionReturn(0);
354 }
355 
356 // ***********************
357 // Zero F_lambda vector
358 
359 ZeroFLmabda::ZeroFLmabda(boost::shared_ptr<ArcLengthCtx> arc_ptr)
360  : arcPtr(arc_ptr) {}
361 
362 MoFEMErrorCode ZeroFLmabda::preProcess() {
363  MoFEMFunctionBegin;
364  switch (snes_ctx) {
365  case CTX_SNESSETFUNCTION: {
366 
367  auto zero_vals = [&](auto v) {
368  MoFEMFunctionBeginHot;
369  int size = problemPtr->getNbLocalDofsRow();
370  int ghosts = problemPtr->getNbGhostDofsRow();
371  double *array;
372  CHKERR VecGetArray(v, &array);
373  for (int i = 0; i != size + ghosts; ++i)
374  array[i] = 0;
375  CHKERR VecRestoreArray(v, &array);
376  MoFEMFunctionReturnHot(0);
377  };
378 
379  Vec l_x_lambda, l_f_lambda;
380  CHKERR VecGhostGetLocalForm(arcPtr->xLambda, &l_x_lambda);
381  CHKERR VecGhostGetLocalForm(arcPtr->F_lambda, &l_f_lambda);
382  CHKERR zero_vals(l_x_lambda);
383  CHKERR zero_vals(l_f_lambda);
384  CHKERR VecGhostRestoreLocalForm(arcPtr->xLambda, &l_x_lambda);
385  CHKERR VecGhostRestoreLocalForm(arcPtr->F_lambda, &l_f_lambda);
386 
387  } break;
388  default:
389  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
390  "Lambda can be zeroed ONLY when the right hand side is evaluated.");
391  }
392  MoFEMFunctionReturn(0);
393 }
394 
395 AssembleFlambda::AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
396  boost::shared_ptr<DirichletDisplacementBc> bc)
397  : arcPtr(arc_ptr) {
398  bCs.push_back(bc);
399 }
400 
401 MoFEMErrorCode AssembleFlambda::preProcess() {
402  MoFEMFunctionBeginHot;
403  MoFEMFunctionReturnHot(0);
404 }
405 MoFEMErrorCode AssembleFlambda::operator()() {
406  MoFEMFunctionBeginHot;
407  MoFEMFunctionReturnHot(0);
408 }
409 
410 MoFEMErrorCode AssembleFlambda::postProcess() {
411  ArcFunctionBegin;
412  switch (snes_ctx) {
413  case CTX_SNESSETFUNCTION: {
414 
415  CHKERR VecAssemblyBegin(arcPtr->F_lambda);
416  CHKERR VecAssemblyEnd(arcPtr->F_lambda);
417  CHKERR VecAssemblyBegin(snes_f);
418  CHKERR VecAssemblyEnd(snes_f);
419  CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES, SCATTER_REVERSE);
420  CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES, SCATTER_REVERSE);
421  CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
422  CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
423 
424  auto set_bc = [&](auto l_snes_f, auto l_f_lambda) {
425  MoFEMFunctionBegin;
426  if (!bCs.empty()) {
427  double *f_array, *f_lambda_array;
428  CHKERR VecGetArray(l_snes_f, &f_array);
429  CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
430  for (auto &bc : bCs) {
431  for (auto idx : bc->dofsIndices) {
432  auto weak_dof = problemPtr->getRowDofsByPetscGlobalDofIdx(idx);
433  if (auto shared_dof = weak_dof.lock()) {
434  f_array[shared_dof->getPetscLocalDofIdx()] = 0;
435  f_lambda_array[shared_dof->getPetscLocalDofIdx()] = 0;
436  }
437  }
438  }
439  CHKERR VecRestoreArray(l_snes_f, &f_array);
440  CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
441  }
442  MoFEMFunctionReturn(0);
443  };
444 
445  auto add_f_lambda = [&](auto l_snes_f, auto l_f_lambda) {
446  MoFEMFunctionBegin;
447  int size = problemPtr->getNbLocalDofsRow();
448  int ghosts = problemPtr->getNbGhostDofsRow();
449  double lambda = arcPtr->getFieldData();
450  int local_lambda_idx = arcPtr->getPetscLocalDofIdx();
451  double *f_array, *f_lambda_array;
452  CHKERR VecGetArray(l_snes_f, &f_array);
453  CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
454  for (int i = 0; i != size; ++i) {
455  f_array[i] += lambda * f_lambda_array[i];
456  }
457  CHKERR VecRestoreArray(l_snes_f, &f_array);
458  CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
459  MoFEMFunctionReturn(0);
460  };
461 
462  auto zero_ghost = [&](auto l_snes_f, auto l_f_lambda) {
463  MoFEMFunctionBegin;
464  int size = problemPtr->getNbLocalDofsRow();
465  int ghosts = problemPtr->getNbGhostDofsRow();
466  double lambda = arcPtr->getFieldData();
467  int local_lambda_idx = arcPtr->getPetscLocalDofIdx();
468  double *f_array, *f_lambda_array;
469  CHKERR VecGetArray(l_snes_f, &f_array);
470  CHKERR VecGetArray(l_f_lambda, &f_lambda_array);
471  f_lambda_array[local_lambda_idx] = 0;
472  for (int i = size; i != size + ghosts; ++i) {
473  f_array[i] = 0;
474  f_lambda_array[i] = 0;
475  }
476  CHKERR VecRestoreArray(l_snes_f, &f_array);
477  CHKERR VecRestoreArray(l_f_lambda, &f_lambda_array);
478  MoFEMFunctionReturn(0);
479  };
480 
481  Vec l_snes_f, l_f_lambda;
482  CHKERR VecGhostGetLocalForm(snes_f, &l_snes_f);
483  CHKERR VecGhostGetLocalForm(arcPtr->F_lambda, &l_f_lambda);
484  CHKERR add_f_lambda(l_snes_f, l_f_lambda);
485  CHKERR set_bc(l_snes_f, l_f_lambda);
486  CHKERR zero_ghost(l_snes_f, l_f_lambda);
487 
488  CHKERR VecGhostRestoreLocalForm(snes_f, &l_snes_f);
489  CHKERR VecGhostRestoreLocalForm(arcPtr->F_lambda, &l_f_lambda);
490 
491  double snes_fnorm, snes_xnorm;
492  CHKERR VecNorm(snes_f, NORM_2, &snes_fnorm);
493  CHKERR VecNorm(snes_x, NORM_2, &snes_xnorm);
494  CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
495 
496  MOFEM_LOG_C("WORLD", Sev::inform, "\tF_lambda2 = %6.4g lambda = %6.4g",
497  arcPtr->F_lambda2, arcPtr->getFieldData());
498  MOFEM_LOG_C("WORLD", Sev::verbose,
499  "\tsnes_f norm = %6.4e snes_x norm = %6.4g", snes_fnorm,
500  snes_xnorm);
501 
502  if (!boost::math::isfinite(snes_fnorm)) {
503  CHKERR arcPtr->mField.getInterface<Tools>()->checkVectorForNotANumber(
504  problemPtr, ROW, snes_f);
505  }
506 
507  } break;
508  default:
509  SETERRQ(
510  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511  "Lambda can be assembled only when the right hand side is evaluated.");
512 
513  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Impossible case");
514  }
515  MoFEMFunctionReturn(0);
516 }
517 
518 // ************************
519 // Simple arc-length method
520 
521 SimpleArcLengthControl::SimpleArcLengthControl(
522  boost::shared_ptr<ArcLengthCtx> &arc_ptr, const bool assemble)
523  : FEMethod(), arcPtr(arc_ptr), aSsemble(assemble) {}
524 
525 SimpleArcLengthControl::~SimpleArcLengthControl() {}
526 
527 MoFEMErrorCode SimpleArcLengthControl::preProcess() {
528  MoFEMFunctionBegin;
529  switch (snes_ctx) {
530  case CTX_SNESSETFUNCTION: {
531  if (aSsemble) {
532  CHKERR VecAssemblyBegin(snes_f);
533  CHKERR VecAssemblyEnd(snes_f);
534  }
535  CHKERR calculateDxAndDlambda(snes_x);
536  CHKERR calculateDb();
537  } break;
538  case CTX_SNESSETJACOBIAN: {
539  if (aSsemble) {
540  CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
541  CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
542  }
543  } break;
544  default:
545  break;
546  }
547  MoFEMFunctionReturn(0);
548 }
549 
550 MoFEMErrorCode SimpleArcLengthControl::operator()() {
551  MoFEMFunctionBegin;
552  switch (snes_ctx) {
553  case CTX_SNESSETFUNCTION: {
554  arcPtr->res_lambda = calculateLambdaInt() - arcPtr->s;
555  CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
556  arcPtr->res_lambda, ADD_VALUES);
557  } break;
558  case CTX_SNESSETJACOBIAN: {
559  arcPtr->dIag = arcPtr->beta;
560  CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
561  arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
562  } break;
563  default:
564  break;
565  }
566  MoFEMFunctionReturn(0);
567 }
568 
569 MoFEMErrorCode SimpleArcLengthControl::postProcess() {
570  MoFEMFunctionBegin;
571  switch (snes_ctx) {
572  case CTX_SNESSETFUNCTION: {
573  if (aSsemble) {
574  CHKERR VecAssemblyBegin(snes_f);
575  CHKERR VecAssemblyEnd(snes_f);
576  }
577  } break;
578  case CTX_SNESSETJACOBIAN: {
579  if (aSsemble) {
580  CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
581  CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
582  }
583  CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
584  SCATTER_FORWARD);
585  CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES, SCATTER_FORWARD);
586  } break;
587  default:
588  break;
589  }
590  MoFEMFunctionReturn(0);
591 }
592 
593 double SimpleArcLengthControl::calculateLambdaInt() {
594  return arcPtr->beta * arcPtr->dLambda;
595 }
596 
597 MoFEMErrorCode SimpleArcLengthControl::calculateDb() {
598  MoFEMFunctionBegin;
599  CHKERR VecZeroEntries(arcPtr->db);
600  CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
601  CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
602  MoFEMFunctionReturn(0);
603 }
604 
605 MoFEMErrorCode SimpleArcLengthControl::calculateDxAndDlambda(Vec x) {
606  ArcFunctionBegin;
607  // Calculate dx
608  CHKERR VecGhostUpdateBegin(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
609  CHKERR VecGhostUpdateEnd(arcPtr->x0, INSERT_VALUES, SCATTER_FORWARD);
610 
611  Vec l_x, l_x0, l_dx;
612  CHKERR VecGhostGetLocalForm(x, &l_x);
613  CHKERR VecGhostGetLocalForm(arcPtr->x0, &l_x0);
614  CHKERR VecGhostGetLocalForm(arcPtr->dx, &l_dx);
615  {
616  double *x_array, *x0_array, *dx_array;
617  CHKERR VecGetArray(l_x, &x_array);
618  CHKERR VecGetArray(l_x0, &x0_array);
619  CHKERR VecGetArray(l_dx, &dx_array);
620  int size =
621  problemPtr->getNbLocalDofsRow() + problemPtr->getNbGhostDofsRow();
622  for (int i = 0; i != size; ++i) {
623  dx_array[i] = x_array[i] - x0_array[i];
624  }
625  CHKERR VecRestoreArray(l_x, &x_array);
626  CHKERR VecRestoreArray(l_x0, &x0_array);
627  CHKERR VecRestoreArray(l_dx, &dx_array);
628  }
629  CHKERR VecGhostRestoreLocalForm(x, &l_x);
630  CHKERR VecGhostRestoreLocalForm(arcPtr->x0, &l_x0);
631  CHKERR VecGhostRestoreLocalForm(arcPtr->dx, &l_dx);
632 
633  // Calculate dlambda
634  if (arcPtr->getPetscLocalDofIdx() != -1) {
635  double *array;
636  CHKERR VecGetArray(arcPtr->dx, &array);
637  arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
638  array[arcPtr->getPetscLocalDofIdx()] = 0;
639  CHKERR VecRestoreArray(arcPtr->dx, &array);
640  }
641  CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
642  SCATTER_FORWARD);
643  CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
644  SCATTER_FORWARD);
645 
646  // Calculate dx2
647  double x_nrm, x0_nrm;
648  CHKERR VecNorm(x, NORM_2, &x_nrm);
649  CHKERR VecNorm(arcPtr->x0, NORM_2, &x0_nrm);
650  CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
651 
652  MOFEM_LOG_C("WORLD", Sev::verbose,
653  "\tx norm = %6.4e x0 norm = %6.4e dx2 = %6.4e", x_nrm, x0_nrm,
654  arcPtr->dx2);
655  MoFEMFunctionReturn(0);
656 }
657 
658 // ***************************
659 // Spherical arc-length control
660 
661 SphericalArcLengthControl::SphericalArcLengthControl(ArcLengthCtx *arc_ptr_raw)
662  : FEMethod(), arcPtrRaw(arc_ptr_raw) {}
663 
664 SphericalArcLengthControl::SphericalArcLengthControl(
665  boost::shared_ptr<ArcLengthCtx> &arc_ptr)
666  : FEMethod(), arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {}
667 
668 SphericalArcLengthControl::~SphericalArcLengthControl() {}
669 
670 MoFEMErrorCode SphericalArcLengthControl::preProcess() {
671  MoFEMFunctionBegin;
672  switch (snes_ctx) {
673  case CTX_SNESSETFUNCTION: {
674  CHKERR calculateDxAndDlambda(snes_x);
675  CHKERR calculateDb();
676  } break;
677  case CTX_SNESSETJACOBIAN: {
678  } break;
679  default:
680  break;
681  }
682  MoFEMFunctionReturn(0);
683 }
684 
685 double SphericalArcLengthControl::calculateLambdaInt() {
686  return arcPtrRaw->alpha * arcPtrRaw->dx2 + pow(arcPtrRaw->dLambda, 2) *
687  pow(arcPtrRaw->beta, 2) *
688  arcPtrRaw->F_lambda2;
689 }
690 
691 MoFEMErrorCode SphericalArcLengthControl::calculateDb() {
692  MoFEMFunctionBegin;
693  CHKERR VecCopy(arcPtrRaw->dx, arcPtrRaw->db);
694  CHKERR VecScale(arcPtrRaw->db, 2 * arcPtrRaw->alpha);
695  MoFEMFunctionReturn(0);
696 }
697 
698 MoFEMErrorCode SphericalArcLengthControl::operator()() {
699  ArcFunctionBegin;
700  switch (snes_ctx) {
701  case CTX_SNESSETFUNCTION: {
702  arcPtrRaw->res_lambda = calculateLambdaInt() - pow(arcPtrRaw->s, 2);
703  CHKERR VecSetValue(snes_f, arcPtrRaw->getPetscGlobalDofIdx(),
704  arcPtrRaw->res_lambda, ADD_VALUES);
705  MOFEM_LOG_C("WORLD", Sev::verbose, "\tres_lambda = %6.4e\n",
706  arcPtrRaw->res_lambda);
707  } break;
708  case CTX_SNESSETJACOBIAN: {
709  arcPtrRaw->dIag =
710  2 * arcPtrRaw->dLambda * pow(arcPtrRaw->beta, 2) * arcPtrRaw->F_lambda2;
711  CHKERR MatSetValue(snes_B, arcPtrRaw->getPetscGlobalDofIdx(),
712  arcPtrRaw->getPetscGlobalDofIdx(), 1, ADD_VALUES);
713  } break;
714  default:
715  break;
716  }
717  MoFEMFunctionReturn(0);
718 }
719 
720 MoFEMErrorCode SphericalArcLengthControl::postProcess() {
721  ArcFunctionBegin;
722  switch (snes_ctx) {
723  case CTX_SNESSETFUNCTION: {
724  MOFEM_LOG_C("WORLD", Sev::verbose, "\tlambda = %6.4e\n",
725  arcPtrRaw->getFieldData());
726  } break;
727  case CTX_SNESSETJACOBIAN: {
728  CHKERR VecGhostUpdateBegin(arcPtrRaw->ghostDiag, INSERT_VALUES,
729  SCATTER_FORWARD);
730  CHKERR VecGhostUpdateEnd(arcPtrRaw->ghostDiag, INSERT_VALUES,
731  SCATTER_FORWARD);
732  CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
733  CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
734  MOFEM_LOG_C("WORLD", Sev::verbose, "\tdiag = %6.4e", arcPtrRaw->dIag);
735  } break;
736  default:
737  break;
738  }
739  MoFEMFunctionReturn(0);
740 }
741 
742 MoFEMErrorCode SphericalArcLengthControl::calculateDxAndDlambda(Vec x) {
743  ArcFunctionBegin;
744  // dx
745  CHKERR VecCopy(x, arcPtrRaw->dx);
746  CHKERR VecAXPY(arcPtrRaw->dx, -1, arcPtrRaw->x0);
747  CHKERR VecGhostUpdateBegin(arcPtrRaw->dx, INSERT_VALUES, SCATTER_FORWARD);
748  CHKERR VecGhostUpdateEnd(arcPtrRaw->dx, INSERT_VALUES, SCATTER_FORWARD);
749  // dlambda
750  if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
751  double *array;
752  CHKERR VecGetArray(arcPtrRaw->dx, &array);
753  arcPtrRaw->dLambda = array[arcPtrRaw->getPetscLocalDofIdx()];
754  array[arcPtrRaw->getPetscLocalDofIdx()] = 0;
755  CHKERR VecRestoreArray(arcPtrRaw->dx, &array);
756  }
757  CHKERR VecGhostUpdateBegin(arcPtrRaw->ghosTdLambda, INSERT_VALUES,
758  SCATTER_FORWARD);
759  CHKERR VecGhostUpdateEnd(arcPtrRaw->ghosTdLambda, INSERT_VALUES,
760  SCATTER_FORWARD);
761  // dx2
762  CHKERR VecDot(arcPtrRaw->dx, arcPtrRaw->dx, &arcPtrRaw->dx2);
763  MOFEM_LOG_C("WORLD", Sev::verbose, "\tdlambda = %6.4e dx2 = %6.4e\n",
764  arcPtrRaw->dLambda, arcPtrRaw->dx2);
765  MoFEMFunctionReturn(0);
766 }
767 
768 MoFEMErrorCode
769 SphericalArcLengthControl::calculateInitDlambda(double *dlambda) {
770  ArcFunctionBegin;
771  *dlambda = std::sqrt(pow(arcPtrRaw->s, 2) /
772  (pow(arcPtrRaw->beta, 2) * arcPtrRaw->F_lambda2));
773  if (!(*dlambda == *dlambda)) {
774  MOFEM_LOG("WORLD", Sev::error)
775  << "s " << arcPtrRaw->s << " " << arcPtrRaw->beta << " "
776  << arcPtrRaw->F_lambda2;
777  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
778  "Increment of lambda is not a number");
779  }
780  MoFEMFunctionReturn(0);
781 }
782 
783 MoFEMErrorCode SphericalArcLengthControl::setDlambdaToX(Vec x, double dlambda) {
784  ArcFunctionBegin;
785  // check if local dof idx is non zero, i.e. that lambda is accessible from
786  // this processor
787  if (arcPtrRaw->getPetscLocalDofIdx() != -1) {
788  double *array;
789  CHKERR VecGetArray(x, &array);
790  double lambda_old = array[arcPtrRaw->getPetscLocalDofIdx()];
791  if (!(dlambda == dlambda)) {
792  MOFEM_LOG("WORLD", Sev::error)
793  << "s " << arcPtrRaw->s << " " << arcPtrRaw->beta << " "
794  << arcPtrRaw->F_lambda2;
795  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
796  "Increment of lambda is not a number");
797  }
798  array[arcPtrRaw->getPetscLocalDofIdx()] = lambda_old + dlambda;
799  MOFEM_LOG_C("WORLD", Sev::verbose, "\tlambda = %6.4e, %6.4e (%6.4e)\n",
800  lambda_old, array[arcPtrRaw->getPetscLocalDofIdx()], dlambda);
801  CHKERR VecRestoreArray(x, &array);
802  }
803  MoFEMFunctionReturn(0);
804 }
ArcLengthCtx::res_lambda
double res_lambda
f_lambda - s
Definition: ArcLengthTools.hpp:82
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
ArcLengthCtx::getFieldData
FieldData & getFieldData()
Get value of load factor.
Definition: ArcLengthTools.hpp:117
ArcLengthCtx::mField
MoFEM::Interface & mField
Definition: ArcLengthTools.hpp:69
SimpleArcLengthControl::postProcess
MoFEMErrorCode postProcess()
Definition: ArcLengthTools.cpp:569
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
SphericalArcLengthControl::calculateLambdaInt
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
Definition: ArcLengthTools.cpp:685
ArcFunctionBegin
#define ArcFunctionBegin
Definition: ArcLengthTools.cpp:13
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
ArcLengthCtx::setS
MoFEMErrorCode setS(double s)
set arc radius
Definition: ArcLengthTools.cpp:22
ArcLengthCtx::db
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
Definition: ArcLengthTools.hpp:85
MoFEM::Tools
Auxiliary tools.
Definition: Tools.hpp:19
PCArcLengthCtx::pC
SmartPetscObj< PC > pC
Definition: ArcLengthTools.hpp:238
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
SimpleArcLengthControl::operator()
MoFEMErrorCode operator()()
Definition: ArcLengthTools.cpp:550
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:65
ArcLengthCtx::dLambda
double dLambda
increment of load factor
Definition: ArcLengthTools.hpp:76
SimpleArcLengthControl::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: ArcLengthTools.hpp:336
SimpleArcLengthControl::aSsemble
const bool aSsemble
Definition: ArcLengthTools.hpp:337
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PCArcLengthCtx::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:240
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
ArcLengthTools.hpp
PCArcLengthCtx::kSP
SmartPetscObj< KSP > kSP
Definition: ArcLengthTools.hpp:237
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEM::NumeredDofEntity::getPart
unsigned int getPart() const
Definition: DofsMultiIndices.hpp:250
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition: ArcLengthTools.cpp:326
ArcLengthMatShell::ArcLengthMatShell
ArcLengthMatShell(Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr, string problem_name)
Definition: ArcLengthTools.cpp:137
SphericalArcLengthControl::~SphericalArcLengthControl
virtual ~SphericalArcLengthControl()
Definition: ArcLengthTools.cpp:668
ZeroFLmabda::ZeroFLmabda
ZeroFLmabda(boost::shared_ptr< ArcLengthCtx > arc_ptr)
Definition: ArcLengthTools.cpp:359
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
ArcLengthCtx::getPetscGlobalDofIdx
DofIdx getPetscGlobalDofIdx()
Get global index of load factor.
Definition: ArcLengthTools.hpp:107
SphericalArcLengthControl::SphericalArcLengthControl
DEPRECATED SphericalArcLengthControl(ArcLengthCtx *arc_ptr_raw)
Definition: ArcLengthTools.cpp:661
ArcLengthCtx::beta
double beta
force scaling factor
Definition: ArcLengthTools.hpp:72
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
ArcLengthCtx::alpha
double alpha
displacement scaling factor
Definition: ArcLengthTools.hpp:73
SimpleArcLengthControl::calculateDxAndDlambda
MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: ArcLengthTools.cpp:605
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ArcLengthCtx::dx2
double dx2
inner_prod(dX,dX)
Definition: ArcLengthTools.hpp:80
SimpleArcLengthControl::SimpleArcLengthControl
SimpleArcLengthControl(boost::shared_ptr< ArcLengthCtx > &arc_ptr, const bool assemble=false)
Definition: ArcLengthTools.cpp:521
ROW
@ ROW
Definition: definitions.h:123
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition: ArcLengthTools.cpp:243
ArcLengthMatShell::setLambda
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
Definition: ArcLengthTools.cpp:143
SphericalArcLengthControl::calculateDxAndDlambda
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: ArcLengthTools.cpp:742
ZeroFLmabda::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: ArcLengthTools.hpp:285
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
ArcLengthMatShell::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:209
ArcLengthCtx::xLambda
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
Definition: ArcLengthTools.hpp:86
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
ArcLengthMatShell::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:207
MoFEM::createPC
auto createPC(MPI_Comm comm)
Definition: PetscSmartObj.hpp:263
SimpleArcLengthControl::~SimpleArcLengthControl
~SimpleArcLengthControl()
Definition: ArcLengthTools.cpp:525
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
PCArcLengthCtx::PCArcLengthCtx
PCArcLengthCtx(Mat shell_Aij, Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr)
Definition: ArcLengthTools.cpp:222
SphericalArcLengthControl::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:376
ArcLengthCtx::arcDofRawPtr
NumeredDofEntity * arcDofRawPtr
Definition: ArcLengthTools.hpp:121
SphericalArcLengthControl::operator()
MoFEMErrorCode operator()()
Definition: ArcLengthTools.cpp:698
SphericalArcLengthControl::postProcess
MoFEMErrorCode postProcess()
Definition: ArcLengthTools.cpp:720
SphericalArcLengthControl::calculateInitDlambda
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
Definition: ArcLengthTools.cpp:769
SimpleArcLengthControl::preProcess
MoFEMErrorCode preProcess()
Definition: ArcLengthTools.cpp:527
ArcLengthCtx::x0
SmartPetscObj< Vec > x0
displacement vector at beginning of step
Definition: ArcLengthTools.hpp:87
ArcLengthCtx::F_lambda2
double F_lambda2
inner_prod(F_lambda,F_lambda);
Definition: ArcLengthTools.hpp:81
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ArcLengthCtx::ArcLengthCtx
ArcLengthCtx(MoFEM::Interface &m_field, const std::string &problem_name, const std::string &field_name="LAMBDA")
Definition: ArcLengthTools.cpp:38
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
ArcLengthCtx::dx
SmartPetscObj< Vec > dx
dx = x-x0
Definition: ArcLengthTools.hpp:88
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ArcLengthCtx::getPetscLocalDofIdx
DofIdx getPetscLocalDofIdx()
Get local index of load factor.
Definition: ArcLengthTools.hpp:113
ArcLengthCtx::dIag
double dIag
diagonal value
Definition: ArcLengthTools.hpp:78
ArcLengthCtx::ghostDiag
SmartPetscObj< Vec > ghostDiag
Definition: ArcLengthTools.hpp:77
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ZeroFLmabda::preProcess
MoFEMErrorCode preProcess()
Definition: ArcLengthTools.cpp:362
SphericalArcLengthControl::setDlambdaToX
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Definition: ArcLengthTools.cpp:783
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
ArcLengthCtx::getPart
int getPart()
Get proc owning lambda dof.
Definition: ArcLengthTools.hpp:121
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
PCArcLengthCtx::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:248
PCArcLengthCtx::shellAij
SmartPetscObj< Mat > shellAij
Definition: ArcLengthTools.hpp:239
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:184
ArcLengthCtx::setAlphaBeta
MoFEMErrorCode setAlphaBeta(double alpha, double beta)
set parameters controlling arc-length equations alpha controls off diagonal therms beta controls diag...
Definition: ArcLengthTools.cpp:29
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
SphericalArcLengthControl::preProcess
MoFEMErrorCode preProcess()
Definition: ArcLengthTools.cpp:670
ArcLengthCtx::ghosTdLambda
SmartPetscObj< Vec > ghosTdLambda
Definition: ArcLengthTools.hpp:75
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::SmartPetscObj< Vec >
SimpleArcLengthControl::calculateDb
MoFEMErrorCode calculateDb()
Calculate db.
Definition: ArcLengthTools.cpp:597
ArcLengthCtx::s
double s
arc length radius
Definition: ArcLengthTools.hpp:71
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
ArcLengthCtx::F_lambda
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
Definition: ArcLengthTools.hpp:83
SphericalArcLengthControl::calculateDb
virtual MoFEMErrorCode calculateDb()
Calculate db.
Definition: ArcLengthTools.cpp:691
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
SimpleArcLengthControl::calculateLambdaInt
double calculateLambdaInt()
Calculate internal lambda.
Definition: ArcLengthTools.cpp:593
Generated by Doxygen 1.8.17 and hosted at University of Glasgow