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();
38 cout <<
"Time taken by function: " <<
duration.count() <<
" us." << endl;
42 using namespace MoFEM;
80 #include <BasicFiniteElements.hpp>
150 Range getEntsOnMeshSkin();
205 no_output = PETSC_FALSE;
206 debug_info = PETSC_FALSE;
209 scale_params = PETSC_FALSE;
211 is_large_strain = PETSC_FALSE;
212 is_fieldsplit_bc_only = PETSC_FALSE;
214 is_lambda_split = PETSC_FALSE;
216 is_neohooke = PETSC_FALSE;
217 with_plasticity = PETSC_FALSE;
219 calculate_reactions = PETSC_FALSE;
222 with_contact = PETSC_FALSE;
224 is_rotating = PETSC_FALSE;
228 is_ho_geometry = PETSC_FALSE;
230 save_restarts = PETSC_FALSE;
231 is_restart = PETSC_FALSE;
235 contact_history = move_history = dirichlet_history =
"-load_history";
236 moab_debug = &mb_post_debug;
246 boost::shared_ptr<OpContactTools::CommonData> common_data_ptr,
247 boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_proc_fe,
251 : dM(dm), mField(m_field), commonDataPtr(common_data_ptr),
252 postProcFe(post_proc_fe), uXScatter(ux_scatter),
253 uYScatter(uy_scatter), uZScatter(uz_scatter) {
254 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
255 volSideFe->getOpPtrVector().push_back(
257 volSideFe->getOpPtrVector().push_back(
260 boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
261 postProcFeSkeleton->generateReferenceElementMesh();
262 postProcFeSkeleton->getOpPtrVector().push_back(
265 postProcFeSkeleton->getOpPtrVector().push_back(
267 "SKELETON_LAMBDA", postProcFeSkeleton->postProcMesh,
268 postProcFeSkeleton->mapGaussPts, commonDataPtr));
270 postProcFeSkin = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
271 CHKERR postProcFeSkin->generateReferenceElementMesh();
272 CHKERR postProcFeSkin->addFieldValuesPostProc(
"U",
"DISPLACEMENT");
274 postProcFeSkin->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
277 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin(
"TAU");
278 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin(
"EP");
282 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin(
"SIGMA");
288 CHKERR TSGetTimeStep(ts, &(commonDataPtr->timeStepSize));
296 auto make_vtk = [&]() {
301 CHKERR postProcFe->writeFile(
"out_plastic_" +
302 boost::lexical_cast<std::string>(ts_step) +
307 auto save_skeleton = [&]() {
311 CHKERR postProcFeSkeleton->writeFile(
312 "out_skeleton_" + boost::lexical_cast<std::string>(ts_step) +
317 auto save_skin = [&]() {
323 CHKERR postProcFeSkin->writeFile(
324 "out_skin_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
328 double max_disp, min_disp;
330 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
332 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
333 INSERT_VALUES, SCATTER_FORWARD);
334 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
335 INSERT_VALUES, SCATTER_FORWARD);
336 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max_disp);
337 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min_disp);
338 PetscPrintf(PETSC_COMM_WORLD,
"%s time %3.4e min %3.4e max %3.4e\n",
339 msg.c_str(), ts_t, min_disp, max_disp);
346 auto &l_reactions = commonDataPtr->bodyReactions;
347 auto &lgauss_pts = commonDataPtr->stateArrayArea;
349 int roller_nb = (*cache).rollerDataVec.size();
350 l_reactions.resize(roller_nb * 3);
351 std::fill(l_reactions.begin(), l_reactions.end(), 0);
356 std::vector<double> gauss_pts(2, 0);
358 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2,
360 PetscObjectComm((PetscObject)dm));
362 "\t TS %d Total contact area (ref): %g / %g", ts_step,
363 gauss_pts[0], gauss_pts[1]);
364 gauss_area = gauss_pts[0];
365 lgauss_pts[0] = lgauss_pts[1] = 0;
367 vector<double> g_reactions(l_reactions.size() * 3, 0);
368 CHKERR MPIU_Allreduce(l_reactions.data(), g_reactions.data(),
369 roller_nb * 3, MPIU_REAL, MPIU_SUM,
370 PetscObjectComm((PetscObject)dm));
371 for (
int dd = 0;
dd != roller_nb; ++
dd)
373 "\t Body %d reactions Time %g Force X: %3.4e Y: "
375 (*cache).rollerDataVec[
dd].iD, ts_t,
376 g_reactions[3 *
dd + 0], g_reactions[3 *
dd + 1],
377 g_reactions[3 *
dd + 2]);
382 std::ostringstream ostrm, ostrm2;
383 ostrm <<
"out_debug_" << ts_step <<
".vtk";
384 ostrm2 <<
"out_debug_" << ts_step <<
".h5m";
390 "PARALLEL=WRITE_PART");
397 auto &vec = commonDataPtr->reactionsVec;
398 CHKERR VecZeroEntries(vec);
399 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
400 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
404 CHKERR mField.
get_moab().get_entities_by_dimension(0, 0, ents,
true);
406 ParallelComm *pcomm =
408 CHKERR pcomm->reduce_tags(commonDataPtr->reactionTag, MPI_SUM, ents);
414 double def_value = 0.;
415 CHKERR mField.
get_moab().tag_clear_data(commonDataPtr->reactionTag,
420 CHKERR VecAssemblyBegin(vec);
421 CHKERR VecAssemblyEnd(vec);
422 CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
423 CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
424 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
425 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
427 const double *react_ptr;
428 CHKERR VecGetArrayRead(vec, &react_ptr);
430 PetscPrintf(PETSC_COMM_WORLD,
431 "Reactions time %3.4e X: %3.4e Y: %3.4e Z: %3.4e \n", ts_t,
432 react_ptr[0], react_ptr[1], react_ptr[2]);
434 CHKERR VecRestoreArrayRead(vec, &react_ptr);
445 auto &m_field = postProcFe->mField;
447 if (m_field.get_comm_rank() == 0) {
448 std::ostringstream ost;
449 for (
auto &rol : (*cache).rollerDataVec) {
452 roller_coords += rol.originCoords;
455 if (!rol.positionDataParamName.empty() &&
457 auto scale_method = boost::make_shared<TimeAccelerogram>(
458 rol.positionDataParamName);
461 roller_coords = rol.BodyDispScaled;
468 ost <<
"out_roller_" << ts_step <<
".vtk";
476 CHKERR VecGhostUpdateBegin(ts_u, INSERT_VALUES, SCATTER_FORWARD);
477 CHKERR VecGhostUpdateEnd(ts_u, INSERT_VALUES, SCATTER_FORWARD);
480 "restart_" + to_string(ts_step) +
"_" + to_string(ts_t) +
".dat";
482 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, rest_name.c_str(),
483 FILE_MODE_WRITE, &viewer);
484 CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
485 VecView(ts_u, viewer);
486 PetscViewerDestroy(&viewer);
489 CHKERR print_max_min(uXScatter,
"Ux");
490 CHKERR print_max_min(uYScatter,
"Uy");
491 CHKERR print_max_min(uZScatter,
"Uz");
496 const double *react_ptr;
497 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
498 if (fabs(react_ptr[0] + 0.12519621572) > 1e-6)
500 "atom test diverged!");
501 CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
507 const double *react_ptr;
508 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
510 if ((0.03663003663 - react_ptr[2]) / 0.03663003663 > 1e-2)
512 "atom test diverged!");
513 CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
520 CHKERR VecMin(std::get<0>(uXScatter), PETSC_NULL, &min_disp);
521 if (fabs(min_disp + 2.6949) > 1e-2)
523 "atom test diverged!");
529 if (fabs(sqrt(gauss_area * 4.0 / M_PI) - 0.4574) / 0.4574 > 4e-2)
531 "atom test diverged!");
564 CHKERR createCommonData();
575 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"",
"none");
577 CHKERR PetscOptionsInt(
"-order",
"approximation order",
"", data.order,
578 &data.order, PETSC_NULL);
579 CHKERR PetscOptionsBool(
"-no_output",
"save no output files",
"",
580 data.no_output, &data.no_output, PETSC_NULL);
581 CHKERR PetscOptionsBool(
"-debug_info",
"print debug info",
"",
582 data.debug_info, &data.debug_info, PETSC_NULL);
583 CHKERR PetscOptionsBool(
"-scale_params",
584 "scale parameters (young modulus = 1)",
"",
585 data.scale_params, &data.scale_params, PETSC_NULL);
586 CHKERR PetscOptionsInt(
"-output_every",
587 "frequncy how often results are dumped on hard disk",
588 "", 1, &data.output_freq, PETSC_NULL);
589 CHKERR PetscOptionsBool(
"-calculate_reactions",
"calculate reactions",
"",
590 data.calculate_reactions, &data.calculate_reactions,
592 CHKERR PetscOptionsInt(
"-reaction_id",
"meshset id for computing reactions",
593 "", data.reaction_id, &data.reaction_id, PETSC_NULL);
594 CHKERR PetscOptionsInt(
"-atom_test",
"number of atom test",
"",
595 data.atom_test_nb, &data.atom_test_nb, PETSC_NULL);
596 CHKERR PetscOptionsBool(
"-is_large_strain",
"is large strains regime",
"",
597 data.is_large_strain, &data.is_large_strain,
600 "-is_fieldsplit_bc_only",
"use fieldsplit for boundary only",
"",
601 data.is_fieldsplit_bc_only, &data.is_fieldsplit_bc_only, PETSC_NULL);
604 "use axiator split approach, with reduced integration",
"",
605 data.is_lambda_split, &data.is_lambda_split, PETSC_NULL);
607 CHKERR PetscOptionsBool(
"-is_neohooke",
"is neohooke material",
"",
608 data.is_neohooke, &data.is_neohooke, PETSC_NULL);
614 "-with_plasticity",
"run calculations with plasticity",
"",
615 data.with_plasticity, &data.with_plasticity, PETSC_NULL);
617 CHKERR PetscOptionsBool(
"-is_rotating",
"is rotating frame used",
"",
618 data.is_rotating, &data.is_rotating, PETSC_NULL);
619 CHKERR PetscOptionsBool(
"-with_contact",
"run calculations with contact",
"",
620 data.with_contact, &data.with_contact, PETSC_NULL);
621 CHKERR PetscOptionsBool(
"-print_rollers",
622 "output file with rollers positions",
"",
623 data.print_rollers, &data.print_rollers, PETSC_NULL);
625 char load_hist_file[255];
626 PetscBool ctg_flag = PETSC_FALSE;
627 CHKERR PetscOptionsString(
"-contact_history",
"load_history.in",
"",
628 "load_history.in", load_hist_file, 255, &ctg_flag);
630 data.contact_history =
"-contact_history";
631 CHKERR PetscOptionsString(
"-move_history",
"load_history.in",
"",
632 "load_history.in", load_hist_file, 255, &ctg_flag);
635 data.move_history =
"-move_history";
637 CHKERR PetscOptionsString(
"-dirichlet_history",
"load_history.in",
"",
638 "load_history.in", load_hist_file, 255, &ctg_flag);
640 data.dirichlet_history =
"-dirichlet_history";
643 CHKERR PetscOptionsString(
"-file_name",
"file name",
"",
"mesh.h5m",
647 char restart_file_name[255];
648 CHKERR PetscOptionsString(
"-restart_file",
"restart file name",
"",
649 "restart.dat", restart_file_name, 255,
651 data.restart_file_str = string(restart_file_name);
653 CHKERR PetscOptionsBool(
"-save_restarts",
654 "save restart files at each iteration",
"",
655 data.save_restarts, &data.save_restarts, PETSC_NULL);
657 ierr = PetscOptionsEnd();
660 if (data.is_large_strain && data.with_plasticity && data.is_neohooke)
662 "Neohooke material is not supported with plasticity.");
672 CHKERR moab.get_entities_by_dimension(0, 3, tets,
false);
674 CHKERR moab.get_adjacencies(tets, 2,
true, faces, moab::Interface::UNION);
676 CHKERR moab.get_adjacencies(tets, 1,
true, edges, moab::Interface::UNION);
682 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
684 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
685 CHKERR moab.create_meshset(MESHSET_SET, part_set);
686 Tag part_tag = pcomm->part_tag();
691 CHKERR moab.get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag, NULL, 1,
692 tagged_sets, moab::Interface::UNION);
694 for (
auto &mit : tagged_sets) {
696 CHKERR moab.tag_get_data(part_tag, &mit, 1, &part);
698 CHKERR moab.get_entities_by_dimension(mit, 3, proc_ents,
true);
699 CHKERR moab.get_entities_by_handle(mit, all_proc_ents,
true);
701 CHKERR moab.get_entities_by_handle(mit, off_proc_ents,
true);
705 Range proc_ents_skin[4];
706 proc_ents_skin[3] = proc_ents;
708 Range all_tets, all_skin;
709 CHKERR m_field.
get_moab().get_entities_by_dimension(0, 3, all_tets,
false);
710 CHKERR skin.find_skin(0, all_tets,
false, all_skin);
711 CHKERR skin.find_skin(0, proc_ents,
false, proc_ents_skin[2]);
713 proc_ents_skin[2] = subtract(proc_ents_skin[2], all_skin);
714 CHKERR moab.get_adjacencies(proc_ents_skin[2], 1,
false, proc_ents_skin[1],
715 moab::Interface::UNION);
716 CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0],
true);
717 for (
int dd = 0;
dd != 3;
dd++)
718 CHKERR moab.add_entities(part_set, proc_ents_skin[
dd]);
719 CHKERR moab.add_entities(part_set, all_proc_ents);
724 CHKERR moab.get_entities_by_handle(0, all_ents);
725 std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
726 CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
727 &*pstat_tag.begin());
731 CHKERR moab.get_entities_by_handle(part_set, ents_to_keep,
false);
732 off_proc_ents = subtract(off_proc_ents, ents_to_keep);
734 CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets,
false);
735 for (
auto m : meshsets) {
736 CHKERR moab.remove_entities(
m, off_proc_ents);
739 CHKERR moab.delete_entities(off_proc_ents);
741 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
748 CHKERR m_field.
get_moab().get_entities_by_dimension(part_set, 3, ents,
true);
750 ents,
simple->getBitRefLevel(),
false);
751 simple->getMeshSet() = part_set;
758 string block_name =
"MAT_ELASTIC";
760 if (!data.with_plasticity)
762 if (it->getName().compare(0, block_name.length(), block_name) == 0) {
763 std::vector<double> block_data;
764 CHKERR it->getAttributes(block_data);
765 const int id = it->getMeshsetId();
770 mField.get_moab().get_entities_by_dimension(meshset, 3, tets,
true);
772 for (
auto ent : tets)
773 data.mapBlockTets[ent] = id;
776 mat_blocks.at(
id).young_modulus = block_data[0];
778 if (block_data.size() > 2)
780 mat_blocks.at(
id).cn_cont = 1. / block_data[0];
782 if (block_data.size() < 2)
784 "provide an appropriate number of entries (2) parameters for "
785 "MAT_ELASTIC block");
804 enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
805 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz",
"bernstein"};
806 PetscInt choice_base_value = AINSWORTH;
808 LASBASETOPT, &choice_base_value, PETSC_NULL);
811 switch (choice_base_value) {
815 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
820 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
825 <<
"Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
832 auto skin_ents = getEntsOnMeshSkin();
833 data.skinEnts = skin_ents;
836 int nb_tets, nb_hexes;
837 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
838 CHKERR mField.get_moab().get_connectivity(&*tets.begin(), 1, nodes);
839 CHKERR mField.get_moab().get_number_entities_by_type(0, MBTET, nb_tets,
true);
840 CHKERR mField.get_moab().get_number_entities_by_type(0, MBHEX, nb_hexes,
843 if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {
844 data.is_ho_geometry = PETSC_TRUE;
858 if (data.with_contact) {
862 CHKERR simple->setFieldOrder(
"SIGMA", data.order - 1, &skin_ents);
864 if (data.with_plasticity) {
866 CHKERR simple->setFieldOrder(
"TAU", std::max(0, data.order - 1));
868 CHKERR simple->setFieldOrder(
"EP", std::max(0, data.order - 1));
871 if (data.is_rotating && data.with_plasticity &&
false) {
874 std::max(0, data.order - 1));
884 simple->getOtherFiniteElements().push_back(
"FORCE_FE");
885 simple->getOtherFiniteElements().push_back(
"PRESSURE_FE");
886 if (data.is_rotating && data.with_plasticity &&
false)
893 if (data.is_ho_geometry) {
895 "MESH_NODE_POSITIONS");
896 CHKERR mField.loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material);
900 if (data.is_rotating && data.with_plasticity &&
false) {
901 CHKERR mField.remove_ents_from_finite_element(
"sFE", skin_ents);
902 CHKERR mField.build_finite_elements(
"sFE");
912 auto get_material_stiffness = [&]() {
919 (*cache).Is(
i,
j,
k,
l) =
930 const double G = (*cache).
young_modulus / (2. * (1. + (*cache).poisson));
932 (*cache).young_modulus / (3. * (1. - 2. * (*cache).poisson));
935 double A = 6. * G / (3. *
K + 4. * G);
939 auto t_D_tmp = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
941 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Axiator);
943 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
945 t_D_axiator(
i,
j,
k,
l) =
946 A * (
K - (2. / 3.) * G) * t_kds(
i,
j) * t_kds(
k,
l);
947 t_D_deviator(
i,
j,
k,
l) = 2 * G * ((t_kds(
i,
k) ^ t_kds(
j,
l)) / 4.);
949 t_D_tmp(
i,
j,
k,
l) = t_D_deviator(
i,
j,
k,
l) + t_D_axiator(
i,
j,
k,
l);
951 (*cache).mtD = *commonDataPtr->mtD;
952 (*cache).scale_constraint = 1.;
960 commonDataPtr = boost::make_shared<OpContactTools::CommonData>();
962 commonDataPtr->mtD = boost::make_shared<MatrixDouble>();
963 commonDataPtr->mtD->resize(36, 1);
964 commonDataPtr->mtD_Axiator = boost::make_shared<MatrixDouble>();
965 commonDataPtr->mtD_Axiator->resize(36, 1);
966 commonDataPtr->mtD_Deviator = boost::make_shared<MatrixDouble>();
967 commonDataPtr->mtD_Deviator->resize(36, 1);
969 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
970 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
971 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
972 commonDataPtr->mPiolaStressPtr = boost::make_shared<MatrixDouble>();
973 commonDataPtr->mPiolaStressPtr->resize(9, 1,
false);
975 commonDataPtr->mDeformationGradient = boost::make_shared<MatrixDouble>();
976 commonDataPtr->matC = boost::make_shared<MatrixDouble>();
977 commonDataPtr->matEigenVal = boost::make_shared<MatrixDouble>();
978 commonDataPtr->matEigenVector = boost::make_shared<MatrixDouble>();
979 commonDataPtr->detF = boost::make_shared<VectorDouble>();
980 commonDataPtr->materialTangent = boost::make_shared<MatrixDouble>();
981 commonDataPtr->dE_dF_mat = boost::make_shared<MatrixDouble>();
982 commonDataPtr->dE_dF_mat->resize(81, 1,
false);
984 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
985 commonDataPtr->contactStressDivergencePtr =
986 boost::make_shared<MatrixDouble>();
987 commonDataPtr->contactNormalPtr = boost::make_shared<MatrixDouble>();
988 commonDataPtr->contactGapPtr = boost::make_shared<VectorDouble>();
989 commonDataPtr->contactGapDiffPtr = boost::make_shared<MatrixDouble>();
990 commonDataPtr->contactNormalDiffPtr = boost::make_shared<MatrixDouble>();
992 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
993 commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
995 commonDataPtr->plasticSurfacePtr = boost::make_shared<VectorDouble>();
996 commonDataPtr->plasticFlowPtr = boost::make_shared<MatrixDouble>();
997 commonDataPtr->plasticTauPtr = boost::make_shared<VectorDouble>();
998 commonDataPtr->plasticStrainPtr = boost::make_shared<MatrixDouble>();
1000 commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
1001 commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
1003 commonDataPtr->plasticTauJumpPtr = boost::make_shared<VectorDouble>();
1004 commonDataPtr->plasticStrainJumpPtr = boost::make_shared<MatrixDouble>();
1005 commonDataPtr->guidingVelocityPtr = boost::make_shared<MatrixDouble>();
1006 commonDataPtr->guidingVelocityPtr->resize(3, 1,
false);
1008 commonDataPtr->plasticGradTauPtr = boost::make_shared<MatrixDouble>();
1009 commonDataPtr->plasticGradStrainPtr = boost::make_shared<MatrixDouble>();
1011 jAc.resize(3, 3,
false);
1012 invJac.resize(3, 3,
false);
1014 if (data.with_plasticity)
1015 CHKERR commonDataPtr->calculateDiffDeviator();
1017 CHKERR get_material_stiffness();
1018 if (data.is_neohooke)
1019 commonDataPtr->isNeoHook =
true;
1026 auto bc_mng = mField.getInterface<
BcManager>();
1030 auto fix_disp = [&](
const std::string blockset_name) {
1033 if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
1035 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
1042 auto remove_ents = [&](
string prob,
const Range &&ents,
const int low,
1048 CHKERR mField.get_moab().get_connectivity(ents, verts,
true);
1049 CHKERR mField.get_moab().get_adjacencies(ents, 1,
false, verts,
1050 moab::Interface::UNION);
1053 CHKERR prb_mng->removeDofsOnEntities(prob,
"U", verts, low, high);
1057 if (data.with_contact) {
1058 Range boundary_ents;
1059 boundary_ents.merge(fix_disp(
"FIX_X"));
1060 boundary_ents.merge(fix_disp(
"FIX_Y"));
1061 boundary_ents.merge(fix_disp(
"FIX_Z"));
1062 boundary_ents.merge(fix_disp(
"FIX_ALL"));
1063 boundary_ents.merge(fix_disp(
"NO_CONTACT"));
1071 CHKERR remove_ents(
simple->getProblemName(), fix_disp(
"FIX_X"), 0, 0);
1072 CHKERR remove_ents(
simple->getProblemName(), fix_disp(
"FIX_Y"), 1, 1);
1073 CHKERR remove_ents(
simple->getProblemName(), fix_disp(
"FIX_Z"), 2, 2);
1074 CHKERR remove_ents(
simple->getProblemName(), fix_disp(
"FIX_ALL"), 0, 2);
1076 std::vector<DataFromBc> bcData;
1077 dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementRemoveDofsBc>(
1078 mField,
"U",
simple->getProblemName(),
"DISPLACEMENT",
true);
1079 dirichlet_bc_ptr->iNitialize();
1082 auto dm =
simple->getDM();
1083 if (data.with_contact && data.with_plasticity) {
1087 }
else if (data.with_plasticity) {
1102 for (
auto mit = neumann_forces.begin(); mit != neumann_forces.end(); mit++) {
1103 if (data.is_ho_geometry)
1106 mit->second->methodsOp.push_back(
1108 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1110 for (
auto mit = nodal_forces.begin(); mit != nodal_forces.end(); mit++) {
1111 if (data.is_ho_geometry)
1114 mit->second->methodsOp.push_back(
1116 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1118 for (
auto mit = edge_forces.begin(); mit != edge_forces.end(); mit++) {
1119 if (data.is_ho_geometry)
1122 mit->second->methodsOp.push_back(
1124 mit->second->methodsOp.push_back(
new LoadScale((*cache).young_modulus_inv));
1128 dirichlet_bc_ptr->methodsOp.push_back(
1140 CHKERR mField.get_fields(&fields_ptr);
1142 auto get_boundary_markers = [&](
string prob, Range bound_ents,
int lo,
1144 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
1146 CHKERR problem_manager->modifyMarkDofs(
1147 prob,
ROW,
"U", lo, hi, ProblemsManager::MarkOP::OR, 1, *marker_ptr);
1150 bound_ents, *marker_ptr);
1155 typedef boost::shared_ptr<std::vector<unsigned char>> MarkerPtr;
1157 std::vector<MarkerPtr> markers_vec;
1158 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
1159 allBoundaryEnts = boost::make_shared<Range>();
1161 auto get_adj = [&](Range &ents) {
1163 CHKERR mField.get_moab().get_connectivity(ents, verts,
true);
1164 CHKERR mField.get_moab().get_adjacencies(ents, 1,
false, verts,
1165 moab::Interface::UNION);
1172 auto add_disp_bc_ops = [&]() {
1175 if (data.with_plasticity)
1179 std::map<char, size_t> xyz_map = {{
'X', 0}, {
'Y', 1}, {
'Z', 2}};
1183 for (
string move_s : {
"MOVE_X",
"MOVE_Y",
"MOVE_Z"}) {
1184 if (it->getName().compare(0, move_s.length(), move_s) == 0) {
1186 xyz_map.at(it->getName()[(it->getName().find(
"MOVE_") + 5)]);
1188 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
1190 std::vector<double> disp_vec;
1191 it->getAttributes(disp_vec);
1192 if (disp_vec.empty())
1195 "provide an appropriate number of params (1) for MOVE block");
1198 disp(ii) = disp_vec[0];
1199 auto verts = get_adj(ents);
1200 allBoundaryEnts->merge(verts);
1202 get_boundary_markers(
simple->getProblemName(), verts, ii, ii);
1203 markers_vec.push_back(markers);
1208 new OpSetBc(
"U",
false, markers));
1211 bcData.back().scaledValues,
1212 bcData.back().ents));
1217 new OpSetBc(
"U",
false, markers));
1220 "U",
"U", [](
double,
double,
double) {
return 1.; },
1221 boost::make_shared<Range>(bcData.back().ents)));
1229 auto merge_boundary_markers = [&](std::vector<MarkerPtr> &mark) {
1230 auto all_markers = boost::make_shared<std::vector<char unsigned>>();
1231 for (
auto &mit : mark) {
1232 all_markers->resize(mit->size(), 0);
1233 for (
int i = 0;
i != mit->size(); ++
i)
1234 (*all_markers)[
i] |= (*mit)[
i];
1239 boundaryMarker = merge_boundary_markers(markers_vec);
1241 auto bc_mng = mField.getInterface<
BcManager>();
1242 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ROTATE_X",
1248 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ROTATE_Z",
1250 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
1251 "ROTATE_ALL",
"U", 0, 3);
1253 auto &bc_map = bc_mng->getBcMapByBlockName();
1254 if (bc_map.size()) {
1255 for (
auto b : bc_map) {
1256 if (std::regex_match(
b.first, std::regex(
"(.*)_ROTATE_(.*)"))) {
1257 boundaryMarker->resize(
b.second->bcMarkers.size(), 0);
1258 for (
int i = 0;
i !=
b.second->bcMarkers.size(); ++
i) {
1261 allBoundaryEnts->merge(
b.second->bcEnts);
1262 (*boundaryMarker)[
i] |=
b.second->bcMarkers[
i];
1268 if (boundaryMarker->empty() && bc_map.size() == 0)
1270 get_boundary_markers(
simple->getProblemName(), Range(), -1, -1);
1273 if (std::regex_match(bc.first, std::regex(
"(.*)_ROTATE_(.*)"))) {
1275 <<
"Set boundary (rotation) residual for " << bc.first;
1276 auto angles = boost::make_shared<VectorDouble>(3);
1277 auto c_coords = boost::make_shared<VectorDouble>(3);
1278 if (bc.second->bcAttributes.size() < 6)
1280 "Wrong size of boundary attributes vector. Set the correct "
1281 "block size attributes (3) angles and (3) coordinates for "
1282 "center of rotation. Size of attributes %d",
1283 bc.second->bcAttributes.size());
1284 std::copy(&bc.second->bcAttributes[0], &bc.second->bcAttributes[3],
1285 angles->data().begin());
1286 std::copy(&bc.second->bcAttributes[3], &bc.second->bcAttributes[6],
1287 c_coords->data().begin());
1288 bcData.push_back(
new DataFromMove(*angles, *bc.second->getBcEntsPtr()));
1291 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1293 new OpRotate(
"U", bcData.back().scaledValues, c_coords,
1294 bc.second->getBcEntsPtr()));
1296 "U", commonDataPtr->mDispPtr,
1297 [](
double,
double,
double) { return 1.; },
1298 bc.second->getBcEntsPtr()));
1302 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
1304 "U",
"U", [](
double,
double,
double) {
return 1.; },
1305 bc.second->getBcEntsPtr()));
1313 auto add_ho_ops = [&](
auto &pipeline) {
1315 auto jac_ptr = boost::make_shared<MatrixDouble>();
1316 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1317 auto det_ptr = boost::make_shared<VectorDouble>();
1319 "MESH_NODE_POSITIONS", jac_ptr));
1333 auto add_domain_base_ops = [&](
auto &pipeline) {
1337 if (data.is_ho_geometry)
1338 CHKERR add_ho_ops(pipeline);
1345 if (data.with_plasticity) {
1347 "EP", commonDataPtr->plasticStrainPtr));
1349 "TAU", commonDataPtr->plasticTauPtr));
1351 "TAU", commonDataPtr->plasticTauDotPtr));
1353 "EP", commonDataPtr->plasticStrainDotPtr));
1355 if (data.is_rotating) {
1357 "EP", commonDataPtr->plasticGradStrainPtr));
1359 "TAU", commonDataPtr->plasticGradTauPtr));
1365 if (data.is_large_strain) {
1366 pipeline.push_back(
new OpLogStrain(
"U", commonDataPtr));
1368 pipeline.push_back(
new OpStrain(
"U", commonDataPtr));
1373 auto add_domain_stress_ops = [&](
auto &pipeline,
auto mat_D_ptr) {
1376 if (data.with_plasticity) {
1377 pipeline.push_back(
new OpPlasticStress(
"U", commonDataPtr, mat_D_ptr));
1383 pipeline.push_back(
new OpStress(
"U", commonDataPtr, mat_D_ptr));
1388 auto add_skeleton_base_ops = [&](
auto &pipeline) {
1389 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
1391 volSideFe->getOpPtrVector().push_back(
1393 volSideFe->getOpPtrVector().push_back(
1400 auto add_domain_ops_lhs = [&](
auto &pipeline,
auto m_D_ptr) {
1403 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1404 if (data.with_plasticity)
1405 pipeline.push_back(
new OpSchurBegin(
"U", commonDataPtr));
1407 if (data.is_large_strain) {
1408 if (data.is_neohooke)
1418 if (data.with_plasticity) {
1419 if (data.is_large_strain)
1421 "U",
"EP", commonDataPtr, m_D_ptr));
1424 "U",
"EP", commonDataPtr, m_D_ptr));
1431 if (data.is_rotating)
1439 auto add_domain_ops_rhs = [&](
auto &pipeline,
auto m_D_ptr) {
1442 auto get_body_force = [
this](
const double,
const double,
const double) {
1446 for (
int i = 0;
i != 3;
i++)
1447 t_source(
i) = (*cache).gravity[
i] * (*cache).density;
1448 const auto time = fe_domain_rhs->ts_t;
1449 t_source(
i) *= time;
1452 auto get_centrifugal_force = [
this](
const double x,
const double y,
1456 t_source(
i) = (*cache).density * (*cache).Omega(
i,
k) *
1457 (*cache).Omega(
k,
j) * t_coords(
j);
1465 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1467 if (data.is_rotating) {
1473 pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
1474 if (data.is_large_strain) {
1475 if (data.is_neohooke)
1480 "U", commonDataPtr, m_D_ptr));
1491 auto add_domain_ops_rhs_constraints = [&](
auto &pipeline) {
1493 if (data.with_plasticity) {
1496 if (data.debug_info)
1502 auto add_domain_ops_lhs_constraints = [&](
auto &pipeline,
auto m_D_ptr) {
1505 if (data.with_plasticity) {
1507 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1509 if (data.is_large_strain) {
1511 "U", commonDataPtr, m_D_ptr));
1532 if (data.with_plasticity) {
1534 "TAU",
"TAU", [](
double,
double,
double) {
return 1.; },
1536 pipeline.push_back(
new OpSchurEnd(
"U", commonDataPtr));
1544 auto add_domain_ops_contact_lhs = [&](
auto &pipeline) {
1545 if (data.with_contact) {
1546 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1549 "SIGMA",
"U", []() {
return 1.; },
true));
1551 "SIGMA",
"U", []() {
return 1.; },
true));
1554 "SIGMA",
"SIGMA", [](
double,
double,
double) {
return 1.; },
1561 auto add_domain_ops_contact_rhs = [&](
auto &pipeline) {
1562 if (data.with_contact) {
1564 "SIGMA", commonDataPtr->contactStressPtr));
1566 "SIGMA", commonDataPtr->contactStressDivergencePtr));
1569 [](
double,
double,
double) {
return 1; }));
1573 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1575 "U", commonDataPtr->contactStressDivergencePtr));
1582 auto add_boundary_base_ops = [&](
auto &pipeline) {
1586 if (data.is_ho_geometry) {
1596 if (data.with_contact) {
1600 "SIGMA", commonDataPtr->contactTractionPtr));
1605 auto add_boundary_ops_lhs = [&](
auto &pipeline) {
1606 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1611 if (data.with_contact) {
1618 "SIGMA",
"SIGMA", [](
double,
double,
double) {
return 1.; },
1629 if (data.is_rotating) {
1630 auto vol_side_fe_ptr = boost::make_shared<MoFEM::DomainSideEle>(mField);
1632 if (data.is_ho_geometry)
1633 CHKERR add_ho_ops(vol_side_fe_ptr->getOpPtrVector());
1635 vol_side_fe_ptr->getOpPtrVector().push_back(
1641 if (data.with_plasticity)
1647 auto add_boundary_ops_rhs = [&](
auto &pipeline) {
1648 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
1649 if (data.debug_info && data.with_contact)
1652 if (data.with_contact)
1659 if (data.is_rotating) {
1660 auto my_vol_side_fe_ptr =
1661 boost::make_shared<MoFEM::DomainSideEle>(mField);
1662 my_vol_side_fe_ptr->getOpPtrVector().push_back(
1664 commonDataPtr->mGradPtr));
1666 my_vol_side_fe_ptr));
1671 auto &tD = commonDataPtr->mtD;
1672 auto &tD_axi = commonDataPtr->mtD_Axiator;
1673 auto &tD_dev = commonDataPtr->mtD_Deviator;
1675 feLambdaLhs = boost::make_shared<DomainEle>(mField);
1676 feLambdaRhs = boost::make_shared<DomainEle>(mField);
1677 auto &lam_pipeline_lhs = feLambdaLhs->getOpPtrVector();
1678 auto &lam_pipeline_rhs = feLambdaRhs->getOpPtrVector();
1683 CHKERR add_disp_bc_ops();
1686 if (data.is_lambda_split) {
1696 CHKERR add_domain_ops_lhs_constraints(
1698 CHKERR add_domain_ops_rhs_constraints(
1701 CHKERR add_domain_base_ops(lam_pipeline_lhs);
1702 CHKERR add_domain_base_ops(lam_pipeline_rhs);
1703 CHKERR add_domain_stress_ops(lam_pipeline_lhs, tD_axi);
1704 CHKERR add_domain_stress_ops(lam_pipeline_rhs, tD_axi);
1705 CHKERR add_domain_ops_lhs(lam_pipeline_lhs, tD_axi);
1706 CHKERR add_domain_ops_rhs(lam_pipeline_rhs, tD_axi);
1722 CHKERR add_domain_ops_lhs_constraints(
1724 CHKERR add_domain_ops_rhs_constraints(
1732 if (data.is_rotating && data.with_plasticity &&
false)
1757 feLambdaLhs->getRuleHook = integration_rule_lambda;
1758 feLambdaRhs->getRuleHook = integration_rule_lambda;
1760 if (data.with_contact) {
1763 CHKERR mField.get_moab().tag_get_handle(
1764 "_ROLLER_NB", 1, MB_TYPE_INTEGER, commonDataPtr->rollerTag,
1765 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1767 data.update_roller = boost::make_shared<BoundaryEle>(mField);
1768 auto &pipeline = data.update_roller->getOpPtrVector();
1769 if (data.is_ho_geometry) {
1778 "SIGMA", commonDataPtr->contactTractionPtr));
1781 if (data.debug_info)
1785 data.update_roller->getRuleHook = integration_rule_bound;
1788 auto create_reactions_element = [&]() {
1790 Range &reaction_ents = commonDataPtr->reactionEnts;
1792 auto get_reactions_meshset_ents = [&]() {
1795 if (it->getMeshsetId() == data.reaction_id) {
1797 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
1799 CHKERR mField.get_moab().get_connectivity(ents, reaction_ents,
true);
1800 CHKERR mField.get_moab().get_adjacencies(
1801 ents, 1,
false, reaction_ents, moab::Interface::UNION);
1802 reaction_ents.merge(ents);
1808 int lents_size, ents_size;
1809 lents_size = reaction_ents.size();
1810 auto dm =
simple->getDM();
1811 CHKERR MPIU_Allreduce(&lents_size, &ents_size, 1, MPIU_INT, MPIU_SUM,
1812 PetscObjectComm((PetscObject)dm));
1813 if (!ents_size && data.reaction_id != -1)
1815 "Warning: Provided meshset (id: %d) for calculating "
1816 "reactions is EMPTY!",
1826 CHKERR get_reactions_meshset_ents();
1828 double defs[] = {0, 0, 0};
1829 CHKERR mField.get_moab().tag_get_handle(
"_REACTION_TAG", 3, MB_TYPE_DOUBLE,
1830 commonDataPtr->reactionTag,
1831 MB_TAG_CREAT | MB_TAG_SPARSE, defs);
1834 data.calc_reactions = boost::make_shared<DomainEle>(mField);
1836 std::vector<int> ghosts{0, 1, 2};
1838 mField.get_comm(), (mField.get_comm_rank() ? 0 : 3), 3,
1839 (mField.get_comm_rank() ? 3 : 0), &*ghosts.begin());
1843 auto &pipeline = data.calc_reactions->getOpPtrVector();
1846 if (data.is_ho_geometry)
1847 CHKERR add_ho_ops(pipeline);
1851 if (data.is_large_strain) {
1852 pipeline.push_back(
new OpLogStrain(
"U", commonDataPtr));
1854 pipeline.push_back(
new OpStrain(
"U", commonDataPtr));
1856 if (data.with_plasticity) {
1858 "EP", commonDataPtr->plasticStrainPtr));
1862 pipeline.push_back(
new OpStress(
"U", commonDataPtr, commonDataPtr->mtD));
1871 if (data.is_large_strain) {
1872 if (data.is_neohooke)
1877 "U", commonDataPtr, commonDataPtr->mtD));
1887 if (data.calculate_reactions)
1888 CHKERR create_reactions_element();
1900 auto preProc = boost::make_shared<FePrePostProcess>(
1901 mField, (*cache).rollerDataVec, bcData, commonDataPtr);
1903 auto create_custom_ts = [&]() {
1904 auto set_dm_section = [&](
auto dm) {
1906 PetscSection section;
1908 simple->getProblemName(), §ion);
1909 CHKERR DMSetDefaultSection(dm, section);
1910 CHKERR DMSetDefaultGlobalSection(dm, section);
1911 CHKERR PetscSectionDestroy(§ion);
1914 auto dm =
simple->getDM();
1916 CHKERR set_dm_section(dm);
1917 boost::shared_ptr<FEMethod>
null;
1918 preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
1921 preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
1924 for (
auto &roller : (*cache).rollerDataVec) {
1925 if (!roller.positionDataParamName.empty())
1926 preProc->methodsOpForRollersPosition.push_back(
1927 boost::shared_ptr<MethodForForceScaling>(
1943 dirichlet_bc_ptr,
null);
1947 auto postproc_bound_el = [&]() {
1950 auto &bmc = *bmc_ptr;
1960 if (data.with_plasticity)
1961 pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
1963 if (pipeline_mng->getBoundaryLhsFE())
1965 pipeline_mng->getBoundaryLhsFE(),
null,
1966 pipeline_mng->getBoundaryLhsFE());
1969 pipeline_mng->getDomainLhsFE(),
null,
null);
1970 if (data.is_lambda_split)
1974 if (pipeline_mng->getSkeletonLhsFE())
1976 pipeline_mng->getSkeletonLhsFE(),
null,
1984 if (pipeline_mng->getDomainRhsFE()) {
1989 if (data.is_lambda_split)
1991 feLambdaRhs,
null,
null);
1993 dirichlet_bc_ptr,
null);
1996 pipeline_mng->getDomainRhsFE(),
null,
1999 if (pipeline_mng->getBoundaryRhsFE())
2001 pipeline_mng->getBoundaryRhsFE(),
null,
2003 if (pipeline_mng->getSkeletonRhsFE())
2005 pipeline_mng->getSkeletonRhsFE(),
null,
2008 for (
auto fit = neumann_forces.begin(); fit != neumann_forces.end();
2011 &fit->second->getLoopFe(), NULL, NULL);
2015 for (
auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
2017 &fit->second->getLoopFe(), NULL, NULL);
2021 for (
auto fit = nodal_forces.begin(); fit != nodal_forces.end();
2024 &fit->second->getLoopFe(), NULL, NULL);
2035 auto print_active_set = [&](std::array<int, 2> &lgauss_pts,
string name,
2038 std::vector<int> gauss_pts(2, 0);
2039 auto dm = mField.getInterface<
Simple>()->getDM();
2040 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
2041 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2043 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\t \t Active %s gauss pts: %d / %d",
2044 name.c_str(), (
int)gauss_pts[0], (
int)gauss_pts[1]);
2045 lgauss_pts[0] = lgauss_pts[1] = 0;
2049 auto postproc_dom = [&]() {
2051 CHKERR print_active_set(commonDataPtr->stateArrayPlast,
"plastic",
2055 auto postproc_bound = [&]() {
2057 CHKERR print_active_set(commonDataPtr->stateArrayCont,
"contact", mField);
2061 if (data.debug_info) {
2062 if (data.with_plasticity)
2063 pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
2065 if (data.with_contact)
2066 pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
2074 auto solver = create_custom_ts();
2076 CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
2078 auto dm =
simple->getDM();
2082 CHKERR TSSetFromOptions(solver);
2083 CHKERR TSSetSolution(solver,
D);
2090 CHKERR TSGetSNES(solver, &snes);
2092 CHKERR SNESGetKSP(snes, &ksp);
2094 CHKERR KSPGetPC(ksp, &pC);
2095 PetscBool is_pcfs = PETSC_FALSE;
2096 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
2099 auto make_is_field_map = [&]() {
2100 PetscSection section;
2101 CHKERR DMGetDefaultSection(dm, §ion);
2104 CHKERR PetscSectionGetNumFields(section, &num_fields);
2109 map<std::string, SmartPetscObj<IS>> is_map;
2110 for (
int ff = 0; ff != num_fields; ff++) {
2112 const char *field_name;
2113 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
2117 is_map[field_name]);
2121 auto is_map = make_is_field_map();
2123 auto set_fieldsplit_on_bc = [&](PC &bc_pc,
string name_prb) {
2126 PetscBool is_bc_field_split;
2127 PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
2128 &is_bc_field_split);
2131 int l_has_bcs_ents = allBoundaryEnts->size();
2133 CHKERR MPIU_Allreduce(&l_has_bcs_ents, &g_has_bcs_ents, 1, MPIU_INT,
2134 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2136 if (is_bc_field_split && g_has_bcs_ents) {
2138 "Applying fieldsplit preconditioner for %d "
2139 "boundary entities.",
2143 CHKERR mField.get_moab().get_connectivity(*allBoundaryEnts, nodes,
2145 allBoundaryEnts->merge(nodes);
2151 name_prb,
ROW,
"U", 0, 3, is_bc, allBoundaryEnts.get());
2154 CHKERR PCFieldSplitSetIS(bc_pc, PETSC_NULL,
2158 CHKERR PCSetFromOptions(bc_pc);
2159 CHKERR PCSetType(bc_pc, PCLU);
2166 auto set_global_mat_and_vec = [&]() {
2170 CHKERR TSSetIFunction(solver, fe->ts_F, PETSC_NULL, PETSC_NULL);
2171 CHKERR TSSetIJacobian(solver, fe->ts_B, fe->ts_B, PETSC_NULL, PETSC_NULL);
2172 CHKERR KSPSetOperators(ksp, fe->ts_B, fe->ts_B);
2177 if (is_pcfs == PETSC_TRUE && !data.is_fieldsplit_bc_only) {
2179 CHKERR set_global_mat_and_vec();
2181 PetscBool is_l0_field_split;
2182 PetscBool is_l1_field_split = PETSC_TRUE;
2184 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l0_field_split);
2186 auto set_fieldsplit_contact = [&]() {
2193 is_map[
"E_IS_SIG"] = sig_data->rowIs;
2194 sigma_ao = sig_data->getSmartRowMap();
2196 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"SIGMA"]);
2197 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"E_IS_SIG"]);
2208 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &ct_ksp);
2210 CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
2211 PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
2212 &is_l1_field_split);
2216 CHKERR PetscFree(ct_ksp);
2221 if (data.with_contact && is_l0_field_split)
2222 CHKERR set_fieldsplit_contact();
2224 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
2225 if (data.with_plasticity && is_l0_field_split && is_l1_field_split) {
2231 is_map[
"E_IS_EP"] = ep_data->rowIs;
2234 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"EP"]);
2237 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"EP"]);
2238 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"E_IS_EP"]);
2240 PCCompositeType pc_type;
2241 CHKERR PCFieldSplitGetType(pC, &pc_type);
2243 if (pc_type == PC_COMPOSITE_SCHUR) {
2245 commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
2247 commonDataPtr->aoSEpEp = sigma_ao;
2249 commonDataPtr->blockMatContainerPtr =
2250 boost::make_shared<BlockMatContainer>();
2252 commonDataPtr->blockMatContainerPtr;
2254 CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
2255 commonDataPtr->SEpEp);
2262 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &tt_ksp);
2264 CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
2265 PetscBool is_tau_field_split;
2266 PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
2267 &is_tau_field_split);
2269 if (is_tau_field_split) {
2275 is_map[
"E_IS_TAU"] = tau_data->rowIs;
2278 CHKERR ep_data->getRowMap(&ep_ao);
2281 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"TAU"]);
2282 CHKERR AOApplicationToPetscIS(ep_ao, is_map[
"TAU"]);
2284 CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map[
"TAU"]);
2285 CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map[
"E_IS_TAU"]);
2287 CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
2288 if (pc_type == PC_COMPOSITE_SCHUR) {
2290 commonDataPtr->aoSTauTau = tau_data->getSmartRowMap();
2292 CHKERR PCFieldSplitSetSchurPre(tau_pc,
2293 PC_FIELDSPLIT_SCHUR_PRE_USER,
2294 commonDataPtr->STauTau);
2300 CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
2302 CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
2304 CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->
getName());
2305 CHKERR PetscFree(bc_ksp);
2308 CHKERR PetscFree(tt_ksp);
2311 }
else if (is_pcfs && data.is_fieldsplit_bc_only) {
2313 char mumps_options[] =
"-fieldsplit_0_mat_mumps_icntl_14 1200 "
2314 "-fieldsplit_0_mat_mumps_icntl_24 1 "
2315 "-fieldsplit_0_mat_mumps_icntl_13 1 "
2316 "-fieldsplit_1_mat_mumps_icntl_14 1200 "
2317 "-fieldsplit_1_mat_mumps_icntl_24 1 "
2318 "-fieldsplit_1_mat_mumps_icntl_13 1";
2319 CHKERR PetscOptionsInsertString(NULL, mumps_options);
2320 CHKERR set_global_mat_and_vec();
2321 CHKERR set_fieldsplit_on_bc(pC,
simple->getProblemName());
2325 CHKERR TSGetSNES(solver, &snes);
2327 CHKERR SNESGetKSP(snes, &ksp);
2329 CHKERR KSPGetPC(ksp, &pCc);
2330 CHKERR PCSetFromOptions(pCc);
2331 CHKERR PCSetType(pCc, PCLU);
2335 auto set_section_monitor = [&]() {
2338 CHKERR TSGetSNES(solver, &snes);
2339 PetscViewerAndFormat *vf;
2340 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2341 PETSC_VIEWER_DEFAULT, &vf);
2344 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
2349 auto create_post_process_element = [&]() {
2351 postProcFe = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
2352 postProcFe->generateReferenceElementMesh();
2354 postProcFe->getOpPtrVector().push_back(
2357 postProcFe->getOpPtrVector().push_back(
2359 if (data.is_large_strain) {
2360 postProcFe->getOpPtrVector().push_back(
2363 postProcFe->getOpPtrVector().push_back(
new OpStrain(
"U", commonDataPtr));
2364 if (data.with_plasticity) {
2366 "TAU", commonDataPtr->plasticTauPtr));
2367 postProcFe->getOpPtrVector().push_back(
2369 "EP", commonDataPtr->plasticStrainPtr));
2370 postProcFe->getOpPtrVector().push_back(
2372 postProcFe->getOpPtrVector().push_back(
2375 postProcFe->getOpPtrVector().push_back(
2376 new OpStress(
"U", commonDataPtr, commonDataPtr->mtD));
2378 if (data.is_large_strain) {
2379 if (data.is_neohooke)
2380 postProcFe->getOpPtrVector().push_back(
2383 postProcFe->getOpPtrVector().push_back(
2385 commonDataPtr->mtD));
2387 if (data.with_contact) {
2388 postProcFe->getOpPtrVector().push_back(
2390 "SIGMA", commonDataPtr->contactStressDivergencePtr));
2391 postProcFe->getOpPtrVector().push_back(
2393 "SIGMA", commonDataPtr->contactStressPtr));
2395 postProcFe->getOpPtrVector().push_back(
2397 postProcFe->mapGaussPts, commonDataPtr));
2399 if (data.is_large_strain)
2400 postProcFe->getOpPtrVector().push_back(
2402 postProcFe->mapGaussPts, commonDataPtr));
2404 postProcFe->getOpPtrVector().push_back(
2406 postProcFe->mapGaussPts, commonDataPtr));
2408 if (data.with_plasticity)
2409 postProcFe->getOpPtrVector().push_back(
2411 postProcFe->mapGaussPts, commonDataPtr));
2413 if (data.calculate_reactions)
2415 postProcFe->getOpPtrVector().push_back(
2417 postProcFe->mapGaussPts, commonDataPtr));
2418 postProcFe->addFieldValuesPostProc(
"U",
"DISPLACEMENT");
2419 if (data.is_ho_geometry)
2420 postProcFe->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
2424 auto scatter_create = [&](
auto coeff) {
2426 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
2427 ROW,
"U", coeff, coeff, is);
2429 CHKERR ISGetLocalSize(is, &loc_size);
2431 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
2433 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
2438 auto set_time_monitor = [&]() {
2440 boost::shared_ptr<MMonitor> monitor_ptr(
2441 new MMonitor(dm, mField, commonDataPtr, postProcFe, uXScatter,
2442 uYScatter, uZScatter));
2443 boost::shared_ptr<ForcesAndSourcesCore>
null;
2445 monitor_ptr,
null,
null);
2449 CHKERR set_section_monitor();
2450 CHKERR create_post_process_element();
2452 if (data.is_restart) {
2454 CHKERR PetscPrintf(PETSC_COMM_WORLD,
2455 "Reading vector in binary from %s file...\n",
2456 data.restart_file_str.c_str());
2458 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
2459 data.restart_file_str.c_str(), FILE_MODE_READ,
2463 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2464 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2467 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
2468 string s = data.restart_file_str;
2469 Tokenizer tok{s, boost::char_separator<char>(
"_")};
2470 auto it = tok.begin();
2471 const int restart_step = std::stoi((++it)->c_str());
2472 data.restart_step = restart_step;
2473 string restart_tt = *(++it);
2474 restart_tt.erase(restart_tt.length() - 4);
2475 const double restart_time = std::atof(restart_tt.c_str());
2476 CHKERR TSSetStepNumber(solver, restart_step);
2477 CHKERR TSSetTime(solver, restart_time);
2480 uXScatter = scatter_create(0);
2481 uYScatter = scatter_create(1);
2482 uZScatter = scatter_create(2);
2483 CHKERR set_time_monitor();
2503 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2504 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2516 auto meshset = mField.getInterface<
Simple>()->getMeshSet();
2517 CHKERR mField.get_moab().get_entities_by_dimension(meshset, 3, body_ents);
2519 Skinner skin(&mField.get_moab());
2521 CHKERR skin.find_skin(0, body_ents,
false, skin_tris);
2522 Range boundary_ents;
2523 ParallelComm *pcomm =
2525 if (pcomm == NULL) {
2526 pcomm =
new ParallelComm(&mField.get_moab(), mField.get_comm());
2528 CHKERR pcomm->filter_pstatus(skin_tris, PSTATUS_SHARED | PSTATUS_MULTISHARED,
2529 PSTATUS_NOT, -1, &boundary_ents);
2531 return boundary_ents;
2544 DMType dm_name =
"DMMOFEM";
2554 simple->getProblemName() =
"Multifield plasticity";
2560 PetscBool is_partitioned = PETSC_FALSE;
2562 &is_partitioned, PETSC_NULL);
2563 if (is_partitioned) {
2580 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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
#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 .
@ 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 ...
auto smartCreateDMVector
Get smart vector from DM.
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
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpSkeletonRhsPipeline()
Get the Op Skeleton Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
#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
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
int main(int argc, char *argv[])
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 >::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
const FTensor::Tensor2< T, Dim, Dim > Vec
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.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
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)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
auto createSmartGhostVector
Create smart ghost vector.
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
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)
DeprecatedCoreInterface Interface
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)
const double D
diffusivity
double duration
impulse duration (ns)
static constexpr int approx_order
FTensor::Index< 'm', 3 > m
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.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
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
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.