35 {
36
38 "-pc_type lu \n"
39 "-pc_factor_mat_solver_type mumps \n"
40 "-mat_mumps_icntl_20 0 \n"
41 "-ksp_atol 1e-10 \n"
42 "-ksp_rtol 1e-10 \n"
43 "-snes_monitor \n"
44 "-snes_type newtonls \n"
45 "-snes_linesearch_type basic \n"
46 "-snes_max_it 100 \n"
47 "-snes_atol 1e-7 \n"
48 "-snes_rtol 1e-7 \n"
49 "-ts_monitor \n"
50 "-ts_type alpha \n";
51
53 if (!
static_cast<bool>(ifstream(
param_file))) {
54 std::ofstream file(
param_file.c_str(), std::ios::ate);
55 if (file.is_open()) {
57 file.close();
58 }
59 }
60
62
63
64 try {
65
68
69 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
70 auto moab_comm_wrap =
71 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
72 if (pcomm == NULL)
73 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
74
75 PetscBool flg = PETSC_TRUE;
79 if (flg != PETSC_TRUE) {
80 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
81 }
82
85 &flg);
86 if (flg != PETSC_TRUE) {
88 }
89
90
91
92 PetscBool is_partitioned = PETSC_FALSE;
94 &is_partitioned, &flg);
95
96 if (is_partitioned == PETSC_TRUE) {
97
98 const char *option;
99 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
100 "PARALLEL_PARTITION;";
102 } else {
103 const char *option;
104 option = "";
106 }
107
108
109 Tag th_step_size, th_step;
110 double def_step_size = 1;
111 CHKERR moab.tag_get_handle(
"_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
112 MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
113 if (
rval == MB_ALREADY_ALLOCATED)
115
116 int def_step = 1;
117 CHKERR moab.tag_get_handle(
"_STEP", 1, MB_TYPE_INTEGER, th_step,
118 MB_TAG_CREAT | MB_TAG_MESH, &def_step);
119 if (
rval == MB_ALREADY_ALLOCATED)
121
122 const void *tag_data_step_size[1];
123 EntityHandle root = moab.get_root_set();
124 CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
125 double &step_size = *(double *)tag_data_step_size[0];
126 const void *tag_data_step[1];
127 CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
128 int &step = *(int *)tag_data_step[0];
129
130 CHKERR PetscPrintf(PETSC_COMM_WORLD,
131 "Start step %D and step_size = %6.4e\n", step,
132 step_size);
133
136
137
140 std::vector<BitRefLevel> bit_levels;
143
144 if (step == 1) {
145
146 problem_bit_level = bit_levels.back();
147
148
150 3);
153
155
156
159
160
163
164
165
167 "MESH_NODE_POSITIONS");
168
169
171 "SPATIAL_POSITION");
174 "SPATIAL_POSITION");
176 "ELASTIC", "LAMBDA");
178 "SPATIAL_POSITION");
180 "ELASTIC", "MESH_NODE_POSITIONS");
182
183
185 "LAMBDA");
187 "LAMBDA");
188
190 "LAMBDA");
191
192
194
195
197 "ELASTIC");
199 "ARC_LENGTH");
200
202 "SPRING");
203
204
206 problem_bit_level);
207
208
212
213
214 {
215
216 EntityHandle no_field_vertex;
217 {
218 const double coords[] = {0, 0, 0};
220 Range range_no_field_vertex;
221 range_no_field_vertex.insert(no_field_vertex);
226 range_no_field_vertex);
227 }
228
229 EntityHandle meshset_fe_arc_length;
230 {
231 CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
232 CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
235 }
236
238 meshset_fe_arc_length, "ARC_LENGTH", false);
239 }
240
241
246
251
256
257
260 "SPATIAL_POSITION");
262 "SPATIAL_POSITION");
264 "SPATIAL_POSITION");
266 "NEUMANN_FE", "MESH_NODE_POSITIONS");
268 "NEUMANN_FE");
271 Range tris;
272 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
274 "NEUMANN_FE");
275 }
278 Range tris;
279 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
281 "NEUMANN_FE");
282 }
283
286 "FORCE_FE");
287 }
288
289
290
291 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
293 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
295
297 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
298 "MESH_NODE_POSITIONS");
299
300 PetscBool linear;
302 &linear);
303
306 CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
307 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
309 false, false);
311 false, false);
313 false, false, false);
314 CHKERR elastic.setOperators(
"SPATIAL_POSITION");
315
316
318 CHKERR post_proc.generateReferenceElementMesh();
320 false);
321 CHKERR post_proc.addFieldValuesPostProc(
"SPATIAL_POSITION");
322 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
323 CHKERR post_proc.addFieldValuesGradientPostProc(
"SPATIAL_POSITION");
324 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
325 elastic.setOfBlocks.begin();
326 for (; sit != elastic.setOfBlocks.end(); sit++) {
328 post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
329 sit->second, post_proc.commonData));
330 }
331
332
334 if (step == 1) {
335
337 "MESH_NODE_POSITIONS");
338 CHKERR m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material, 0);
340 "SPATIAL_POSITION");
342 "SPATIAL_POSITION");
344 1., "MESH_NODE_POSITIONS", "SPATIAL_POSITION");
346 "SPATIAL_POSITION");
348 "SPATIAL_POSITION");
349 }
350
351
353
354
356
359
360 if (is_partitioned) {
361 SETERRQ(PETSC_COMM_SELF, 1,
362 "Not implemented, problem with arc-length force multiplayer");
363 } else {
367 }
369
370
375
377
378
381 "ELASTIC_MECHANICS",
COL, &F);
384 Mat Aij;
386 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
387 &Aij);
388
389 boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
391
395 CHKERR MatGetLocalSize(Aij, &
m, &
n);
396 boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
398
399 Mat ShellAij;
400 CHKERR MatCreateShell(PETSC_COMM_WORLD,
m,
n,
M,
N, mat_ctx.get(),
401 &ShellAij);
402 CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
404
405 ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
406
408
409 Range node_set;
411 EntityHandle meshset = cit->getMeshset();
412 Range nodes;
413 CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes,
true);
415 node_set.merge(nodes);
416 }
417 PetscPrintf(PETSC_COMM_WORLD, "Nb. nodes in load path: %u\n",
418 node_set.size());
419
421
422 double scaled_reference_load = 1;
423 double *scale_lhs = &(arc_ctx->getFieldData());
424 double *scale_rhs = &(scaled_reference_load);
426 m_field, Aij, arc_ctx->F_lambda, scale_lhs, scale_rhs);
428 neumann_forces.getLoopSpatialFe();
429 if (linear) {
432 }
433 fe_neumann.
uSeF =
true;
435 it)) {
437 }
441 }
442
443 boost::shared_ptr<FEMethod> my_dirichlet_bc =
445 m_field,
"SPATIAL_POSITION", Aij,
D, F));
447 &(my_dirichlet_bc->problemPtr));
449 ->iNitialize();
450
452
453 boost::shared_ptr<ArcLengthCtx> arcPtr;
454 Range &nodeSet;
455
457 Range &node_set)
458 : arcPtr(arc_ptr), nodeSet(node_set) {}
459
462
463
467 CHKERR VecGhostUpdateBegin(
snes_f, INSERT_VALUES, SCATTER_FORWARD);
468 CHKERR VecGhostUpdateEnd(
snes_f, INSERT_VALUES, SCATTER_FORWARD);
469 CHKERR VecZeroEntries(arcPtr->F_lambda);
470 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, INSERT_VALUES,
471 SCATTER_FORWARD);
472 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
473 SCATTER_FORWARD);
474 } break;
475 default:
476 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
477 }
478
480 }
481
486
487 CHKERR VecGhostUpdateBegin(
snes_f, ADD_VALUES, SCATTER_REVERSE);
488 CHKERR VecGhostUpdateEnd(
snes_f, ADD_VALUES, SCATTER_REVERSE);
491 } break;
492 default:
493 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
494 }
496 }
497
500 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
502 Range::iterator nit = nodeSet.begin();
503 for (; nit != nodeSet.end(); nit++) {
504 NumeredDofEntityByEnt::iterator dit, hi_dit;
505 dit = numered_dofs_rows->get<
Ent_mi_tag>().lower_bound(*nit);
506 hi_dit = numered_dofs_rows->get<
Ent_mi_tag>().upper_bound(*nit);
507 for (; dit != hi_dit; dit++) {
508 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e -> ", "LAMBDA", 0,
509 arcPtr->getFieldData());
510 PetscPrintf(PETSC_COMM_WORLD, "%s [ %d ] %6.4e\n",
511 dit->get()->getName().c_str(),
512 dit->get()->getDofCoeffIdx(),
513 dit->get()->getFieldData());
514 }
515 }
517 }
518 };
519
520 struct AddLambdaVectorToFInternal :
public FEMethod {
521
522 boost::shared_ptr<ArcLengthCtx> arcPtr;
523 boost::shared_ptr<DirichletSpatialPositionsBc> bC;
524
525 AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
526 boost::shared_ptr<FEMethod> &bc)
527 : arcPtr(arc_ptr),
530
534 }
538 }
541 switch (snes_ctx) {
543
544 CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
545 SCATTER_REVERSE);
546 CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
547 SCATTER_REVERSE);
548 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
549 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
550 for (std::vector<int>::iterator vit = bC->dofsIndices.begin();
551 vit != bC->dofsIndices.end(); vit++) {
552 CHKERR VecSetValue(arcPtr->F_lambda, *vit, 0, INSERT_VALUES);
553 }
554 CHKERR VecAssemblyBegin(arcPtr->F_lambda);
555 CHKERR VecAssemblyEnd(arcPtr->F_lambda);
556 CHKERR VecDot(arcPtr->F_lambda, arcPtr->F_lambda, &arcPtr->F_lambda2);
557 PetscPrintf(PETSC_COMM_WORLD, "\tFlambda2 = %6.4e\n",
558 arcPtr->F_lambda2);
559
560 CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
561 PetscPrintf(PETSC_COMM_WORLD, "\tlambda = %6.4e\n",
562 arcPtr->getFieldData());
563 double fnorm;
564 CHKERR VecNorm(snes_f, NORM_2, &fnorm);
565 PetscPrintf(PETSC_COMM_WORLD, "\tfnorm = %6.4e\n", fnorm);
566 } break;
567 default:
568 SETERRQ(PETSC_COMM_SELF, 1, "not implemented");
569 }
571 }
572 };
573
575 AddLambdaVectorToFInternal assemble_F_lambda(arc_ctx, my_dirichlet_bc);
576
577 SNES snes;
578 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
579 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
581 CHKERR SNESSetJacobian(snes, ShellAij, Aij,
SnesMat, &snes_ctx);
582 CHKERR SNESSetFromOptions(snes);
583
584 PetscReal my_tol;
586 &flg);
587 if (flg == PETSC_TRUE) {
588 PetscReal atol, rtol, stol;
589 PetscInt maxit, maxf;
590 CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
591 atol = my_tol;
592 rtol = atol * 1e2;
593 CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
594 }
595
596 KSP ksp;
597 CHKERR SNESGetKSP(snes, &ksp);
598 PC pc;
599 CHKERR KSPGetPC(ksp, &pc);
600 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
602 CHKERR PCSetType(pc, PCSHELL);
603 CHKERR PCShellSetContext(pc, pc_ctx.get());
606
607 if (flg == PETSC_TRUE) {
608 PetscReal rtol, atol, dtol;
609 PetscInt maxits;
610 CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
611 atol = my_tol * 1e-2;
612 rtol = atol * 1e-2;
613 CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
614 }
615
617 snes_ctx.get_loops_to_do_Rhs();
618 snes_ctx.get_preProcess_to_do_Rhs().push_back(my_dirichlet_bc);
619 snes_ctx.get_preProcess_to_do_Rhs().push_back(&pre_post_method);
620 loops_to_do_Rhs.push_back(
622
623 loops_to_do_Rhs.push_back(
625
626
627 loops_to_do_Rhs.push_back(
629
630
631 boost::ptr_map<std::string, EdgeForce> edge_forces;
632 string fe_name_str = "FORCE_FE";
633 edge_forces.insert(fe_name_str,
new EdgeForce(m_field));
635 it)) {
636 CHKERR edge_forces.at(fe_name_str)
637 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
638 }
639 for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
640 edge_forces.begin();
641 eit != edge_forces.end(); eit++) {
642 loops_to_do_Rhs.push_back(
644 }
645
646
647 boost::ptr_map<std::string, NodalForce> nodal_forces;
648
649 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
651 it)) {
652 CHKERR nodal_forces.at(fe_name_str)
653 .addForce("SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
654 }
655 for (boost::ptr_map<std::string, NodalForce>::iterator fit =
656 nodal_forces.begin();
657 fit != nodal_forces.end(); fit++) {
658 loops_to_do_Rhs.push_back(
660 }
661
662
663 loops_to_do_Rhs.push_back(
665 loops_to_do_Rhs.push_back(
667 snes_ctx.get_postProcess_to_do_Rhs().push_back(&pre_post_method);
668 snes_ctx.get_postProcess_to_do_Rhs().push_back(my_dirichlet_bc);
669
671 snes_ctx.get_loops_to_do_Mat();
672 snes_ctx.get_preProcess_to_do_Mat().push_back(my_dirichlet_bc);
673 loops_to_do_Mat.push_back(
675
676 loops_to_do_Mat.push_back(
678
679 loops_to_do_Mat.push_back(
681 loops_to_do_Mat.push_back(
683 snes_ctx.get_postProcess_to_do_Mat().push_back(my_dirichlet_bc);
684
686 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
687 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
688 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
689
690 PetscScalar step_size_reduction;
692 &step_size_reduction, &flg);
693 if (flg != PETSC_TRUE) {
694 step_size_reduction = 1.;
695 }
696
697 PetscInt max_steps;
699 &flg);
700 if (flg != PETSC_TRUE) {
701 max_steps = 5;
702 }
703
704 int its_d;
706 &flg);
707 if (flg != PETSC_TRUE) {
708 its_d = 4;
709 }
710 PetscScalar max_reduction = 10, min_reduction = 0.1;
712 &max_reduction, &flg);
714 &min_reduction, &flg);
715
716 double gamma = 0.5, reduction = 1;
717
718 if (step == 1) {
719 step_size = step_size_reduction;
720 } else {
721 reduction = step_size_reduction;
722 step++;
723 }
724 double step_size0 = step_size;
725
726 if (step > 1) {
728 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"X0_SPATIAL_POSITION",
COL,
729 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
730 double x0_nrm;
731 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
732 CHKERR PetscPrintf(PETSC_COMM_WORLD,
733 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
734 arc_ctx->dLambda);
735 CHKERR arc_ctx->setAlphaBeta(1, 0);
736 } else {
737 CHKERR arc_ctx->setS(step_size);
738 CHKERR arc_ctx->setAlphaBeta(0, 1);
739 }
740
742
745 CHKERR VecDuplicate(arc_ctx->x0, &x00);
746 bool converged_state = false;
747
748 for (int jj = 0; step < max_steps; step++, jj++) {
749
751 CHKERR VecCopy(arc_ctx->x0, x00);
752
753 if (step == 1) {
754
755 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load Step %D step_size = %6.4e\n",
756 step, step_size);
757 CHKERR arc_ctx->setS(step_size);
758 CHKERR arc_ctx->setAlphaBeta(0, 1);
759 CHKERR VecCopy(
D, arc_ctx->x0);
760 double dlambda;
761 CHKERR arc_method.calculateInitDlambda(&dlambda);
762 CHKERR arc_method.setDlambdaToX(
D, dlambda);
763
764 } else if (step == 2) {
765
766 CHKERR arc_ctx->setAlphaBeta(1, 0);
767 CHKERR arc_method.calculateDxAndDlambda(
D);
768 step_size = std::sqrt(arc_method.calculateLambdaInt());
769 step_size0 = step_size;
770 CHKERR arc_ctx->setS(step_size);
771 double dlambda = arc_ctx->dLambda;
772 double dx_nrm;
773 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
774 CHKERR PetscPrintf(PETSC_COMM_WORLD,
775 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
776 "dx_nrm = %6.4e dx2 = %6.4e\n",
777 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
778 CHKERR VecCopy(
D, arc_ctx->x0);
779 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
780 CHKERR arc_method.setDlambdaToX(
D, dlambda);
781
782 } else {
783
784 if (jj == 0) {
785 step_size0 = step_size;
786 }
787
788 CHKERR arc_method.calculateDxAndDlambda(
D);
789 step_size *= reduction;
790 if (step_size > max_reduction * step_size0) {
791 step_size = max_reduction * step_size0;
792 } else if (step_size < min_reduction * step_size0) {
793 step_size = min_reduction * step_size0;
794 }
795 CHKERR arc_ctx->setS(step_size);
796 double dlambda = reduction * arc_ctx->dLambda;
797 double dx_nrm;
798 CHKERR VecScale(arc_ctx->dx, reduction);
799 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
800 CHKERR PetscPrintf(PETSC_COMM_WORLD,
801 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
802 "dx_nrm = %6.4e dx2 = %6.4e\n",
803 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
804 CHKERR VecCopy(
D, arc_ctx->x0);
805 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
806 CHKERR arc_method.setDlambdaToX(
D, dlambda);
807 }
808
809 CHKERR SNESSolve(snes, PETSC_NULL,
D);
810 int its;
811 CHKERR SNESGetIterationNumber(snes, &its);
812 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"number of Newton iterations = %D\n",
813 its);
814
815 SNESConvergedReason reason;
816 CHKERR SNESGetConvergedReason(snes, &reason);
817 if (reason < 0) {
818
820 CHKERR VecCopy(x00, arc_ctx->x0);
821
822 double x0_nrm;
823 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
824 CHKERR PetscPrintf(PETSC_COMM_WORLD,
825 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
826 arc_ctx->dLambda);
827 CHKERR arc_ctx->setAlphaBeta(1, 0);
828
829 reduction = 0.1;
830 converged_state = false;
831
832 continue;
833
834 } else {
835
836 if (step > 1 && converged_state) {
837
838 reduction = pow((double)its_d / (double)(its + 1), gamma);
839 if (step_size >= max_reduction * step_size0 && reduction > 1) {
840 reduction = 1;
841 } else if (step_size <= min_reduction * step_size0 && reduction < 1) {
842 reduction = 1;
843 }
844 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"reduction step_size = %6.4e\n",
845 reduction);
846 }
847
848
850 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
852 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"X0_SPATIAL_POSITION",
COL,
853 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
854 converged_state = true;
855 }
856
857 if (step % 1 == 0) {
858
859
860
861
862
863
864
865
866
867
869 post_proc);
870 std::ostringstream o1;
871 o1 << "out_" << step << ".h5m";
872 CHKERR post_proc.writeFile(o1.str().c_str());
873 }
874
875 CHKERR pre_post_method.potsProcessLoadPath();
876 }
877
880
881
885 CHKERR MatDestroy(&ShellAij);
886 CHKERR SNESDestroy(&snes);
887 }
889
891
892 return 0;
893}
const std::string default_options
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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 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 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 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 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 modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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 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 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
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
const FTensor::Tensor2< T, Dim, Dim > Vec
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)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
FTensor::Index< 'm', 3 > m
Store variables for ArcLength analysis.
shell matrix for arc-length method
PetscErrorCode preProcess()
function is run at the beginning of loop
AssembleRhsVectors(ArcLengthCtx *arc_ptr)
PetscErrorCode postProcess()
function is run at the end of loop
Set Dirichlet boundary conditions on spatial displacements.
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
const Problem * problemPtr
raw pointer to problem
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
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
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 refernce 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.