DGProjection.cpp

Go to the documentation of this file.

/**
 * @file DGProjection.cpp
 */
 
namespace MoFEM {
 
static auto dg_save_range(moab::Interface &moab, const std::string name,
                          const Range r, std::vector<Tag> tags = {}) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  if (r.size()) {
    CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1,
                           tags.data(), tags.size());
  } else {
    MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
  }
  MoFEMFunctionReturn(0);
};
 
struct DGSetIntegration {
 
  ForcesAndSourcesCore::RuleHookFun funRule = [](int, int, int p) {
    return 2 * (p + 1);
  };
 
  DGSetIntegration(boost::shared_ptr<Range> edges_ptr,
                              ForcesAndSourcesCore::RuleHookFun fun_rule)
      : edgesPtr(edges_ptr), funRule(fun_rule) {};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 3;
    constexpr int numEdges = 3;
    constexpr int refinementLevels = 1;
 
    auto &m_field = fe_raw_ptr->mField;
    auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
    auto fe_handle = fe_ptr->getFEEntityHandle();
    auto type = type_from_handle(fe_handle);
    if (type != MBTRI) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "only MBTRI is implemented");
    }
 
 
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = funRule(order_row, order_col, order_data);
      if (rule < QUAD_2D_TABLE_SIZE) {
        if (QUAD_2D_TABLE[rule]->dim != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
        }
        if (QUAD_2D_TABLE[rule]->order < rule) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
        fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                    &fe_ptr->gaussPts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                    &fe_ptr->gaussPts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                    &fe_ptr->gaussPts(2, 0), 1);
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
                    1);
        data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
        std::copy(
            Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
            data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (edgesPtr) {
 
      auto get_edges = [&]() {
        std::bitset<numEdges> edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (edgesPtr->find(edge) != edgesPtr->end()) {
            edges.set(ee);
          } else {
            edges.reset(ee);
          }
        }
        return edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        auto *shape_ptr =
            data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().data();
        CHKERR Tools::shapeFunMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
                                    &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto refine_quadrature = [&]() {
        MoFEMFunctionBegin;
 
        const int max_level = refinementLevels;
 
        moab::Core moab_ref;
        double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
        EntityHandle nodes[numNodes];
        for (int nn = 0; nn != numNodes; nn++)
          CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
        EntityHandle tri;
        CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
        MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
        MoFEM::Interface &m_field_ref = mofem_ref_core;
        {
          Range tris(tri, tri);
          Range edges;
          CHKERR m_field_ref.get_moab().get_adjacencies(tris, 1, true, edges,
                                                        moab::Interface::UNION);
          CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
              tris, BitRefLevel().set(0), false, VERBOSE);
        }
 
        auto edges = get_edges();
 
        EntityHandle meshset;
        CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
        for (int ee = 0; ee != numEdges; ee++) {
          if (edges[ee]) {
            EntityHandle ent;
            CHKERR moab_ref.side_element(tri, 1, ee, ent);
            CHKERR moab_ref.add_entities(meshset, &ent, 1);
          }
        }
 
        // refine mesh
        auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
        for (int ll = 0; ll != max_level; ll++) {
          Range ref_edges;
          CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ref_edges,
                                               true);
          CHKERR m_field_ref.getInterface<BitRefManager>()
              ->filterEntitiesByRefLevel(BitRefLevel().set(ll),
                                         BitRefLevel().set(), ref_edges);
 
          Range tris;
          CHKERR m_field_ref.getInterface<BitRefManager>()
              ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                             BitRefLevel().set(), MBTRI, tris);
          CHKERR m_ref->addVerticesInTheMiddleOfEdges(
              ref_edges, BitRefLevel().set(ll + 1));
          CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
          CHKERR m_field_ref.getInterface<BitRefManager>()
              ->updateMeshsetByEntitiesChildren(
                  meshset, BitRefLevel().set(ll + 1), meshset, MBEDGE, true);
        }
 
        // get ref coords
        Range tris;
        CHKERR m_field_ref.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                           BitRefLevel().set(), MBTRI, tris);
 
        if (debug) {
          CHKERR dg_save_range(moab_ref, "ref_tris.vtk", tris);
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
        }
 
        MatrixDouble ref_coords(tris.size(), 9, false);
        int tt = 0;
        for (Range::iterator tit = tris.begin(); tit != tris.end();
             tit++, tt++) {
          int num_nodes;
          const EntityHandle *conn;
          CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
          CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
        }
 
        auto &data = fe_ptr->dataOnElement[H1];
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
        MatrixDouble &shape_n = data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
        int gg = 0;
        for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
          double *tri_coords = &ref_coords(tt, 0);
          FTensor::Tensor1<double, 3> t_normal;
          CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
          auto det = t_normal.l2();
          for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
            for (int dd = 0; dd != 2; dd++) {
              ref_gauss_pts(dd, gg) = shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
                                      shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
                                      shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
            }
            ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
          }
        }
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR refine_quadrature();
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  struct Fe : public ForcesAndSourcesCore {
    using ForcesAndSourcesCore::dataOnElement;
 
  private:
    using ForcesAndSourcesCore::ForcesAndSourcesCore;
  };
 
  boost::shared_ptr<Range> edgesPtr;
 
};
 
template <>
MoFEMErrorCode
setDGSetIntegrationPoints<2>(ForcesAndSourcesCore::GaussHookFun &set_hook,
                             ForcesAndSourcesCore::RuleHookFun rule,
                             boost::shared_ptr<Range> edges_ptr) {
  MoFEMFunctionBegin;
  set_hook = DGSetIntegration(edges_ptr, rule);
  MoFEMFunctionReturn(0);
}
 
OpDGProjectionMassMatrix::OpDGProjectionMassMatrix(
    int order, boost::shared_ptr<MatrixDouble> mass_ptr,
    boost::shared_ptr<EntitiesFieldData> data_l2,
    const FieldApproximationBase b, const FieldSpace s, int verb, Sev sev)
    : OpBaseDerivativesBase(mass_ptr, data_l2, b, s, verb, sev),
      baseOrder(order) {}
 
MoFEMErrorCode
OpDGProjectionMassMatrix::doWork(int side, EntityType type,
                                 EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  if (sPace != L2) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Space should be set to L2");
  }
 
  CHKERR calculateBase(
 
      [this]() { return this->baseOrder; }
 
  );
  CHKERR calculateMass();
 
  MoFEMFunctionReturn(0);
}
 
OpDGProjectionCoefficients::OpDGProjectionCoefficients(
    boost::shared_ptr<MatrixDouble> data_ptr,
    boost::shared_ptr<MatrixDouble> coeffs_ptr,
    boost::shared_ptr<MatrixDouble> mass_ptr,
    boost::shared_ptr<EntitiesFieldData> data_l2,
    const FieldApproximationBase b, const FieldSpace s,
    const LogManager::SeverityLevel sev)
    : OpBaseDerivativesBase(mass_ptr, data_l2, b, s, VERBOSE, sev),
      dataPtr(data_ptr), coeffsPtr(coeffs_ptr) {}
 
MoFEMErrorCode
OpDGProjectionCoefficients::doWork(int side, EntityType type,
                                   EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  auto nb_integration_pts = getGaussPts().size2();
 
  const auto fe_type = getFEType();
  constexpr auto side_number = 0;
  auto &ent_data = dataL2->dataOnEntities[fe_type][side_number];
  auto nb_base_functions = ent_data.getN(base).size2();
 
  auto nb_data_coeffs = dataPtr->size1();
  auto &rhs = *coeffsPtr;
  auto &data = *dataPtr;
 
  rhs.resize(nb_base_functions, nb_data_coeffs, false);
  rhs.clear();
 
  using Tensor0 = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>;
  MatrixDouble trans_data = trans(data);
 
  // assemble rhs
  auto t_base = ent_data.getFTensor0N(base, 0, 0);
  auto t_w = getFTensor0IntegrationWeight();
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    Tensor0 t_rhs(&*rhs.data().begin());
    for (auto bb = 0; bb != nb_base_functions; ++bb) {
      double alpha = t_w * t_base;
      Tensor0 t_data(&trans_data(gg, 0));
      for (auto cc = 0; cc != nb_data_coeffs; ++cc) {
        t_rhs += alpha * t_data;
        ++t_data;
        ++t_rhs;
      }
      ++t_base;
    }
    ++t_w;
  }
 
  // solve for coefficients
  for (auto cc = 0; cc != nb_data_coeffs; ++cc) {
    ublas::matrix_column<MatrixDouble> mc(rhs, cc);
    cholesky_solve(*baseMassPtr, mc, ublas::lower());
  }
 
  MoFEMFunctionReturn(0);
}
 
OpDGProjectionEvaluation::OpDGProjectionEvaluation(
    boost::shared_ptr<MatrixDouble> data_ptr,
    boost::shared_ptr<MatrixDouble> coeffs_ptr,
    boost::shared_ptr<EntitiesFieldData> data_l2,
    const FieldApproximationBase b, const FieldSpace s,
    const LogManager::SeverityLevel sev)
    : OpDGProjectionCoefficients(data_ptr, coeffs_ptr,
                                 boost::make_shared<MatrixDouble>(), data_l2, b,
                                 s, sev) {}
 
MoFEMErrorCode
OpDGProjectionEvaluation::doWork(int side, EntityType type,
                                 EntitiesFieldData::EntData &data) {
 
  MoFEMFunctionBegin;
 
  if (sPace != L2) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Space should be set to L2");
  }
 
  auto fe_type = getFEType();
  constexpr auto side_number = 0;
  auto order = dataL2->dataOnEntities[fe_type][side_number].getOrder();
  CHKERR calculateBase(
 
      [order]() { return order; }
 
  );
 
  auto &ent_data = dataL2->dataOnEntities[fe_type][side_number];
  const auto nb_base_functions = ent_data.getN(base).size2();
  auto nb_data_coeffs = dataPtr->size1();
  if (!nb_data_coeffs) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Number of data coefficients should be set");
  }
  auto &data = *dataPtr;
  auto &coeffs = *coeffsPtr;
 
  auto nb_integration_pts = getGaussPts().size2();
  data.resize(nb_integration_pts, nb_data_coeffs);
  data.clear();
 
  using Tensor0 = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>;
 
  auto t_base = ent_data.getFTensor0N(base, 0, 0);
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    Tensor0 t_coeffs(&*coeffs.data().begin());
    for (auto bb = 0; bb != nb_base_functions; ++bb) {
      double alpha = t_base;
      Tensor0 t_data(&data(gg, 0));
      for (auto cc = 0; cc != nb_data_coeffs; ++cc) {
        t_data += alpha * t_coeffs;
        ++t_data;
        ++t_coeffs;
      }
      ++t_base;
    }
  }
 
  *dataPtr = trans(*dataPtr);  
 
  MoFEMFunctionReturn(0);
}
} // namespace MoFEM