EdgeFlippingOps.hpp

 
 
/** \file EdgeFlippingOps.hpp
 * \example EdgeFlippingOps.hpp
 */
 
namespace EdgeFlippingOps {
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpLhsSymTensorLeastSquaresProjImpl;
template <int DIM, IntegrationType I>
using OpLhsSymTensorLeastSquaresProj =
    OpLhsSymTensorLeastSquaresProjImpl<DIM, I, AssemblyDomainEleOp>;
 
/**
 * @brief Filter for scalar fields
 *
 */
template <AssemblyType A, IntegrationType I>
struct filterScalarSolution : public AssemblyDomainEleOp {
 
  filterScalarSolution(
      const std::string field_name, boost::shared_ptr<VectorDouble> old_sol_ptr,
      boost::shared_ptr<VectorDouble> dot_new_sol_ptr,
      boost::shared_ptr<VectorDouble> new_sol_ptr,
      boost::shared_ptr<MatrixDouble> grad_new_sol_ptr,
      ScalarFun scalar_function = [](double, double,
                                     double) constexpr { return 1; })
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        oldSolPtr(old_sol_ptr), dotNewSolPtr(dot_new_sol_ptr),
        newSolPtr(new_sol_ptr), gradNewSolPtr(grad_new_sol_ptr),
        scalarFunction(scalar_function) {}
 
  MoFEMErrorCode doWork(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_old_sol = getFTensor0FromVec(*oldSolPtr);
    auto t_new_sol = getFTensor0FromVec(*newSolPtr);
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double alpha =
          t_w * scalarFunction(t_coords(0), t_coords(1), t_coords(2));
 
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        t_old_sol += alpha * t_old_sol;
 
        if (t_w > 0.0)
          continue;
        if (t_w > 0.0 && t_old_sol < fabs(1e-12))
          continue;
        if (t_w < fabs(1e-12))
          continue;
        ++t_base;
        ++t_diff_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_old_sol;
      ++t_new_sol;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFun scalarFunction;
  boost::shared_ptr<VectorDouble> dotNewSolPtr;
  boost::shared_ptr<VectorDouble> newSolPtr;
  boost::shared_ptr<MatrixDouble> gradNewSolPtr;
  boost::shared_ptr<VectorDouble> oldSolPtr;
};
 
/**
 * @brief Rhs for mapping scalar fields with least squares
 *
 */
template <AssemblyType A, IntegrationType I, int DIM>
struct OpRhsScalarLeastSquaresProj : public AssemblyDomainEleOp {
 
  OpRhsScalarLeastSquaresProj(
      const std::string field_name, boost::shared_ptr<VectorDouble> old_sol_ptr,
      boost::shared_ptr<VectorDouble> dot_new_sol_ptr,
      boost::shared_ptr<VectorDouble> new_sol_ptr,
      boost::shared_ptr<MatrixDouble> grad_new_sol_ptr,
      MoFEM::Interface &m_field,
      ScalarFun scalar_function = [](double, double,
                                     double) constexpr { return 1; })
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        oldSolPtr(old_sol_ptr), dotNewSolPtr(dot_new_sol_ptr),
        newSolPtr(new_sol_ptr), gradNewSolPtr(grad_new_sol_ptr),
        mField(m_field), scalarFunction(scalar_function) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_old_sol = getFTensor0FromVec(*oldSolPtr);
    auto t_new_sol = getFTensor0FromVec(*newSolPtr);
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double alpha =
          t_w * vol * scalarFunction(t_coords(0), t_coords(1), t_coords(2));
 
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        locF[bb] += (t_base * alpha) * (t_old_sol);
        ++t_base;
        ++t_diff_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_old_sol;
      ++t_new_sol;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFun scalarFunction;
  boost::shared_ptr<VectorDouble> dotNewSolPtr;
  boost::shared_ptr<VectorDouble> newSolPtr;
  boost::shared_ptr<MatrixDouble> gradNewSolPtr;
  boost::shared_ptr<VectorDouble> oldSolPtr;
  MoFEM::Interface &mField;
};
 
/**
 * @brief Rhs for mapping vector fields with a scalar basis with least squares
 *
 */
template <AssemblyType A, IntegrationType I, int DIM>
struct OpRhsH1VectorLeastSquaresProj : public AssemblyDomainEleOp {
 
  OpRhsH1VectorLeastSquaresProj(
      const std::string field_name, boost::shared_ptr<MatrixDouble> old_sol_ptr,
      boost::shared_ptr<MatrixDouble> dot_new_sol_ptr,
      boost::shared_ptr<MatrixDouble> new_sol_ptr,
      boost::shared_ptr<MatrixDouble> grad_new_sol_ptr,
      MoFEM::Interface &m_field,
      ScalarFun scalar_function = [](double, double,
                                     double) constexpr { return 1; })
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        oldSolPtr(old_sol_ptr), dotNewSolPtr(dot_new_sol_ptr),
        newSolPtr(new_sol_ptr), gradNewSolPtr(grad_new_sol_ptr),
        mField(m_field), scalarFunction(scalar_function) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(DIM, i);
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_old_sol = getFTensor1FromMat<SPACE_DIM>(*oldSolPtr);
    auto t_new_sol = getFTensor1FromMat<SPACE_DIM>(*newSolPtr);
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double alpha =
          t_w * vol * scalarFunction(t_coords(0), t_coords(1), t_coords(2));
 
      auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(locF);
 
      int bb = 0;
      for (; bb != nbRows / SPACE_DIM; ++bb) {
        t_nf(i) += (t_base * alpha) * (t_old_sol(i));
        ++t_base;
        ++t_diff_base;
        ++t_nf;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_old_sol;
      ++t_new_sol;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFun scalarFunction;
  boost::shared_ptr<MatrixDouble> dotNewSolPtr;
  boost::shared_ptr<MatrixDouble> newSolPtr;
  boost::shared_ptr<MatrixDouble> gradNewSolPtr;
  boost::shared_ptr<MatrixDouble> oldSolPtr;
  MoFEM::Interface &mField;
};
 
/**
 * @brief Rhs for mapping symmetric tensor fields with least squares
 *
 */
template <AssemblyType A, IntegrationType I, int DIM>
struct OpRhsSymTensorLeastSquaresProj : public AssemblyDomainEleOp {
 
  OpRhsSymTensorLeastSquaresProj(
      const std::string field_name, boost::shared_ptr<MatrixDouble> old_sol_ptr,
      boost::shared_ptr<MatrixDouble> dot_new_sol_ptr,
      boost::shared_ptr<MatrixDouble> new_sol_ptr,
      boost::shared_ptr<MatrixDouble> grad_new_sol_ptr,
      MoFEM::Interface &m_field,
      ScalarFun scalar_function = [](double, double,
                                     double) constexpr { return 1; })
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        oldSolPtr(old_sol_ptr), dotNewSolPtr(dot_new_sol_ptr),
        newSolPtr(new_sol_ptr), gradNewSolPtr(grad_new_sol_ptr),
        mField(m_field), scalarFunction(scalar_function) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, j);
    constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
    FTENSOR_INDEX(size_symm, L);
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_L = PlasticOps::symm_L_tensor(FTensor::Number<DIM>());
    auto &nf = AssemblyDomainEleOp::locF;
 
    auto t_old_sol = getFTensor2SymmetricFromMat<DIM>(*oldSolPtr);
    auto t_new_sol = getFTensor2SymmetricFromMat<DIM>(*newSolPtr);
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double alpha =
          t_w * vol * scalarFunction(t_coords(0), t_coords(1), t_coords(2));
 
      FTensor::Tensor1<double, size_symm> t_old_sol_L;
      t_old_sol_L(L) = alpha * t_old_sol(i, j) * t_L(i, j, L);
 
      auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
      int bb = 0;
      for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
        t_nf(L) += t_base * (t_old_sol_L(L));
        ++t_base;
        ++t_diff_base;
        ++t_nf;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_old_sol;
      ++t_new_sol;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFun scalarFunction;
  boost::shared_ptr<MatrixDouble> dotNewSolPtr;
  boost::shared_ptr<MatrixDouble> newSolPtr;
  boost::shared_ptr<MatrixDouble> gradNewSolPtr;
  boost::shared_ptr<MatrixDouble> oldSolPtr;
  MoFEM::Interface &mField;
};
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpLhsSymTensorLeastSquaresProjImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpLhsSymTensorLeastSquaresProjImpl(const std::string row_field_name,
                                     const std::string col_field_name);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
 
private:
};
 
template <int DIM, typename AssemblyDomainEleOp>
OpLhsSymTensorLeastSquaresProjImpl<DIM, GAUSS, AssemblyDomainEleOp>::
    OpLhsSymTensorLeastSquaresProjImpl(const std::string row_field_name,
                                       const std::string col_field_name)
    : AssemblyDomainEleOp(row_field_name, col_field_name,
                          AssemblyDomainEleOp::OPROWCOL) {
  AssemblyDomainEleOp::sYmm = false;
}
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpLhsSymTensorLeastSquaresProjImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(DIM, i);
  FTENSOR_INDEX(DIM, j);
  FTENSOR_INDEX(DIM, k);
  FTENSOR_INDEX(DIM, l);
  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
  FTENSOR_INDEX(size_symm, L);
  FTENSOR_INDEX(size_symm, O);
 
  auto &locMat = AssemblyDomainEleOp::locMat;
 
  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
  const auto nb_row_base_functions = row_data.getN().size2();
 
  auto t_L = PlasticOps::symm_L_tensor(FTensor::Number<DIM>());
 
  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
  auto t_row_base = row_data.getFTensor0N();
  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
    double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
    ++t_w;
 
    FTensor::Tensor2<double, size_symm, size_symm> t_res_mat;
    t_res_mat(O, L) = alpha * (t_L(i, j, O) * t_L(i, j, L));
 
    size_t rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
      auto t_mat = PlasticOps::get_mat_tensor_sym_dtensor_sym(
          rr, locMat, FTensor::Number<SPACE_DIM>());
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
        t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
        ++t_mat;
        ++t_col_base;
      }
 
      ++t_row_base;
    }
 
    for (; rr < nb_row_base_functions; ++rr)
      ++t_row_base;
  }
 
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Rhs for mapping vector fields with an Hdiv basis with least squares
 *
 */
template <AssemblyType A, IntegrationType I, int DIM>
struct OpRhsHdivLeastSquaresProj : public AssemblyDomainEleOp {
 
  OpRhsHdivLeastSquaresProj(
      const std::string field_name, boost::shared_ptr<MatrixDouble> old_sol_ptr,
      boost::shared_ptr<MatrixDouble> dot_new_sol_ptr,
      boost::shared_ptr<MatrixDouble> new_sol_ptr,
      boost::shared_ptr<MatrixDouble> grad_new_sol_ptr,
      ScalarFun resistance_function, MoFEM::Interface &m_field,
      ScalarFun scalar_function = [](double, double,
                                     double) constexpr { return 1; })
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        oldSolPtr(old_sol_ptr), dotNewSolPtr(dot_new_sol_ptr),
        newSolPtr(new_sol_ptr), gradNewSolPtr(grad_new_sol_ptr),
        resistanceFunction(resistance_function), mField(m_field),
        scalarFunction(scalar_function) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(DIM, i);
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor1N<3>();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_old_sol = getFTensor1FromMat<SPACE_DIM>(*oldSolPtr);
    auto t_new_sol = getFTensor1FromMat<SPACE_DIM>(*newSolPtr);
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double alpha =
          t_w * vol * scalarFunction(t_coords(0), t_coords(1), t_coords(2));
 
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        locF[bb] += (t_base(i) * alpha) *
                    resistanceFunction(t_coords(0), t_coords(1), t_coords(2)) *
                    (t_old_sol(i));
        ++t_base;
        ++t_diff_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_old_sol;
      ++t_new_sol;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFun scalarFunction;
  boost::shared_ptr<MatrixDouble> dotNewSolPtr;
  boost::shared_ptr<MatrixDouble> newSolPtr;
  boost::shared_ptr<MatrixDouble> gradNewSolPtr;
  boost::shared_ptr<MatrixDouble> oldSolPtr;
  ScalarFun resistanceFunction;
  MoFEM::Interface &mField;
};
 
/**
 * @brief Rhs for testing EP mapping with initial conditions
 *
 */
template <AssemblyType A, IntegrationType I, int DIM>
struct OpRhsSetInitEP : public AssemblyDomainEleOp {
 
  OpRhsSetInitEP(const std::string field_name,
                 boost::shared_ptr<MatrixDouble> dot_EP_ptr,
                 boost::shared_ptr<MatrixDouble> EP_ptr,
                 boost::shared_ptr<MatrixDouble> grad_EP_ptr,
                 boost::shared_ptr<double> initial_EP_ptr,
                 boost::shared_ptr<double> peak_EP_ptr)
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        dotEPPtr(dot_EP_ptr), EPPtr(EP_ptr), gradEPPtr(grad_EP_ptr),
        initialEPPtr(initial_EP_ptr), peakEPPtr(peak_EP_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(DIM, i);
    FTENSOR_INDEX(DIM, j);
    constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
    FTENSOR_INDEX(size_symm, L);
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    auto t_EP = getFTensor2SymmetricFromMat<SPACE_DIM>(*EPPtr);
    // auto t_grad_temp = getFTensor1FromMat<SPACE_DIM>(*gradTPtr);
 
    const double alpha = t_w * vol;
 
    FTensor::Tensor2_symmetric<double, DIM> t_set_EP;
 
    auto t_L = PlasticOps::symm_L_tensor(FTensor::Number<DIM>());
    auto &nf = AssemblyDomainEleOp::locF;
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      t_set_EP(0, 0) = init_T(*initialEPPtr, *peakEPPtr, t_coords(0),
                              t_coords(1), t_coords(2));
      t_set_EP(1, 1) = t_set_EP(0, 0) - 0.01;
      t_set_EP(0, 1) = t_set_EP(0, 0) - 0.02;
 
      // if (DIM == 3) {
      //   t_set_EP(2, 2) = t_set_EP(0, 0) - 3.0;
      //   t_set_EP(0, 2) = t_set_EP(0, 0) - 4.0;
      //   t_set_EP(1, 2) = t_set_EP(0, 0) - 5.0;
      // }
 
      const double alpha = t_w * vol;
 
      FTensor::Tensor1<double, size_symm> t_set_EP_L;
      t_set_EP_L(L) = alpha * t_set_EP(i, j) * t_L(i, j, L);
 
      auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
      int bb = 0;
      for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
        t_nf(L) -= t_base * (t_set_EP_L(L));
        ++t_base;
        ++t_diff_base;
        ++t_nf;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_EP;
      // ++t_grad_g;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dotEPPtr;
  boost::shared_ptr<MatrixDouble> EPPtr;
  boost::shared_ptr<MatrixDouble> gradEPPtr;
  boost::shared_ptr<MatrixDouble> gradQPtr;
  boost::shared_ptr<double> initialEPPtr;
  boost::shared_ptr<double> peakEPPtr;
};
 
} // namespace EdgeFlippingOps