23int main(
int argc, 
char *argv[]) {
 
   25  const string default_options = 
"-ksp_type fgmres \n" 
   27                                 "-pc_factor_mat_solver_type mumps \n" 
   28                                 "-mat_mumps_icntl_20 0 \n" 
   32                                 "-snes_type newtonls \n" 
   33                                 "-snes_linesearch_type basic \n" 
   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);
 
   44      file << default_options;
 
   54    moab::Core mb_instance;
 
   55    moab::Interface &moab = mb_instance;
 
   57    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, 
MYPCOMM_INDEX);
 
   59        boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, 
false);
 
   61      pcomm = 
new ParallelComm(&moab, moab_comm_wrap->get_comm());
 
   63    PetscBool flg = PETSC_TRUE;
 
   67    if (flg != PETSC_TRUE) {
 
   68      SETERRQ(PETSC_COMM_SELF, 1, 
"*** ERROR -my_file (MESH FILE NEEDED)");
 
   74    if (flg != PETSC_TRUE) {
 
   80    PetscBool is_partitioned = PETSC_FALSE;
 
   82                               &is_partitioned, &flg);
 
   84    if (is_partitioned == PETSC_TRUE) {
 
   87      option = 
"PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=" 
   88               "PARALLEL_PARTITION;";
 
   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)
 
  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)
 
  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];
 
  118    CHKERR PetscPrintf(PETSC_COMM_WORLD,
 
  119                       "Start step %D and step_size = %6.4e\n", step,
 
  128    std::vector<BitRefLevel> bit_levels;
 
  134      problem_bit_level = bit_levels.back();
 
  155                                             "MESH_NODE_POSITIONS");
 
  164          "ELASTIC", 
"LAMBDA"); 
 
  168          "ELASTIC", 
"MESH_NODE_POSITIONS");
 
  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);
 
  219          CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
 
  220          CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
 
  226            meshset_fe_arc_length, 
"ARC_LENGTH", 
false);
 
  254          "NEUMANN_FE", 
"MESH_NODE_POSITIONS");
 
  260        CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, 
true);
 
  267        CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, 
true);
 
  279    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
 
  281    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
 
  285        m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, 
"SPATIAL_POSITION",
 
  286        "MESH_NODE_POSITIONS");
 
  301                    false, 
false, 
false);
 
  312    std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
 
  325                                                        "MESH_NODE_POSITIONS");
 
  326      CHKERR m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material, 0);
 
  332          1., 
"MESH_NODE_POSITIONS", 
"SPATIAL_POSITION");
 
  348    if (is_partitioned) {
 
  349      SETERRQ(PETSC_COMM_SELF, 1,
 
  350              "Not implemented, problem with arc-length force multiplayer");
 
  369        "ELASTIC_MECHANICS", 
COL, &
F);
 
  374        ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"ELASTIC_MECHANICS",
 
  377    boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
 
  381    CHKERR MatGetSize(Aij, &M, &
N);
 
  383    CHKERR MatGetLocalSize(Aij, &
m, &
n);
 
  384    boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
 
  388    CHKERR MatCreateShell(PETSC_COMM_WORLD, 
m, 
n, M, 
N, mat_ctx.get(),
 
  390    CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
 
  393    ArcLengthSnesCtx snes_ctx(m_field, 
"ELASTIC_MECHANICS", arc_ctx);
 
  399      CHKERR moab.get_entities_by_type(meshset, MBVERTEX, nodes, 
true);
 
  401      node_set.merge(nodes);
 
  403    PetscPrintf(PETSC_COMM_WORLD, 
"Nb. nodes in load path: %u\n",
 
  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);
 
  419    fe_neumann.
uSeF = 
true;
 
  429    boost::shared_ptr<FEMethod> my_dirichlet_bc =
 
  431            m_field, 
"SPATIAL_POSITION", Aij, 
D, 
F));
 
  433                               &(my_dirichlet_bc->problemPtr));
 
  437    struct AssembleRhsVectors : 
public FEMethod {
 
  439      boost::shared_ptr<ArcLengthCtx> arcPtr;
 
  442      AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
 
  444          : arcPtr(arc_ptr), nodeSet(node_set) {}
 
  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,
 
  458          CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, INSERT_VALUES,
 
  462          SETERRQ(PETSC_COMM_SELF, 1, 
"not implemented");
 
  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);
 
  479          SETERRQ(PETSC_COMM_SELF, 1, 
"not implemented");
 
  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());
 
  506    struct AddLambdaVectorToFInternal : 
public FEMethod {
 
  508      boost::shared_ptr<ArcLengthCtx> arcPtr;
 
  509      boost::shared_ptr<DirichletSpatialPositionsBc> bC;
 
  511      AddLambdaVectorToFInternal(boost::shared_ptr<ArcLengthCtx> &arc_ptr,
 
  512                                 boost::shared_ptr<FEMethod> &bc)
 
  514            bC(boost::shared_ptr<DirichletSpatialPositionsBc>(
 
  530          CHKERR VecGhostUpdateBegin(arcPtr->F_lambda, ADD_VALUES,
 
  532          CHKERR VecGhostUpdateEnd(arcPtr->F_lambda, ADD_VALUES,
 
  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);
 
  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",
 
  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());
 
  552          CHKERR VecNorm(snes_f, NORM_2, &fnorm);
 
  553          PetscPrintf(PETSC_COMM_WORLD, 
"\tfnorm = %6.4e\n", fnorm);
 
  556          SETERRQ(PETSC_COMM_SELF, 1, 
"not implemented");
 
  562    AssembleRhsVectors pre_post_method(arc_ctx, node_set);
 
  563    AddLambdaVectorToFInternal assemble_F_lambda(arc_ctx, my_dirichlet_bc);
 
  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);
 
  577    if (flg == PETSC_TRUE) {
 
  578      PetscReal atol, rtol, stol;
 
  579      PetscInt maxit, maxf;
 
  580      CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
 
  583      CHKERR SNESSetTolerances(snes, atol, rtol, stol, maxit, maxf);
 
  587    CHKERR SNESGetKSP(snes, &ksp);
 
  589    CHKERR KSPGetPC(ksp, &pc);
 
  590    boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
 
  592    CHKERR PCSetType(pc, PCSHELL);
 
  593    CHKERR PCShellSetContext(pc, pc_ctx.get());
 
  597    if (flg == PETSC_TRUE) {
 
  598      PetscReal rtol, atol, dtol;
 
  600      CHKERR KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
 
  601      atol = my_tol * 1e-2;
 
  603      CHKERR KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
 
  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(
 
  613    loops_to_do_Rhs.push_back(
 
  617    loops_to_do_Rhs.push_back(
 
  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));
 
  626      CHKERR edge_forces.at(fe_name_str)
 
  627          .addForce(
"SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
 
  629    for (boost::ptr_map<std::string, EdgeForce>::iterator eit =
 
  631         eit != edge_forces.end(); eit++) {
 
  632      loops_to_do_Rhs.push_back(
 
  637    boost::ptr_map<std::string, NodalForce> nodal_forces;
 
  639    nodal_forces.insert(fe_name_str, 
new NodalForce(m_field));
 
  642      CHKERR nodal_forces.at(fe_name_str)
 
  643          .addForce(
"SPATIAL_POSITION", arc_ctx->F_lambda, it->getMeshsetId());
 
  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(
 
  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);
 
  661        snes_ctx.getSetOperators();
 
  662    snes_ctx.getPreProcSetOperators().push_back(my_dirichlet_bc);
 
  663    loops_to_do_Mat.push_back(
 
  666    loops_to_do_Mat.push_back(
 
  669    loops_to_do_Mat.push_back(
 
  671    loops_to_do_Mat.push_back(
 
  673    snes_ctx.getPostProcSetOperators().push_back(my_dirichlet_bc);
 
  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);
 
  680    PetscScalar step_size_reduction;
 
  682                               &step_size_reduction, &flg);
 
  683    if (flg != PETSC_TRUE) {
 
  684      step_size_reduction = 1.;
 
  690    if (flg != PETSC_TRUE) {
 
  697    if (flg != PETSC_TRUE) {
 
  700    PetscScalar max_reduction = 10, min_reduction = 0.1;
 
  702                               &max_reduction, &flg);
 
  704                               &min_reduction, &flg);
 
  706    double gamma = 0.5, reduction = 1;
 
  709      step_size = step_size_reduction;
 
  711      reduction = step_size_reduction;
 
  714    double step_size0 = step_size;
 
  718          "ELASTIC_MECHANICS", 
"SPATIAL_POSITION", 
"X0_SPATIAL_POSITION", 
COL,
 
  719          arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
 
  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,
 
  725      CHKERR arc_ctx->setAlphaBeta(1, 0);
 
  727      CHKERR arc_ctx->setS(step_size);
 
  728      CHKERR arc_ctx->setAlphaBeta(0, 1);
 
  735    CHKERR VecDuplicate(arc_ctx->x0, &x00);
 
  736    bool converged_state = 
false;
 
  738    for (
int jj = 0; step < max_steps; step++, jj++) {
 
  741      CHKERR VecCopy(arc_ctx->x0, x00);
 
  745        CHKERR PetscPrintf(PETSC_COMM_WORLD, 
"Load Step %D step_size = %6.4e\n",
 
  747        CHKERR arc_ctx->setS(step_size);
 
  748        CHKERR arc_ctx->setAlphaBeta(0, 1);
 
  749        CHKERR VecCopy(
D, arc_ctx->x0);
 
  754      } 
else if (step == 2) {
 
  756        CHKERR arc_ctx->setAlphaBeta(1, 0);
 
  759        step_size0 = step_size;
 
  760        CHKERR arc_ctx->setS(step_size);
 
  761        double dlambda = arc_ctx->dLambda;
 
  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);
 
  775          step_size0 = step_size;
 
  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;
 
  785        CHKERR arc_ctx->setS(step_size);
 
  786        double dlambda = reduction * arc_ctx->dLambda;
 
  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);
 
  799      CHKERR SNESSolve(snes, PETSC_NULLPTR, 
D);
 
  801      CHKERR SNESGetIterationNumber(snes, &its);
 
  802      CHKERR PetscPrintf(PETSC_COMM_WORLD, 
"number of Newton iterations = %D\n",
 
  805      SNESConvergedReason reason;
 
  806      CHKERR SNESGetConvergedReason(snes, &reason);
 
  810        CHKERR VecCopy(x00, arc_ctx->x0);
 
  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,
 
  817        CHKERR arc_ctx->setAlphaBeta(1, 0);
 
  820        converged_state = 
false;
 
  826        if (step > 1 && converged_state) {
 
  828          reduction = pow((
double)its_d / (
double)(its + 1), gamma);
 
  829          if (step_size >= max_reduction * step_size0 && reduction > 1) {
 
  831          } 
else if (step_size <= min_reduction * step_size0 && reduction < 1) {
 
  834          CHKERR PetscPrintf(PETSC_COMM_WORLD, 
"reduction step_size = %6.4e\n",
 
  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;
 
  860        std::ostringstream o1;
 
  861        o1 << 
"out_" << step << 
".h5m";
 
  865      CHKERR pre_post_method.potsProcessLoadPath();
 
  875    CHKERR MatDestroy(&ShellAij);
 
  876    CHKERR SNESDestroy(&snes);