MatAxisymmetric.hpp

Go to the documentation of this file.

#ifndef MAT_AXISYMMETRIC_HPP
#define MAT_AXISYMMETRIC_HPP
 
namespace MatOps {
 
struct OpCalculateAxisymmetricGradient
    : public ForcesAndSourcesCore::UserDataOperator {
  OpCalculateAxisymmetricGradient(
      boost::shared_ptr<MatrixDouble> disp_ptr,
      boost::shared_ptr<MatrixDouble> grad_ptr,
      boost::shared_ptr<MatrixDouble> axisymmetric_grad_ptr)
      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        dispPtr(disp_ptr), gradPtr(grad_ptr),
        axisymmetricGradPtr(axisymmetric_grad_ptr) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) {
    MoFEMFunctionBegin;
    const int nb_gauss_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_disp = MatrixSizeHelper<GetFTensor1FromMatType<2, -1, DL>, DL>::get(
        *dispPtr, nb_gauss_pts)();
    auto t_grad =
        MatrixSizeHelper<GetFTensor2FromMatType<2, 2, -1, DL>, DL>::get(
            *gradPtr, nb_gauss_pts)();
    auto t_axisymmetric_grad =
        MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
            *axisymmetricGradPtr, nb_gauss_pts)();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    FTENSOR_INDEX(2, j);
    FTENSOR_INDEX(2, I);
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, J);
    FTensor::Number<0> N0;
    FTensor::Number<2> N2;
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
      const double radius = t_coords(0);
      if (radius <= 0)
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Axisymmetric integration point is on or below the axis");
      t_axisymmetric_grad(i, J) = 0;
      t_axisymmetric_grad(j, I) = t_grad(j, I);
      t_axisymmetric_grad(N2, N2) = t_disp(N0) / radius;
      ++t_disp;
      ++t_grad;
      ++t_axisymmetric_grad;
      ++t_coords;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dispPtr;
  boost::shared_ptr<MatrixDouble> gradPtr;
  boost::shared_ptr<MatrixDouble> axisymmetricGradPtr;
};
 
struct AxisymmetricAssemblyData {
  boost::shared_ptr<MatrixDouble> inPlaneStress =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> inPlaneTangent =
      boost::make_shared<MatrixDouble>();
};
 
struct OpCalculateAxisymmetricStress
    : public ForcesAndSourcesCore::UserDataOperator {
  OpCalculateAxisymmetricStress(
      boost::shared_ptr<MatrixDouble> stress_ptr,
      boost::shared_ptr<AxisymmetricAssemblyData> data_ptr)
      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        stressPtr(stress_ptr), dataPtr(data_ptr) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) {
    MoFEMFunctionBegin;
    const int nb_gauss_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_stress =
        MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::get(
            *stressPtr, nb_gauss_pts)();
    auto t_in_plane_stress =
        MatrixSizeHelper<GetFTensor2FromMatType<2, 2, -1, DL>, DL>::size(
            *dataPtr->inPlaneStress, nb_gauss_pts)();
    FTENSOR_INDEX(2, i);
    FTENSOR_INDEX(2, J);
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
      t_in_plane_stress(i, J) = t_stress(i, J);
      ++t_stress;
      ++t_in_plane_stress;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> stressPtr;
  boost::shared_ptr<AxisymmetricAssemblyData> dataPtr;
};
 
struct OpCalculateAxisymmetricTangent
    : public ForcesAndSourcesCore::UserDataOperator {
  OpCalculateAxisymmetricTangent(
      boost::shared_ptr<MatrixDouble> tangent_ptr,
      boost::shared_ptr<AxisymmetricAssemblyData> data_ptr)
      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        tangentPtr(tangent_ptr), dataPtr(data_ptr) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) {
    MoFEMFunctionBegin;
    const int nb_gauss_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_tangent =
        MatrixSizeHelper<GetFTensor4FromMatType<3, 3, 3, 3, -1, DL>, DL>::get(
            *tangentPtr, nb_gauss_pts)();
    auto t_in_plane_tangent =
        MatrixSizeHelper<GetFTensor4FromMatType<2, 2, 2, 2, -1, DL>, DL>::size(
            *dataPtr->inPlaneTangent, nb_gauss_pts)();
    FTENSOR_INDEX(2, i);
    FTENSOR_INDEX(2, J);
    FTENSOR_INDEX(2, k);
    FTENSOR_INDEX(2, L);
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
      t_in_plane_tangent(i, J, k, L) = t_tangent(i, J, k, L);
      ++t_tangent;
      ++t_in_plane_tangent;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> tangentPtr;
  boost::shared_ptr<AxisymmetricAssemblyData> dataPtr;
};
 
template <AssemblyType A, typename DomainEleOp>
using OpAxisymmetricBase =
    typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
 
template <AssemblyType A, typename DomainEleOp>
struct OpAxisymmetricRhs : public OpAxisymmetricBase<A, DomainEleOp> {
  using OpBase = OpAxisymmetricBase<A, DomainEleOp>;
 
  OpAxisymmetricRhs(const std::string field_name,
                    boost::shared_ptr<MatrixDouble> full_stress_ptr)
      : OpBase(field_name, field_name, DomainEleOp::OPROW),
        fullStressPtr(full_stress_ptr) {}
 
protected:
  boost::shared_ptr<MatrixDouble> fullStressPtr;
 
  MoFEMErrorCode iNtegrate(EntData &row_data) override {
    MoFEMFunctionBegin;
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_full_stress =
        MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::get(
            *fullStressPtr, OpBase::nbIntegrationPts)();
    auto t_w = OpBase::getFTensor0IntegrationWeight();
    auto t_row_base = row_data.getFTensor0N();
    FTensor::Number<0> N0;
    FTensor::Number<2> N2;
    for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
      auto t_nf = OpBase::template getNf<2>();
      const double alpha = 2 * M_PI * t_w * OpBase::getMeasure();
      int rr = 0;
      for (; rr != OpBase::nbRows / 2; ++rr) {
        t_nf(N0) += alpha * t_row_base * t_full_stress(N2, N2);
        ++t_nf;
        ++t_row_base;
      }
      for (; rr < OpBase::nbRowBaseFunctions; ++rr)
        ++t_row_base;
      ++t_full_stress;
      ++t_w;
    }
    MoFEMFunctionReturn(0);
  }
};
 
template <AssemblyType A, typename DomainEleOp>
struct OpAxisymmetricLhs : public OpAxisymmetricBase<A, DomainEleOp> {
  using OpBase = OpAxisymmetricBase<A, DomainEleOp>;
 
  OpAxisymmetricLhs(const std::string field_name,
                    boost::shared_ptr<MatrixDouble> full_tangent_ptr)
      : OpBase(field_name, field_name, DomainEleOp::OPROWCOL),
        fullTangentPtr(full_tangent_ptr) {}
 
protected:
  boost::shared_ptr<MatrixDouble> fullTangentPtr;
 
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) override {
    MoFEMFunctionBegin;
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_D =
        MatrixSizeHelper<GetFTensor4FromMatType<3, 3, 3, 3, -1, DL>, DL>::get(
            *fullTangentPtr, OpBase::nbIntegrationPts)();
    auto t_w = OpBase::getFTensor0IntegrationWeight();
    auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
    auto t_row_base = row_data.getFTensor0N();
    auto t_row_diff = row_data.getFTensor1DiffN<2>();
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, J);
    FTENSOR_INDEX(3, k);
    FTENSOR_INDEX(3, L);
    FTENSOR_INDEX(2, j);
    FTENSOR_INDEX(2, I);
    FTENSOR_INDEX(2, l);
    FTensor::Number<0> N0;
    FTensor::Number<2> N2;
    FTensor::Tensor2<double, 3, 3> t_hoop;
    t_hoop(i, J) = 0;
    t_hoop(N2, N2) = 1;
    for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
      const double radius = t_coords(0);
      const double alpha = 2 * M_PI * t_w * OpBase::getMeasure();
      int rr = 0;
      for (; rr != OpBase::nbRows / 2; ++rr) {
        auto t_m = OpBase::template getLocMat<2>(2 * rr);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff = col_data.getFTensor1DiffN<2>(gg, 0);
        for (int cc = 0; cc != OpBase::nbCols / 2; ++cc) {
          t_m(j, N0) += alpha * t_col_base * t_row_diff(I) *
                        (t_D(j, I, k, L) * t_hoop(k, L));
          t_m(N0, l) += alpha * t_row_base * t_hoop(i, J) *
                        (t_D(i, J, l, I) * t_col_diff(I));
          t_m(N0, N0) += alpha * t_row_base * t_col_base / radius *
                          t_hoop(i, J) *
                          (t_D(i, J, k, L) * t_hoop(k, L));
          ++t_m;
          ++t_col_base;
          ++t_col_diff;
        }
        ++t_row_base;
        ++t_row_diff;
      }
      for (; rr < OpBase::nbRowBaseFunctions; ++rr) {
        ++t_row_base;
        ++t_row_diff;
      }
      ++t_D;
      ++t_w;
      ++t_coords;
    }
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace MatOps
 
#endif // MAT_AXISYMMETRIC_HPP