Smoother.hpp

Go to the documentation of this file.

#ifndef __SMOOTHER_HPP__
#define __SMOOTHER_HPP__
 
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
 
struct Smoother {
 
  static MoFEMErrorCode extractMeshEdges(MoFEM::Interface &m_field,
                                            Range &output_edges,
                                            Range &output_vertices) {
    MoFEMFunctionBegin;
 
    Range fineFeatureTris;
    FTENSOR_INDEX(3, i);
    Range vols;
    CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, vols);
    if (vols.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "no 3d entities found in the mesh");
    }
 
    Skinner skinner(&m_field.get_moab());
    Range skin_tris;
    CHKERR skinner.find_skin(0, vols, false, skin_tris);
    if (skin_tris.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "no skin triangles found in the mesh");
    }
 
    auto tri_normal = [&](const EntityHandle tri,
                          FTensor::Tensor1<double, 3> &normal) {
      MoFEMFunctionBeginHot;
      CHKERR m_field.getInterface<Tools>()->getTriNormal(tri, &normal(0));
      normal.normalize();
      MoFEMFunctionReturnHot(0);
    };
 
    const double fine_planar_cos = std::cos(1e-6); //FIXME: parameter?
 
    Range skin_edges;
    CHKERR m_field.get_moab().get_adjacencies(skin_tris, 1, true, skin_edges,
                                              moab::Interface::UNION);
 
    for (auto edge : skin_edges) {
      Range adj_tris;
      CHKERR m_field.get_moab().get_adjacencies(&edge, 1, 2, false, adj_tris,
                                                moab::Interface::UNION);
      adj_tris = intersect(adj_tris, skin_tris);
      if (adj_tris.size() != 2) {
        continue;
      }
 
      FTensor::Tensor1<double, 3> n0;
      FTensor::Tensor1<double, 3> n1;
      CHKERR tri_normal(adj_tris[0], n0);
      CHKERR tri_normal(adj_tris[1], n1);
      const double dot = n0(i) * n1(i);
      if (std::abs(dot) <= fine_planar_cos) {
        fineFeatureTris.insert(adj_tris[0]);
        fineFeatureTris.insert(adj_tris[1]);
        output_edges.insert(edge);
      }
    }
 
    std::map<EntityHandle, std::vector<EntityHandle>> vert_edges;
    for (auto edge : output_edges) {
      const EntityHandle *conn = nullptr;
      int nconn = 0;
      CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
      if (nconn != 2)
        continue;
      vert_edges[conn[0]].push_back(edge);
      vert_edges[conn[1]].push_back(edge);
    }
 
    Range straight_edges;
    for (auto edge : output_edges) {
      const EntityHandle *conn = nullptr;
      int nconn = 0;
      CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
      if (nconn != 2)
        continue;
 
      bool keep_edge = true;
      for (int vv = 0; vv < 2; ++vv) {
        const auto v = conn[vv];
        const auto it = vert_edges.find(v);
        if (it == vert_edges.end()) {
          keep_edge = false;
          break;
        }
        const auto &inc_edges = it->second;
        if (inc_edges.size() > 2) {
          keep_edge = false;
          break;
        }
        if (inc_edges.size() == 2) {
          const auto other_edge =
              (inc_edges[0] == edge) ? inc_edges[1] : inc_edges[0];
          const EntityHandle *conn_other = nullptr;
          int nconn_other = 0;
          CHKERR m_field.get_moab().get_connectivity(other_edge, conn_other,
                                                     nconn_other, false);
          if (nconn_other != 2) {
            keep_edge = false;
            break;
          }
          const EntityHandle v_other =
              (conn_other[0] == v) ? conn_other[1] : conn_other[0];
 
          double pv[3], pa[3], pb[3];
          CHKERR m_field.get_moab().get_coords(&v, 1, pv);
          EntityHandle va = conn[1 - vv];
          CHKERR m_field.get_moab().get_coords(&va, 1, pa);
          CHKERR m_field.get_moab().get_coords(&v_other, 1, pb);
 
          FTensor::Tensor1<double, 3> a;
          FTensor::Tensor1<double, 3> b;
          FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pv(
              &pv[0], &pv[1], &pv[2]);
          FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pa(
              &pa[0], &pa[1], &pa[2]);
          FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pb(
              &pb[0], &pb[1], &pb[2]);
          a(i) = t_pa(i) - t_pv(i);
          b(i) = t_pb(i) - t_pv(i);
          const double na = a.l2();
          const double nb = b.l2();
          if (na == 0 || nb == 0) {
            keep_edge = false;
            break;
          }
          const double dot = (a(i) * b(i)) / (na * nb);
          if (std::abs(dot) < fine_planar_cos) {
            keep_edge = false;
            break;
          }
        }
      }
 
      if (keep_edge) {
        straight_edges.insert(edge);
      }
    }
 
    output_edges = straight_edges;
    output_vertices.clear();
    if (!output_edges.empty()) {
      std::map<EntityHandle, int> vert_count;
      for (auto edge : output_edges) {
        const EntityHandle *conn = nullptr;
        int nconn = 0;
        CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
        if (nconn != 2)
          continue;
        vert_count[conn[0]]++;
        vert_count[conn[1]]++;
      }
      Range filtered_edges;
      for (auto edge : output_edges) {
        const EntityHandle *conn = nullptr;
        int nconn = 0;
        CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
        if (nconn != 2)
          continue;
        const int d0 = vert_count[conn[0]];
        const int d1 = vert_count[conn[1]];
        if (!(d0 == 1 && d1 == 1)) {
          filtered_edges.insert(edge);
          for (int vv = 0; vv < 2; ++vv) {
            const auto v = conn[vv];
            auto it = vert_count.find(v);
            if (it != vert_count.end() && it->second == 1) {
              output_vertices.insert(v);
            }
          }
        }
      }
      output_edges = filtered_edges;
    }
 
    if (!fineFeatureTris.empty())
      CHKERR CutMeshInterface::SaveData(m_field.get_moab())(
          "tri_features_test.vtk", fineFeatureTris);
    if (!output_edges.empty())
      CHKERR CutMeshInterface::SaveData(m_field.get_moab())(
          "edges_test.vtk", output_edges);
    if (!output_vertices.empty())
      CHKERR CutMeshInterface::SaveData(m_field.get_moab())(
          "vertices_test.vtk", output_vertices);
 
    MoFEMFunctionReturn(0);
  }
 
  using DomainEleOp = MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator;
  using LinearForm = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
      GAUSS>;
  using BiLinearForm =
      FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
  using OpInternalForcePiola = LinearForm::OpGradTimesTensor<1, 3, 3>;
  using OpKPiola = BiLinearForm::OpGradTensorGrad<1, 3, 3, 1>;
 
  struct CommonData {
    boost::shared_ptr<MatrixDouble> matGradPtr;
    boost::shared_ptr<MatrixDouble> matFirstPiolaStressPtr;
    boost::shared_ptr<MatrixDouble> matTangentPtr;
    std::string spatialPositions;
    std::string meshPositions;
    std::vector<MatrixDouble> sTress;
 
    CommonData()
        : matGradPtr(new MatrixDouble()),
          matFirstPiolaStressPtr(new MatrixDouble()),
          matTangentPtr(new MatrixDouble()) {}
  };
 
  template <typename TYPE>
  struct CalculatePK1;
 
  struct BlockData {
    Range tEts;
    boost::shared_ptr<CalculatePK1<adouble>>
        materialAdoublePtr;
    boost::shared_ptr<CalculatePK1<double>>
        materialDoublePtr;
  };
 
  std::map<int, BlockData> setOfBlocks;
 
  CommonData commonData;
 
  template <typename TYPE>
  struct CalculatePK1 : public DomainEleOp {
 
    CalculatePK1(const std::string field_name)
        : DomainEleOp(field_name, DomainEleOp::OPROW),
          t_F(resizeAndSet(F)), t_invF(resizeAndSet(invF)),
          t_P(resizeAndSet(P)) {}
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    MatrixBoundedArray<TYPE, 9> F, invF, P;
    FTensor::Tensor2<FTensor::PackPtr<TYPE *, 0>, 3, 3> t_F, t_invF, t_P;
 
    int gG = 0;                ///< Gauss point number
    CommonData *commonDataPtr = nullptr; ///< common data shared between entities (f.e.
                               ///< field values at Gauss pts.)
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
        *opPtr = nullptr; ///< pointer to finite element tetrahedral operator
 
    virtual MoFEMErrorCode calculateP_PiolaKirchhoffI() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    }
 
  protected:
    inline static auto resizeAndSet(MatrixBoundedArray<TYPE, 9> &m) {
      m.resize(3, 3, false);
      return getFTensor2FromArray3by3(m, FTensor::Number<0>(), 0);
    }
  };
 
  // struct MyVolumeFE : public MoFEM::VolumeElementForcesAndSourcesCore {
 
  //   std::string meshPositionsFieldName = "NONE";
  //   int addToRule = 0;
 
  //   MyVolumeFE(MoFEM::Interface &m_field)
  //       : MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
  // };
 
  using MyVolumeFE = MoFEM::VolumeElementForcesAndSourcesCore;
  boost::shared_ptr<MyVolumeFE> feRhsPtr;
  boost::shared_ptr<MyVolumeFE> feLhsPtr;
 
  MyVolumeFE &feRhs; ///< calculate right hand side for tetrahedral elements
  MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
  MyVolumeFE &feLhs; ///< calculate left hand side for tetrahedral elements
  MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
 
  Smoother(MoFEM::Interface &m_field)
      : feRhsPtr(new MyVolumeFE(m_field)),
        feLhsPtr(new MyVolumeFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr) {}
 
 
  struct OpJacobianSmoother
      : public DomainEleOp {
 
    BlockData &dAta;
    CommonData &commonData;
    const int tag;
    const bool computeJacobian;
    VectorDouble activeVars;
    VectorDouble results;
    MatrixDouble jacMat;
 
    OpJacobianSmoother(const std::string field_name, BlockData &data,
                       CommonData &common_data, int tag, bool jacobian)
        : DomainEleOp(field_name, UserDataOperator::OPROW), dAta(data),
          commonData(common_data), tag(tag), computeJacobian(jacobian) {
      std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
 
      const auto nb_gauss_pts = getGaussPts().size2();
      auto &mat_grad = *commonData.matGradPtr;
      if (mat_grad.size2() != nb_gauss_pts)
        MoFEMFunctionReturnHot(0);
 
      auto t_grad = getFTensor2FromMat<3, 3>(mat_grad);
 
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
      FTensor::Index<'l', 3> l;
 
      if (!computeJacobian) {
        auto &mat_p = *commonData.matFirstPiolaStressPtr;
        mat_p.resize(9, nb_gauss_pts, false);
        mat_p.clear();
        auto t_p = getFTensor2FromMat<3, 3>(mat_p);
        if (!dAta.materialDoublePtr) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Material pointer not set");
        }
        for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
          dAta.materialDoublePtr->opPtr = this;
          dAta.materialDoublePtr->gG = gg;
          for (int ii = 0; ii < 3; ++ii) {
            for (int jj = 0; jj < 3; ++jj) {
              dAta.materialDoublePtr->F(ii, jj) = t_grad(ii, jj);
            }
          }
 
          CHKERR dAta.materialDoublePtr->calculateP_PiolaKirchhoffI();
          for (int ii = 0; ii < 3; ++ii)
            for (int jj = 0; jj < 3; ++jj)
              t_p(ii, jj) = dAta.materialDoublePtr->P(ii, jj);
 
          ++t_grad;
          ++t_p;
        }
        MoFEMFunctionReturnHot(0);
      }
 
      if (!dAta.materialAdoublePtr) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "AD material pointer not set");
      }
 
      MatrixDouble &mat_tangent = *commonData.matTangentPtr;
      mat_tangent.resize(81, nb_gauss_pts, false);
      mat_tangent.clear();
      auto t_dP_dF = getFTensor4FromMat<3, 3, 3, 3, 1>(mat_tangent);
 
      activeVars.resize(9, false);
      results.resize(9, false);
      jacMat.resize(9, 9, false);
 
      bool tag_recorded = false;
      for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
        int idx = 0;
        for (int ii = 0; ii < 3; ++ii) {
          for (int jj = 0; jj < 3; ++jj) {
            activeVars[idx++] = t_grad(ii, jj);
          }
        }
 
        if (!tag_recorded) {
          trace_on(tag);
          idx = 0;
          dAta.materialAdoublePtr->opPtr = this;
          dAta.materialAdoublePtr->gG = gg;
          for (int ii = 0; ii < 3; ++ii) {
            for (int jj = 0; jj < 3; ++jj) {
              dAta.materialAdoublePtr->F(ii, jj) <<= activeVars[idx++];
            }
          }
          CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI();
          idx = 0;
          for (int ii = 0; ii < 3; ++ii) {
            for (int jj = 0; jj < 3; ++jj) {
              dAta.materialAdoublePtr->P(ii, jj) >>= results[idx++];
            }
          }
          trace_off();
          tag_recorded = true;
        }
 
        double *jac_ptr[9];
        for (int rr = 0; rr < 9; ++rr)
          jac_ptr[rr] = &jacMat(rr, 0);
 
        int r = ::jacobian(tag, 9, 9, &activeVars[0], jac_ptr);
        if (r < 0)
          SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                  "ADOL-C function evaluation with error");
 
        for (int ii = 0; ii < 3; ++ii) {
          for (int jj = 0; jj < 3; ++jj) {
            for (int kk = 0; kk < 3; ++kk) {
              for (int ll = 0; ll < 3; ++ll) {
                t_dP_dF(ii, jj, kk, ll) =
                    jacMat(ii * 3 + jj, kk * 3 + ll);
              }
            }
          }
        }
 
        ++t_grad;
        ++t_dP_dF;
      }
 
      MoFEMFunctionReturn(0);
    }
  }; 
 
};
 
#endif //__SMOOTHER_HPP__