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 const char *list_continuity[] = {
"continuous",
"discontinuous"};
175 PetscInt choice_continuity_value =
CONTINUOUS;
187 enum spaces { H1SPACE, L2SPACE, HCURLSPACE, HDIVSPACE, LASBASETSPACE };
188 const char *list_spaces[] = {
"h1",
"l2",
"hcurl",
"hdiv"};
189 PetscInt choice_space_value = H1SPACE;
191 LASBASETSPACE, &choice_space_value, &flg);
192 if (flg != PETSC_TRUE)
195 if (choice_space_value == H1SPACE)
197 else if (choice_space_value == L2SPACE)
199 else if (choice_space_value == HCURLSPACE)
201 else if (choice_space_value == HDIVSPACE)
205 CHKERR simpleInterface->addDomainField(
"U", space, base, 1);
207 CHKERR simpleInterface->addDomainBrokenField(
"U", space, base, 1);
212 CHKERR simpleInterface->setUp();
214 auto bc_mng = mField.getInterface<
BcManager>();
215 CHKERR bc_mng->removeSideDOFs(simpleInterface->getProblemName(),
"ZERO_FLUX",
251 auto post_proc_fe = boost::make_shared<MyPostProc>(mField);
252 post_proc_fe->generateReferenceElementMesh();
256 if (space ==
HCURL) {
257 auto jac_ptr = boost::make_shared<MatrixDouble>();
258 post_proc_fe->getOpPtrVector().push_back(
261 post_proc_fe->getOpPtrVector().push_back(
274 auto u_ptr = boost::make_shared<VectorDouble>();
275 post_proc_fe->getOpPtrVector().push_back(
277 post_proc_fe->getOpPtrVector().push_back(
281 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
302 post_proc_fe->getOpPtrVector(), {space});
303 auto u_ptr = boost::make_shared<MatrixDouble>();
304 post_proc_fe->getOpPtrVector().push_back(
307 post_proc_fe->getOpPtrVector().push_back(
311 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
329 auto scale_tag_val = [&]() {
331 auto &post_proc_mesh = post_proc_fe->getPostProcMesh();
333 CHKERR post_proc_mesh.get_entities_by_type(0, MBVERTEX, nodes);
335 CHKERR post_proc_mesh.tag_get_handle(
"U",
th);
337 CHKERR post_proc_mesh.tag_get_length(
th, length);
338 std::vector<double> data(nodes.size() * length);
339 CHKERR post_proc_mesh.tag_get_data(
th, nodes, &*data.begin());
341 for (
int i = 0;
i != nodes.size(); ++
i) {
343 for (
int d = 0;
d != length; ++
d)
344 v += pow(data[length *
i +
d], 2);
346 max_v = std::max(max_v,
v);
350 CHKERR post_proc_mesh.tag_set_data(
th, nodes, &*data.begin());
354 auto prb_ptr = mField.get_problem(simpleInterface->getProblemName());
357 auto dofs_ptr = prb_ptr->getNumeredRowDofsPtr();
359 for (
auto dof_ptr : (*dofs_ptr)) {
360 MOFEM_LOG(
"PLOTBASE", Sev::verbose) << *dof_ptr;
361 auto &val =
const_cast<double &
>(dof_ptr->getFieldData());
365 CHKERR post_proc_fe->writeFile(
366 "out_base_dof_" + boost::lexical_cast<std::string>(nb) +
".h5m");
367 CHKERR post_proc_fe->getPostProcMesh().delete_mesh();
380 int main(
int argc,
char *argv[]) {
388 DMType dm_name =
"DMMOFEM";
393 auto core_log = logging::core::get();
395 LogManager::createSink(LogManager::getStrmWorld(),
"PLOTBASE"));
396 LogManager::setLog(
"PLOTBASE");
426 char ref_mesh_file_name[255];
429 strcpy(ref_mesh_file_name,
"ref_mesh2d.h5m");
431 strcpy(ref_mesh_file_name,
"ref_mesh3d.h5m");
434 "Dimension not implemented");
438 CHKERR moab_ref.load_file(ref_mesh_file_name, 0,
"");
446 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
447 CHKERR moab_ref.add_entities(meshset, elems);
448 CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
449 CHKERR moab_ref.delete_entities(&meshset, 1);
453 CHKERR moab_ref.get_connectivity(elems, elem_nodes,
false);
456 std::map<EntityHandle, int> nodes_pts_map;
459 gaussPts.resize(
SPACE_DIM + 1, elem_nodes.size(),
false);
461 Range::iterator nit = elem_nodes.begin();
462 for (
int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
464 CHKERR moab_ref.get_coords(&*nit, 1, coords);
465 for (
auto d : {0, 1, 2})
466 gaussPts(
d, gg) = coords[
d];
467 nodes_pts_map[*nit] = gg;
479 Range::iterator tit = elems.begin();
480 for (
int tt = 0; tit != elems.end(); ++tit, ++tt) {
483 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
484 for (
int nn = 0; nn != num_nodes; ++nn) {
485 refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
495 const int num_nodes = gaussPts.size2();
499 switch (numeredEntFiniteElementPtr->getEntType()) {
503 &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
507 for (
int gg = 0; gg != num_nodes; gg++) {
508 double ksi = gaussPts(0, gg);
509 double eta = gaussPts(1, gg);
519 &gaussPts(0, 0), &gaussPts(1, 0),
520 &gaussPts(2, 0), num_nodes);
524 for (
int gg = 0; gg != num_nodes; gg++) {
525 double ksi = gaussPts(0, gg);
526 double eta = gaussPts(1, gg);
527 double zeta = gaussPts(2, gg);
540 "Not implemented element type");
546 ReadUtilIface *iface;
547 CHKERR getPostProcMesh().query_interface(iface);
549 std::vector<double *> arrays;
554 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
556 mapGaussPts.resize(gaussPts.size2());
557 for (
int gg = 0; gg != num_nodes; ++gg)
558 mapGaussPts[gg] = startv + gg;
561 int def_in_the_loop = -1;
562 CHKERR getPostProcMesh().tag_get_handle(
"NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
563 th, MB_TAG_CREAT | MB_TAG_SPARSE,
569 const int num_nodes_on_ele =
refEleMap.size2();
576 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
579 CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
583 "Dimension not implemented");
587 for (
unsigned int tt = 0; tt !=
refEleMap.size1(); ++tt) {
588 for (
int nn = 0; nn != num_nodes_on_ele; ++nn)
589 conn[num_nodes_on_ele * tt + nn] = mapGaussPts[
refEleMap(tt, nn)];
594 CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
596 auto physical_elements =
Range(starte, starte + num_el - 1);
597 CHKERR getPostProcMesh().tag_clear_data(
th, physical_elements, &(nInTheLoop));
599 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
603 mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes,
true);
604 coords.resize(3 * fe_num_nodes,
false);
605 CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
614 arrays[0], arrays[1], arrays[2]);
615 const double *t_coords_ele_x = &coords[0];
616 const double *t_coords_ele_y = &coords[1];
617 const double *t_coords_ele_z = &coords[2];
618 for (
int gg = 0; gg != num_nodes; ++gg) {
620 t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
622 for (
int nn = 0; nn != fe_num_nodes; ++nn) {
623 t_coords(
i) += t_n * t_ele_coords(
i);
624 for (
auto ii : {0, 1, 2})
625 if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
638 ParallelComm *pcomm_post_proc_mesh =
640 if (pcomm_post_proc_mesh != NULL)
641 delete pcomm_post_proc_mesh;
648 auto resolve_shared_ents = [&]() {
651 ParallelComm *pcomm_post_proc_mesh =
652 ParallelComm::get_pcomm(&(getPostProcMesh()),
MYPCOMM_INDEX);
653 if (pcomm_post_proc_mesh == NULL) {
656 pcomm_post_proc_mesh =
new ParallelComm(
657 &(getPostProcMesh()),
661 CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
666 CHKERR resolve_shared_ents();