25#include <boost/tokenizer.hpp>
32 chrono::high_resolution_clock::time_point
start;
33 chrono::high_resolution_clock::time_point
stop;
36 stop = chrono::high_resolution_clock::now();
37 auto duration = chrono::duration_cast<chrono::microseconds>(
stop -
start);
38 cout <<
"Time taken by function: " << duration.count() <<
" us." << endl;
250 boost::shared_ptr<OpContactTools::CommonData> common_data_ptr,
251 boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_proc_fe,
264 boost::make_shared<PostProcFaceOnRefinedMesh>(
mField);
298 for (
auto &roller : (*cache).rollerDataVec) {
299 if (roller.methodOpForRollerPosition) {
300 roller.BodyDispScaled.clear();
302 this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
304 if (roller.methodOpForRollerDirection) {
305 roller.BodyDirectionScaled.clear();
307 this, roller.methodOpForRollerDirection,
308 roller.BodyDirectionScaled);
320 auto make_vtk = [&]() {
326 boost::lexical_cast<std::string>(
ts_step) +
331 auto save_skeleton = [&]() {
336 "out_skeleton_" + boost::lexical_cast<std::string>(
ts_step) +
341 auto save_skin = [&]() {
348 "out_skin_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
352 double max_disp, min_disp;
354 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
356 CHKERR VecScatterBegin(std::get<1>(tuple),
ts_u, std::get<0>(tuple),
357 INSERT_VALUES, SCATTER_FORWARD);
358 CHKERR VecScatterEnd(std::get<1>(tuple),
ts_u, std::get<0>(tuple),
359 INSERT_VALUES, SCATTER_FORWARD);
360 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max_disp);
361 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min_disp);
362 PetscPrintf(PETSC_COMM_WORLD,
"%s time %3.4e min %3.4e max %3.4e\n",
363 msg.c_str(),
ts_t, min_disp, max_disp);
373 int roller_nb = (*cache).rollerDataVec.size();
374 l_reactions.resize(roller_nb * 3);
375 std::fill(l_reactions.begin(), l_reactions.end(), 0);
380 std::vector<double> gauss_pts(2, 0);
382 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2,
384 PetscObjectComm((PetscObject)dm));
386 "\t TS %d Total contact area (ref): %g / %g",
ts_step,
387 gauss_pts[0], gauss_pts[1]);
388 gauss_area = gauss_pts[0];
389 lgauss_pts[0] = lgauss_pts[1] = 0;
391 vector<double> g_reactions(l_reactions.size() * 3, 0);
392 CHKERR MPIU_Allreduce(l_reactions.data(), g_reactions.data(),
393 roller_nb * 3, MPIU_REAL, MPIU_SUM,
394 PetscObjectComm((PetscObject)dm));
395 for (
int dd = 0;
dd != roller_nb; ++
dd)
397 "\t Body %d reactions Time %g Force X: %3.4e Y: "
399 (*cache).rollerDataVec[
dd].iD,
ts_t,
400 g_reactions[3 *
dd + 0], g_reactions[3 *
dd + 1],
401 g_reactions[3 *
dd + 2]);
406 std::ostringstream ostrm, ostrm2;
407 ostrm <<
"out_debug_" <<
ts_step <<
".vtk";
408 ostrm2 <<
"out_debug_" <<
ts_step <<
".h5m";
414 "PARALLEL=WRITE_PART");
422 CHKERR VecZeroEntries(vec);
423 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
424 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
430 ParallelComm *pcomm =
438 double def_value = 0.;
444 CHKERR VecAssemblyBegin(vec);
445 CHKERR VecAssemblyEnd(vec);
446 CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
447 CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
448 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
449 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
451 const double *react_ptr;
452 CHKERR VecGetArrayRead(vec, &react_ptr);
454 PetscPrintf(PETSC_COMM_WORLD,
455 "Reactions time %3.4e X: %3.4e Y: %3.4e Z: %3.4e \n",
ts_t,
456 react_ptr[0], react_ptr[1], react_ptr[2]);
458 CHKERR VecRestoreArrayRead(vec, &react_ptr);
471 if (m_field.get_comm_rank() == 0) {
472 std::ostringstream ost;
473 for (
auto &rol : (*cache).rollerDataVec) {
476 roller_coords += rol.originCoords;
482 ost <<
"out_roller_" <<
ts_step <<
".vtk";
490 CHKERR VecGhostUpdateBegin(
ts_u, INSERT_VALUES, SCATTER_FORWARD);
491 CHKERR VecGhostUpdateEnd(
ts_u, INSERT_VALUES, SCATTER_FORWARD);
494 "restart_" + to_string(
ts_step) +
"_" + to_string(
ts_t) +
".dat";
496 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, rest_name.c_str(),
497 FILE_MODE_WRITE, &viewer);
498 CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
499 VecView(
ts_u, viewer);
500 PetscViewerDestroy(&viewer);
510 const double *react_ptr;
512 if (fabs(react_ptr[0] + 0.12519621572) > 1e-6)
514 "atom test diverged!");
521 const double *react_ptr;
524 if ((0.03663003663 - react_ptr[2]) / 0.03663003663 > 1e-2)
526 "atom test diverged!");
535 if (fabs((min_disp + 2.70) / 2.70) > 5e-2)
537 "atom test diverged!");
543 const double *react_ptr;
548 if (abs(0.056090728 - (react_ptr[2] * 4.0)) / 0.056090728 > 5.0e-2)
550 "atom test diverged!");
554 const double *react_ptr;
557 if (abs(19030.7 + react_ptr[2]) / 19030.7 > 1e-3)
559 "atom test diverged!");
605 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"",
"none");
607 CHKERR PetscOptionsInt(
"-order",
"approximation order",
"",
data.
order,
609 CHKERR PetscOptionsBool(
"-no_output",
"save no output files",
"",
611 CHKERR PetscOptionsBool(
"-debug_info",
"print debug info",
"",
613 CHKERR PetscOptionsBool(
"-scale_params",
614 "scale parameters (young modulus = 1)",
"",
616 CHKERR PetscOptionsInt(
"-output_every",
617 "frequncy how often results are dumped on hard disk",
619 CHKERR PetscOptionsBool(
"-calculate_reactions",
"calculate reactions",
"",
622 CHKERR PetscOptionsInt(
"-reaction_id",
"meshset id for computing reactions",
624 CHKERR PetscOptionsInt(
"-atom_test",
"number of atom test",
"",
626 CHKERR PetscOptionsBool(
"-is_large_strain",
"is large strains regime",
"",
630 "-is_fieldsplit_bc_only",
"use fieldsplit for boundary only",
"",
633 "-is_reduced_integration",
634 "use axiator split approach, with reduced integration",
"",
636 CHKERR PetscOptionsBool(
"-test_user_base_tau",
"test_user_base_for_tau",
"",
639 CHKERR PetscOptionsBool(
"-is_alm",
640 "use Augmented Lagrangian method for enforcing the "
641 "KKT condtions for plasticity",
644 CHKERR PetscOptionsBool(
"-is_neohooke",
"is neohooke material",
"",
651 "-with_plasticity",
"run calculations with plasticity",
"",
654 CHKERR PetscOptionsBool(
"-is_rotating",
"is rotating frame used",
"",
656 CHKERR PetscOptionsBool(
"-with_contact",
"run calculations with contact",
"",
658 CHKERR PetscOptionsBool(
"-print_rollers",
659 "output file with rollers positions",
"",
662 char load_hist_file[255];
663 PetscBool ctg_flag = PETSC_FALSE;
664 CHKERR PetscOptionsString(
"-contact_history",
"load_history.in",
"",
665 "load_history.in", load_hist_file, 255, &ctg_flag);
668 CHKERR PetscOptionsString(
"-move_history",
"load_history.in",
"",
669 "load_history.in", load_hist_file, 255, &ctg_flag);
674 CHKERR PetscOptionsString(
"-dirichlet_history",
"load_history.in",
"",
675 "load_history.in", load_hist_file, 255, &ctg_flag);
680 CHKERR PetscOptionsString(
"-file_name",
"file name",
"",
"mesh.h5m",
684 char restart_file_name[255];
685 CHKERR PetscOptionsString(
"-restart_file",
"restart file name",
"",
686 "restart.dat", restart_file_name, 255,
690 CHKERR PetscOptionsBool(
"-save_restarts",
691 "save restart files at each iteration",
"",
694 ierr = PetscOptionsEnd();
699 "Neohooke material is not supported with plasticity.");
709 CHKERR moab.get_entities_by_dimension(0, 3, tets,
false);
711 CHKERR moab.get_adjacencies(tets, 2,
true, faces, moab::Interface::UNION);
713 CHKERR moab.get_adjacencies(tets, 1,
true, edges, moab::Interface::UNION);
719 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
721 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
722 CHKERR moab.create_meshset(MESHSET_SET, part_set);
723 Tag part_tag = pcomm->part_tag();
728 CHKERR moab.get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag, NULL, 1,
729 tagged_sets, moab::Interface::UNION);
731 for (
auto &mit : tagged_sets) {
733 CHKERR moab.tag_get_data(part_tag, &mit, 1, &part);
735 CHKERR moab.get_entities_by_dimension(mit, 3, proc_ents,
true);
736 CHKERR moab.get_entities_by_handle(mit, all_proc_ents,
true);
738 CHKERR moab.get_entities_by_handle(mit, off_proc_ents,
true);
742 Range proc_ents_skin[4];
743 proc_ents_skin[3] = proc_ents;
745 Range all_tets, all_skin;
746 CHKERR m_field.
get_moab().get_entities_by_dimension(0, 3, all_tets,
false);
747 CHKERR skin.find_skin(0, all_tets,
false, all_skin);
748 CHKERR skin.find_skin(0, proc_ents,
false, proc_ents_skin[2]);
750 proc_ents_skin[2] = subtract(proc_ents_skin[2], all_skin);
751 CHKERR moab.get_adjacencies(proc_ents_skin[2], 1,
false, proc_ents_skin[1],
752 moab::Interface::UNION);
753 CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0],
true);
754 for (
int dd = 0;
dd != 3;
dd++)
755 CHKERR moab.add_entities(part_set, proc_ents_skin[
dd]);
756 CHKERR moab.add_entities(part_set, all_proc_ents);
761 CHKERR moab.get_entities_by_handle(0, all_ents);
762 std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
763 CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
764 &*pstat_tag.begin());
768 CHKERR moab.get_entities_by_handle(part_set, ents_to_keep,
false);
769 off_proc_ents = subtract(off_proc_ents, ents_to_keep);
771 CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets,
false);
772 for (
auto m : meshsets) {
773 CHKERR moab.remove_entities(
m, off_proc_ents);
776 CHKERR moab.delete_entities(off_proc_ents);
778 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
785 CHKERR m_field.
get_moab().get_entities_by_dimension(part_set, 3, ents,
true);
787 ents,
simple->getBitRefLevel(),
false);
788 simple->getMeshSet() = part_set;
795 string block_name =
"MAT_ELASTIC";
799 if (it->getName().compare(0, block_name.length(), block_name) == 0) {
800 std::vector<double> block_data;
801 CHKERR it->getAttributes(block_data);
802 const int id = it->getMeshsetId();
807 mField.
get_moab().get_entities_by_dimension(meshset, 3, tets,
true);
809 for (
auto ent : tets)
813 mat_blocks.at(
id).young_modulus = block_data[0];
815 if (block_data.size() > 2)
817 mat_blocks.at(
id).cn_cont = 1. / block_data[0];
819 if (block_data.size() < 2)
821 "provide an appropriate number of entries (2) parameters for "
822 "MAT_ELASTIC block");
838 auto get_base = [
this](
bool has_hexes =
false) {
844 enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
845 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz",
847 PetscInt choice_base_value = has_hexes ? DEMKOWICZ : AINSWORTH;
849 choice_base_value = DEMKOWICZ;
851 LASBASETOPT, &choice_base_value, PETSC_NULL);
853 switch (choice_base_value) {
857 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
862 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
867 <<
"Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
880 int nb_tets, nb_hexes;
887 if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {
898 CHKERR get_base(nb_hexes);
922 pipeline_mng->setDomainRhsIntegrationRule(
926 for (
string name : {
"TAU",
"EP"}) {
931 auto field_order_table =
932 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
935 auto get_tau_bubble_order_zero = [](
int p) {
return 0; };
937 auto get_tau_bubble_order_tet = [](
int p) ->
int {
938 return std::min(
p + 1, 2);
941 field_order_table[MBVERTEX] = get_tau_bubble_order_zero;
942 field_order_table[MBEDGE] = get_tau_bubble_order_zero;
943 field_order_table[MBTRI] = get_tau_bubble_order_zero;
944 field_order_table[MBTET] = get_tau_bubble_order_tet;
945 const_cast<Field *
>(field_ptr)->rebuildDofsOrderMap();
969 simple->getOtherFiniteElements().push_back(
"FORCE_FE");
970 simple->getOtherFiniteElements().push_back(
"PRESSURE_FE");
985 "MESH_NODE_POSITIONS");
1004 auto get_material_stiffness = [&]() {
1011 (*cache).Is(
i,
j,
k,
l) =
1021 cache = &mit.second;
1022 const double G = (*cache).
young_modulus / (2. * (1. + (*cache).poisson));
1024 (*cache).young_modulus / (3. * (1. - 2. * (*cache).poisson));
1027 double A = 6. * G / (3. *
K + 4. * G);
1031 auto t_D_tmp = getFTensor4DdgFromMat<3, 3, 0>(*
commonDataPtr->mtD);
1033 getFTensor4DdgFromMat<3, 3, 0>(*
commonDataPtr->mtD_Axiator);
1035 getFTensor4DdgFromMat<3, 3, 0>(*
commonDataPtr->mtD_Deviator);
1037 t_D_axiator(
i,
j,
k,
l) =
1038 A * (
K - (2. / 3.) * G) * t_kds(
i,
j) * t_kds(
k,
l);
1039 t_D_deviator(
i,
j,
k,
l) = 2 * G * ((t_kds(
i,
k) ^ t_kds(
j,
l)) / 4.);
1041 t_D_tmp(
i,
j,
k,
l) = t_D_deviator(
i,
j,
k,
l) + t_D_axiator(
i,
j,
k,
l);
1044 (*cache).scale_constraint = 1.;
1052 commonDataPtr = boost::make_shared<OpContactTools::CommonData>();
1056 commonDataPtr->mtD_Axiator = boost::make_shared<MatrixDouble>();
1058 commonDataPtr->mtD_Deviator = boost::make_shared<MatrixDouble>();
1061 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
1062 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
1063 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
1064 commonDataPtr->mPiolaStressPtr = boost::make_shared<MatrixDouble>();
1067 commonDataPtr->mDeformationGradient = boost::make_shared<MatrixDouble>();
1069 commonDataPtr->matEigenVal = boost::make_shared<MatrixDouble>();
1070 commonDataPtr->matEigenVector = boost::make_shared<MatrixDouble>();
1072 commonDataPtr->materialTangent = boost::make_shared<MatrixDouble>();
1073 commonDataPtr->dE_dF_mat = boost::make_shared<MatrixDouble>();
1076 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
1078 boost::make_shared<MatrixDouble>();
1079 commonDataPtr->contactNormalPtr = boost::make_shared<MatrixDouble>();
1080 commonDataPtr->contactGapPtr = boost::make_shared<VectorDouble>();
1081 commonDataPtr->contactGapDiffPtr = boost::make_shared<MatrixDouble>();
1082 commonDataPtr->contactNormalDiffPtr = boost::make_shared<MatrixDouble>();
1084 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
1085 commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
1087 commonDataPtr->plasticSurfacePtr = boost::make_shared<VectorDouble>();
1088 commonDataPtr->plasticFlowPtr = boost::make_shared<MatrixDouble>();
1089 commonDataPtr->plasticTauPtr = boost::make_shared<VectorDouble>();
1090 commonDataPtr->plasticStrainPtr = boost::make_shared<MatrixDouble>();
1092 commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
1093 commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
1095 commonDataPtr->plasticTauJumpPtr = boost::make_shared<VectorDouble>();
1096 commonDataPtr->plasticStrainJumpPtr = boost::make_shared<MatrixDouble>();
1097 commonDataPtr->guidingVelocityPtr = boost::make_shared<MatrixDouble>();
1100 commonDataPtr->plasticGradTauPtr = boost::make_shared<MatrixDouble>();
1101 commonDataPtr->plasticGradStrainPtr = boost::make_shared<MatrixDouble>();
1103 jAc.resize(3, 3,
false);
1104 invJac.resize(3, 3,
false);
1109 CHKERR get_material_stiffness();
1122 auto get_ents = [&](
const std::string blockset_name) {
1125 if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
1134 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
1136 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
1138 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
1140 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
1144 Range boundary_ents;
1146 if (bc_mng->checkBlock(bc,
"FIX_"))
1147 boundary_ents.merge(*bc.second->getBcEntsPtr());
1148 boundary_ents.merge(get_ents(
"NO_CONTACT"));
1156 std::vector<DataFromBc>
bcData;
1158 mField,
"U",
simple->getProblemName(),
"DISPLACEMENT",
true);
1161 auto dm =
simple->getDM();
1189 mit->second->methodsOp.push_back(
1191 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1197 mit->second->methodsOp.push_back(
1199 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1205 mit->second->methodsOp.push_back(
1207 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1226 auto get_boundary_markers = [&](
string prob,
Range bound_ents,
int lo,
1228 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
1230 CHKERR problem_manager->modifyMarkDofs(
1231 prob,
ROW,
"U", lo, hi, ProblemsManager::MarkOP::OR, 1, *marker_ptr);
1234 bound_ents, *marker_ptr);
1239 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"MOVE_X",
"U",
1241 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"MOVE_Y",
"U",
1243 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"MOVE_Z",
"U",
1246 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ROTATE_X",
1252 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ROTATE_Z",
1254 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ROTATE_ALL",
1257 auto &bc_map = bc_mng->getBcMapByBlockName();
1259 bc_mng->getMergedBlocksMarker(vector<string>{
"MOVE_",
"ROTATE_"});
1261 auto add_disp_bc_ops = [&]() {
1263 std::map<char, size_t> xyz_map = {{
'X', 0}, {
'Y', 1}, {
'Z', 2}};
1270 if (bc_mng->checkBlock(bc,
"MOVE_")) {
1272 int idx = xyz_map.at(bc.first[(bc.first.find(
"MOVE_") + 5)]);
1274 auto disp_vec = bc.second->bcAttributes;
1275 if (disp_vec.empty())
1277 "provide an appropriate number of params (1) for MOVE block");
1279 disp(idx) = disp_vec[0];
1284 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1291 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1293 "U",
"U", [](
double,
double,
double) {
return 1.; },
1294 boost::make_shared<Range>(
bcData.back().ents)));
1300 if (std::regex_match(bc.first, std::regex(
"(.*)_ROTATE_(.*)"))) {
1302 <<
"Set boundary (rotation) residual for " << bc.first;
1303 auto angles = boost::make_shared<VectorDouble>(3);
1304 auto c_coords = boost::make_shared<VectorDouble>(3);
1305 if (bc.second->bcAttributes.size() < 6)
1307 "Wrong size of boundary attributes vector. Set the correct "
1308 "block size attributes (3) angles and (3) coordinates for "
1309 "center of rotation. Size of attributes %d",
1310 bc.second->bcAttributes.size());
1311 std::copy(&bc.second->bcAttributes[0], &bc.second->bcAttributes[3],
1312 angles->data().begin());
1313 std::copy(&bc.second->bcAttributes[3], &bc.second->bcAttributes[6],
1314 c_coords->data().begin());
1318 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1321 bc.second->getBcEntsPtr()));
1324 [](
double,
double,
double) { return 1.; },
1325 bc.second->getBcEntsPtr()));
1329 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1331 "U",
"U", [](
double,
double,
double) {
return 1.; },
1332 bc.second->getBcEntsPtr()));
1340 auto add_ho_ops = [&](
auto &pipeline) {
1342 auto jac_ptr = boost::make_shared<MatrixDouble>();
1343 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1344 auto det_ptr = boost::make_shared<VectorDouble>();
1346 "MESH_NODE_POSITIONS", jac_ptr));
1360 auto add_domain_base_ops = [&](
auto &pipeline) {
1365 CHKERR add_ho_ops(pipeline);
1400 auto add_domain_stress_ops = [&](
auto &pipeline,
auto mat_D_ptr) {
1415 auto add_skeleton_base_ops = [&](
auto &pipeline) {
1427 auto add_domain_ops_lhs = [&](
auto &pipeline,
auto m_D_ptr) {
1466 auto add_domain_ops_rhs = [&](
auto &pipeline,
auto m_D_ptr) {
1473 for (
int i = 0;
i != 3;
i++)
1474 t_source(
i) = (*cache).gravity[
i] * (*cache).density;
1475 const auto time = fe_domain_rhs->ts_t;
1476 t_source(
i) *= time;
1479 auto get_centrifugal_force = [
this](
const double x,
const double y,
1483 t_source(
i) = (*cache).density * (*cache).Omega(
i,
k) *
1484 (*cache).Omega(
k,
j) * t_coords(
j);
1500 pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
1518 auto add_domain_ops_rhs_constraints = [&](
auto &pipeline) {
1529 auto add_domain_ops_lhs_constraints = [&](
auto &pipeline,
auto m_D_ptr) {
1561 "TAU",
"TAU", [](
double,
double,
double) {
return 1.; },
1571 auto add_domain_ops_contact_lhs = [&](
auto &pipeline) {
1576 "SIGMA",
"U", []() {
return 1.; },
true));
1578 "SIGMA",
"U", []() {
return 1.; },
true));
1581 "SIGMA",
"SIGMA", [](
double,
double,
double) {
return 1.; },
1588 auto add_domain_ops_contact_rhs = [&](
auto &pipeline) {
1596 [](
double,
double,
double) {
return 1; }));
1609 auto add_boundary_base_ops = [&](
auto &pipeline) {
1632 auto add_boundary_ops_lhs = [&](
auto &pipeline) {
1645 "SIGMA",
"SIGMA", [](
double,
double,
double) {
return 1.; },
1656 auto vol_side_fe_ptr = boost::make_shared<MoFEM::DomainSideEle>(
mField);
1659 CHKERR add_ho_ops(vol_side_fe_ptr->getOpPtrVector());
1661 vol_side_fe_ptr->getOpPtrVector().push_back(
1673 auto add_boundary_ops_rhs = [&](
auto &pipeline) {
1686 auto my_vol_side_fe_ptr =
1687 boost::make_shared<MoFEM::DomainSideEle>(
mField);
1688 my_vol_side_fe_ptr->getOpPtrVector().push_back(
1692 my_vol_side_fe_ptr));
1711 auto &lam_pipeline_lhs =
feLambdaLhs->getOpPtrVector();
1712 auto &lam_pipeline_rhs =
feLambdaRhs->getOpPtrVector();
1717 CHKERR add_disp_bc_ops();
1730 CHKERR add_domain_ops_lhs_constraints(
1732 CHKERR add_domain_ops_rhs_constraints(
1735 CHKERR add_domain_base_ops(lam_pipeline_lhs);
1736 CHKERR add_domain_base_ops(lam_pipeline_rhs);
1737 CHKERR add_domain_stress_ops(lam_pipeline_lhs, tD_axi);
1738 CHKERR add_domain_stress_ops(lam_pipeline_rhs, tD_axi);
1739 CHKERR add_domain_ops_lhs(lam_pipeline_lhs, tD_axi);
1740 CHKERR add_domain_ops_rhs(lam_pipeline_rhs, tD_axi);
1756 CHKERR add_domain_ops_lhs_constraints(
1758 CHKERR add_domain_ops_rhs_constraints(
1774 auto integration_rule_bound = [](int, int,
int approx_order) {
1786 auto integration_rule_lambda = [](int, int,
int approx_order) {
1796 feLambdaLhs->getRuleHook = integration_rule_lambda;
1797 feLambdaRhs->getRuleHook = integration_rule_lambda;
1804 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1827 auto create_reactions_element = [&]() {
1831 auto get_reactions_meshset_ents = [&]() {
1840 ents, 1,
false, reaction_ents, moab::Interface::UNION);
1841 reaction_ents.merge(ents);
1847 int lents_size, ents_size;
1848 lents_size = reaction_ents.size();
1849 auto dm =
simple->getDM();
1850 CHKERR MPIU_Allreduce(&lents_size, &ents_size, 1, MPIU_INT, MPIU_SUM,
1851 PetscObjectComm((PetscObject)dm));
1854 "WORLD", Sev::warning,
1855 "Warning: Provided meshset (-reaction_id: %d) for calculating "
1856 "reactions is EMPTY! Calculating reactions for all FIX_ blocks.",
1859 if (bc_mng->checkBlock(bc,
"FIX_"))
1860 reaction_ents.merge(*bc.second->getBcEntsPtr());
1868 CHKERR get_reactions_meshset_ents();
1870 double defs[] = {0, 0, 0};
1873 MB_TAG_CREAT | MB_TAG_SPARSE, defs);
1878 std::vector<int> ghosts{0, 1, 2};
1889 CHKERR add_ho_ops(pipeline);
1930 CHKERR create_reactions_element();
1945 auto create_custom_ts = [&]() {
1946 auto set_dm_section = [&](
auto dm) {
1948 PetscSection section;
1950 simple->getProblemName(), §ion);
1951 CHKERR DMSetDefaultSection(dm, section);
1952 CHKERR DMSetDefaultGlobalSection(dm, section);
1953 CHKERR PetscSectionDestroy(§ion);
1956 auto dm =
simple->getDM();
1958 CHKERR set_dm_section(dm);
1959 boost::shared_ptr<FEMethod> null;
1960 preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
1963 preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
1993 auto postproc_bound_el = [&]() {
1996 auto &bmc = *bmc_ptr;
2007 pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
2009 if (pipeline_mng->getBoundaryLhsFE())
2011 pipeline_mng->getBoundaryLhsFE(), null,
2012 pipeline_mng->getBoundaryLhsFE());
2015 pipeline_mng->getDomainLhsFE(), null, null);
2020 if (pipeline_mng->getSkeletonLhsFE())
2022 pipeline_mng->getSkeletonLhsFE(), null,
2030 if (pipeline_mng->getDomainRhsFE()) {
2042 pipeline_mng->getDomainRhsFE(), null,
2045 if (pipeline_mng->getBoundaryRhsFE())
2047 pipeline_mng->getBoundaryRhsFE(), null,
2049 if (pipeline_mng->getSkeletonRhsFE())
2051 pipeline_mng->getSkeletonRhsFE(), null,
2057 &fit->second->getLoopFe(), NULL, NULL);
2063 &fit->second->getLoopFe(), NULL, NULL);
2070 &fit->second->getLoopFe(), NULL, NULL);
2081 auto print_active_set = [&](std::array<int, 2> &lgauss_pts,
string name,
2084 std::vector<int> gauss_pts(2, 0);
2086 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
2087 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2089 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\t \t Active %s gauss pts: %d / %d",
2090 name.c_str(), (
int)gauss_pts[0], (
int)gauss_pts[1]);
2091 lgauss_pts[0] = lgauss_pts[1] = 0;
2095 auto postproc_dom = [&]() {
2101 auto postproc_bound = [&]() {
2109 pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
2112 pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
2120 auto solver = create_custom_ts();
2122 CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
2124 auto dm =
simple->getDM();
2128 CHKERR TSSetFromOptions(solver);
2129 CHKERR TSSetSolution(solver,
D);
2136 CHKERR TSGetSNES(solver, &snes);
2138 CHKERR SNESGetKSP(snes, &ksp);
2140 CHKERR KSPGetPC(ksp, &pC);
2141 PetscBool is_pcfs = PETSC_FALSE;
2142 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
2145 auto make_is_field_map = [&]() {
2146 PetscSection section;
2147 CHKERR DMGetDefaultSection(dm, §ion);
2150 CHKERR PetscSectionGetNumFields(section, &num_fields);
2155 map<std::string, SmartPetscObj<IS>> is_map;
2156 for (
int ff = 0; ff != num_fields; ff++) {
2167 auto is_map = make_is_field_map();
2169 auto set_fieldsplit_on_bc = [&](PC &bc_pc,
string name_prb) {
2172 PetscBool is_bc_field_split;
2173 PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
2174 &is_bc_field_split);
2176 auto block_prefix =
simple->getProblemName();
2180 bc_mng->
getBlockIS(block_prefix,
"MOVE_Z",
"U", name_prb, 2, 2);
2181 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_Y",
"U", name_prb, 1,
2183 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_X",
"U", name_prb, 0,
2185 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_X",
"U", name_prb, 1,
2187 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_Z",
"U", name_prb, 0,
2189 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_ALL",
"U", name_prb,
2191 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_ALL",
"U", name_prb, 0,
2195 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
2197 if (is_bc_field_split && is_all_bc_size) {
2199 <<
"Fieldsplit preconditioner for boundary block size: "
2201 CHKERR PCFieldSplitSetIS(bc_pc, NULL,
2227 auto set_global_mat_and_vec = [&]() {
2231 CHKERR TSSetIFunction(solver, fe->ts_F, PETSC_NULL, PETSC_NULL);
2232 CHKERR TSSetIJacobian(solver, fe->ts_B, fe->ts_B, PETSC_NULL, PETSC_NULL);
2233 CHKERR KSPSetOperators(ksp, fe->ts_B, fe->ts_B);
2240 CHKERR set_global_mat_and_vec();
2242 PetscBool is_l0_field_split = PETSC_TRUE;
2243 PetscBool is_l1_field_split = PETSC_FALSE;
2247 auto set_fieldsplit_contact = [&]() {
2254 is_map[
"E_IS_SIG"] = sig_data->rowIs;
2255 sigma_ao = sig_data->getSmartRowMap();
2257 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"SIGMA"]);
2258 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"E_IS_SIG"]);
2269 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &ct_ksp);
2271 CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
2272 PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
2273 &is_l1_field_split);
2277 CHKERR PetscFree(ct_ksp);
2285 CHKERR set_fieldsplit_contact();
2287 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
2294 is_map[
"E_IS_EP"] = ep_data->rowIs;
2297 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"EP"]);
2300 CHKERR PCFieldSplitSetIS(pC,
"ep", is_map[
"EP"]);
2301 CHKERR PCFieldSplitSetIS(pC,
"etau", is_map[
"E_IS_EP"]);
2303 PCCompositeType pc_type;
2304 CHKERR PCFieldSplitGetType(pC, &pc_type);
2306 if (pc_type == PC_COMPOSITE_SCHUR) {
2313 boost::make_shared<BlockMatContainer>();
2317 CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
2325 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &tt_ksp);
2327 CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
2328 PetscBool is_tau_field_split;
2329 PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
2330 &is_tau_field_split);
2332 if (is_tau_field_split) {
2338 is_map[
"E_IS_TAU"] = tau_data->rowIs;
2341 CHKERR ep_data->getRowMap(&ep_ao);
2344 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"TAU"]);
2345 CHKERR AOApplicationToPetscIS(ep_ao, is_map[
"TAU"]);
2347 CHKERR PCFieldSplitSetIS(tau_pc,
"tau", is_map[
"TAU"]);
2348 CHKERR PCFieldSplitSetIS(tau_pc,
"elastic", is_map[
"E_IS_TAU"]);
2350 CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
2351 if (pc_type == PC_COMPOSITE_SCHUR) {
2355 CHKERR PCFieldSplitSetSchurPre(tau_pc,
2356 PC_FIELDSPLIT_SCHUR_PRE_USER,
2363 CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
2365 CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
2367 CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->
getName());
2368 CHKERR PetscFree(bc_ksp);
2371 CHKERR PetscFree(tt_ksp);
2376 char mumps_options[] =
"-fieldsplit_0_mat_mumps_icntl_14 1200 "
2377 "-fieldsplit_0_mat_mumps_icntl_24 1 "
2378 "-fieldsplit_0_mat_mumps_icntl_20 0 "
2379 "-fieldsplit_0_mat_mumps_icntl_13 1 "
2380 "-fieldsplit_1_mat_mumps_icntl_14 1200 "
2381 "-fieldsplit_1_mat_mumps_icntl_24 1 "
2382 "-fieldsplit_1_mat_mumps_icntl_20 0 "
2383 "-fieldsplit_1_mat_mumps_icntl_13 1";
2384 CHKERR PetscOptionsInsertString(NULL, mumps_options);
2385 CHKERR set_global_mat_and_vec();
2386 CHKERR set_fieldsplit_on_bc(pC,
simple->getProblemName());
2402 auto set_section_monitor = [&]() {
2405 CHKERR TSGetSNES(solver, &snes);
2406 PetscViewerAndFormat *vf;
2407 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2408 PETSC_VIEWER_DEFAULT, &vf);
2411 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
2416 auto create_post_process_element = [&]() {
2487 postProcFe->addFieldValuesPostProc(
"U",
"DISPLACEMENT");
2489 postProcFe->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
2493 auto scatter_create = [&](
auto coeff) {
2495 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
2496 ROW,
"U", coeff, coeff, is);
2498 CHKERR ISGetLocalSize(is, &loc_size);
2502 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
2507 auto set_time_monitor = [&]() {
2509 boost::shared_ptr<MMonitor> monitor_ptr(
2512 boost::shared_ptr<ForcesAndSourcesCore> null;
2514 monitor_ptr, null, null);
2518 CHKERR set_section_monitor();
2519 CHKERR create_post_process_element();
2523 CHKERR PetscPrintf(PETSC_COMM_WORLD,
2524 "Reading vector in binary from %s file...\n",
2527 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
2532 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2536 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
2538 Tokenizer tok{s, boost::char_separator<char>(
"_")};
2539 auto it = tok.begin();
2540 const int restart_step = std::stoi((++it)->c_str());
2542 string restart_tt = *(++it);
2543 restart_tt.erase(restart_tt.length() - 4);
2544 const double restart_time = std::atof(restart_tt.c_str());
2545 CHKERR TSSetStepNumber(solver, restart_step);
2546 CHKERR TSSetTime(solver, restart_time);
2552 CHKERR set_time_monitor();
2572 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2573 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2590 CHKERR skin.find_skin(0, body_ents,
false, skin_tris);
2591 Range boundary_ents;
2592 ParallelComm *pcomm =
2594 if (pcomm == NULL) {
2597 CHKERR pcomm->filter_pstatus(skin_tris, PSTATUS_SHARED | PSTATUS_MULTISHARED,
2598 PSTATUS_NOT, -1, &boundary_ents);
2600 return boundary_ents;
2613 DMType dm_name =
"DMMOFEM";
2616 moab::Core mb_instance;
2617 moab::Interface &moab = mb_instance;
2623 simple->getProblemName() =
"Multifield plasticity";
2629 PetscBool is_partitioned = PETSC_FALSE;
2631 &is_partitioned, PETSC_NULL);
2632 if (is_partitioned) {
2649 CHKERR mit.second.scaleParameters();
multi_index_container< boost::shared_ptr< Field >, indexed_by< hashed_unique< tag< BitFieldId_mi_tag >, const_mem_fun< Field, const BitFieldId &, &Field::getId >, HashBit< BitFieldId >, EqBit< BitFieldId > >, ordered_unique< tag< Meshset_mi_tag >, member< Field, EntityHandle, &Field::meshSet > >, ordered_unique< tag< FieldName_mi_tag >, const_mem_fun< Field, boost::string_ref, &Field::getNameRef > >, ordered_non_unique< tag< BitFieldId_space_mi_tag >, const_mem_fun< Field, FieldSpace, &Field::getSpace > > > > Field_multiIndex
Field_multiIndex for Field.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpSkeletonRhsPipeline()
Get the Op Skeleton Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
std::map< int, BlockParamData > mat_blocks
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixTensorTimesGradU< 3 > OpMiXLambdaGradURhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpCentrifugalForce2
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhsNoFS
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, 3, 3, 1 > OpLogStrainMatrixLhs
constexpr bool TEST_H1_SPACE
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, 3, 3, 0 > OpStiffnessMatrixLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
Tensors class implemented by Walter Landry.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_contact)
DEPRECATED auto smartCreateDMVector(DM dm)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
auto createTS(MPI_Comm comm)
DEPRECATED auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
MoFEMErrorCode set_plastic_tau_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic_tau)
MoFEMErrorCode set_plastic_ep_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic)
constexpr auto field_name
static constexpr int approx_order
MoFEMErrorCode getOptionsFromCommandLine()
chrono::high_resolution_clock::time_point stop
chrono::high_resolution_clock::time_point start
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Simple interface for fast problem set-up.
SmartPetscObj< IS > getBlockIS(const std::string block_prefix, const std::string block_name, const std::string field_name, const std::string problem_name, int lo, int hi, SmartPetscObj< IS > is_expand=SmartPetscObj< IS >())
Get block IS.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode remove_ents_from_finite_element(const std::string name, const EntityHandle meshset, const EntityType type, int verb=DEFAULT_VERBOSITY)=0
remove entities from given refinement level to finite element database
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
structure for User Loop Methods on finite elements
Provide data structure for (tensor) field approximation.
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Calculate HO coordinates at gauss points.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for symmetric tensorial field rank 2.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
const std::string getProblemName() const
Get the Problem Name.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.