#include "tutorials/cor-12_cohesive_interface/src/InterfaceGapArcLengthControl.hpp"

Inheritance diagram for CohesiveElement::ArcLengthIntElemFEMethod:

Collaboration diagram for CohesiveElement::ArcLengthIntElemFEMethod:

Public Member Functions
	ArcLengthIntElemFEMethod (moab::Interface &moab, boost::shared_ptr< ArcLengthCtx > arcptr)

	~ArcLengthIntElemFEMethod ()

MoFEMErrorCode	remove_damaged_prisms_nodes ()
	remove nodes of prims which are fully damaged

MoFEMErrorCode	preProcess ()

MoFEMErrorCode	calculate_lambda_int (double &_lambda_int_)

virtual MoFEMErrorCode	calculate_db ()

MoFEMErrorCode	operator() ()

MoFEMErrorCode	postProcess ()

MoFEMErrorCode	calculate_dx_and_dlambda (Vec &x)

MoFEMErrorCode	calculate_init_dlambda (double *dlambda)

MoFEMErrorCode	set_dlambda_to_x (Vec &x, double dlambda)

Public Attributes
moab::Interface &	mOab

boost::shared_ptr< ArcLengthCtx >	arcPtr

Vec	GhostLambdaInt

Range	Faces3

Range	Faces4

Range	Edges3

Range	Edges4

Range	Nodes3

Range	Nodes4

Tag	thDamagedPrism

double	lambda_int

Detailed Description

Definition at line 10 of file InterfaceGapArcLengthControl.hpp.

Constructor & Destructor Documentation

◆ ArcLengthIntElemFEMethod()

CohesiveElement::ArcLengthIntElemFEMethod::ArcLengthIntElemFEMethod	(	moab::Interface &	moab,
		boost::shared_ptr< ArcLengthCtx >	arcptr
	)

inline

Definition at line 22 of file InterfaceGapArcLengthControl.hpp.

      : FEMethod(), mOab(moab), arcPtr(arcptr) {
    PetscInt ghosts[1] = {0};
    ParallelComm *pcomm = ParallelComm::get_pcomm(&mOab, MYPCOMM_INDEX);
    if (pcomm->rank() == 0) {
      CHKERRABORT(PETSC_COMM_WORLD, VecCreateGhost(PETSC_COMM_WORLD, 1, 1, 0,
                                                   ghosts, &GhostLambdaInt));
    } else {
      CHKERRABORT(PETSC_COMM_WORLD, VecCreateGhost(PETSC_COMM_WORLD, 0, 1, 1,
                                                   ghosts, &GhostLambdaInt));
    }
    Range prisms;
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_entities_by_type(0, MBPRISM, prisms, false));
    for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
      EntityHandle f3, f4;
      CHKERRABORT(PETSC_COMM_WORLD, mOab.side_element(*pit, 2, 3, f3));
      CHKERRABORT(PETSC_COMM_WORLD, mOab.side_element(*pit, 2, 4, f4));
      Faces3.insert(f3);
      Faces4.insert(f4);
    }
 
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_adjacencies(Faces3, 1, false, Edges3));
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_adjacencies(Faces4, 1, false, Edges4));
    CHKERRABORT(PETSC_COMM_WORLD, mOab.get_connectivity(Faces3, Nodes3, true));
    CHKERRABORT(PETSC_COMM_WORLD, mOab.get_connectivity(Faces4, Nodes4, true));
    //Faces3.insert(Edges3.begin(),Edges3.end());
    Faces3.insert(Nodes3.begin(), Nodes3.end());
    //Faces4.insert(Edges4.begin(),Edges4.end());
    Faces4.insert(Nodes4.begin(), Nodes4.end());
 
    double def_damaged = 0;
    CHKERRABORT(
        PETSC_COMM_WORLD,
        mOab.tag_get_handle("DAMAGED_PRISM", 1, MB_TYPE_INTEGER, thDamagedPrism,
                            MB_TAG_CREAT | MB_TAG_SPARSE, &def_damaged));
  }

◆ ~ArcLengthIntElemFEMethod()

CohesiveElement::ArcLengthIntElemFEMethod::~ArcLengthIntElemFEMethod ( )

inline

Definition at line 62 of file InterfaceGapArcLengthControl.hpp.

                              {
    CHKERRABORT(PETSC_COMM_WORLD, VecDestroy(&GhostLambdaInt));
  }

Member Function Documentation

◆ calculate_db()

virtual MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::calculate_db ( )

inlinevirtual

Definition at line 149 of file InterfaceGapArcLengthControl.hpp.

                                        {
    MoFEMFunctionBeginHot;
    CHKERR VecZeroEntries(arcPtr->db);
    CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
    NumeredDofEntityByLocalIdx::iterator dit, hi_dit;
    dit = problemPtr->getNumeredRowDofsPtr()
              ->get<PetscLocalIdx_mi_tag>()
              .lower_bound(0);
    hi_dit = problemPtr->getNumeredRowDofsPtr()
                 ->get<PetscLocalIdx_mi_tag>()
                 .upper_bound(problemPtr->getNbLocalDofsRow() +
                              problemPtr->getNbGhostDofsRow());
    double *array;
    CHKERR VecGetArray(arcPtr->db, &array);
    for (; dit != hi_dit; dit++) {
      if (dit->get()->getEntType() != MBVERTEX) {
        array[dit->get()->getPetscLocalDofIdx()] = 0;
        continue;
      }
      if (Nodes3.find(dit->get()->getEnt()) != Nodes3.end()) {
        array[dit->get()->getPetscLocalDofIdx()] = +arcPtr->alpha;
      } else if (Nodes4.find(dit->get()->getEnt()) != Nodes4.end()) {
        array[dit->get()->getPetscLocalDofIdx()] = -arcPtr->alpha;
      } else
        array[dit->get()->getPetscLocalDofIdx()] = 0;
    }
    CHKERR VecRestoreArray(arcPtr->db, &array);
    MoFEMFunctionReturnHot(0);
  }

◆ calculate_dx_and_dlambda()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::calculate_dx_and_dlambda ( Vec & x )

inline

Definition at line 225 of file InterfaceGapArcLengthControl.hpp.

                                                  {
    MoFEMFunctionBeginHot;
    //dx
    CHKERR VecCopy(x, arcPtr->dx);
    CHKERR VecAXPY(arcPtr->dx, -1, arcPtr->x0);
    //if LAMBDA dof is on this partition
    if (arcPtr->getPetscLocalDofIdx() != -1) {
      double *array;
      CHKERR VecGetArray(arcPtr->dx, &array);
      arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
      array[arcPtr->getPetscLocalDofIdx()] = 0;
      CHKERR VecRestoreArray(arcPtr->dx, &array);
    }
    //brodcast dlambda
    CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
                             SCATTER_FORWARD);
    //calculate dx2 (dot product)
    CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                "\tdlambda = %6.4e dx2 = %6.4e", arcPtr->dLambda,
                arcPtr->dx2);
    MoFEMFunctionReturnHot(0);
  }

◆ calculate_init_dlambda()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::calculate_init_dlambda ( double * dlambda )

inline

Definition at line 251 of file InterfaceGapArcLengthControl.hpp.

                                                         {
    MoFEMFunctionBeginHot;
 
    *dlambda = arcPtr->s / (arcPtr->beta * sqrt(arcPtr->F_lambda2));
    MOFEM_LOG_C(
        "ARC_LENGTH", Sev::inform,
        "\tInit dlambda = %6.4e s = %6.4e beta = %6.4e F_lambda2 = %6.4e",
        *dlambda, arcPtr->s, arcPtr->beta, arcPtr->F_lambda2);
    double a = *dlambda;
    if (a - a != 0) {
      std::ostringstream sss;
      sss << "s " << arcPtr->s << " " << arcPtr->beta << " "
          << arcPtr->F_lambda2;
      SETERRQ(PETSC_COMM_SELF, 1, sss.str().c_str());
    }
 
    MoFEMFunctionReturnHot(0);
  }

◆ calculate_lambda_int()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::calculate_lambda_int ( double & _lambda_int_ )

inline

Definition at line 106 of file InterfaceGapArcLengthControl.hpp.

                                                            {
    MoFEMFunctionBeginHot;
    ParallelComm *pcomm = ParallelComm::get_pcomm(&mOab, MYPCOMM_INDEX);
    NumeredDofEntityByLocalIdx::iterator dit, hi_dit;
    dit = problemPtr->getNumeredRowDofsPtr()
              ->get<PetscLocalIdx_mi_tag>()
              .lower_bound(0);
    hi_dit = problemPtr->getNumeredRowDofsPtr()
                 ->get<PetscLocalIdx_mi_tag>()
                 .upper_bound(problemPtr->getNbLocalDofsRow());
    double *array;
    double *array_int_lambda;
    CHKERR VecZeroEntries(GhostLambdaInt);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGetArray(arcPtr->dx, &array);
    CHKERR VecGetArray(GhostLambdaInt, &array_int_lambda);
    array_int_lambda[0] = 0;
    for (; dit != hi_dit; dit++) {
      if (dit->get()->getEntType() != MBVERTEX)
        continue;
      if (pcomm->rank() != dit->get()->getPart())
        continue;
      if (Nodes3.find(dit->get()->getEnt()) != Nodes3.end()) {
        array_int_lambda[0] += array[dit->get()->getPetscLocalDofIdx()];
      }
      if (Nodes4.find(dit->get()->getEnt()) != Nodes4.end()) {
        array_int_lambda[0] -= array[dit->get()->getPetscLocalDofIdx()];
      }
    }
    CHKERR VecRestoreArray(arcPtr->dx, &array);
    CHKERR VecRestoreArray(GhostLambdaInt, &array_int_lambda);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGetArray(GhostLambdaInt, &array_int_lambda);
    _lambda_int_ = arcPtr->alpha * array_int_lambda[0] +
                   arcPtr->dLambda * arcPtr->beta * sqrt(arcPtr->F_lambda2);
    CHKERR VecRestoreArray(GhostLambdaInt, &array_int_lambda);
    MoFEMFunctionReturnHot(0);
  }

◆ operator()()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::operator() ( )

inline

Definition at line 180 of file InterfaceGapArcLengthControl.hpp.

                              {
    MoFEMFunctionBeginHot;
 
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      //calculate residual for arc length row
      arcPtr->res_lambda = lambda_int - arcPtr->s;
      CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
                         arcPtr->res_lambda, ADD_VALUES);
      MOFEM_LOG_C(
          "ARC_LENGTH", Sev::inform,
          "\tres_lambda = %6.4e lambda_int = %6.4e s = %6.4e",
          arcPtr->res_lambda, lambda_int, arcPtr->s);
    } break;
    case CTX_SNESSETJACOBIAN: {
      //calculate diagonal therm
      arcPtr->dIag = arcPtr->beta * sqrt(arcPtr->F_lambda2);
      CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
                         arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
    } break;
    default:
      break;
    }
 
    MoFEMFunctionReturnHot(0);
  }

◆ postProcess()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::postProcess ( )

inline

Definition at line 207 of file InterfaceGapArcLengthControl.hpp.

                               {
    MoFEMFunctionBeginHot;
    switch (snes_ctx) {
    case CTX_SNESSETJACOBIAN: {
      CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES,
                               SCATTER_FORWARD);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tdiag = %6.4e", arcPtr->dIag);
      CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
      CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
    } break;
    default:
      break;
    }
    MoFEMFunctionReturnHot(0);
  }

◆ preProcess()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::preProcess ( )

inline

Definition at line 92 of file InterfaceGapArcLengthControl.hpp.

                              {
    MoFEMFunctionBeginHot;
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      CHKERR calculate_dx_and_dlambda(snes_x);
      CHKERR calculate_db();
      CHKERR calculate_lambda_int(lambda_int);
    } break;
    default:
      break;
    }
    MoFEMFunctionReturnHot(0);
  }

◆ remove_damaged_prisms_nodes()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::remove_damaged_prisms_nodes ( )

inline

remove nodes of prims which are fully damaged

Definition at line 69 of file InterfaceGapArcLengthControl.hpp.

                                               {
    MoFEMFunctionBeginHot;
    Range prisms;
    CHKERR mOab.get_entities_by_type(0, MBPRISM, prisms, false);
    std::vector<int> is_prism_damaged(prisms.size());
    CHKERR mOab.tag_get_data(thDamagedPrism, prisms,
                             &*is_prism_damaged.begin());
    Range::iterator pit = prisms.begin();
    std::vector<int>::iterator vit = is_prism_damaged.begin();
    for (; pit != prisms.end(); pit++, vit++) {
      if (*vit > 0) {
        Range nodes;
        CHKERR mOab.get_connectivity(&*pit, 1, nodes, true);
        for (Range::iterator nit = nodes.begin(); nit != nodes.end(); nit++) {
          Faces3.erase(*nit);
          Faces4.erase(*nit);
        }
      }
    }
    MoFEMFunctionReturnHot(0);
  }

◆ set_dlambda_to_x()

MoFEMErrorCode CohesiveElement::ArcLengthIntElemFEMethod::set_dlambda_to_x	(	Vec &	x,
		double	dlambda
	)

inline

Definition at line 270 of file InterfaceGapArcLengthControl.hpp.

                                                          {
    MoFEMFunctionBeginHot;
 
    if (arcPtr->getPetscLocalDofIdx() != -1) {
      double *array;
      CHKERR VecGetArray(x, &array);
      double lambda_old = array[arcPtr->getPetscLocalDofIdx()];
      if (!(dlambda == dlambda)) {
        std::ostringstream sss;
        sss << "s " << arcPtr->s << " " << arcPtr->beta << " "
            << arcPtr->F_lambda2;
        SETERRQ(PETSC_COMM_SELF, 1, sss.str().c_str());
      }
      array[arcPtr->getPetscLocalDofIdx()] = lambda_old + dlambda;
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "\tlambda = %6.4e, %6.4e (%6.4e)", lambda_old,
                  array[arcPtr->getPetscLocalDofIdx()], dlambda);
      CHKERR VecRestoreArray(x, &array);
    }
 
    MoFEMFunctionReturnHot(0);
  }