v0.14.0
Loading...
Searching...
No Matches
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>
10using 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
22MoFEMErrorCode 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
29MoFEMErrorCode 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
38ArcLengthCtx::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
133ArcLengthMatShell::ArcLengthMatShell(Mat aij, ArcLengthCtx *arc_ptr_raw,
134 string problem_name)
135 : Aij(aij, true), problemName(problem_name), arcPtrRaw(arc_ptr_raw) {}
136
137ArcLengthMatShell::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
143MoFEMErrorCode 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
184MoFEMErrorCode 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
203PCArcLengthCtx::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
212PCArcLengthCtx::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
222PCArcLengthCtx::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
233PCArcLengthCtx::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
243MoFEMErrorCode 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
326MoFEMErrorCode 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
359ZeroFLmabda::ZeroFLmabda(boost::shared_ptr<ArcLengthCtx> arc_ptr)
360 : arcPtr(arc_ptr) {}
361
362MoFEMErrorCode 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
395AssembleFlambda::AssembleFlambda(boost::shared_ptr<ArcLengthCtx> arc_ptr,
396 boost::shared_ptr<DirichletDisplacementBc> bc)
397 : arcPtr(arc_ptr) {
398 bCs.push_back(bc);
399}
400
401MoFEMErrorCode AssembleFlambda::preProcess() {
402 MoFEMFunctionBeginHot;
403 MoFEMFunctionReturnHot(0);
404}
405MoFEMErrorCode AssembleFlambda::operator()() {
406 MoFEMFunctionBeginHot;
407 MoFEMFunctionReturnHot(0);
408}
409
410MoFEMErrorCode 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
521SimpleArcLengthControl::SimpleArcLengthControl(
522 boost::shared_ptr<ArcLengthCtx> &arc_ptr, const bool assemble)
523 : FEMethod(), arcPtr(arc_ptr), aSsemble(assemble) {}
524
525SimpleArcLengthControl::~SimpleArcLengthControl() {}
526
527MoFEMErrorCode 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
550MoFEMErrorCode 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
569MoFEMErrorCode 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
593double SimpleArcLengthControl::calculateLambdaInt() {
594 return arcPtr->beta * arcPtr->dLambda;
595}
596
597MoFEMErrorCode 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
605MoFEMErrorCode 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
661SphericalArcLengthControl::SphericalArcLengthControl(ArcLengthCtx *arc_ptr_raw)
662 : FEMethod(), arcPtrRaw(arc_ptr_raw) {}
663
664SphericalArcLengthControl::SphericalArcLengthControl(
665 boost::shared_ptr<ArcLengthCtx> &arc_ptr)
666 : FEMethod(), arcPtrRaw(arc_ptr.get()), arcPtr(arc_ptr) {}
667
668SphericalArcLengthControl::~SphericalArcLengthControl() {}
669
670MoFEMErrorCode 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
685double SphericalArcLengthControl::calculateLambdaInt() {
686 return arcPtrRaw->alpha * arcPtrRaw->dx2 + pow(arcPtrRaw->dLambda, 2) *
687 pow(arcPtrRaw->beta, 2) *
688 arcPtrRaw->F_lambda2;
689}
690
691MoFEMErrorCode SphericalArcLengthControl::calculateDb() {
692 MoFEMFunctionBegin;
693 CHKERR VecCopy(arcPtrRaw->dx, arcPtrRaw->db);
694 CHKERR VecScale(arcPtrRaw->db, 2 * arcPtrRaw->alpha);
695 MoFEMFunctionReturn(0);
696}
697
698MoFEMErrorCode 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
720MoFEMErrorCode 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
742MoFEMErrorCode 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
768MoFEMErrorCode
769SphericalArcLengthControl::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
783MoFEMErrorCode 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}
PCApplyArcLength
MoFEMErrorCode PCApplyArcLength(PC pc, Vec pc_f, Vec pc_x)
Definition: ArcLengthTools.cpp:243
ArcLengthMatMultShellOp
MoFEMErrorCode ArcLengthMatMultShellOp(Mat A, Vec x, Vec f)
Definition: ArcLengthTools.cpp:184
PCSetupArcLength
MoFEMErrorCode PCSetupArcLength(PC pc)
Definition: ArcLengthTools.cpp:326
ArcFunctionBegin
#define ArcFunctionBegin
Definition: ArcLengthTools.cpp:13
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
ROW
@ ROW
Definition: definitions.h:123
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
MoFEM::createPC
auto createPC(MPI_Comm comm)
Definition: PetscSmartObj.hpp:263
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ArcLengthCtx
Store variables for ArcLength analysis.
Definition: ArcLengthTools.hpp:65
ArcLengthCtx::db
SmartPetscObj< Vec > db
db derivative of f(dx*dx), i.e. db = d[ f(dx*dx) ]/dx
Definition: ArcLengthTools.hpp:85
ArcLengthCtx::s
double s
arc length radius
Definition: ArcLengthTools.hpp:71
ArcLengthCtx::res_lambda
double res_lambda
f_lambda - s
Definition: ArcLengthTools.hpp:82
ArcLengthCtx::getPetscLocalDofIdx
DofIdx getPetscLocalDofIdx()
Get local index of load factor.
Definition: ArcLengthTools.hpp:113
ArcLengthCtx::dx
SmartPetscObj< Vec > dx
dx = x-x0
Definition: ArcLengthTools.hpp:88
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
ArcLengthCtx::ArcLengthCtx
ArcLengthCtx(MoFEM::Interface &m_field, const std::string &problem_name, const std::string &field_name="LAMBDA")
Definition: ArcLengthTools.cpp:38
ArcLengthCtx::alpha
double alpha
displacement scaling factor
Definition: ArcLengthTools.hpp:73
ArcLengthCtx::F_lambda
SmartPetscObj< Vec > F_lambda
F_lambda reference load vector.
Definition: ArcLengthTools.hpp:83
ArcLengthCtx::setS
MoFEMErrorCode setS(double s)
set arc radius
Definition: ArcLengthTools.cpp:22
ArcLengthCtx::mField
MoFEM::Interface & mField
Definition: ArcLengthTools.hpp:69
ArcLengthCtx::F_lambda2
double F_lambda2
inner_prod(F_lambda,F_lambda);
Definition: ArcLengthTools.hpp:81
ArcLengthCtx::dIag
double dIag
diagonal value
Definition: ArcLengthTools.hpp:78
ArcLengthCtx::dx2
double dx2
inner_prod(dX,dX)
Definition: ArcLengthTools.hpp:80
ArcLengthCtx::getFieldData
FieldData & getFieldData()
Get value of load factor.
Definition: ArcLengthTools.hpp:117
ArcLengthCtx::ghostDiag
SmartPetscObj< Vec > ghostDiag
Definition: ArcLengthTools.hpp:77
ArcLengthCtx::x0
SmartPetscObj< Vec > x0
displacement vector at beginning of step
Definition: ArcLengthTools.hpp:87
ArcLengthCtx::dLambda
double dLambda
increment of load factor
Definition: ArcLengthTools.hpp:76
ArcLengthCtx::beta
double beta
force scaling factor
Definition: ArcLengthTools.hpp:72
ArcLengthCtx::ghosTdLambda
SmartPetscObj< Vec > ghosTdLambda
Definition: ArcLengthTools.hpp:75
ArcLengthCtx::getPart
int getPart()
Get proc owning lambda dof.
Definition: ArcLengthTools.hpp:121
ArcLengthCtx::arcDofRawPtr
NumeredDofEntity * arcDofRawPtr
Definition: ArcLengthTools.hpp:124
ArcLengthCtx::getPetscGlobalDofIdx
DofIdx getPetscGlobalDofIdx()
Get global index of load factor.
Definition: ArcLengthTools.hpp:107
ArcLengthCtx::xLambda
SmartPetscObj< Vec > xLambda
solution of eq. K*xLambda = F_lambda
Definition: ArcLengthTools.hpp:86
ArcLengthMatShell
shell matrix for arc-length method
Definition: ArcLengthTools.hpp:205
ArcLengthMatShell::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:209
ArcLengthMatShell::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:207
ArcLengthMatShell::ArcLengthMatShell
ArcLengthMatShell(Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr, string problem_name)
Definition: ArcLengthTools.cpp:137
ArcLengthMatShell::setLambda
MoFEMErrorCode setLambda(Vec ksp_x, double *lambda, ScatterMode scattermode)
Definition: ArcLengthTools.cpp:143
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
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
MoFEM::NumeredDofEntity::getPart
unsigned int getPart() const
Definition: DofsMultiIndices.hpp:250
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition: PetscSmartObj.hpp:79
MoFEM::Tools
Auxiliary tools.
Definition: Tools.hpp:19
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
PCArcLengthCtx
structure for Arc Length pre-conditioner
Definition: ArcLengthTools.hpp:235
PCArcLengthCtx::arcPtrRaw
ArcLengthCtx * arcPtrRaw
Definition: ArcLengthTools.hpp:248
PCArcLengthCtx::Aij
SmartPetscObj< Mat > Aij
Definition: ArcLengthTools.hpp:240
PCArcLengthCtx::shellAij
SmartPetscObj< Mat > shellAij
Definition: ArcLengthTools.hpp:239
PCArcLengthCtx::PCArcLengthCtx
PCArcLengthCtx(Mat shell_Aij, Mat aij, boost::shared_ptr< ArcLengthCtx > arc_ptr)
Definition: ArcLengthTools.cpp:222
PCArcLengthCtx::kSP
SmartPetscObj< KSP > kSP
Definition: ArcLengthTools.hpp:237
PCArcLengthCtx::pC
SmartPetscObj< PC > pC
Definition: ArcLengthTools.hpp:238
SimpleArcLengthControl::preProcess
MoFEMErrorCode preProcess()
Definition: ArcLengthTools.cpp:527
SimpleArcLengthControl::calculateDb
MoFEMErrorCode calculateDb()
Calculate db.
Definition: ArcLengthTools.cpp:597
SimpleArcLengthControl::operator()
MoFEMErrorCode operator()()
Definition: ArcLengthTools.cpp:550
SimpleArcLengthControl::SimpleArcLengthControl
SimpleArcLengthControl(boost::shared_ptr< ArcLengthCtx > &arc_ptr, const bool assemble=false)
Definition: ArcLengthTools.cpp:521
SimpleArcLengthControl::calculateLambdaInt
double calculateLambdaInt()
Calculate internal lambda.
Definition: ArcLengthTools.cpp:593
SimpleArcLengthControl::calculateDxAndDlambda
MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: ArcLengthTools.cpp:605
SimpleArcLengthControl::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: ArcLengthTools.hpp:336
SimpleArcLengthControl::postProcess
MoFEMErrorCode postProcess()
Definition: ArcLengthTools.cpp:569
SphericalArcLengthControl::postProcess
MoFEMErrorCode postProcess()
Definition: ArcLengthTools.cpp:720
SphericalArcLengthControl::calculateDxAndDlambda
virtual MoFEMErrorCode calculateDxAndDlambda(Vec x)
Definition: ArcLengthTools.cpp:742
SphericalArcLengthControl::calculateDb
virtual MoFEMErrorCode calculateDb()
Calculate db.
Definition: ArcLengthTools.cpp:691
SphericalArcLengthControl::SphericalArcLengthControl
DEPRECATED SphericalArcLengthControl(ArcLengthCtx *arc_ptr_raw)
Definition: ArcLengthTools.cpp:661
SphericalArcLengthControl::calculateLambdaInt
virtual double calculateLambdaInt()
Calculate f_lambda(dx,lambda)
Definition: ArcLengthTools.cpp:685
SphericalArcLengthControl::calculateInitDlambda
virtual MoFEMErrorCode calculateInitDlambda(double *dlambda)
Definition: ArcLengthTools.cpp:769
SphericalArcLengthControl::setDlambdaToX
virtual MoFEMErrorCode setDlambdaToX(Vec x, double dlambda)
Definition: ArcLengthTools.cpp:783
ZeroFLmabda::ZeroFLmabda
ZeroFLmabda(boost::shared_ptr< ArcLengthCtx > arc_ptr)
Definition: ArcLengthTools.cpp:359
ZeroFLmabda::arcPtr
boost::shared_ptr< ArcLengthCtx > arcPtr
Definition: ArcLengthTools.hpp:285
ZeroFLmabda::preProcess
MoFEMErrorCode preProcess()
Definition: ArcLengthTools.cpp:362
Generated by Doxygen 1.9.5 and hosted at University of Glasgow