HenckyOps Namespace Reference

Classes
struct	CommonData
struct	HenkyIntegrators
struct	isEq
struct	OpCalculateEigenValsImpl
struct	OpCalculateEigenValsImpl< DIM, GAUSS, DomainEleOp >
struct	OpCalculateHenckyPlasticStressImpl
struct	OpCalculateHenckyPlasticStressImpl< DIM, GAUSS, DomainEleOp >
struct	OpCalculateHenckyStressImpl
struct	OpCalculateHenckyStressImpl< DIM, GAUSS, DomainEleOp >
struct	OpCalculateLogC_dCImpl
struct	OpCalculateLogC_dCImpl< DIM, GAUSS, DomainEleOp >
struct	OpCalculateLogCImpl
struct	OpCalculateLogCImpl< DIM, GAUSS, DomainEleOp >
struct	OpCalculatePiolaStressImpl
struct	OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp >
struct	OpHenckyTangentImpl
struct	OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp >

Functions
auto	is_eq (const double &a, const double &b)
template<int DIM>
auto	get_uniq_nb (double *ptr)
template<int DIM>
auto	sort_eigen_vals (FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
template<int DIM>
MoFEMErrorCode	addMatBlockOps (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
template<int DIM, IntegrationType I, typename DomainEleOp >
auto	commonDataFactory (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode	opFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode	opFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode	opFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode	opFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)

Variables
constexpr double	eps = std::numeric_limits<float>::epsilon()
auto	f = [](double v) { return 0.5 * std::log(v); }
auto	d_f = [](double v) { return 0.5 / v; }
auto	dd_f = [](double v) { return -0.5 / (v * v); }

Function Documentation

addMatBlockOps()

template<int DIM>

MoFEMErrorCode HenckyOps::addMatBlockOps	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		std::string	block_name,
		boost::shared_ptr< MatrixDouble >	mat_D_Ptr,
		Sev	sev,
		double	scale = 1
	)

[Calculate elasticity tensor]

[Calculate elasticity tensor]

Definition at line 597 of file HenckyOps.hpp.

                                 {
  MoFEMFunctionBegin;
 
  struct OpMatBlocks : public DomainEleOp {
    OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
                double bulk_modulus_K, double shear_modulus_G,
                MoFEM::Interface &m_field, Sev sev,
                std::vector<const CubitMeshSets *> meshset_vec_ptr,
                double scale)
        : DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
          bulkModulusKDefault(bulk_modulus_K),
          shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
      std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
      CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                        "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      for (auto &b : blockData) {
 
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
                            b.shearModulusG * scaleYoungModulus);
          MoFEMFunctionReturnHot(0);
        }
      }
 
      CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
                        shearModulusGDefault * scaleYoungModulus);
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<MatrixDouble> matDPtr;
    const double scaleYoungModulus;
 
    struct BlockData {
      double bulkModulusK;
      double shearModulusG;
      Range blockEnts;
    };
 
    double bulkModulusKDefault;
    double shearModulusGDefault;
    std::vector<BlockData> blockData;
 
    MoFEMErrorCode
    extractBlockData(MoFEM::Interface &m_field,
                     std::vector<const CubitMeshSets *> meshset_vec_ptr,
                     Sev sev) {
      MoFEMFunctionBegin;
 
      for (auto m : meshset_vec_ptr) {
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
        std::vector<double> block_data;
        CHKERR m->getAttributes(block_data);
        if (block_data.size() != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Expected that block has two attribute");
        }
        auto get_block_ents = [&]() {
          Range ents;
          CHKERR
          m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
          return ents;
        };
 
        double young_modulus = block_data[0];
        double poisson_ratio = block_data[1];
        double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
        double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
            << "E = " << young_modulus << " nu = " << poisson_ratio;
 
        blockData.push_back(
            {bulk_modulus_K, shear_modulus_G, get_block_ents()});
      }
      MOFEM_LOG_CHANNEL("WORLD");
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
                              double bulk_modulus_K, double shear_modulus_G) {
      MoFEMFunctionBegin;
      //! [Calculate elasticity tensor]
      auto set_material_stiffness = [&]() {
        FTensor::Index<'i', DIM> i;
        FTensor::Index<'j', DIM> j;
        FTensor::Index<'k', DIM> k;
        FTensor::Index<'l', DIM> l;
        constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
        double A = (DIM == 2)
                       ? 2 * shear_modulus_G /
                             (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
                       : 1;
        auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
        t_D(i, j, k, l) =
            2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
            A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
                t_kd(k, l);
      };
      //! [Calculate elasticity tensor]
      constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
      mat_D_ptr->resize(size_symm * size_symm, 1);
      set_material_stiffness();
      MoFEMFunctionReturn(0);
    }
  };
 
  double E = young_modulus;
  double nu = poisson_ratio;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
  CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
                            PETSC_NULL);
  CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
                            PETSC_NULL);
  ierr = PetscOptionsEnd();
 
  double bulk_modulus_K = E / (3 * (1 - 2 * nu));
  double shear_modulus_G = E / (2 * (1 + nu));
  pip.push_back(new OpMatBlocks(
      field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
 
      // Get blockset using regular expression
      m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
          (boost::format("%s(.*)") % block_name).str()
 
              )),
      scale
 
      ));
 
  MoFEMFunctionReturn(0);
}

commonDataFactory()

template<int DIM, IntegrationType I, typename DomainEleOp >

auto HenckyOps::commonDataFactory	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		std::string	block_name,
		Sev	sev,
		double	scale = 1
	)

Definition at line 743 of file HenckyOps.hpp.

                                                                           {
 
  auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
  common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
  common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
 
  CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, field_name, block_name,
                                        common_ptr->matDPtr, sev, scale),
                    "addMatBlockOps");
 
  using H = HenkyIntegrators<DomainEleOp>;
 
  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
      field_name, common_ptr->matGradPtr));
  pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
      field_name, common_ptr));
  pip.push_back(
      new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
  pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
      field_name, common_ptr));
  pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I>(
      field_name, common_ptr));
  pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
      field_name, common_ptr));
 
  return common_ptr;
}

get_uniq_nb()

template<int DIM>

auto HenckyOps::get_uniq_nb ( double *  ptr )

inline

Definition at line 32 of file HenckyOps.hpp.

                                                        {
  std::array<double, DIM> tmp;
  std::copy(ptr, &ptr[DIM], tmp.begin());
  std::sort(tmp.begin(), tmp.end());
  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
};

is_eq()

auto HenckyOps::is_eq ( const double &  a, const double &  b  )

inline

Definition at line 28 of file HenckyOps.hpp.

                                                    {
  return isEq::check(a, b);
};

opFactoryDomainLhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >

MoFEMErrorCode HenckyOps::opFactoryDomainLhs	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		boost::shared_ptr< HenckyOps::CommonData >	common_ptr,
		Sev	sev
	)

Definition at line 808 of file HenckyOps.hpp.

             {
  MoFEMFunctionBegin;
 
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template BiLinearForm<I>;
  using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
 
  using H = HenkyIntegrators<DomainEleOp>;
  pip.push_back(
      new typename H::template OpHenckyTangent<DIM, I>(field_name, common_ptr));
  pip.push_back(
      new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
 
  MoFEMFunctionReturn(0);
}

opFactoryDomainLhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >

MoFEMErrorCode HenckyOps::opFactoryDomainLhs	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		std::string	block_name,
		Sev	sev,
		double	scale = 1
	)

Definition at line 829 of file HenckyOps.hpp.

                                                                           {
  MoFEMFunctionBegin;
 
  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
      m_field, pip, field_name, block_name, sev, scale);
  CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
                                                    common_ptr, sev);
 
  MoFEMFunctionReturn(0);
}

opFactoryDomainRhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >

MoFEMErrorCode HenckyOps::opFactoryDomainRhs	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		boost::shared_ptr< HenckyOps::CommonData >	common_ptr,
		Sev	sev
	)

Definition at line 775 of file HenckyOps.hpp.

             {
  MoFEMFunctionBegin;
 
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template LinearForm<I>;
  using OpInternalForcePiola =
      typename B::template OpGradTimesTensor<1, DIM, DIM>;
  pip.push_back(
      new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
 
  MoFEMFunctionReturn(0);
}

opFactoryDomainRhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >

MoFEMErrorCode HenckyOps::opFactoryDomainRhs	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	field_name,
		std::string	block_name,
		Sev	sev,
		double	scale = 1
	)

Definition at line 793 of file HenckyOps.hpp.

                                                                           {
  MoFEMFunctionBegin;
 
  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
      m_field, pip, field_name, block_name, sev, scale);
  CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
                                                    common_ptr, sev);
 
  MoFEMFunctionReturn(0);
}

sort_eigen_vals()

template<int DIM>

auto HenckyOps::sort_eigen_vals ( FTensor::Tensor1< double, DIM > &  eig, FTensor::Tensor2< double, DIM, DIM > &  eigen_vec  )

inline

Definition at line 41 of file HenckyOps.hpp.

                                                                         {
 
  isEq::absMax =
      std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
 
  int i = 0, j = 1, k = 2;
 
  if (is_eq(eig(0), eig(1))) {
    i = 0;
    j = 2;
    k = 1;
  } else if (is_eq(eig(0), eig(2))) {
    i = 0;
    j = 1;
    k = 2;
  } else if (is_eq(eig(1), eig(2))) {
    i = 1;
    j = 0;
    k = 2;
  }
 
  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
      eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
 
      eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
 
      eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
 
  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
 
  {
    FTensor::Index<'i', DIM> i;
    FTensor::Index<'j', DIM> j;
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  }
};

Variable Documentation

d_f

auto HenckyOps::d_f = [](double v) { return 0.5 / v; }

Definition at line 16 of file HenckyOps.hpp.

dd_f

auto HenckyOps::dd_f = [](double v) { return -0.5 / (v * v); }

Definition at line 17 of file HenckyOps.hpp.

eps

constexpr double HenckyOps::eps = std::numeric_limits<float>::epsilon()

constexpr

Definition at line 13 of file HenckyOps.hpp.

f

auto HenckyOps::f = [](double v) { return 0.5 * std::log(v); }

Definition at line 15 of file HenckyOps.hpp.