23 {
24
25 const string default_options = "-ksp_type fgmres \n"
26 "-pc_type lu \n"
27 "-pc_factor_mat_solver_type mumps \n"
28 "-mat_mumps_icntl_20 0 \n"
29 "-ksp_atol 1e-10 \n"
30 "-ksp_rtol 1e-10 \n"
31 "-snes_monitor \n"
32 "-snes_type newtonls \n"
33 "-snes_linesearch_type basic \n"
34 "-snes_max_it 100 \n"
35 "-snes_atol 1e-7 \n"
36 "-snes_rtol 1e-7 \n"
37 "-ts_monitor \n"
38 "-ts_type alpha \n";
39
40 string param_file = "param_file.petsc";
41 if (!static_cast<bool>(ifstream(param_file))) {
42 std::ofstream file(param_file.c_str(), std::ios::ate);
43 if (file.is_open()) {
44 file << default_options;
45 file.close();
46 }
47 }
48
50
51
52 try {
53
54 moab::Core mb_instance;
55 moab::Interface &moab = mb_instance;
56
57 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
58 auto moab_comm_wrap =
59 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
60 if (pcomm == NULL)
61 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
62
63 PetscBool flg = PETSC_TRUE;
67 if (flg != PETSC_TRUE) {
68 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
69 }
70
73 &flg);
74 if (flg != PETSC_TRUE) {
76 }
77
78
79
80 PetscBool is_partitioned = PETSC_FALSE;
82 &is_partitioned, &flg);
83
84 if (is_partitioned == PETSC_TRUE) {
85
86 const char *option;
87 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
88 "PARALLEL_PARTITION;";
90 } else {
91 const char *option;
92 option = "";
94 }
95
96
97 Tag th_step_size, th_step;
98 double def_step_size = 1;
99 CHKERR moab.tag_get_handle(
"_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
100 MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
101 if (
rval == MB_ALREADY_ALLOCATED)
103
104 int def_step = 1;
105 CHKERR moab.tag_get_handle(
"_STEP", 1, MB_TYPE_INTEGER, th_step,
106 MB_TAG_CREAT | MB_TAG_MESH, &def_step);
107 if (
rval == MB_ALREADY_ALLOCATED)
109
110 const void *tag_data_step_size[1];
112 CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
113 double &step_size = *(double *)tag_data_step_size[0];
114 const void *tag_data_step[1];
115 CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
116 int &step = *(int *)tag_data_step[0];
117
118 CHKERR PetscPrintf(PETSC_COMM_WORLD,
119 "Start step %D and step_size = %6.4e\n", step,
120 step_size);
121
124
125
128 std::vector<BitRefLevel> bit_levels;
131
132 if (step == 1) {
133
134 problem_bit_level = bit_levels.back();
135
136
138 3);
141
143
144
147
148
151
152
153
155 "MESH_NODE_POSITIONS");
156
157
159 "SPATIAL_POSITION");
162 "SPATIAL_POSITION");
164 "ELASTIC", "LAMBDA");
166 "SPATIAL_POSITION");
168 "ELASTIC", "MESH_NODE_POSITIONS");
170
171
173 "LAMBDA");
175 "LAMBDA");
176
178 "LAMBDA");
179
180
182
183
185 "ELASTIC");
187 "ARC_LENGTH");
188
190 "SPRING");
191
192
194 problem_bit_level);
195
196
200
201
202 {
203
205 {
206 const double coords[] = {0, 0, 0};
208 Range range_no_field_vertex;
209 range_no_field_vertex.insert(no_field_vertex);
214 range_no_field_vertex);
215 }
216
218 {
219 CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
220 CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
223 }
224
226 meshset_fe_arc_length, "ARC_LENGTH", false);
227 }
228
229
234
239
244
245
248 "SPATIAL_POSITION");
250 "SPATIAL_POSITION");
252 "SPATIAL_POSITION");
254 "NEUMANN_FE", "MESH_NODE_POSITIONS");
256 "NEUMANN_FE");
260 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
262 "NEUMANN_FE");
263 }
267 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
269 "NEUMANN_FE");
270 }
271
274 "FORCE_FE");
275 }
276
277
278
279 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
281 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
283
285 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
286 "MESH_NODE_POSITIONS");
287
288 PetscBool linear;
290 &linear);
291
294 CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
295 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
297 false, false);
299 false, false);
301 false, false, false);
302 CHKERR elastic.setOperators(
"SPATIAL_POSITION");
303
304
306 CHKERR post_proc.generateReferenceElementMesh();
308 false);
309 CHKERR post_proc.addFieldValuesPostProc(
"SPATIAL_POSITION");
310 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
311 CHKERR post_proc.addFieldValuesGradientPostProc(
"SPATIAL_POSITION");
312 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
313 elastic.setOfBlocks.begin();
314 for (; sit != elastic.setOfBlocks.end(); sit++) {
316 post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
317 sit->second, post_proc.commonData));
318 }
319
320
322 if (step == 1) {
323
325 "MESH_NODE_POSITIONS");
326 CHKERR m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material, 0);
328 "SPATIAL_POSITION");
330 "SPATIAL_POSITION");
332 1., "MESH_NODE_POSITIONS", "SPATIAL_POSITION");
334 "SPATIAL_POSITION");
336 "SPATIAL_POSITION");
337 }
338
339
341
342
344
347
348 if (is_partitioned) {
349 SETERRQ(PETSC_COMM_SELF, 1,
350 "Not implemented, problem with arc-length force multiplayer");
351 } else {
355 }
357
358
363
365
366
369 "ELASTIC_MECHANICS",
COL, &
F);
372 Mat Aij;
374 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
375 &Aij);
376
377 boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
379
381 CHKERR MatGetSize(Aij, &M, &
N);
383 CHKERR MatGetLocalSize(Aij, &
m, &
n);
384 boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
386
387 Mat ShellAij;
388 CHKERR MatCreateShell(PETSC_COMM_WORLD,
m,
n, M,
N, mat_ctx.get(),
389 &ShellAij);
390 CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
392
393 ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
394
399 CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes,
true);
401 node_set.merge(nodes);
402 }
403 PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
404 node_set.size());
405
407
408 double scaled_reference_load = 1;
409 double *scale_lhs = &(arc_ctx->getFieldData());
410 double *scale_rhs = &(scaled_reference_load);
412 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
414 neumann_forces.getLoopSpatialFe();
415 if (linear) {
418 }
419 fe_neumann.
uSeF =
true;
421 it)) {
423 }
427 }
428
429 boost::shared_ptr<FEMethod> my_dirichlet_bc =
431 m_field,
"SPATIAL_POSITION", Aij,
D,
F));
433 &(my_dirichlet_bc->problemPtr));
436
437 struct AssembleRhsVectors :
public FEMethod {
438
439 boost::shared_ptr<ArcLengthCtx> arcPtr;
441
442 AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
444 : arcPtr(arc_ptr), nodeSet(node_set) {}
445
448
449
450 switch (snes_ctx) {
452 CHKERR VecZeroEntries(snes_f);
453 CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
454 CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
455 CHKERR VecZeroEntries(arcPtr->F_lambda);
456 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, INSERT_VALUES,
457 SCATTER_FORWARD);
458 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
459 SCATTER_FORWARD);
460 } break;
461 default:
462 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
463 }
464
466 }
467
470 switch (snes_ctx) {
472
473 CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
474 CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
475 CHKERR VecAssemblyBegin(snes_f);
476 CHKERR VecAssemblyEnd(snes_f);
477 } break;
478 default:
479 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
480 }
482 }
483
486 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
487 problemPtr->getNumeredRowDofsPtr();
488 Range::iterator nit = nodeSet.begin();
489 for (; nit != nodeSet.end(); nit++) {
490 NumeredDofEntityByEnt::iterator dit, hi_dit;
491 dit = numered_dofs_rows->get<
Ent_mi_tag>().lower_bound(*nit);
492 hi_dit = numered_dofs_rows->get<
Ent_mi_tag>().upper_bound(*nit);
493 for (; dit != hi_dit; dit++) {
494 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ", "LAMBDA", 0,
495 arcPtr->getFieldData());
496 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e\n",
497 dit->get()->getName().c_str(),
498 dit->get()->getDofCoeffIdx(),
499 dit->get()->getFieldData());
500 }
501 }
503 }
504 };
505
506 struct AddLambdaVectorToFInternal :
public FEMethod {
507
508 boost::shared_ptr<ArcLengthCtx> arcPtr;
509 boost::shared_ptr<DirichletSpatialPositionsBc> bC;
510
511 AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
512 boost::shared_ptr<FEMethod> &bc)
513 : arcPtr(arc_ptr),
516
520 }
524 }
527 switch (snes_ctx) {
529
530 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
531 SCATTER_REVERSE);
532 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
533 SCATTER_REVERSE);
534 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
535 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
536 for (std::vector<int>::iterator vit = bC->dofsIndices.begin();
537 vit != bC->dofsIndices.end(); vit++) {
538 CHKERR VecSetValue(arcPtr->F_lambda, *vit, 0, INSERT_VALUES);
539 }
540 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
541 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
542 CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
543 PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n",
544 arcPtr->F_lambda2);
545
546 CHKERR VecAssemblyBegin(snes_f);
547 CHKERR VecAssemblyEnd(snes_f);
548 CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
549 PetscPrintf(PETSC_COMM_WORLD, "\tlambda = %6.4e\n",
550 arcPtr->getFieldData());
551 double fnorm;
552 CHKERR VecNorm(snes_f, NORM_2, &fnorm);
553 PetscPrintf(PETSC_COMM_WORLD, "\tfnorm = %6.4e\n", fnorm);
554 } break;
555 default:
556 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
557 }
559 }
560 };
561
562 AssembleRhsVectors pre_post_method(arc_ctx, node_set);
563 AddLambdaVectorToFInternal assemble_F_lambda(arc_ctx, my_dirichlet_bc);
564
565 SNES snes;
566 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
567 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
569 CHKERR SNESSetJacobian(snes, ShellAij, Aij,
SnesMat, &snes_ctx);
570 CHKERR SNESSetFromOptions(snes);
571
573
574 PetscReal my_tol;
576 &flg);
577 if (flg == PETSC_TRUE) {
578 PetscReal atol, rtol, stol;
579 PetscInt maxit, maxf;
580 CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
581 atol = my_tol;
582 rtol = atol * 1e2;
583 CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
584 }
585
586 KSP ksp;
587 CHKERR SNESGetKSP(snes, &ksp);
588 PC pc;
589 CHKERR KSPGetPC(ksp, &pc);
590 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
592 CHKERR PCSetType(pc, PCSHELL);
593 CHKERR PCShellSetContext(pc, pc_ctx.get());
596
597 if (flg == PETSC_TRUE) {
598 PetscReal rtol, atol, dtol;
599 PetscInt maxits;
600 CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
601 atol = my_tol * 1e-2;
602 rtol = atol * 1e-2;
603 CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
604 }
605
607 snes_ctx.getComputeRhs();
608 snes_ctx.getPreProcComputeRhs().push_back(my_dirichlet_bc);
609 snes_ctx.getPreProcComputeRhs().push_back(&pre_post_method);
610 loops_to_do_Rhs.push_back(
612
613 loops_to_do_Rhs.push_back(
615
616
617 loops_to_do_Rhs.push_back(
619
620
621 boost::ptr_map<std::string, EdgeForce> edge_forces;
622 string fe_name_str = "FORCE_FE";
623 edge_forces.insert(fe_name_str,
new EdgeForce(m_field));
625 it)) {
626 CHKERR edge_forces.at(fe_name_str)
627 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
628 }
629 for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
630 edge_forces.begin();
631 eit != edge_forces.end(); eit++) {
632 loops_to_do_Rhs.push_back(
634 }
635
636
637 boost::ptr_map<std::string, NodalForce> nodal_forces;
638
639 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
641 it)) {
642 CHKERR nodal_forces.at(fe_name_str)
643 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
644 }
645 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
646 nodal_forces.begin();
647 fit != nodal_forces.end(); fit++) {
648 loops_to_do_Rhs.push_back(
650 }
651
652
653 loops_to_do_Rhs.push_back(
655 loops_to_do_Rhs.push_back(
657 snes_ctx.getPostProcComputeRhs().push_back(&pre_post_method);
658 snes_ctx.getPostProcComputeRhs().push_back(my_dirichlet_bc);
659
661 snes_ctx.getSetOperators();
662 snes_ctx.getPreProcSetOperators().push_back(my_dirichlet_bc);
663 loops_to_do_Mat.push_back(
665
666 loops_to_do_Mat.push_back(
668
669 loops_to_do_Mat.push_back(
671 loops_to_do_Mat.push_back(
673 snes_ctx.getPostProcSetOperators().push_back(my_dirichlet_bc);
674
676 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
677 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
678 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
679
680 PetscScalar step_size_reduction;
682 &step_size_reduction, &flg);
683 if (flg != PETSC_TRUE) {
684 step_size_reduction = 1.;
685 }
686
687 PetscInt max_steps;
689 &flg);
690 if (flg != PETSC_TRUE) {
691 max_steps = 5;
692 }
693
694 int its_d;
696 &flg);
697 if (flg != PETSC_TRUE) {
698 its_d = 4;
699 }
700 PetscScalar max_reduction = 10, min_reduction = 0.1;
702 &max_reduction, &flg);
704 &min_reduction, &flg);
705
706 double gamma = 0.5, reduction = 1;
707
708 if (step == 1) {
709 step_size = step_size_reduction;
710 } else {
711 reduction = step_size_reduction;
712 step++;
713 }
714 double step_size0 = step_size;
715
716 if (step > 1) {
718 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"X0_SPATIAL_POSITION",
COL,
719 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
720 double x0_nrm;
721 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
722 CHKERR PetscPrintf(PETSC_COMM_WORLD,
723 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
724 arc_ctx->dLambda);
725 CHKERR arc_ctx->setAlphaBeta(1, 0);
726 } else {
727 CHKERR arc_ctx->setS(step_size);
728 CHKERR arc_ctx->setAlphaBeta(0, 1);
729 }
730
732
733 Vec D0, x00;
735 CHKERR VecDuplicate(arc_ctx->x0, &x00);
736 bool converged_state = false;
737
738 for (int jj = 0; step < max_steps; step++, jj++) {
739
741 CHKERR VecCopy(arc_ctx->x0, x00);
742
743 if (step == 1) {
744
745 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load Step %D step_size = %6.4e\n",
746 step, step_size);
747 CHKERR arc_ctx->setS(step_size);
748 CHKERR arc_ctx->setAlphaBeta(0, 1);
749 CHKERR VecCopy(
D, arc_ctx->x0);
750 double dlambda;
751 CHKERR arc_method.calculateInitDlambda(&dlambda);
752 CHKERR arc_method.setDlambdaToX(
D, dlambda);
753
754 } else if (step == 2) {
755
756 CHKERR arc_ctx->setAlphaBeta(1, 0);
757 CHKERR arc_method.calculateDxAndDlambda(
D);
758 step_size = std::sqrt(arc_method.calculateLambdaInt());
759 step_size0 = step_size;
760 CHKERR arc_ctx->setS(step_size);
761 double dlambda = arc_ctx->dLambda;
762 double dx_nrm;
763 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
764 CHKERR PetscPrintf(PETSC_COMM_WORLD,
765 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
766 "dx_nrm = %6.4e dx2 = %6.4e\n",
767 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
768 CHKERR VecCopy(
D, arc_ctx->x0);
769 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
770 CHKERR arc_method.setDlambdaToX(
D, dlambda);
771
772 } else {
773
774 if (jj == 0) {
775 step_size0 = step_size;
776 }
777
778 CHKERR arc_method.calculateDxAndDlambda(
D);
779 step_size *= reduction;
780 if (step_size > max_reduction * step_size0) {
781 step_size = max_reduction * step_size0;
782 } else if (step_size < min_reduction * step_size0) {
783 step_size = min_reduction * step_size0;
784 }
785 CHKERR arc_ctx->setS(step_size);
786 double dlambda = reduction * arc_ctx->dLambda;
787 double dx_nrm;
788 CHKERR VecScale(arc_ctx->dx, reduction);
789 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
790 CHKERR PetscPrintf(PETSC_COMM_WORLD,
791 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
792 "dx_nrm = %6.4e dx2 = %6.4e\n",
793 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
794 CHKERR VecCopy(
D, arc_ctx->x0);
795 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
796 CHKERR arc_method.setDlambdaToX(
D, dlambda);
797 }
798
799 CHKERR SNESSolve(snes, PETSC_NULLPTR,
D);
800 int its;
801 CHKERR SNESGetIterationNumber(snes, &its);
802 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"number of Newton iterations = %D\n",
803 its);
804
805 SNESConvergedReason reason;
806 CHKERR SNESGetConvergedReason(snes, &reason);
807 if (reason < 0) {
808
810 CHKERR VecCopy(x00, arc_ctx->x0);
811
812 double x0_nrm;
813 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
814 CHKERR PetscPrintf(PETSC_COMM_WORLD,
815 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
816 arc_ctx->dLambda);
817 CHKERR arc_ctx->setAlphaBeta(1, 0);
818
819 reduction = 0.1;
820 converged_state = false;
821
822 continue;
823
824 } else {
825
826 if (step > 1 && converged_state) {
827
828 reduction = pow((double)its_d / (double)(its + 1), gamma);
829 if (step_size >= max_reduction * step_size0 && reduction > 1) {
830 reduction = 1;
831 } else if (step_size <= min_reduction * step_size0 && reduction < 1) {
832 reduction = 1;
833 }
834 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"reduction step_size = %6.4e\n",
835 reduction);
836 }
837
838
840 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
842 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"X0_SPATIAL_POSITION",
COL,
843 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
844 converged_state = true;
845 }
846
847 if (step % 1 == 0) {
848
849
850
851
852
853
854
855
856
857
859 post_proc);
860 std::ostringstream o1;
861 o1 << "out_" << step << ".h5m";
862 CHKERR post_proc.writeFile(o1.str().c_str());
863 }
864
865 CHKERR pre_post_method.potsProcessLoadPath();
866 }
867
870
871
875 CHKERR MatDestroy(&ShellAij);
876 CHKERR SNESDestroy(&snes);
877 }
879
881
882 return 0;
883}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOFIELD
scalar or vector of scalars describe (no true field)
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_MESHSET(const EntityHandle meshset, const std::string &name, const bool recursive=false)=0
add MESHSET element to finite element database given by name
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
const double n
refractive index of diffusive medium
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
FTensor::Index< 'M', 3 > M
FTensor::Index< 'm', 3 > m
Store variables for ArcLength analysis.
shell matrix for arc-length method
Set Dirichlet boundary conditions on spatial displacements.
MoFEMErrorCode iNitialize()
Force on edges and lines.
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual moab::Interface & get_moab()=0
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
structure for User Loop Methods on finite elements
Matrix manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEM::FEMethodsSequence FEMethodsSequence
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
MoFEMErrorCode addPressure(int ms_id)
MoFEMErrorCode addForce(int ms_id)
NonLinear surface pressure element (obsolete implementation)
structure grouping operators and data used for calculation of nonlinear elastic element
structure for Arc Length pre-conditioner
Implementation of spherical arc-length method.