MatVolumeLengthQuality.cpp

Go to the documentation of this file.

/**
 * @file MatVolumeLengthQuality.cpp
 * @brief ADOL-C implementation of the volume-length quality material
 * @version 0.1
 * @date 2026-04-03
 *
 * @copyright Copyright (c) 2026
 *
 */
 
namespace MatOps {
 
static const char *VolumeLengthQualityTypeNames[] = {
    "QUALITY", "BARRIER_AND_QUALITY", "BARRIER_AND_CHANGE_QUALITY",
    "BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME"};
 
template <int DIM>
struct MatVolumeLengthQuality : public MatElasticImpl<DIM> {
  using MatElasticImpl<DIM>::MatElasticImpl;
 
  using A = MatElasticImpl<DIM>;
 
  static inline constexpr std::array<std::array<int, 2>, 6> edgeNodes = {
      std::array<int, 2>{0, 1}, std::array<int, 2>{0, 2},
      std::array<int, 2>{0, 3}, std::array<int, 2>{1, 2},
      std::array<int, 2>{1, 3}, std::array<int, 2>{2, 3}};
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    MOFEM_LOG_CHANNEL("WORLD");
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "",
                      "Get VolumeLengthQuality material options", "none");
    CHKERR PetscOptionsEList(
        "-volume_length_type", "Volume length quality type", "",
        VolumeLengthQualityTypeNames, LASTOP_VOLUMELENGTHQUALITYTYPE,
        VolumeLengthQualityTypeNames[tYpe], &tYpe, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-volume_length_alpha",
                              "volume length alpha parameter", "", aLpha,
                              &aLpha, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-volume_length_gamma",
                              "volume length parameter (barrier)", "", gAmma,
                              &gAmma, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "VolumeLengthQuality")
        << " type = " << VolumeLengthQualityTypeNames[tYpe]
        << " alpha = " << aLpha << " gamma = " << gAmma;
 
    mFieldPtr = m_field_ptr;
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode setParams(FEMethod *fe_ptr, int gg) override {
    MoFEMFunctionBegin;
    (void)gg;
    const auto ent = fe_ptr->getFEEntityHandle();
 
    if (mFieldPtr == nullptr) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "MoFEM interface pointer not set");
    }
 
    const EntityHandle *conn = nullptr;
    int num_nodes = 0;
    CHKERR mFieldPtr->get_moab().get_connectivity(ent, conn, num_nodes, true);
    if (PetscUnlikely(num_nodes != 4)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "VolumeLengthQuality is implemented for 4-node tetrahedra only");
    }
 
    CHKERR mFieldPtr->get_moab().get_coords(conn, num_nodes, coordsTet.data());
    set_param_vec(A::tAg, coordsTet.size(), coordsTet.data());
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
    static_assert(DIM == 3,
                  "MatVolumeLengthQuality is implemented for 3D only");
 
    A::matOpsDataPtr->insertCommonData("grad", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->insertActiveData("F", MatrixDouble());
    A::matOpsDataPtr->insertDependentData("P", MatrixDouble());
    A::matOpsDataPtr->insertDependentDerivativesData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
 
    auto t_F = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getActiveDataPtr("F")->data().data());
    auto t_P = getFTensor2FromPtr<DIM, DIM>(
        A::matOpsDataPtr->getDependentDataPtr("P")->data().data());
 
    auto t_kd = FTensor::Kronecker_Delta<int>();
 
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, j);
    FTENSOR_INDEX(DIM, I);
    FTENSOR_INDEX(DIM, J);
 
    // That assumes that nodal positions are approximated.
    t_F(i, J) = t_kd(i, J);
 
    FTensor::Tensor2<adouble, DIM, DIM> ta_F;
    FTensor::Tensor2<adouble, DIM, DIM> ta_P;
    FTensor::Tensor2<adouble, DIM, DIM> ta_invF;
    FTensor::Tensor2<adouble, DIM, DIM> dX_dChiT;
    FTensor::Tensor2<adouble, DIM, DIM> Q;
    FTensor::Tensor1<adouble, DIM> delta_chi;
    FTensor::Tensor1<adouble, DIM> delta_X;
 
    adouble det_aF;
 
    trace_on(A::tAg);
 
    std::array<adouble, 12> refCoords;
    for (auto dd = 0; dd != 12; ++dd) {
      refCoords[dd] = mkparam(coordsTet[dd]);
    }
 
    ta_F(i, J) <<= t_F(i, J);
    det_aF = determinantTensor(ta_F);
    CHKERR invertTensor(ta_F, det_aF, ta_invF);
 
    const auto x00 = refCoords[0];
    const auto x01 = refCoords[1];
    const auto x02 = refCoords[2];
    const auto x10 = refCoords[3];
    const auto x11 = refCoords[4];
    const auto x12 = refCoords[5];
    const auto x20 = refCoords[6];
    const auto x21 = refCoords[7];
    const auto x22 = refCoords[8];
    const auto x30 = refCoords[9];
    const auto x31 = refCoords[10];
    const auto x32 = refCoords[11];
 
    const auto ax = x10 - x00;
    const auto ay = x11 - x01;
    const auto az = x12 - x02;
 
    const auto bx = x20 - x00;
    const auto by = x21 - x01;
    const auto bz = x22 - x02;
 
    const auto cx = x30 - x00;
    const auto cy = x31 - x01;
    const auto cz = x32 - x02;
 
    const adouble signed_volume =
        (ax * (by * cz - bz * cy) - ay * (bx * cz - bz * cx) +
         az * (bx * cy - by * cx)) /
        6.0;
 
    dX_dChiT(i, J) = 0.0;
 
    adouble lrms_squared = 0.0;
    adouble lrms_squared0 = 0.0;
 
    for (auto ee = 0; ee != 6; ++ee) {
      for (auto dd = 0; dd != DIM; ++dd) {
        delta_chi(dd) =
            refCoords[3 * edgeNodes[ee][0] + dd] -
            refCoords[3 * edgeNodes[ee][1] + dd];
      }
 
      delta_X(i) = ta_F(i, J) * delta_chi(J);
 
      lrms_squared += (1. / 6.) * delta_X(i) * delta_X(i);
      lrms_squared0 += (1. / 6.) * delta_chi(I) * delta_chi(I);
      dX_dChiT(i, I) += delta_X(i) * delta_chi(I);
    }
 
    Q(i, j) = ta_invF(j, i) - 0.5 * dX_dChiT(i, j) / lrms_squared;
 
    const adouble lrms_cubed0 = lrms_squared0 * sqrt(lrms_squared0);
    const adouble b =
        det_aF / (lrms_squared * sqrt(lrms_squared) / lrms_cubed0);
    const adouble current_volume = det_aF * signed_volume;
    const adouble q =
        6. * sqrt(2.) * current_volume / (lrms_squared * sqrt(lrms_squared));
 
    switch (tYpe) {
    case QUALITY:
      ta_P(i, J) = q * Q(i, J);
      break;
    case BARRIER_AND_QUALITY: {
      const adouble t_mp = q / (1.0 - gAmma) - 1.0 / (q - gAmma);
      ta_P(i, J) = t_mp * Q(i, J);
      break;
    }
    case BARRIER_AND_CHANGE_QUALITY: {
      const adouble t_mp = b / (1.0 - gAmma) - 1.0 / (b - gAmma);
      ta_P(i, J) = t_mp * Q(i, J);
      break;
    }
    case BARRIER_AND_CHANGE_QUALITY_SCALED_BY_VOLUME: {
      adouble t_mp = current_volume;
      t_mp *= b / (1.0 - gAmma) - 1.0 / (b - gAmma);
      ta_P(i, J) = t_mp * Q(i, J);
      break;
    }
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unknown volume-length-quality type");
    }
 
    ta_P(i, J) *= aLpha / signed_volume;
    ta_P(i, J) >>= t_P(i, J);
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  int tYpe = BARRIER_AND_QUALITY;
  double aLpha = 1.0;
  double gAmma = 0.0;
  std::array<double, 12> coordsTet = {
 
      0.0, 0.0, 0.0,
 
      1.0, 0.0, 0.0,
 
      0.0, 1.0, 0.0,
 
      0.0, 0.0, 1.0};
  MoFEM::Interface *mFieldPtr = nullptr;
};
 
template <>
boost::shared_ptr<PhysicalEquations>
createMatOpsPhysicalEquationsPtr<ELASTICITY::VOLUMELENGTHQUALITY, MODEL_3D>(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag) {
  return boost::make_shared<MatVolumeLengthQuality<3>>(mat_ops_data_ptr, tag);
}
 
} // namespace MatOps