520 {
522
523#ifndef NDEBUG
526 << "Refinement for quad is not implemented";
527#endif
528
529 moab::Core core_ref;
530 moab::Interface &moab_ref = core_ref;
531
532 auto create_reference_element = [&moab_ref]() {
534 constexpr double base_coords[] = {
535
536 0, 0,
537 0,
538
539 1, 0,
540 0,
541
542 1, 1,
543 0,
544
545 0, 1,
546 0
547
548 };
550 for (int nn = 0; nn < 4; nn++)
551 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
553 CHKERR moab_ref.create_element(MBQUAD, nodes, 4, quad);
555 };
556
557 auto add_ho_nodes = [&]() {
560 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
562 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
563 CHKERR moab_ref.add_entities(meshset, quads);
564 CHKERR moab_ref.convert_entities(meshset,
true,
true,
true);
565 CHKERR moab_ref.delete_entities(&meshset, 1);
567 };
568
569 auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
572 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
574 CHKERR moab_ref.get_connectivity(quads, quads_nodes,
false);
576 gauss_pts.resize(quads_nodes.size(), 4, false);
577 size_t gg = 0;
578 for (auto node : quads_nodes) {
579 CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
580 little_map[node] = gg;
581 ++gg;
582 }
583 gauss_pts = trans(gauss_pts);
585 };
586
587 auto set_ref_quads = [&](std::map<EntityHandle, int> &little_map) {
590 CHKERR moab_ref.get_entities_by_type(0, MBQUAD, quads,
true);
591 size_t hh = 0;
593 for (auto quad : quads) {
595 int num_nodes;
596 CHKERR moab_ref.get_connectivity(quad, conn, num_nodes,
false);
597 if (ref_quads.size2() != num_nodes) {
598 ref_quads.resize(quads.size(), num_nodes);
599 }
600 for (int nn = 0; nn != num_nodes; ++nn) {
601 ref_quads(hh, nn) = little_map[conn[nn]];
602 }
603 ++hh;
604 }
606 };
607
608 auto set_shape_functions = [&]() {
612 const auto nb_gauss_pts = gauss_pts.size2();
613 shape_functions.resize(nb_gauss_pts, 4);
614 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
615 const double ksi = gauss_pts(0, gg);
616 const double zeta = gauss_pts(1, gg);
621 }
623 };
624
628
629 CHKERR create_reference_element();
632 std::map<EntityHandle, int> little_map;
633 CHKERR set_gauss_pts(little_map);
634 CHKERR set_ref_quads(little_map);
635 CHKERR set_shape_functions();
636
638}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
std::vector< ublas::matrix< int > > levelRef
std::vector< MatrixDouble > levelShapeFunctions
std::vector< MatrixDouble > levelGaussPtsOnRefMesh
double zeta
Viscous hardening.