EshelbianPlasticity::OpTractionBc Struct Reference

#include <users_modules/eshelbian_plasticit/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::OpTractionBc:

Collaboration diagram for EshelbianPlasticity::OpTractionBc:

Public Member Functions
	OpTractionBc (std::string row_field, FeTractionBc &bc_fe)
MoFEMErrorCode	doWork (int side, EntityType type, EntData &data)

Protected Attributes
FeTractionBc &	bcFe
MatrixDouble	matPP
MatrixDouble	vecPv

Detailed Description

Definition at line 612 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

OpTractionBc()

EshelbianPlasticity::OpTractionBc::OpTractionBc ( std::string  row_field, FeTractionBc &  bc_fe  )

inline

Definition at line 613 of file EshelbianPlasticity.hpp.

      : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
        bcFe(bc_fe) {}

Member Function Documentation

doWork()

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

Definition at line 1083 of file EshelbianOperators.cpp.

                                                                            {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  auto t_normal = getFTensor1Normal();
  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
  FTensor::Tensor1<double, 3> t_unit_normal;
  t_unit_normal(i) = t_normal(i) / nrm2;
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = data.getN().size1();
  int nb_base_functions = data.getN().size2() / 3;
  double ts_t = getFEMethod()->ts_t;
 
  auto integrate_matrix = [&]() {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
    matPP.clear();
 
    auto t_row_base_fun = data.getFTensor1N<3>();
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double a = t_row_base_fun(i) * t_unit_normal(i);
        auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
        for (int cc = 0; cc != nb_dofs / 3; ++cc) {
          const double b = t_col_base_fun(i) * t_unit_normal(i);
          matPP(rr, cc) += t_w * a * b;
          ++t_col_base_fun;
        }
        ++t_row_base_fun;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
 
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto integrate_rhs = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    vecPv.resize(3, nb_dofs / 3, false);
    vecPv.clear();
 
    auto t_row_base_fun = data.getFTensor1N<3>();
    double ts_t = getFEMethod()->ts_t;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            vecPv(dd, rr) += t * bc.vals[dd];
        ++t_row_base_fun;
      }
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
      ++t_w;
    }
    MoFEMFunctionReturn(0);
  };
 
  auto integrate_rhs_cook = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    vecPv.resize(3, nb_dofs / 3, false);
    vecPv.clear();
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base_fun = data.getFTensor1N<3>();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      auto calc_tau = [](double y) {
        y -= 44;
        y /= (60 - 44);
        return -y * (y - 1) / 0.25;
      };
 
      const double tau = calc_tau(t_coords(1));
 
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            vecPv(dd, rr) += t * tau * bc.vals[dd];
        ++t_row_base_fun;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
      ++t_w;
      ++t_coords;
    }
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcFe.bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs) {
 
        CHKERR integrate_matrix();
        if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
          CHKERR integrate_rhs_cook(bc);
        else
          CHKERR integrate_rhs(bc);
 
        const auto info =
            lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
                         nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
        if (info > 0)
          SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                   "The leading minor of order %i is "
                   "not positive; definite;\nthe "
                   "solution could not be computed",
                   info);
 
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            for (int rr = 0; rr != nb_dofs / 3; ++rr)
              data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

bcFe

FeTractionBc& EshelbianPlasticity::OpTractionBc::bcFe

protected

Definition at line 619 of file EshelbianPlasticity.hpp.

matPP

MatrixDouble EshelbianPlasticity::OpTractionBc::matPP

protected

Definition at line 620 of file EshelbianPlasticity.hpp.

vecPv

MatrixDouble EshelbianPlasticity::OpTractionBc::vecPv

protected

Definition at line 621 of file EshelbianPlasticity.hpp.

---

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

• users_modules/eshelbian_plasticit/src/EshelbianPlasticity.hpp
• users_modules/eshelbian_plasticit/src/impl/EshelbianOperators.cpp