405 using namespace MoFEM;
409 static char help[] =
"-my_block_config set block data\n"
423 using namespace boost::numeric;
426 #include <boost/program_options.hpp>
428 namespace po = boost::program_options;
432 int main(
int argc,
char *argv[]) {
441 PetscBool flg_block_config, flg_file;
443 char block_config_file[255];
444 PetscBool flg_order_force;
446 PetscInt order_force = 2;
447 PetscBool flg_eps_u, flg_eps_rho, flg_eps_l;
449 double eps_rho = 1e-3;
451 PetscBool is_curl = PETSC_TRUE;
452 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
453 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
455 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
459 "default approximation order for force approximation",
"", order_force,
460 &order_force, &flg_order_force);
461 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
462 "",
"block_conf.in", block_config_file, 255,
464 CHKERR PetscOptionsReal(
"-my_eps_rho",
"foce optimality parameter ",
"",
465 eps_rho, &eps_rho, &flg_eps_rho);
466 CHKERR PetscOptionsReal(
"-my_eps_u",
"displacement optimality parameter ",
467 "", eps_u, &eps_u, &flg_eps_u);
468 CHKERR PetscOptionsReal(
"-my_eps_l",
"displacement optimality parameter ",
469 "", eps_l, &eps_l, &flg_eps_l);
470 CHKERR PetscOptionsBool(
"-my_curl",
"use Hcurl space to approximate forces",
471 "", is_curl, &is_curl, NULL);
472 ierr = PetscOptionsEnd();
476 if (flg_file != PETSC_TRUE) {
478 "*** ERROR -my_file (MESH FILE NEEDED)");
481 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
483 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
487 option =
"PARALLEL=READ_PART;"
488 "PARALLEL_RESOLVE_SHARED_ENTS;"
489 "PARTITION=PARALLEL_PARTITION;";
510 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
516 std::vector<Range> setOrderToEnts(10);
520 Range set_order_ents;
521 std::map<int, BlockOptionData> block_data;
522 if (flg_block_config) {
523 double read_eps_u, read_eps_rho, read_eps_l;
525 ifstream ini_file(block_config_file);
526 if (!ini_file.is_open()) {
527 SETERRQ(PETSC_COMM_SELF, 1,
528 "*** -my_block_config does not exist ***");
531 po::variables_map vm;
532 po::options_description config_file_options;
533 config_file_options.add_options()(
534 "eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
535 "eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
536 "eps_l", po::value<double>(&read_eps_l)->default_value(-1));
538 std::ostringstream str_order;
539 str_order <<
"block_" << it->getMeshsetId() <<
".displacement_order";
540 config_file_options.add_options()(
541 str_order.str().c_str(),
542 po::value<int>(&block_data[it->getMeshsetId()].oRder)
543 ->default_value(
order));
544 std::ostringstream str_cond;
545 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
546 config_file_options.add_options()(
547 str_cond.str().c_str(),
548 po::value<double>(&block_data[it->getMeshsetId()].yOung)
549 ->default_value(-1));
550 std::ostringstream str_capa;
551 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
552 config_file_options.add_options()(
553 str_capa.str().c_str(),
554 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
555 ->default_value(-2));
557 po::parsed_options parsed =
558 parse_config_file(ini_file, config_file_options,
true);
562 if (block_data[it->getMeshsetId()].oRder == -1)
564 if (block_data[it->getMeshsetId()].oRder ==
order)
566 PetscPrintf(PETSC_COMM_WORLD,
"Set block %d order to %d\n",
567 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
569 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
573 CHKERR moab.get_connectivity(block_ents, nodes,
true);
574 Range ents_to_set_order, ents3d;
575 CHKERR moab.get_adjacencies(nodes, 3,
false, ents3d,
576 moab::Interface::UNION);
577 CHKERR moab.get_adjacencies(ents3d, 2,
false, ents_to_set_order,
578 moab::Interface::UNION);
579 CHKERR moab.get_adjacencies(ents3d, 1,
false, ents_to_set_order,
580 moab::Interface::UNION);
581 ents_to_set_order = subtract(
582 ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
583 ents_to_set_order = subtract(
584 ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
585 set_order_ents.merge(ents3d);
586 set_order_ents.merge(ents_to_set_order);
587 setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
591 std::vector<std::string> additional_parameters;
592 additional_parameters =
593 collect_unrecognized(parsed.options, po::include_positional);
594 for (std::vector<std::string>::iterator vit =
595 additional_parameters.begin();
596 vit != additional_parameters.end(); vit++) {
597 CHKERR PetscPrintf(PETSC_COMM_WORLD,
598 "** WARNING Unrecognized option %s\n",
601 }
catch (
const std::exception &ex) {
602 std::ostringstream ss;
603 ss << ex.what() << std::endl;
606 if (read_eps_u > 0) {
609 if (read_eps_rho > 0) {
610 eps_rho = read_eps_rho;
612 if (read_eps_l > 0) {
617 PetscPrintf(PETSC_COMM_WORLD,
"epsU = %6.4e epsRho = %6.4e\n", eps_u,
643 Range ents_1st_layer;
647 ents_1st_layer,
true);
649 vertex_to_fix,
false);
652 if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
654 "Should be one vertex only, but is %d", vertex_to_fix.size());
659 ents_1st_layer.subset_by_type(MBTRI), MBTRI,
"RHO");
660 Range ents_2nd_layer;
664 ents_2nd_layer,
true);
668 for (
int oo = 2; oo != setOrderToEnts.size(); oo++) {
669 if (setOrderToEnts[oo].size() > 0) {
676 const int through_thickness_order = 2;
679 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
681 CHKERR moab.get_adjacencies(ents3d, 2,
false, ents,
682 moab::Interface::UNION);
683 CHKERR moab.get_adjacencies(ents3d, 1,
false, ents,
684 moab::Interface::UNION);
687 CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
690 CHKERR moab.get_adjacencies(prisms, 2,
false, quads,
691 moab::Interface::UNION);
693 prism_tris = quads.subset_by_type(MBTRI);
694 quads = subtract(quads, prism_tris);
696 CHKERR moab.get_adjacencies(quads, 1,
false, quads_edges,
697 moab::Interface::UNION);
698 Range prism_tris_edges;
699 CHKERR moab.get_adjacencies(prism_tris, 1,
false, prism_tris_edges,
700 moab::Interface::UNION);
701 quads_edges = subtract(quads_edges, prism_tris_edges);
703 prisms.merge(quads_edges);
707 ents = subtract(ents, set_order_ents);
708 ents = subtract(ents, prisms);
718 through_thickness_order);
733 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>());
734 boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>());
747 Range ents_2nd_layer;
751 elastic.
setOfBlocks[2].tEts, MBPRISM,
"ELASTIC_PRISM");
787 int id = it->getMeshsetId();
788 if (block_data[
id].yOung > 0) {
790 CHKERR PetscPrintf(PETSC_COMM_WORLD,
791 "Block %d set Young modulus %3.4g\n",
id,
794 if (block_data[
id].pOisson >= -1) {
795 elastic.
setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
796 CHKERR PetscPrintf(PETSC_COMM_WORLD,
797 "Block %d set Poisson ratio %3.4g\n",
id,
827 "DISPLACEMENTS_PENALTY");
864 DMType dm_name =
"MOFEM";
870 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
871 CHKERR DMSetType(dm_control, dm_name);
875 CHKERR DMSetFromOptions(dm_control);
886 CHKERR DMSetUp(dm_control);
892 ublas::matrix<Mat> nested_matrices(2, 2);
893 ublas::vector<IS> nested_is_rows(2);
894 ublas::vector<IS> nested_is_cols(2);
895 for (
int i = 0;
i != 2;
i++) {
896 nested_is_rows[
i] = PETSC_NULL;
897 nested_is_cols[
i] = PETSC_NULL;
898 for (
int j = 0;
j != 2;
j++) {
899 nested_matrices(
i,
j) = PETSC_NULL;
903 ublas::matrix<Mat> sub_nested_matrices(2, 2);
904 ublas::vector<IS> sub_nested_is_rows(2);
905 ublas::vector<IS> sub_nested_is_cols(2);
906 for (
int i = 0;
i != 2;
i++) {
907 sub_nested_is_rows[
i] = PETSC_NULL;
908 sub_nested_is_cols[
i] = PETSC_NULL;
909 for (
int j = 0;
j != 2;
j++) {
910 sub_nested_matrices(
i,
j) = PETSC_NULL;
914 DM dm_sub_volume_control;
916 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
917 CHKERR DMSetType(dm_sub_volume_control, dm_name);
928 CHKERR DMSetUp(dm_sub_volume_control);
932 boost::shared_ptr<Problem::SubProblemData> sub_data =
935 CHKERR sub_data->getRowIs(&nested_is_rows[0]);
936 CHKERR sub_data->getColIs(&nested_is_cols[0]);
938 nested_matrices(0, 0) = PETSC_NULL;
942 DM dm_sub_sub_elastic;
944 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
945 CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
956 CHKERR DMSetUp(dm_sub_sub_elastic);
960 CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
961 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
962 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
963 CHKERR MatZeroEntries(Kuu);
964 CHKERR VecZeroEntries(Du);
965 CHKERR VecZeroEntries(Fu);
970 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
971 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
977 dirihlet_bc_ptr.get());
980 elastic.getLoopFeRhs().snes_f = Fu;
981 fat_prism_rhs.
snes_f = Fu;
983 &elastic.getLoopFeRhs());
987 elastic.getLoopFeLhs().snes_B = Kuu;
988 fat_prism_lhs.
snes_B = Kuu;
990 &elastic.getLoopFeLhs());
995 dirihlet_bc_ptr.get());
996 CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
997 CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
998 CHKERR VecAssemblyBegin(Fu);
999 CHKERR VecAssemblyEnd(Fu);
1006 boost::shared_ptr<Problem::SubProblemData> sub_data =
1009 CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
1010 CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1011 sub_nested_matrices(0, 0) = Kuu;
1014 ->isCreateFromProblemFieldToOtherProblemField(
1015 "ELASTIC_PROB",
"U",
ROW,
"SUB_CONTROL_PROB",
"UPSILON",
ROW,
1016 PETSC_NULL, &isUpsilon);
1017 sub_nested_is_rows[1] = isUpsilon;
1018 sub_nested_is_cols[1] = isUpsilon;
1019 sub_nested_matrices(1, 1) = Kuu;
1020 PetscObjectReference((PetscObject)Kuu);
1021 PetscObjectReference((PetscObject)isUpsilon);
1025 cerr <<
"Kuu" << endl;
1026 MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
1031 CHKERR DMDestroy(&dm_sub_sub_elastic);
1035 DM dm_sub_disp_penalty;
1038 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
1039 CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
1049 CHKERR DMSetUp(dm_sub_disp_penalty);
1052 CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
1053 CHKERR MatZeroEntries(S);
1058 "DISPLACEMENTS_PENALTY", &face_element);
1059 CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
1060 CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
1064 cerr <<
"S" << endl;
1065 MatView(S, PETSC_VIEWER_DRAW_WORLD);
1072 boost::shared_ptr<Problem::SubProblemData> sub_data =
1076 sub_nested_matrices(1, 0) = S;
1078 CHKERR DMDestroy(&dm_sub_disp_penalty);
1083 DM dm_sub_force_penalty;
1086 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
1087 CHKERR DMSetType(dm_sub_force_penalty, dm_name);
1096 CHKERR DMSetUp(dm_sub_force_penalty);
1099 CHKERR DMCreateMatrix(dm_sub_force_penalty, &
D);
1105 new OpCalculateInvJacForFace(inv_jac));
1120 CHKERR MatAssemblyBegin(
D, MAT_FINAL_ASSEMBLY);
1121 CHKERR MatAssemblyEnd(
D, MAT_FINAL_ASSEMBLY);
1128 if (is_curl == PETSC_FALSE) {
1129 int nb_dofs_to_fix = 0;
1130 int index_to_fix = 0;
1131 if (!vertex_to_fix.empty()) {
1132 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1139 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1140 cerr << *dof_ptr << endl;
1144 CHKERR MatZeroRowsColumns(
D, nb_dofs_to_fix, &index_to_fix,
1145 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1147 std::vector<int> dofs_to_fix;
1148 for (
auto p_eit = edges_to_fix.pair_begin();
1149 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1152 auto lo = row_dofs->lower_bound(
1156 bit_number, p_eit->second));
1157 for (; lo != hi; ++lo)
1158 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1160 CHKERR MatZeroRowsColumns(
D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1161 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1166 cerr <<
"D" << endl;
1167 MatView(
D, PETSC_VIEWER_DRAW_WORLD);
1172 boost::shared_ptr<Problem::SubProblemData> sub_data =
1174 CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1175 CHKERR sub_data->getColIs(&nested_is_cols[1]);
1176 nested_matrices(1, 1) =
D;
1178 CHKERR DMDestroy(&dm_sub_force_penalty);
1185 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1186 CHKERR DMSetType(dm_sub_force, dm_name);
1195 CHKERR DMSetUp(dm_sub_force);
1198 CHKERR DMCreateMatrix(dm_sub_force, &UB);
1199 CHKERR MatZeroEntries(UB);
1201 CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1202 CHKERR MatZeroEntries(UPSILONB);
1206 new OpCalculateInvJacForFace(inv_jac));
1225 CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1226 CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1227 CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1228 CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1235 if (is_curl == PETSC_FALSE) {
1236 int nb_dofs_to_fix = 0;
1237 int index_to_fix = 0;
1238 if (!vertex_to_fix.empty()) {
1239 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1246 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1247 cerr << *dof_ptr << endl;
1251 CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1253 CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1254 PETSC_NULL, PETSC_NULL);
1256 std::vector<int> dofs_to_fix;
1257 for (
auto p_eit = edges_to_fix.pair_begin();
1258 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1261 auto lo = row_dofs->lower_bound(
1265 bit_number, p_eit->second));
1266 for (; lo != hi; ++lo)
1267 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1269 CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1270 PETSC_NULL, PETSC_NULL);
1271 CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1272 0, PETSC_NULL, PETSC_NULL);
1276 CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1281 cerr <<
"UBT" << endl;
1282 MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1287 boost::shared_ptr<Problem::SubProblemData> sub_data =
1291 nested_matrices(0, 1) = UBT;
1294 cerr <<
"UPSILONB" << endl;
1295 MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1302 nested_matrices(1, 0) = UPSILONB;
1304 CHKERR DMDestroy(&dm_sub_force);
1308 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1309 &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1311 nested_matrices(0, 0) = SubA;
1313 CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1314 CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1317 cerr <<
"Nested SubA" << endl;
1318 MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1322 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1323 &nested_is_cols[0], &nested_matrices(0, 0), &
A);
1325 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
1326 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
1329 cerr <<
"Nested A" << endl;
1330 MatView(
A, PETSC_VIEWER_STDOUT_WORLD);
1334 CHKERR DMCreateGlobalVector(dm_control, &
D);
1335 CHKERR DMCreateGlobalVector(dm_control, &
F);
1346 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
1347 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
1355 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1356 CHKERR KSPSetDM(solver, dm_control);
1357 CHKERR KSPSetFromOptions(solver);
1358 CHKERR KSPSetOperators(solver,
A,
A);
1359 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1360 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1361 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1363 CHKERR KSPGetPC(solver, &pc);
1364 CHKERR PCSetType(pc, PCFIELDSPLIT);
1365 PetscBool is_pcfs = PETSC_FALSE;
1366 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1369 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1370 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1371 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1375 CHKERR PCFieldSplitGetSubKSP(pc, &
n, &sub_ksp);
1378 CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1379 CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1380 CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1381 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1382 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1383 CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1387 CHKERR PCSetUp(sub_pc_0);
1391 "Only works with pre-conditioner PCFIELDSPLIT");
1399 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1400 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1405 CHKERR VecView(
D, PETSC_VIEWER_DRAW_WORLD);
1411 for (
int i = 0;
i != 2;
i++) {
1412 if (sub_nested_is_rows[
i]) {
1413 CHKERR ISDestroy(&sub_nested_is_rows[
i]);
1415 if (sub_nested_is_cols[
i]) {
1416 CHKERR ISDestroy(&sub_nested_is_cols[
i]);
1418 for (
int j = 0;
j != 2;
j++) {
1419 if (sub_nested_matrices(
i,
j)) {
1420 CHKERR MatDestroy(&sub_nested_matrices(
i,
j));
1424 for (
int i = 0;
i != 2;
i++) {
1425 if (nested_is_rows[
i]) {
1426 CHKERR ISDestroy(&nested_is_rows[
i]);
1428 if (nested_is_cols[
i]) {
1429 CHKERR ISDestroy(&nested_is_cols[
i]);
1431 for (
int j = 0;
j != 2;
j++) {
1432 if (nested_matrices(
i,
j)) {
1433 CHKERR MatDestroy(&nested_matrices(
i,
j));
1438 CHKERR MatDestroy(&SubA);
1443 CHKERR DMDestroy(&dm_sub_volume_control);
1457 elastic.getLoopFeEnergy().eNergy = 0;
1459 &elastic.getLoopFeEnergy());
1460 PetscPrintf(PETSC_COMM_WORLD,
"Elastic energy %6.4e\n",
1461 elastic.getLoopFeEnergy().eNergy);
1468 new OpCalculateInvJacForFace(inv_jac));
1487 CHKERR DMDestroy(&dm_control);