432int main(
int argc,
char *argv[]) {
438 moab::Core mb_instance;
439 moab::Interface &moab = mb_instance;
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;";
500 CHKERR mmanager_ptr->printForceSet();
502 CHKERR mmanager_ptr->printMaterialsSet();
509 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
515 std::vector<Range> setOrderToEnts(10);
519 Range set_order_ents;
520 std::map<int, BlockOptionData> block_data;
521 if (flg_block_config) {
522 double read_eps_u, read_eps_rho, read_eps_l;
524 ifstream ini_file(block_config_file);
525 if (!ini_file.is_open()) {
526 SETERRQ(PETSC_COMM_SELF, 1,
527 "*** -my_block_config does not exist ***");
530 po::variables_map vm;
531 po::options_description config_file_options;
532 config_file_options.add_options()(
533 "eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
534 "eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
535 "eps_l", po::value<double>(&read_eps_l)->default_value(-1));
537 std::ostringstream str_order;
538 str_order <<
"block_" << it->getMeshsetId() <<
".displacement_order";
539 config_file_options.add_options()(
540 str_order.str().c_str(),
541 po::value<int>(&block_data[it->getMeshsetId()].oRder)
542 ->default_value(
order));
543 std::ostringstream str_cond;
544 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
545 config_file_options.add_options()(
546 str_cond.str().c_str(),
547 po::value<double>(&block_data[it->getMeshsetId()].yOung)
548 ->default_value(-1));
549 std::ostringstream str_capa;
550 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
551 config_file_options.add_options()(
552 str_capa.str().c_str(),
553 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
554 ->default_value(-2));
556 po::parsed_options parsed =
557 parse_config_file(ini_file, config_file_options,
true);
561 if (block_data[it->getMeshsetId()].oRder == -1)
563 if (block_data[it->getMeshsetId()].oRder ==
order)
565 PetscPrintf(PETSC_COMM_WORLD,
"Set block %d order to %d\n",
566 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
568 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
572 CHKERR moab.get_connectivity(block_ents, nodes,
true);
573 Range ents_to_set_order, ents3d;
574 CHKERR moab.get_adjacencies(nodes, 3,
false, ents3d,
575 moab::Interface::UNION);
576 CHKERR moab.get_adjacencies(ents3d, 2,
false, ents_to_set_order,
577 moab::Interface::UNION);
578 CHKERR moab.get_adjacencies(ents3d, 1,
false, ents_to_set_order,
579 moab::Interface::UNION);
580 ents_to_set_order = subtract(
581 ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
582 ents_to_set_order = subtract(
583 ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
584 set_order_ents.merge(ents3d);
585 set_order_ents.merge(ents_to_set_order);
586 setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
590 std::vector<std::string> additional_parameters;
591 additional_parameters =
592 collect_unrecognized(parsed.options, po::include_positional);
593 for (std::vector<std::string>::iterator vit =
594 additional_parameters.begin();
595 vit != additional_parameters.end(); vit++) {
596 CHKERR PetscPrintf(PETSC_COMM_WORLD,
597 "** WARNING Unrecognized option %s\n",
600 }
catch (
const std::exception &ex) {
601 std::ostringstream ss;
602 ss << ex.what() << std::endl;
605 if (read_eps_u > 0) {
608 if (read_eps_rho > 0) {
609 eps_rho = read_eps_rho;
611 if (read_eps_l > 0) {
616 PetscPrintf(PETSC_COMM_WORLD,
"epsU = %6.4e epsRho = %6.4e\n", eps_u,
642 Range ents_1st_layer;
644 if (mmanager_ptr->checkMeshset(202,
SIDESET)) {
646 ents_1st_layer,
true);
648 vertex_to_fix,
false);
649 CHKERR mmanager_ptr->getEntitiesByDimension(202,
SIDESET, 1, edges_to_fix,
651 if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
653 "Should be one vertex only, but is %d", vertex_to_fix.size());
658 ents_1st_layer.subset_by_type(MBTRI), MBTRI,
"RHO");
659 Range ents_2nd_layer;
661 if (mmanager_ptr->checkMeshset(101,
SIDESET)) {
663 ents_2nd_layer,
true);
667 for (
int oo = 2; oo != setOrderToEnts.size(); oo++) {
668 if (setOrderToEnts[oo].size() > 0) {
675 const int through_thickness_order = 2;
678 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
680 CHKERR moab.get_adjacencies(ents3d, 2,
false, ents,
681 moab::Interface::UNION);
682 CHKERR moab.get_adjacencies(ents3d, 1,
false, ents,
683 moab::Interface::UNION);
686 CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
689 CHKERR moab.get_adjacencies(prisms, 2,
false, quads,
690 moab::Interface::UNION);
692 prism_tris = quads.subset_by_type(MBTRI);
693 quads = subtract(quads, prism_tris);
695 CHKERR moab.get_adjacencies(quads, 1,
false, quads_edges,
696 moab::Interface::UNION);
697 Range prism_tris_edges;
698 CHKERR moab.get_adjacencies(prism_tris, 1,
false, prism_tris_edges,
699 moab::Interface::UNION);
700 quads_edges = subtract(quads_edges, prism_tris_edges);
702 prisms.merge(quads_edges);
706 ents = subtract(ents, set_order_ents);
707 ents = subtract(ents, prisms);
717 through_thickness_order);
732 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>());
733 boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>());
746 Range ents_2nd_layer;
750 elastic.
setOfBlocks[2].tEts, MBPRISM,
"ELASTIC_PRISM");
752 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
786 int id = it->getMeshsetId();
787 if (block_data[
id].yOung > 0) {
789 CHKERR PetscPrintf(PETSC_COMM_WORLD,
790 "Block %d set Young modulus %3.4g\n",
id,
793 if (block_data[
id].pOisson >= -1) {
794 elastic.
setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
795 CHKERR PetscPrintf(PETSC_COMM_WORLD,
796 "Block %d set Poisson ratio %3.4g\n",
id,
826 "DISPLACEMENTS_PENALTY");
863 DMType dm_name =
"MOFEM";
869 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
870 CHKERR DMSetType(dm_control, dm_name);
874 CHKERR DMSetFromOptions(dm_control);
885 CHKERR DMSetUp(dm_control);
890 ublas::matrix<Mat> nested_matrices(2, 2);
891 ublas::vector<IS> nested_is_rows(2);
892 ublas::vector<IS> nested_is_cols(2);
893 for (
int i = 0;
i != 2;
i++) {
894 nested_is_rows[
i] = PETSC_NULL;
895 nested_is_cols[
i] = PETSC_NULL;
896 for (
int j = 0;
j != 2;
j++) {
897 nested_matrices(
i,
j) = PETSC_NULL;
901 ublas::matrix<Mat> sub_nested_matrices(2, 2);
902 ublas::vector<IS> sub_nested_is_rows(2);
903 ublas::vector<IS> sub_nested_is_cols(2);
904 for (
int i = 0;
i != 2;
i++) {
905 sub_nested_is_rows[
i] = PETSC_NULL;
906 sub_nested_is_cols[
i] = PETSC_NULL;
907 for (
int j = 0;
j != 2;
j++) {
908 sub_nested_matrices(
i,
j) = PETSC_NULL;
912 DM dm_sub_volume_control;
914 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
915 CHKERR DMSetType(dm_sub_volume_control, dm_name);
926 CHKERR DMSetUp(dm_sub_volume_control);
930 boost::shared_ptr<Problem::SubProblemData> sub_data =
933 CHKERR sub_data->getRowIs(&nested_is_rows[0]);
934 CHKERR sub_data->getColIs(&nested_is_cols[0]);
936 nested_matrices(0, 0) = PETSC_NULL;
940 DM dm_sub_sub_elastic;
942 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
943 CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
954 CHKERR DMSetUp(dm_sub_sub_elastic);
958 CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
959 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
960 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
961 CHKERR MatZeroEntries(Kuu);
962 CHKERR VecZeroEntries(Du);
963 CHKERR VecZeroEntries(Fu);
968 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
969 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
971 dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
972 dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
975 dirihlet_bc_ptr.get());
979 fat_prism_rhs.
snes_f = Fu;
986 fat_prism_lhs.
snes_B = Kuu;
993 dirihlet_bc_ptr.get());
994 CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
995 CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
996 CHKERR VecAssemblyBegin(Fu);
997 CHKERR VecAssemblyEnd(Fu);
1004 boost::shared_ptr<Problem::SubProblemData> sub_data =
1007 CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
1008 CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1009 sub_nested_matrices(0, 0) = Kuu;
1012 ->isCreateFromProblemFieldToOtherProblemField(
1013 "ELASTIC_PROB",
"U",
ROW,
"SUB_CONTROL_PROB",
"UPSILON",
ROW,
1014 PETSC_NULL, &isUpsilon);
1015 sub_nested_is_rows[1] = isUpsilon;
1016 sub_nested_is_cols[1] = isUpsilon;
1017 sub_nested_matrices(1, 1) = Kuu;
1018 PetscObjectReference((PetscObject)Kuu);
1019 PetscObjectReference((PetscObject)isUpsilon);
1023 cerr <<
"Kuu" << endl;
1024 MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
1029 CHKERR DMDestroy(&dm_sub_sub_elastic);
1033 DM dm_sub_disp_penalty;
1036 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
1037 CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
1047 CHKERR DMSetUp(dm_sub_disp_penalty);
1050 CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
1051 CHKERR MatZeroEntries(S);
1056 "DISPLACEMENTS_PENALTY", &face_element);
1057 CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
1058 CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
1062 cerr <<
"S" << endl;
1063 MatView(S, PETSC_VIEWER_DRAW_WORLD);
1070 boost::shared_ptr<Problem::SubProblemData> sub_data =
1074 sub_nested_matrices(1, 0) = S;
1076 CHKERR DMDestroy(&dm_sub_disp_penalty);
1081 DM dm_sub_force_penalty;
1084 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
1085 CHKERR DMSetType(dm_sub_force_penalty, dm_name);
1094 CHKERR DMSetUp(dm_sub_force_penalty);
1097 CHKERR DMCreateMatrix(dm_sub_force_penalty, &
D);
1100 auto det_ptr = boost::make_shared<VectorDouble>();
1101 auto jac_ptr = boost::make_shared<MatrixDouble>();
1102 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1126 CHKERR MatAssemblyBegin(
D, MAT_FINAL_ASSEMBLY);
1127 CHKERR MatAssemblyEnd(
D, MAT_FINAL_ASSEMBLY);
1134 if (is_curl == PETSC_FALSE) {
1135 int nb_dofs_to_fix = 0;
1136 int index_to_fix = 0;
1137 if (!vertex_to_fix.empty()) {
1138 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1145 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1146 cerr << *dof_ptr << endl;
1150 CHKERR MatZeroRowsColumns(
D, nb_dofs_to_fix, &index_to_fix,
1151 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1153 std::vector<int> dofs_to_fix;
1154 for (
auto p_eit = edges_to_fix.pair_begin();
1155 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1158 auto lo = row_dofs->lower_bound(
1162 bit_number, p_eit->second));
1163 for (; lo != hi; ++lo)
1165 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1167 CHKERR MatZeroRowsColumns(
D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1168 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1173 cerr <<
"D" << endl;
1174 MatView(
D, PETSC_VIEWER_DRAW_WORLD);
1179 boost::shared_ptr<Problem::SubProblemData> sub_data =
1181 CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1182 CHKERR sub_data->getColIs(&nested_is_cols[1]);
1183 nested_matrices(1, 1) =
D;
1185 CHKERR DMDestroy(&dm_sub_force_penalty);
1192 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1193 CHKERR DMSetType(dm_sub_force, dm_name);
1202 CHKERR DMSetUp(dm_sub_force);
1205 CHKERR DMCreateMatrix(dm_sub_force, &UB);
1206 CHKERR MatZeroEntries(UB);
1208 CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1209 CHKERR MatZeroEntries(UPSILONB);
1212 auto det_ptr = boost::make_shared<VectorDouble>();
1213 auto jac_ptr = boost::make_shared<MatrixDouble>();
1214 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1237 CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1238 CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1239 CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1240 CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1247 if (is_curl == PETSC_FALSE) {
1248 int nb_dofs_to_fix = 0;
1249 int index_to_fix = 0;
1250 if (!vertex_to_fix.empty()) {
1251 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1258 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1259 cerr << *dof_ptr << endl;
1263 CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1265 CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1266 PETSC_NULL, PETSC_NULL);
1268 std::vector<int> dofs_to_fix;
1269 for (
auto p_eit = edges_to_fix.pair_begin();
1270 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1273 auto lo = row_dofs->lower_bound(
1277 bit_number, p_eit->second));
1278 for (; lo != hi; ++lo)
1280 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1282 CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1283 PETSC_NULL, PETSC_NULL);
1284 CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1285 0, PETSC_NULL, PETSC_NULL);
1289 CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1294 cerr <<
"UBT" << endl;
1295 MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1300 boost::shared_ptr<Problem::SubProblemData> sub_data =
1304 nested_matrices(0, 1) = UBT;
1307 cerr <<
"UPSILONB" << endl;
1308 MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1315 nested_matrices(1, 0) = UPSILONB;
1317 CHKERR DMDestroy(&dm_sub_force);
1321 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1322 &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1324 nested_matrices(0, 0) = SubA;
1326 CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1327 CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1330 cerr <<
"Nested SubA" << endl;
1331 MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1335 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1336 &nested_is_cols[0], &nested_matrices(0, 0), &A);
1338 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1339 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1342 cerr <<
"Nested A" << endl;
1343 MatView(A, PETSC_VIEWER_STDOUT_WORLD);
1347 CHKERR DMCreateGlobalVector(dm_control, &
D);
1348 CHKERR DMCreateGlobalVector(dm_control, &
F);
1359 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
1360 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
1368 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1369 CHKERR KSPSetDM(solver, dm_control);
1370 CHKERR KSPSetFromOptions(solver);
1371 CHKERR KSPSetOperators(solver, A, A);
1372 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1373 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1374 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1376 CHKERR KSPGetPC(solver, &pc);
1377 CHKERR PCSetType(pc, PCFIELDSPLIT);
1378 PetscBool is_pcfs = PETSC_FALSE;
1379 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1381 CHKERR PCSetOperators(pc, A, A);
1382 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1383 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1384 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1388 CHKERR PCFieldSplitGetSubKSP(pc, &
n, &sub_ksp);
1391 CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1392 CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1393 CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1394 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1395 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1396 CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1400 CHKERR PCSetUp(sub_pc_0);
1404 "Only works with pre-conditioner PCFIELDSPLIT");
1412 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1413 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1418 CHKERR VecView(
D, PETSC_VIEWER_DRAW_WORLD);
1424 for (
int i = 0;
i != 2;
i++) {
1425 if (sub_nested_is_rows[
i]) {
1426 CHKERR ISDestroy(&sub_nested_is_rows[
i]);
1428 if (sub_nested_is_cols[
i]) {
1429 CHKERR ISDestroy(&sub_nested_is_cols[
i]);
1431 for (
int j = 0;
j != 2;
j++) {
1432 if (sub_nested_matrices(
i,
j)) {
1433 CHKERR MatDestroy(&sub_nested_matrices(
i,
j));
1437 for (
int i = 0;
i != 2;
i++) {
1438 if (nested_is_rows[
i]) {
1439 CHKERR ISDestroy(&nested_is_rows[
i]);
1441 if (nested_is_cols[
i]) {
1442 CHKERR ISDestroy(&nested_is_cols[
i]);
1444 for (
int j = 0;
j != 2;
j++) {
1445 if (nested_matrices(
i,
j)) {
1446 CHKERR MatDestroy(&nested_matrices(
i,
j));
1451 CHKERR MatDestroy(&SubA);
1456 CHKERR DMDestroy(&dm_sub_volume_control);
1473 PetscPrintf(PETSC_COMM_WORLD,
"Elastic energy %6.4e\n",
1480 auto det_ptr = boost::make_shared<VectorDouble>();
1481 auto jac_ptr = boost::make_shared<MatrixDouble>();
1482 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1503 CHKERR DMDestroy(&dm_control);