EshelbianPlasticity::OpTopologicalObjectivePythonImpl Struct Reference

Inheritance diagram for EshelbianPlasticity::OpTopologicalObjectivePythonImpl:

Collaboration diagram for EshelbianPlasticity::OpTopologicalObjectivePythonImpl:

Public Types
using	OP = ForcesAndSourcesCore::UserDataOperator

Public Member Functions
	OpTopologicalObjectivePythonImpl (boost::shared_ptr< DataAtIntegrationPts > data_at_pts_ptr, boost::shared_ptr< TopologicalData > topo_p, boost::shared_ptr< ObjectiveFunctionData > python_ptr, const ObjectiveModelType eval_energy_model=PYTHON_MODEL)
MoFEMErrorCode	doWork (int side, EntityType type, EntData &data)

Private Attributes
ObjectiveModelType	evalEnergyModel = PYTHON_MODEL
boost::shared_ptr< DataAtIntegrationPts >	dataAtPts
boost::shared_ptr< TopologicalData >	topoData
boost::shared_ptr< ObjectiveFunctionData >	pythonPtr

Detailed Description

Definition at line 338 of file EshelbianTopologicalDerivativeOperators.cpp.

Member Typedef Documentation

OP

using EshelbianPlasticity::OpTopologicalObjectivePythonImpl::OP = ForcesAndSourcesCore::UserDataOperator

Definition at line 340 of file EshelbianTopologicalDerivativeOperators.cpp.

Constructor & Destructor Documentation

OpTopologicalObjectivePythonImpl()

EshelbianPlasticity::OpTopologicalObjectivePythonImpl::OpTopologicalObjectivePythonImpl ( boost::shared_ptr< DataAtIntegrationPts >  data_at_pts_ptr, boost::shared_ptr< TopologicalData >  topo_p, boost::shared_ptr< ObjectiveFunctionData >  python_ptr, const ObjectiveModelType  eval_energy_model = PYTHON_MODEL  )

inline

Definition at line 342 of file EshelbianTopologicalDerivativeOperators.cpp.

      : OP(NOSPACE, OP::OPSPACE), dataAtPts(data_at_pts_ptr), topoData(topo_p),
        pythonPtr(python_ptr), evalEnergyModel(eval_energy_model) {}

Member Function Documentation

doWork()

MoFEMErrorCode EshelbianPlasticity::OpTopologicalObjectivePythonImpl::doWork	(	int	side,
		EntityType	type,
		EntData &	data
	)

Definition at line 362 of file EshelbianTopologicalDerivativeOperators.cpp.

                                                                       {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!dataAtPts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "DataAtIntegrationPts pointer is null");
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_gauss_pts = getGaussPts().size2();
  if (!nb_gauss_pts)
    MoFEMFunctionReturnHot(0);
 
  auto stress_full_ptr = boost::make_shared<MatrixDouble>();
  auto get_stress_full =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(*stress_full_ptr, nb_gauss_pts);
  stress_full_ptr->clear();
  auto strain_full_ptr = boost::make_shared<MatrixDouble>();
  auto get_strain_full =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(*strain_full_ptr, nb_gauss_pts);
  strain_full_ptr->clear();
 
  auto t_stress = get_stress_full();
  auto t_strain = get_strain_full();
 
  auto t_biot = dataAtPts->getFTensorAdjointPdstretch(nb_gauss_pts);
  auto t_u = dataAtPts->getFTensorStretch(nb_gauss_pts);
 
  auto next = [&]() {
    ++t_stress;
    ++t_strain;
    ++t_biot;
    ++t_u;
  };
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
 
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
    // we have to handle all variants, that will render how the Jacobian
    // gradient is evaluated.
    t_stress(i, j) = t_biot(i, j);
    t_strain(i, j) = t_u(i, j);
    next();
  }
 
  auto evaluate_python_objective = [&]() {
    MoFEMFunctionBegin;
    if (!pythonPtr)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "ObjectiveFunctionData pointer is null");
 
    auto &coords = OP::getCoordsAtGaussPts();
    CHKERR pythonPtr->evalInteriorObjectiveFunction(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDStrainAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientU(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDDisplacementAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDStressAtPts(), false);
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_energy_of_hencky_model = [&]() {
    MoFEMFunctionBegin;
 
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDDisplacementAtPts(), nb_gauss_pts);
    topoData->getObjDDisplacementAtPts()->clear();
    MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                     DL>::size(*topoData->getObjDStressAtPts(), nb_gauss_pts);
    topoData->getObjDStressAtPts()->clear();
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
 
    auto eval_evergy = [&](auto &&t_D) {
      auto get_obj =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
              *topoData->getObjAtPts(), nb_gauss_pts);
      auto get_dstrain_obj =
          MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                           DL>::size(*topoData->getObjDStrainAtPts(),
                                     nb_gauss_pts);
      auto t_obj = get_obj();
      auto t_dstrain_obj = get_dstrain_obj();
      auto t_log_u = dataAtPts->getFTensorLogStretchTotal(nb_gauss_pts);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        t_obj(0) = 0.5 * (t_log_u(i, j) * t_D(i, j, k, l) * t_log_u(k, l));
        t_dstrain_obj(i, j) = t_D(i, j, k, l) * t_log_u(k, l);
        ++t_log_u;
        ++t_obj;
        ++t_dstrain_obj;
        ++t_D;
      }
    };
 
    if (!EshelbianCore::hasNonHomogeneousMaterialBlock) {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matD));
    } else {
      eval_evergy(getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matD));
    }
 
    topoData->getObjDRotationAtPts()->clear();
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_energy_of_hencky_model_nostreach = [&]() {
    MoFEMFunctionBegin;
 
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDDisplacementAtPts(), nb_gauss_pts);
    topoData->getObjDDisplacementAtPts()->clear();
    MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                     DL>::size(*topoData->getObjDStrainAtPts(), nb_gauss_pts);
    topoData->getObjDStrainAtPts()->clear();
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
    topoData->getObjDRotationAtPts()->clear();
 
    auto eval_evergy = [&](auto &&t_inv_D) {
      auto get_obj =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
              *topoData->getObjAtPts(), nb_gauss_pts);
      auto get_dstress_obj =
          MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                           DL>::size(*topoData->getObjDStressAtPts(),
                                     nb_gauss_pts);
      auto t_obj = get_obj();
      auto t_dstress_obj = get_dstress_obj();
      auto t_stress = dataAtPts->getFTensorAdjointPdstretch(nb_gauss_pts);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        t_obj(0) =
            0.5 * (t_stress(i, j) * t_inv_D(i, j, k, l) * t_stress(k, l));
        t_dstress_obj(i, j) = t_inv_D(i, j, k, l) * t_stress(k, l);
        ++t_stress;
        ++t_obj;
        ++t_dstress_obj;
        ++t_inv_D;
      }
    };
 
    if (!EshelbianCore::hasNonHomogeneousMaterialBlock &&
        EshelbianCore::stretchHandling !=
            EshelbianCore::NO_STREACH_NONLINEAR) {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD));
    } else {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matInvD));
    }
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_biot_stress = [&]() {
    MoFEMFunctionBegin;
    // Python differentiates the objective with respect to Biot stress. The
    // material equations need the corresponding Piola and rotation derivatives.
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
    topoData->getObjDRotationAtPts()->clear();
 
    auto t_obj_dbiot = topoData->getFTensorObjDStress(nb_gauss_pts);
    auto t_obj_domega = topoData->getFTensorObjDRotation(nb_gauss_pts);
    auto t_R = dataAtPts->getFTensorRotMat(nb_gauss_pts);
    auto t_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
    auto t_grad_h1 = dataAtPts->getFTensorSmallWGradH1(nb_gauss_pts);
    auto t_omega = dataAtPts->getFTensorRotAxis(nb_gauss_pts);
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
    FTENSOR_INDEX(SPACE_DIM, m);
    FTENSOR_INDEX(SPACE_DIM, o);
    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_dJ_dbiot;
      t_dJ_dbiot(l, o) = t_obj_dbiot(l, o);
 
      switch (EshelbianCore::gradApproximator) {
      case SMALL_ROT:
        // In the linear/small formulation Python sees the Piola field itself:
        // B = P. There is no rotation or H1-gradient pullback to apply.
        t_obj_dbiot(i, k) = t_dJ_dbiot(i, k);
        t_obj_domega(m) = 0;
        break;
      case NO_H1_CONFIGURATION: {
        FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_diff_R;
        switch (EshelbianCore::rotSelector) {
        case LARGE_ROT:
          t_diff_R(i, l, m) =
              LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, l, m);
          break;
        default:
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "rotationSelector not handled");
        }
 
        // Python sees B = R^T P.
        t_obj_dbiot(i, k) = t_R(i, l) * t_dJ_dbiot(l, k);
        t_obj_domega(m) = t_dJ_dbiot(l, k) * t_diff_R(i, l, m) * t_P(i, k);
      } break;
      case LARGE_ROT:
      case MODERATE_ROT: {
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_h1;
        t_h1(o, k) = t_kd(o, k) + t_grad_h1(o, k);
 
        FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_diff_R;
        switch (EshelbianCore::rotSelector) {
        case SMALL_ROT:
          t_diff_R(i, l, m) = levi_civita(i, l, m);
          break;
        case LARGE_ROT:
          t_diff_R(i, l, m) =
              LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, l, m);
          break;
        default:
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "rotationSelector not handled");
        }
 
        // Python sees B = R^T P H1^T.
        t_obj_dbiot(i, k) = t_R(i, l) * (t_dJ_dbiot(l, o) * t_h1(o, k));
        t_obj_domega(m) =
            t_dJ_dbiot(l, o) * (t_diff_R(i, l, m) * t_P(i, k)) * t_h1(o, k);
      } break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "gradApproximator not handled");
      }
 
      ++t_obj_dbiot;
      ++t_obj_domega;
      ++t_R;
      ++t_P;
      ++t_grad_h1;
      ++t_omega;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_stretch = [&]() {
    MoFEMFunctionBegin;
    // Python sees the physical stretch tensor. The stretch field stores log
    // stretch, so convert dJ/dU to dJ/dlogU once here.
    auto t_obj_dstretch = topoData->getFTensorObjDStrain(nb_gauss_pts);
    auto t_diff_stretch = dataAtPts->getFTensorDiffStretch(nb_gauss_pts);
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_dJ_dstretch;
      t_dJ_dstretch(i, j) = t_obj_dstretch(i, j);
 
      t_obj_dstretch(k, l) = t_dJ_dstretch(i, j) * t_diff_stretch(i, j, k, l);
 
      ++t_obj_dstretch;
      ++t_diff_stretch;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_stretch_to_stress_for_no_stretch = [&](auto t_inv_D) {
    MoFEMFunctionBegin;
    // In no-stretch mode, log stretch is computed from the stress field. Fold
    // the objective stretch derivative into the stress derivative buffer.
    auto t_obj_dstress = topoData->getFTensorObjDStress(nb_gauss_pts);
    auto t_obj_dstretch = topoData->getFTensorObjDStrain(nb_gauss_pts);
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_dstretch_dstress;
      t_dstretch_dstress(i, j) =
          ((t_obj_dstretch(k, l) || t_obj_dstretch(l, k)) / 2.) *
          t_inv_D(k, l, i, j);
 
      t_obj_dstress(i, j) += t_dstretch_dstress(i, j);
 
      ++t_obj_dstress;
      ++t_obj_dstretch;
      ++t_inv_D;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  switch (evalEnergyModel) {
  case PYTHON_MODEL:
    CHKERR evaluate_python_objective();
    CHKERR conversion_of_biot_stress();
    if (EshelbianCore::isNoStretch()) {
      if (!EshelbianCore::hasNonHomogeneousMaterialBlock &&
          EshelbianCore::stretchHandling !=
              EshelbianCore::NO_STREACH_NONLINEAR) {
        CHKERR conversion_of_stretch_to_stress_for_no_stretch(
            getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD));
      } else {
        CHKERR conversion_of_stretch_to_stress_for_no_stretch(
            getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matInvD));
      }
    } else {
      CHKERR conversion_of_stretch();
    }
    break;
  case HENCKY_MODEL:
    if (EshelbianCore::isNoStretch()) {
      CHKERR evaluate_energy_of_hencky_model_nostreach();
      CHKERR conversion_of_biot_stress();
    } else {
      CHKERR evaluate_energy_of_hencky_model();
    }
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Objective model type not handled");
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianPlasticity::OpTopologicalObjectivePythonImpl::dataAtPts

private

Definition at line 357 of file EshelbianTopologicalDerivativeOperators.cpp.

evalEnergyModel

ObjectiveModelType EshelbianPlasticity::OpTopologicalObjectivePythonImpl::evalEnergyModel = PYTHON_MODEL

private

Definition at line 356 of file EshelbianTopologicalDerivativeOperators.cpp.

pythonPtr

boost::shared_ptr<ObjectiveFunctionData> EshelbianPlasticity::OpTopologicalObjectivePythonImpl::pythonPtr

private

Definition at line 359 of file EshelbianTopologicalDerivativeOperators.cpp.

topoData

boost::shared_ptr<TopologicalData> EshelbianPlasticity::OpTopologicalObjectivePythonImpl::topoData

private

Definition at line 358 of file EshelbianTopologicalDerivativeOperators.cpp.

---

The documentation for this struct was generated from the following file:

• users_modules/eshelbian_plasticity/src/impl/EshelbianTopologicalDerivativeOperators.cpp