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