12 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
18 for (
int nn = 0; nn < 4; nn++) {
19 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
22 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
24 MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
30 const int max_level = 4;
31 for (
int ll = 0; ll != max_level; ll++) {
34 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
38 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
43 CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges,
50 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
56 CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
57 CHKERR moab_ref.add_entities(meshset, tets);
58 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
59 CHKERR moab_ref.delete_entities(&meshset, 1);
63 CHKERR moab_ref.get_connectivity(tets, elem_nodes,
false);
65 const int nb_gauss_pts = elem_nodes.size();
68 Range::iterator nit = elem_nodes.begin();
69 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
70 CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
72 gauss_pts = trans(gauss_pts);
75 shape_fun.resize(nb_gauss_pts, 4);
77 &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
78 double diff_shape_fun[12];
90 p, &shape_fun(0, 0), diff_shape_fun, t_phi, nb_gauss_pts);
94 auto get_tag = [&](
int d) {
95 double def_val[] = {0, 0, 0};
96 std::string tag_name =
"B_" + boost::lexical_cast<std::string>(ll) +
97 "_D" + boost::lexical_cast<std::string>(
d);
99 CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE,
th,
100 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
103 Tag th0 = get_tag(0);
104 Tag th1 = get_tag(1);
105 Tag th2 = get_tag(2);
108 for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
111 auto save_tag = [&](Tag
th,
const int d,
int idx) {
113 double data[] = {
phi(gg, idx + 3 *
d + 0),
phi(gg, idx + 3 *
d + 1),
114 phi(gg, idx + 3 *
d + 2)};
115 CHKERR moab_ref.tag_set_data(
th, &*nit, 1, data);
118 const int idx = 9 * ll;
119 CHKERR save_tag(th0, 0, idx);
120 CHKERR save_tag(th1, 1, idx);
121 CHKERR save_tag(th2, 2, idx);
127 CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
128 CHKERR moab_ref.add_entities(meshset, tets);
129 CHKERR moab_ref.write_file(file_name.c_str(),
"VTK",
"", &meshset, 1);