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));
  435        ->iNitialize();
  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 
  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 .
#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.
MoFEMErrorCode printForceSet() const
Print meshsets with force boundary conditions.
#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 printMaterialsSet() const
Print meshsets with material properties.
MoFEMErrorCode printDisplacementSet() const
Print meshsets with displacement boundary conditions.
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
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)
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.
Force on edges and lines.
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
virtual MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
virtual MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
virtual MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
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.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEM::FEMethodsSequence FEMethodsSequence
@ CTX_SNESSETFUNCTION
Setting up nonlinear function evaluation.
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.