10 using namespace MoFEM;
12 static char help[] =
"...\n\n";
77 CHKERR setIntegrationRules();
79 CHKERR boundaryCondition();
92 PetscBool load_file = PETSC_FALSE;
96 if (load_file == PETSC_FALSE) {
98 auto &moab = mField.get_moab();
103 double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
105 for (
int nn = 0; nn < 4; nn++) {
106 CHKERR moab.create_vertex(&tet_coords[3 * nn], nodes[nn]);
109 CHKERR moab.create_element(MBTET, nodes, 4, tet);
111 for (
auto d : {1, 2})
112 CHKERR moab.get_adjacencies(&tet, 1,
d,
true, adj);
118 double tri_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
120 for (
int nn = 0; nn < 3; nn++) {
121 CHKERR moab.create_vertex(&tri_coords[3 * nn], nodes[nn]);
124 CHKERR moab.create_element(MBTRI, nodes, 3, tri);
126 CHKERR moab.get_adjacencies(&tri, 1, 1,
true, adj);
129 CHKERR mField.rebuild_database();
130 CHKERR mField.getInterface(simpleInterface);
135 0,
SPACE_DIM, simpleInterface->getBitRefLevel());
139 CHKERR mField.getInterface(simpleInterface);
140 CHKERR simpleInterface->getOptions();
141 CHKERR simpleInterface->loadFile();
154 enum bases { AINSWORTH, AINSWORTH_LOBATTO, DEMKOWICZ, BERNSTEIN, LASBASETOP };
155 const char *list_bases[] = {
"ainsworth",
"ainsworth_lobatto",
"demkowicz",
159 PetscInt choice_base_value = AINSWORTH;
161 &choice_base_value, &flg);
162 if (flg != PETSC_TRUE)
165 if (choice_base_value == AINSWORTH)
167 if (choice_base_value == AINSWORTH_LOBATTO)
169 else if (choice_base_value == DEMKOWICZ)
171 else if (choice_base_value == BERNSTEIN)
174 enum spaces { H1SPACE, L2SPACE, HCURLSPACE, HDIVSPACE, LASBASETSPACE };
175 const char *list_spaces[] = {
"h1",
"l2",
"hcurl",
"hdiv"};
176 PetscInt choice_space_value = H1SPACE;
178 LASBASETSPACE, &choice_space_value, &flg);
179 if (flg != PETSC_TRUE)
182 if (choice_space_value == H1SPACE)
184 else if (choice_space_value == L2SPACE)
186 else if (choice_space_value == HCURLSPACE)
188 else if (choice_space_value == HDIVSPACE)
191 CHKERR simpleInterface->addDomainField(
"U", space, base, 1);
196 CHKERR simpleInterface->setUp();
231 auto post_proc_fe = boost::make_shared<MyPostProc>(mField);
232 post_proc_fe->generateReferenceElementMesh();
236 if (space ==
HCURL) {
237 auto jac_ptr = boost::make_shared<MatrixDouble>();
238 post_proc_fe->getOpPtrVector().push_back(
241 post_proc_fe->getOpPtrVector().push_back(
254 auto u_ptr = boost::make_shared<VectorDouble>();
255 post_proc_fe->getOpPtrVector().push_back(
257 post_proc_fe->getOpPtrVector().push_back(
261 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
281 auto u_ptr = boost::make_shared<MatrixDouble>();
282 post_proc_fe->getOpPtrVector().push_back(
285 post_proc_fe->getOpPtrVector().push_back(
289 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
307 auto scale_tag_val = [&]() {
309 auto &post_proc_mesh = post_proc_fe->getPostProcMesh();
311 CHKERR post_proc_mesh.get_entities_by_type(0, MBVERTEX, nodes);
313 CHKERR post_proc_mesh.tag_get_handle(
"U",
th);
315 CHKERR post_proc_mesh.tag_get_length(
th, length);
316 std::vector<double> data(nodes.size() * length);
317 CHKERR post_proc_mesh.tag_get_data(
th, nodes, &*data.begin());
319 for (
int i = 0;
i != nodes.size(); ++
i) {
321 for (
int d = 0;
d != length; ++
d)
322 v += pow(data[length *
i +
d], 2);
324 max_v = std::max(max_v,
v);
328 CHKERR post_proc_mesh.tag_set_data(
th, nodes, &*data.begin());
333 auto dofs_ptr = mField.get_dofs();
335 for (
auto dof_ptr : (*dofs_ptr)) {
336 MOFEM_LOG(
"PLOTBASE", Sev::verbose) << *dof_ptr;
337 auto &val =
const_cast<double &
>(dof_ptr->getFieldData());
341 CHKERR post_proc_fe->writeFile(
342 "out_base_dof_" + boost::lexical_cast<std::string>(nb) +
".h5m");
343 CHKERR post_proc_fe->getPostProcMesh().delete_mesh();
356 int main(
int argc,
char *argv[]) {
364 DMType dm_name =
"DMMOFEM";
369 auto core_log = logging::core::get();
371 LogManager::createSink(LogManager::getStrmWorld(),
"PLOTBASE"));
372 LogManager::setLog(
"PLOTBASE");
402 char ref_mesh_file_name[255];
405 strcpy(ref_mesh_file_name,
"ref_mesh2d.h5m");
407 strcpy(ref_mesh_file_name,
"ref_mesh3d.h5m");
410 "Dimension not implemented");
414 CHKERR moab_ref.load_file(ref_mesh_file_name, 0,
"");
422 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
423 CHKERR moab_ref.add_entities(meshset, elems);
424 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
425 CHKERR moab_ref.delete_entities(&meshset, 1);
429 CHKERR moab_ref.get_connectivity(elems, elem_nodes,
false);
432 std::map<EntityHandle, int> nodes_pts_map;
435 gaussPts.resize(
SPACE_DIM + 1, elem_nodes.size(),
false);
437 Range::iterator nit = elem_nodes.begin();
438 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
440 CHKERR moab_ref.get_coords(&*nit, 1, coords);
441 for (
auto d : {0, 1, 2})
442 gaussPts(
d, gg) = coords[
d];
443 nodes_pts_map[*nit] = gg;
455 Range::iterator tit = elems.begin();
456 for (
int tt = 0; tit != elems.end(); ++tit, ++tt) {
459 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
460 for (
int nn = 0; nn != num_nodes; ++nn) {
461 refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
471 const int num_nodes = gaussPts.size2();
475 switch (numeredEntFiniteElementPtr->getEntType()) {
479 &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
483 for (
int gg = 0; gg != num_nodes; gg++) {
484 double ksi = gaussPts(0, gg);
485 double eta = gaussPts(1, gg);
495 &gaussPts(0, 0), &gaussPts(1, 0),
496 &gaussPts(2, 0), num_nodes);
500 for (
int gg = 0; gg != num_nodes; gg++) {
501 double ksi = gaussPts(0, gg);
502 double eta = gaussPts(1, gg);
503 double zeta = gaussPts(2, gg);
516 "Not implemented element type");
522 ReadUtilIface *iface;
523 CHKERR getPostProcMesh().query_interface(iface);
525 std::vector<double *> arrays;
530 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
532 mapGaussPts.resize(gaussPts.size2());
533 for (
int gg = 0; gg != num_nodes; ++gg)
534 mapGaussPts[gg] = startv + gg;
537 int def_in_the_loop = -1;
538 CHKERR getPostProcMesh().tag_get_handle(
"NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
539 th, MB_TAG_CREAT | MB_TAG_SPARSE,
545 const int num_nodes_on_ele =
refEleMap.size2();
552 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
555 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
559 "Dimension not implemented");
563 for (
unsigned int tt = 0; tt !=
refEleMap.size1(); ++tt) {
564 for (
int nn = 0; nn != num_nodes_on_ele; ++nn)
565 conn[num_nodes_on_ele * tt + nn] = mapGaussPts[
refEleMap(tt, nn)];
570 CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
572 auto physical_elements =
Range(starte, starte + num_el - 1);
573 CHKERR getPostProcMesh().tag_clear_data(
th, physical_elements, &(nInTheLoop));
575 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
579 mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes,
true);
580 coords.resize(3 * fe_num_nodes,
false);
581 CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
590 arrays[0], arrays[1], arrays[2]);
591 const double *t_coords_ele_x = &coords[0];
592 const double *t_coords_ele_y = &coords[1];
593 const double *t_coords_ele_z = &coords[2];
594 for (
int gg = 0; gg != num_nodes; ++gg) {
596 t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
598 for (
int nn = 0; nn != fe_num_nodes; ++nn) {
599 t_coords(
i) += t_n * t_ele_coords(
i);
600 for (
auto ii : {0, 1, 2})
601 if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
614 ParallelComm *pcomm_post_proc_mesh =
616 if (pcomm_post_proc_mesh != NULL)
617 delete pcomm_post_proc_mesh;
624 auto resolve_shared_ents = [&]() {
627 ParallelComm *pcomm_post_proc_mesh =
628 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
629 if (pcomm_post_proc_mesh == NULL) {
632 pcomm_post_proc_mesh =
new ParallelComm(
633 &(getPostProcMesh()),
637 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
642 CHKERR resolve_shared_ents();