vtk_cgg_base.cpp

Go to the documentation of this file.

 
#include <MoFEM.hpp>
#include <h1_hdiv_hcurl_l2.h>
using namespace MoFEM;
#include <CGGTonsorialBubbleBase.hpp>
 
using namespace FTensor;
 
MoFEMErrorCode VTK_cgg_bubble_base_MBTET(const string file_name) {
  MoFEMFunctionBegin;
 
  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
 
  moab::Core core_ref;
  moab::Interface &moab_ref = core_ref;
 
  EntityHandle nodes[4];
  for (int nn = 0; nn < 4; nn++) {
    CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
  }
  EntityHandle tet;
  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
 
  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
  MoFEM::Interface &m_field_ref = m_core_ref;
 
  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
      0, 3, BitRefLevel().set(0));
 
  const int max_level = 4;
  for (int ll = 0; ll != max_level; ll++) {
    Range edges;
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                       BitRefLevel().set(), MBEDGE, edges);
    Range tets;
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                       BitRefLevel(ll).set(), MBTET, tets);
    // refine mesh
    MeshRefinement *m_ref;
    CHKERR m_field_ref.getInterface(m_ref);
    CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
                                                BitRefLevel().set(ll + 1));
    CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
  }
 
  Range tets;
  CHKERR m_field_ref.getInterface<BitRefManager>()
      ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                     BitRefLevel().set(max_level), MBTET, tets);
 
  // Use 10 node tets to print base
  if (1) {
    EntityHandle meshset;
    CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
    CHKERR moab_ref.add_entities(meshset, tets);
    CHKERR moab_ref.convert_entities(meshset, true, false, false);
    CHKERR moab_ref.delete_entities(&meshset, 1);
  }
 
  Range elem_nodes;
  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
 
  const int nb_gauss_pts = elem_nodes.size();
  MatrixDouble gauss_pts(nb_gauss_pts, 4);
  gauss_pts.clear();
  Range::iterator nit = elem_nodes.begin();
  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
    CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
  }
  gauss_pts = trans(gauss_pts);
 
  MatrixDouble shape_fun;
  shape_fun.resize(nb_gauss_pts, 4);
  CHKERR ShapeMBTET(&*shape_fun.data().begin(), &gauss_pts(0, 0),
                    &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
  double diff_shape_fun[12];
  CHKERR ShapeDiffMBTET(diff_shape_fun);
 
  int p = 3;
  MatrixDouble phi(nb_gauss_pts, 9 * NBVOLUMETET_CCG_BUBBLE(p));
  Tensor2<PackPtr<double *, 9>, 3, 3> t_phi(&phi(0, 0), &phi(0, 1), &phi(0, 2),
 
                                            &phi(0, 3), &phi(0, 4), &phi(0, 5),
 
                                            &phi(0, 6), &phi(0, 7), &phi(0, 8));
 
  CHKERR EshelbianPlasticity::CGG_BubbleBase_MBTET(
      p, &shape_fun(0, 0), diff_shape_fun, t_phi, nb_gauss_pts);
 
  for (int ll = 0; ll != NBVOLUMETET_CCG_BUBBLE(p); ++ll) {
 
    auto get_tag = [&](int d) {
      double def_val[] = {0, 0, 0};
      std::string tag_name = "B_" + boost::lexical_cast<std::string>(ll) +
                              "_D" + boost::lexical_cast<std::string>(d);
      Tag th;
      CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
                                     MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
      return th;
    };
    Tag th0 = get_tag(0);
    Tag th1 = get_tag(1);
    Tag th2 = get_tag(2);
 
    int gg = 0;
    for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
         nit++, gg++) {
 
      auto save_tag = [&](Tag th, const int d, int idx) {
        MoFEMFunctionBegin;
        double data[] = {phi(gg, idx + 3 * d + 0), phi(gg, idx + 3 * d + 1),
                         phi(gg, idx + 3 * d + 2)};
        CHKERR moab_ref.tag_set_data(th, &*nit, 1, data);
        MoFEMFunctionReturn(0);
      };
      const int idx = 9 * ll;
      CHKERR save_tag(th0, 0, idx);
      CHKERR save_tag(th1, 1, idx);
      CHKERR save_tag(th2, 2, idx);
 
    }
  }
 
  EntityHandle meshset;
  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
  CHKERR moab_ref.add_entities(meshset, tets);
  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &meshset, 1);
 
  MoFEMFunctionReturn(0);
}
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
    CHKERR VTK_cgg_bubble_base_MBTET("out_vtk_cgg_bubble_base_on_tet.vtk");
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}