v0.14.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
ContactPlasticity Struct Reference
Collaboration diagram for ContactPlasticity:
[legend]

Classes

struct  MMonitor
 
struct  ProblemData
 

Public Member Functions

 ContactPlasticity (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
MoFEMErrorCode loadMeshBlockSets ()
 

Static Public Member Functions

static MoFEMErrorCode getAnalysisParameters ()
 
static MoFEMErrorCode myMeshPartition (MoFEM::Interface &m_field)
 

Public Attributes

FieldApproximationBase base
 

Static Public Attributes

static ProblemData data
 

Private Member Functions

MoFEMErrorCode setUP ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode OPs ()
 
MoFEMErrorCode bC ()
 
MoFEMErrorCode tsSolve ()
 
MoFEMErrorCode postProcess ()
 
MoFEMErrorCode checkResults ()
 
Range getEntsOnMeshSkin ()
 

Private Attributes

MoFEM::InterfacemField
 
MatrixDouble invJac
 
MatrixDouble jAc
 
boost::shared_ptr< OpContactTools::CommonDatacommonDataPtr
 
boost::shared_ptr< PostProcVolumeOnRefinedMeshpostProcFe
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
boost::ptr_map< std::string, NeumannForcesSurfaceneumann_forces
 
boost::ptr_map< std::string, NodalForcenodal_forces
 
boost::ptr_map< std::string, EdgeForceedge_forces
 
boost::shared_ptr< DirichletDisplacementRemoveDofsBcdirichlet_bc_ptr
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
boost::ptr_vector< DataFromMovebcData
 
boost::shared_ptr< DomainSideElevolSideFe
 
boost::shared_ptr< DomainElefeLambdaRhs
 
boost::shared_ptr< DomainElefeLambdaLhs
 
SmartPetscObj< DM > dmEP
 
SmartPetscObj< DM > dmTAU
 
SmartPetscObj< DM > dmContact
 

Detailed Description

Definition at line 107 of file multifield_plasticity.cpp.

Constructor & Destructor Documentation

◆ ContactPlasticity()

ContactPlasticity::ContactPlasticity ( MoFEM::Interface m_field)
inline

Definition at line 109 of file multifield_plasticity.cpp.

109: mField(m_field) {}
MoFEM::Interface & mField

Member Function Documentation

◆ bC()

MoFEMErrorCode ContactPlasticity::bC ( )
private

Definition at line 1116 of file multifield_plasticity.cpp.

1116 {
1118 auto bc_mng = mField.getInterface<BcManager>();
1119 auto prb_mng = mField.getInterface<ProblemsManager>();
1120 auto simple = mField.getInterface<Simple>();
1121
1122 auto get_ents = [&](const std::string blockset_name) {
1123 Range fix_ents;
1125 if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
1126 0) {
1127 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
1128 true);
1129 }
1130 }
1131 return fix_ents;
1132 };
1133
1134 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
1135 "U", 0, 0);
1136 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
1137 "U", 1, 1);
1138 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
1139 "U", 2, 2);
1140 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
1141 "U", 0, 2);
1142
1143 if (data.with_contact) {
1144 Range boundary_ents;
1145 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
1146 if (bc_mng->checkBlock(bc, "FIX_"))
1147 boundary_ents.merge(*bc.second->getBcEntsPtr());
1148 boundary_ents.merge(get_ents("NO_CONTACT"));
1149 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1150 boundary_ents);
1151 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1152 mField.getInterface<Simple>()->getProblemName(), "SIGMA", boundary_ents,
1153 0, 2);
1154 }
1155
1156 std::vector<DataFromBc> bcData;
1157 dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementRemoveDofsBc>(
1158 mField, "U", simple->getProblemName(), "DISPLACEMENT", true);
1159 dirichlet_bc_ptr->iNitialize();
1160
1161 auto dm = simple->getDM();
1162 // FIXME: contact with fieldsplit does not work in parallel, requires further
1163 // investigation
1164 // TODO: will come back into this in the future
1165 if (data.with_contact && data.with_plasticity && false) {
1169 } else if (data.with_plasticity) {
1172 }
1173 // CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1174
1175 // forces and pressures on surface
1177 PETSC_NULL, "U");
1178 // nodal forces
1180 // edge forces
1182
1183 // FIXME: add higher order operators !!
1184 // TODO: add higher order operators (check nonlinear_dynamics)
1185 for (auto mit = neumann_forces.begin(); mit != neumann_forces.end(); mit++) {
1186 if (data.is_ho_geometry)
1187 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1188 false, false);
1189 mit->second->methodsOp.push_back(
1190 new TimeForceScale("-load_history", false));
1191 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1192 }
1193 for (auto mit = nodal_forces.begin(); mit != nodal_forces.end(); mit++) {
1194 if (data.is_ho_geometry)
1195 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1196 false, false);
1197 mit->second->methodsOp.push_back(
1198 new TimeForceScale("-load_history", false));
1199 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1200 }
1201 for (auto mit = edge_forces.begin(); mit != edge_forces.end(); mit++) {
1202 if (data.is_ho_geometry)
1203 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1204 false, false);
1205 mit->second->methodsOp.push_back(
1206 new TimeForceScale("-load_history", false));
1207 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1208 }
1209
1210 // dirichlet_bc_ptr->dIag = 1;
1211 dirichlet_bc_ptr->methodsOp.push_back(
1213
1215}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#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.
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)
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)
boost::ptr_vector< DataFromMove > bcData
boost::ptr_map< std::string, EdgeForce > edge_forces
boost::ptr_map< std::string, NeumannForcesSurface > neumann_forces
SmartPetscObj< DM > dmContact
boost::shared_ptr< DirichletDisplacementRemoveDofsBc > dirichlet_bc_ptr
static ProblemData data
boost::ptr_map< std::string, NodalForce > nodal_forces
SmartPetscObj< DM > dmEP
SmartPetscObj< DM > dmTAU
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:348
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.

◆ checkResults()

MoFEMErrorCode ContactPlasticity::checkResults ( )
private

Definition at line 2581 of file multifield_plasticity.cpp.

2581{ return 0; }

◆ createCommonData()

MoFEMErrorCode ContactPlasticity::createCommonData ( )
private

Definition at line 1001 of file multifield_plasticity.cpp.

1001 {
1003
1004 auto get_material_stiffness = [&]() {
1006
1007 constexpr auto t_kd = Kronecker_Delta<int>();
1008 constexpr auto t_kds = Kronecker_Delta_symmetric<int>();
1009
1010 // auto &t_D = commonDataPtr->tD;
1011 (*cache).Is(i, j, k, l) =
1012 (0.5 * t_kd(i, k) * t_kd(j, l)) + (0.5 * t_kd(i, l) * t_kd(j, k));
1013
1014 // Tensor4<double, 3, 3, 3, 3> II;
1015 // II(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
1016
1017 auto Is_ddg = tensor4_to_ddg((*cache).Is);
1018
1019 for (auto &mit : mat_blocks) {
1020
1021 cache = &mit.second;
1022 const double G = (*cache).young_modulus / (2. * (1. + (*cache).poisson));
1023 const double K =
1024 (*cache).young_modulus / (3. * (1. - 2. * (*cache).poisson));
1025
1026 // for plane stress
1027 double A = 6. * G / (3. * K + 4. * G);
1028 // is_plane_strain and 3D
1029 A = 1.;
1030
1031 auto t_D_tmp = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
1032 auto t_D_axiator =
1033 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Axiator);
1034 auto t_D_deviator =
1035 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
1036
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.);
1040
1041 t_D_tmp(i, j, k, l) = t_D_deviator(i, j, k, l) + t_D_axiator(i, j, k, l);
1042
1043 (*cache).mtD = *commonDataPtr->mtD;
1044 (*cache).scale_constraint = 1.;
1045 }
1046
1047 // commonDataPtr->isDualBase = data.is_dual_base;
1048
1050 };
1051
1052 commonDataPtr = boost::make_shared<OpContactTools::CommonData>();
1053
1054 commonDataPtr->mtD = boost::make_shared<MatrixDouble>();
1055 commonDataPtr->mtD->resize(36, 1);
1056 commonDataPtr->mtD_Axiator = boost::make_shared<MatrixDouble>();
1057 commonDataPtr->mtD_Axiator->resize(36, 1);
1058 commonDataPtr->mtD_Deviator = boost::make_shared<MatrixDouble>();
1059 commonDataPtr->mtD_Deviator->resize(36, 1);
1060
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>();
1065 commonDataPtr->mPiolaStressPtr->resize(9, 1, false);
1066
1067 commonDataPtr->mDeformationGradient = boost::make_shared<MatrixDouble>();
1068 commonDataPtr->matC = boost::make_shared<MatrixDouble>();
1069 commonDataPtr->matEigenVal = boost::make_shared<MatrixDouble>();
1070 commonDataPtr->matEigenVector = boost::make_shared<MatrixDouble>();
1071 commonDataPtr->detF = boost::make_shared<VectorDouble>();
1072 commonDataPtr->materialTangent = boost::make_shared<MatrixDouble>();
1073 commonDataPtr->dE_dF_mat = boost::make_shared<MatrixDouble>();
1074 commonDataPtr->dE_dF_mat->resize(81, 1, false);
1075
1076 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
1077 commonDataPtr->contactStressDivergencePtr =
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>();
1083
1084 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
1085 commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
1086
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>();
1091
1092 commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
1093 commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
1094
1095 commonDataPtr->plasticTauJumpPtr = boost::make_shared<VectorDouble>();
1096 commonDataPtr->plasticStrainJumpPtr = boost::make_shared<MatrixDouble>();
1097 commonDataPtr->guidingVelocityPtr = boost::make_shared<MatrixDouble>();
1098 commonDataPtr->guidingVelocityPtr->resize(3, 1, false);
1099
1100 commonDataPtr->plasticGradTauPtr = boost::make_shared<MatrixDouble>();
1101 commonDataPtr->plasticGradStrainPtr = boost::make_shared<MatrixDouble>();
1102
1103 jAc.resize(3, 3, false);
1104 invJac.resize(3, 3, false);
1105
1107 CHKERR commonDataPtr->calculateDiffDeviator();
1108
1109 CHKERR get_material_stiffness();
1110 if (data.is_neohooke)
1111 commonDataPtr->isNeoHook = true;
1112
1114}
Kronecker Delta class symmetric.
Kronecker Delta class.
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
std::map< int, BlockParamData > mat_blocks
BlockParamData * cache
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
constexpr AssemblyType A
double young_modulus
boost::shared_ptr< OpContactTools::CommonData > commonDataPtr

◆ getAnalysisParameters()

MoFEMErrorCode ContactPlasticity::getAnalysisParameters ( )
static

Definition at line 602 of file multifield_plasticity.cpp.

602 {
604
605 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
606 // Set approximation order
607 CHKERR PetscOptionsInt("-order", "approximation order", "", data.order,
608 &data.order, PETSC_NULL);
609 CHKERR PetscOptionsBool("-no_output", "save no output files", "",
610 data.no_output, &data.no_output, PETSC_NULL);
611 CHKERR PetscOptionsBool("-debug_info", "print debug info", "",
612 data.debug_info, &data.debug_info, PETSC_NULL);
613 CHKERR PetscOptionsBool("-scale_params",
614 "scale parameters (young modulus = 1)", "",
615 data.scale_params, &data.scale_params, PETSC_NULL);
616 CHKERR PetscOptionsInt("-output_every",
617 "frequncy how often results are dumped on hard disk",
618 "", 1, &data.output_freq, PETSC_NULL);
619 CHKERR PetscOptionsBool("-calculate_reactions", "calculate reactions", "",
621 PETSC_NULL);
622 CHKERR PetscOptionsInt("-reaction_id", "meshset id for computing reactions",
623 "", data.reaction_id, &data.reaction_id, PETSC_NULL);
624 CHKERR PetscOptionsInt("-atom_test", "number of atom test", "",
625 data.atom_test_nb, &data.atom_test_nb, PETSC_NULL);
626 CHKERR PetscOptionsBool("-is_large_strain", "is large strains regime", "",
628 PETSC_NULL);
629 CHKERR PetscOptionsBool(
630 "-is_fieldsplit_bc_only", "use fieldsplit for boundary only", "",
632 CHKERR PetscOptionsBool(
633 "-is_reduced_integration",
634 "use axiator split approach, with reduced integration", "",
636 CHKERR PetscOptionsBool("-test_user_base_tau", "test_user_base_for_tau", "",
638 PETSC_NULL);
639 CHKERR PetscOptionsBool("-is_alm",
640 "use Augmented Lagrangian method for enforcing the "
641 "KKT condtions for plasticity",
642 "", data.is_alm, &data.is_alm, PETSC_NULL);
643
644 CHKERR PetscOptionsBool("-is_neohooke", "is neohooke material", "",
645 data.is_neohooke, &data.is_neohooke, PETSC_NULL);
646 // CHKERR PetscOptionsBool(
647 // "-use_fieldsplit_for_bc", "use fieldsplit for boundary conditions", "",
648 // data.use_fieldsplit_for_bc, &data.use_fieldsplit_for_bc, PETSC_NULL);
649
650 CHKERR PetscOptionsBool(
651 "-with_plasticity", "run calculations with plasticity", "",
653
654 CHKERR PetscOptionsBool("-is_rotating", "is rotating frame used", "",
655 data.is_rotating, &data.is_rotating, PETSC_NULL);
656 CHKERR PetscOptionsBool("-with_contact", "run calculations with contact", "",
657 data.with_contact, &data.with_contact, PETSC_NULL);
658 CHKERR PetscOptionsBool("-print_rollers",
659 "output file with rollers positions", "",
660 data.print_rollers, &data.print_rollers, PETSC_NULL);
661
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);
666 if (ctg_flag)
667 data.contact_history = "-contact_history";
668 CHKERR PetscOptionsString("-move_history", "load_history.in", "",
669 "load_history.in", load_hist_file, 255, &ctg_flag);
670
671 if (ctg_flag)
672 data.move_history = "-move_history";
673
674 CHKERR PetscOptionsString("-dirichlet_history", "load_history.in", "",
675 "load_history.in", load_hist_file, 255, &ctg_flag);
676 if (ctg_flag)
677 data.dirichlet_history = "-dirichlet_history";
678
679 char mesh_file_name[255];
680 CHKERR PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
681 mesh_file_name, 255, PETSC_NULL);
683
684 char restart_file_name[255];
685 CHKERR PetscOptionsString("-restart_file", "restart file name", "",
686 "restart.dat", restart_file_name, 255,
688 data.restart_file_str = string(restart_file_name);
689
690 CHKERR PetscOptionsBool("-save_restarts",
691 "save restart files at each iteration", "",
692 data.save_restarts, &data.save_restarts, PETSC_NULL);
693
694 ierr = PetscOptionsEnd();
695 CHKERRQ(ierr);
696
698 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
699 "Neohooke material is not supported with plasticity.");
700
702}
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76

◆ getEntsOnMeshSkin()

Range ContactPlasticity::getEntsOnMeshSkin ( )
private

Definition at line 2583 of file multifield_plasticity.cpp.

2583 {
2584 Range body_ents;
2585 auto meshset = mField.getInterface<Simple>()->getMeshSet();
2586 CHKERR mField.get_moab().get_entities_by_dimension(meshset, 3, body_ents);
2587 // cout << "body_ents for skin: " << body_ents.size() << endl;
2588 Skinner skin(&mField.get_moab());
2589 Range skin_tris;
2590 CHKERR skin.find_skin(0, body_ents, false, skin_tris);
2591 Range boundary_ents;
2592 ParallelComm *pcomm =
2593 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2594 if (pcomm == NULL) {
2595 pcomm = new ParallelComm(&mField.get_moab(), mField.get_comm());
2596 }
2597 CHKERR pcomm->filter_pstatus(skin_tris, PSTATUS_SHARED | PSTATUS_MULTISHARED,
2598 PSTATUS_NOT, -1, &boundary_ents);
2599
2600 return boundary_ents;
2601};
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
virtual MPI_Comm & get_comm() const =0

◆ loadMeshBlockSets()

MoFEMErrorCode ContactPlasticity::loadMeshBlockSets ( )

Definition at line 793 of file multifield_plasticity.cpp.

793 {
795 string block_name = "MAT_ELASTIC";
796 // TODO: implement blocks for plasticity
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();
803 EntityHandle meshset = it->getMeshset();
804 Range tets;
805 // on a particular part, the meshset can be empty
806 rval =
807 mField.get_moab().get_entities_by_dimension(meshset, 3, tets, true);
808 rval = MB_SUCCESS;
809 for (auto ent : tets)
810 data.mapBlockTets[ent] = id;
811
813 mat_blocks.at(id).young_modulus = block_data[0];
814 mat_blocks.at(id).poisson = block_data[1];
815 if (block_data.size() > 2)
816 mat_blocks.at(id).density = block_data[2];
817 mat_blocks.at(id).cn_cont = 1. / block_data[0];
818
819 if (block_data.size() < 2)
820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "provide an appropriate number of entries (2) parameters for "
822 "MAT_ELASTIC block");
823 }
824 }
825
826 if (mat_blocks.size() < 2)
828
829 cache = &mat_blocks.begin()->second;
830
832}
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
std::map< EntityHandle, int > mapBlockTets

◆ myMeshPartition()

MoFEMErrorCode ContactPlasticity::myMeshPartition ( MoFEM::Interface m_field)
static

Definition at line 704 of file multifield_plasticity.cpp.

704 {
706
707 auto &moab = m_field.get_moab();
708 Range tets;
709 CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
710 Range faces;
711 CHKERR moab.get_adjacencies(tets, 2, true, faces, moab::Interface::UNION);
712 Range edges;
713 CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
714
715 CHKERR m_field.getInterface<ProblemsManager>()->partitionMesh(
716 tets, 3, 0, m_field.get_comm_size());
717
718 EntityHandle part_set;
719 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
720 if (pcomm == NULL)
721 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
722 CHKERR moab.create_meshset(MESHSET_SET, part_set);
723 Tag part_tag = pcomm->part_tag();
724 Range proc_ents;
725 Range all_proc_ents;
726 Range off_proc_ents;
727 Range tagged_sets;
728 CHKERR moab.get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag, NULL, 1,
729 tagged_sets, moab::Interface::UNION);
730
731 for (auto &mit : tagged_sets) {
732 int part;
733 CHKERR moab.tag_get_data(part_tag, &mit, 1, &part);
734 if (part == m_field.get_comm_rank()) {
735 CHKERR moab.get_entities_by_dimension(mit, 3, proc_ents, true);
736 CHKERR moab.get_entities_by_handle(mit, all_proc_ents, true);
737 } else
738 CHKERR moab.get_entities_by_handle(mit, off_proc_ents, true);
739 }
740
741 Skinner skin(&moab);
742 Range proc_ents_skin[4];
743 proc_ents_skin[3] = proc_ents;
744
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]);
749
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);
757
758 // not sure if necessary...
759 {
760 Range all_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());
765 }
766
767 Range ents_to_keep;
768 CHKERR moab.get_entities_by_handle(part_set, ents_to_keep, false);
769 off_proc_ents = subtract(off_proc_ents, ents_to_keep);
770 Range meshsets;
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);
774 // CHKERR moab.add_entities(part_set, &m, 1);
775 }
776 CHKERR moab.delete_entities(off_proc_ents);
777
778 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
779
780 BitRefLevel bit_level0;
781 bit_level0.set(0);
782 Simple *simple = m_field.getInterface<Simple>();
783
784 Range ents;
785 CHKERR m_field.get_moab().get_entities_by_dimension(part_set, 3, ents, true);
786 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
787 ents, simple->getBitRefLevel(), false);
788 simple->getMeshSet() = part_set;
789
791}
FTensor::Index< 'm', SPACE_DIM > m
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)
Definition: ddTensor0.hpp:33
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual int get_comm_rank() const =0

◆ OPs()

MoFEMErrorCode ContactPlasticity::OPs ( )
private

Definition at line 1217 of file multifield_plasticity.cpp.

1217 {
1221 auto problem_manager = mField.getInterface<ProblemsManager>();
1222 auto bc_mng = mField.getInterface<BcManager>();
1223 const Field_multiIndex *fields_ptr;
1224 CHKERR mField.get_fields(&fields_ptr);
1225
1226 auto get_boundary_markers = [&](string prob, Range bound_ents, int lo,
1227 int hi) {
1228 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
1229
1230 CHKERR problem_manager->modifyMarkDofs(
1231 prob, ROW, "U", lo, hi, ProblemsManager::MarkOP::OR, 1, *marker_ptr);
1232
1233 CHKERR problem_manager->markDofs(prob, ROW, ProblemsManager::AND,
1234 bound_ents, *marker_ptr);
1235
1236 return marker_ptr;
1237 };
1238
1239 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_X", "U",
1240 0, 0);
1241 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Y", "U",
1242 1, 1);
1243 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Z", "U",
1244 2, 2);
1245
1246 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_X",
1247 "U", 1, 2);
1248 // TODO: implement this
1249 // CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(),
1250 // "ROTATE_Y",
1251 // "U", 0, 0);
1252 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_Z",
1253 "U", 0, 1);
1254 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_ALL",
1255 "U", 0, 3);
1256
1257 auto &bc_map = bc_mng->getBcMapByBlockName();
1259 bc_mng->getMergedBlocksMarker(vector<string>{"MOVE_", "ROTATE_"});
1260
1261 auto add_disp_bc_ops = [&]() {
1263 std::map<char, size_t> xyz_map = {{'X', 0}, {'Y', 1}, {'Z', 2}};
1264
1266 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1268
1269 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
1270 if (bc_mng->checkBlock(bc, "MOVE_")) {
1271
1272 int idx = xyz_map.at(bc.first[(bc.first.find("MOVE_") + 5)]);
1273
1274 auto disp_vec = bc.second->bcAttributes;
1275 if (disp_vec.empty())
1276 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1277 "provide an appropriate number of params (1) for MOVE block");
1278 VectorDouble disp({0, 0, 0});
1279 disp(idx) = disp_vec[0];
1280
1281 bcData.push_back(new DataFromMove(disp, *bc.second->getBcEntsPtr()));
1282 // rhs
1283 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1284 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1285 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1286 new OpEssentialRHS_U("U", commonDataPtr, bcData.back().scaledValues,
1287 bcData.back().ents));
1288 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("U"));
1289 // lhs
1290 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1291 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1292 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpSpringLhsNoFS(
1293 "U", "U", [](double, double, double) { return 1.; },
1294 boost::make_shared<Range>(bcData.back().ents)));
1295 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("U"));
1296 }
1297 }
1298
1299 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
1300 if (std::regex_match(bc.first, std::regex("(.*)_ROTATE_(.*)"))) {
1301 MOFEM_LOG("WORLD", Sev::inform)
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)
1306 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
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());
1315 bcData.push_back(new DataFromMove(*angles, *bc.second->getBcEntsPtr()));
1316 // rhs
1317 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1318 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1319 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1320 new OpRotate("U", bcData.back().scaledValues, c_coords,
1321 bc.second->getBcEntsPtr()));
1322 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpSpringRhs(
1323 "U", commonDataPtr->mDispPtr,
1324 [](double, double, double) { return 1.; },
1325 bc.second->getBcEntsPtr()));
1326 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("U"));
1327 // lhs
1328 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1329 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1330 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpSpringLhsNoFS(
1331 "U", "U", [](double, double, double) { return 1.; },
1332 bc.second->getBcEntsPtr()));
1333 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("U"));
1334 }
1335 }
1336
1338 };
1339
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>();
1345 pipeline.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1346 "MESH_NODE_POSITIONS", jac_ptr));
1347 // pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
1348 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
1349 pipeline.push_back(
1350 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
1351 pipeline.push_back(new OpSetHOWeights(det_ptr));
1352 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
1353 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
1354 pipeline.push_back(new OpSetHOWeights(det_ptr));
1355 // pipeline.push_back(new OpMakeHighOrderGeometryWeightsOnVolume());
1356 // pipeline.push_back(new OpSetHOGeometryCoordsOnVolume());
1358 };
1359
1360 auto add_domain_base_ops = [&](auto &pipeline) {
1362 pipeline.push_back(
1364 if (data.is_ho_geometry)
1365 CHKERR add_ho_ops(pipeline);
1366
1367 pipeline.push_back(
1369 pipeline.push_back(
1371
1372 if (data.with_plasticity) {
1373 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1374 "EP", commonDataPtr->plasticStrainPtr));
1375 pipeline.push_back(new OpCalculateScalarFieldValues(
1376 "TAU", commonDataPtr->plasticTauPtr));
1377 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
1378 "TAU", commonDataPtr->plasticTauDotPtr));
1379 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<3>(
1380 "EP", commonDataPtr->plasticStrainDotPtr));
1381
1382 if (data.is_rotating) {
1384 "EP", commonDataPtr->plasticGradStrainPtr));
1385 pipeline.push_back(new OpCalculateScalarFieldGradient<3>(
1386 "TAU", commonDataPtr->plasticGradTauPtr));
1387 pipeline.push_back(
1389 }
1390 }
1391
1392 if (data.is_large_strain) {
1393 pipeline.push_back(new OpLogStrain("U", commonDataPtr));
1394 } else
1395 pipeline.push_back(new OpStrain("U", commonDataPtr));
1396
1398 };
1399
1400 auto add_domain_stress_ops = [&](auto &pipeline, auto mat_D_ptr) {
1402
1403 if (data.with_plasticity) {
1404 pipeline.push_back(new OpPlasticStress("U", commonDataPtr, mat_D_ptr));
1405
1406 pipeline.push_back(
1408
1409 } else
1410 pipeline.push_back(new OpStress("U", commonDataPtr, mat_D_ptr));
1411
1413 };
1414
1415 auto add_skeleton_base_ops = [&](auto &pipeline) {
1416 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
1417
1418 volSideFe->getOpPtrVector().push_back(
1420 volSideFe->getOpPtrVector().push_back(
1422 // pipeline.push_back(new OpCalculateJumpOnSkeleton("SKELETON_LAMBDA",
1423 // commonDataPtr,
1424 // volSideFe));
1425 };
1426
1427 auto add_domain_ops_lhs = [&](auto &pipeline, auto m_D_ptr) {
1429
1430 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1432 pipeline.push_back(new OpSchurBegin("U", commonDataPtr));
1433
1434 if (data.is_large_strain) {
1435 if (data.is_neohooke)
1436 pipeline.push_back(
1438 else
1439 pipeline.push_back(
1441 pipeline.push_back(
1442 new OpLogStrainMatrixLhs("U", "U", commonDataPtr->materialTangent));
1443 } else
1444 pipeline.push_back(new OpStiffnessMatrixLhs("U", "U", m_D_ptr));
1445 if (data.with_plasticity) {
1447 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<true>(
1448 "U", "EP", commonDataPtr, m_D_ptr));
1449 else
1450 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<false>(
1451 "U", "EP", commonDataPtr, m_D_ptr));
1452 // FIXME: dummy one for testing
1453 // TODO: this can be deleted
1454 // pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dTAU(
1455 // "U", "TAU", commonDataPtr));
1456 }
1457
1458 if (data.is_rotating)
1459 pipeline.push_back(new OpRotatingFrameLhs("U", "U", commonDataPtr));
1460
1461 pipeline.push_back(new OpUnSetBc("U"));
1462
1464 };
1465
1466 auto add_domain_ops_rhs = [&](auto &pipeline, auto m_D_ptr) {
1468
1469 auto get_body_force = [this](const double, const double, const double) {
1470 auto *pipeline_mng = mField.getInterface<PipelineManager>();
1471 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
1472 Tensor1<double, 3> t_source;
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;
1477 return t_source;
1478 };
1479 auto get_centrifugal_force = [this](const double x, const double y,
1480 const double z) {
1481 Tensor1<double, 3> t_coords(x, y, z);
1482 Tensor1<double, 3> t_source;
1483 t_source(i) = (*cache).density * (*cache).Omega(i, k) *
1484 (*cache).Omega(k, j) * t_coords(j);
1485 // auto *pipeline_mng = mField.getInterface<PipelineManager>();
1486 // auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
1487 // const auto time = fe_domain_rhs->ts_t;
1488 // t_source(i) *= time;
1489 return t_source;
1490 };
1491
1492 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1493
1494 if (data.is_rotating) {
1495 pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
1496 // pipeline.push_back(new OpCentrifugalForce2("U",
1497 // get_centrifugal_force));
1498 }
1499
1500 pipeline.push_back(new OpBodyForce("U", get_body_force));
1501 if (data.is_large_strain) {
1502 if (data.is_neohooke)
1503 pipeline.push_back(
1505 else
1506 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1507 "U", commonDataPtr, m_D_ptr));
1508 pipeline.push_back(
1510 } else
1511 pipeline.push_back(
1513
1514 pipeline.push_back(new OpUnSetBc("U"));
1516 };
1517
1518 auto add_domain_ops_rhs_constraints = [&](auto &pipeline) {
1520 if (data.with_plasticity) {
1521 pipeline.push_back(new OpCalculatePlasticFlowRhs("EP", commonDataPtr));
1522 pipeline.push_back(new OpCalculateConstraintRhs("TAU", commonDataPtr));
1523 if (data.debug_info)
1524 pipeline.push_back(new OpGetGaussPtsPlasticState("U", commonDataPtr));
1525 }
1527 };
1528
1529 auto add_domain_ops_lhs_constraints = [&](auto &pipeline, auto m_D_ptr) {
1531
1532 if (data.with_plasticity) {
1533
1534 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1535
1536 if (data.is_large_strain) {
1537 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1538 "U", commonDataPtr, m_D_ptr));
1539 pipeline.push_back(
1541 pipeline.push_back(
1543 } else {
1544 pipeline.push_back(
1546 pipeline.push_back(
1548 }
1549
1550 pipeline.push_back(
1552 pipeline.push_back(
1554 pipeline.push_back(
1556 pipeline.push_back(
1557 new OpCalculateConstraintLhs_dTAU("TAU", "TAU", commonDataPtr));
1558
1559 if (data.with_plasticity) {
1560 pipeline.push_back(new OpSchurPreconditionMassTau(
1561 "TAU", "TAU", [](double, double, double) { return 1.; },
1562 commonDataPtr));
1563 pipeline.push_back(new OpSchurEnd("U", commonDataPtr));
1564 }
1565 pipeline.push_back(new OpUnSetBc("U"));
1566 }
1567
1569 };
1570
1571 auto add_domain_ops_contact_lhs = [&](auto &pipeline) {
1572 if (data.with_contact) {
1573 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1574
1575 pipeline.push_back(new OpMixDivULhs(
1576 "SIGMA", "U", []() { return 1.; }, true));
1577 pipeline.push_back(new OpLambdaGraULhs(
1578 "SIGMA", "U", []() { return 1.; }, true));
1579
1580 pipeline.push_back(new OpSchurPreconditionMassContactVol(
1581 "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
1582 commonDataPtr));
1583
1584 pipeline.push_back(new OpUnSetBc("U"));
1585 }
1586 };
1587
1588 auto add_domain_ops_contact_rhs = [&](auto &pipeline) {
1589 if (data.with_contact) {
1590 pipeline.push_back(new OpCalculateHVecTensorField<3, 3>(
1591 "SIGMA", commonDataPtr->contactStressPtr));
1592 pipeline.push_back(new OpCalculateHVecTensorDivergence<3, 3>(
1593 "SIGMA", commonDataPtr->contactStressDivergencePtr));
1594 pipeline.push_back(
1595 new OpMixDivURhs("SIGMA", commonDataPtr->mDispPtr,
1596 [](double, double, double) { return 1; }));
1597 pipeline.push_back(
1598 new OpMiXLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
1599
1600 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1601 pipeline.push_back(new OpMixUTimesDivLambdaRhs(
1602 "U", commonDataPtr->contactStressDivergencePtr));
1603 pipeline.push_back(
1604 new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
1605 pipeline.push_back(new OpUnSetBc("U"));
1606 }
1607 };
1608
1609 auto add_boundary_base_ops = [&](auto &pipeline) {
1611 pipeline.push_back(
1613 if (data.is_ho_geometry) {
1614 pipeline.push_back(new OpSetHOWeightsOnFace());
1615 pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
1616 pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
1617 }
1618 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1619 // pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
1620
1621 pipeline.push_back(
1623 if (data.with_contact) {
1624 pipeline.push_back(
1626 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1627 "SIGMA", commonDataPtr->contactTractionPtr));
1628 }
1630 };
1631
1632 auto add_boundary_ops_lhs = [&](auto &pipeline) {
1633 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1634
1635 // if (data.with_plasticity)
1636 // pipeline.push_back(new OpSchurBeginBoundary("U", commonDataPtr));
1637
1638 if (data.with_contact) {
1639 pipeline.push_back(
1640 new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
1641 pipeline.push_back(new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA",
1642 commonDataPtr));
1643 // FIXME: precondition matrix
1644 pipeline.push_back(new OpSchurPreconditionMassContactBound(
1645 "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
1646 commonDataPtr));
1647 }
1648
1649 // pipeline.push_back(new OpSpringLhs(
1650 // "U", "U",
1651 // [this](double, double, double) { return (*cache).spring_stiffness;
1652 // }
1653 // ));
1654
1655 if (data.is_rotating) {
1656 auto vol_side_fe_ptr = boost::make_shared<MoFEM::DomainSideEle>(mField);
1657
1658 if (data.is_ho_geometry)
1659 CHKERR add_ho_ops(vol_side_fe_ptr->getOpPtrVector());
1660
1661 vol_side_fe_ptr->getOpPtrVector().push_back(
1663
1664 pipeline.push_back(new OpRotatingFrameBoundaryLhs("U", "U", commonDataPtr,
1665 vol_side_fe_ptr));
1666 }
1668 pipeline.push_back(new OpSchurEndBoundary("U", commonDataPtr));
1669
1670 pipeline.push_back(new OpUnSetBc("U"));
1671 };
1672
1673 auto add_boundary_ops_rhs = [&](auto &pipeline) {
1674 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1676 pipeline.push_back(new OpGetGaussPtsContactState("SIGMA", commonDataPtr));
1677
1678 if (data.with_contact)
1679 pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
1680 // pipeline.push_back(new OpSpringRhs(
1681 // "U", commonDataPtr->mDispPtr,
1682 // [this](double, double, double) { return (*cache).spring_stiffness;
1683 // }));
1684
1685 if (data.is_rotating) {
1686 auto my_vol_side_fe_ptr =
1687 boost::make_shared<MoFEM::DomainSideEle>(mField);
1688 my_vol_side_fe_ptr->getOpPtrVector().push_back(
1690 commonDataPtr->mGradPtr));
1691 pipeline.push_back(new OpRotatingFrameBoundaryRhs("U", commonDataPtr,
1692 my_vol_side_fe_ptr));
1693 }
1694 pipeline.push_back(new OpUnSetBc("U"));
1695 };
1696
1697 auto &tD = commonDataPtr->mtD;
1698 auto &tD_axi = commonDataPtr->mtD_Axiator;
1699 auto &tD_dev = commonDataPtr->mtD_Deviator;
1700
1701 feLambdaLhs = boost::make_shared<DomainEle>(mField);
1702 feLambdaRhs = boost::make_shared<DomainEle>(mField);
1703
1705 feLambdaLhs->getUserPolynomialBase() =
1706 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
1707 feLambdaRhs->getUserPolynomialBase() =
1708 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
1709 }
1710
1711 auto &lam_pipeline_lhs = feLambdaLhs->getOpPtrVector();
1712 auto &lam_pipeline_rhs = feLambdaRhs->getOpPtrVector();
1713
1714 // boundary
1715 CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
1716 CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
1717 CHKERR add_disp_bc_ops();
1718
1719 // pipeline for mechanical
1721
1722 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
1723 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
1724 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
1725 tD_dev);
1726 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
1727 tD_dev);
1728 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1729 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD_dev);
1730 CHKERR add_domain_ops_lhs_constraints(
1731 pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1732 CHKERR add_domain_ops_rhs_constraints(
1733 pipeline_mng->getOpDomainRhsPipeline());
1734
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);
1741
1742 } else {
1743
1744 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
1745 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
1746 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(), tD);
1747 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(), tD);
1748 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD);
1749 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD);
1750
1751 // pipeline for constraints (only mu-part)
1752 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
1753 tD_dev);
1754 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
1755 tD_dev);
1756 CHKERR add_domain_ops_lhs_constraints(
1757 pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1758 CHKERR add_domain_ops_rhs_constraints(
1759 pipeline_mng->getOpDomainRhsPipeline());
1760 }
1761
1762 // contact integration volume
1763 add_domain_ops_contact_lhs(pipeline_mng->getOpDomainLhsPipeline());
1764 add_domain_ops_contact_rhs(pipeline_mng->getOpDomainRhsPipeline());
1765
1766 if (data.is_rotating && data.with_plasticity && false)
1767 add_skeleton_base_ops(pipeline_mng->getOpSkeletonRhsPipeline());
1768
1769 // contact integration on boundary
1770 add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
1771 add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
1772
1773 // setting integration rules
1774 auto integration_rule_bound = [](int, int, int approx_order) {
1775 return 2 * approx_order + 1;
1776 // return 3 * approx_order;
1777 };
1778
1779 auto integration_rule = [this](int, int, int approx_order) {
1781 return 2 * approx_order - 1;
1782 // return 3 * approx_order;
1783 return 2 * (approx_order - 1);
1784 };
1785
1786 auto integration_rule_lambda = [](int, int, int approx_order) {
1787 // return 3 * approx_order;
1788 return 1;
1789 };
1790
1793 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bound);
1794 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bound);
1795
1796 feLambdaLhs->getRuleHook = integration_rule_lambda;
1797 feLambdaRhs->getRuleHook = integration_rule_lambda;
1798
1799 if (data.with_contact) {
1800 // store data about the roller
1801 int def_val = 0;
1802 CHKERR mField.get_moab().tag_get_handle(
1803 "_ROLLER_NB", 1, MB_TYPE_INTEGER, commonDataPtr->rollerTag,
1804 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1805
1806 data.update_roller = boost::make_shared<BoundaryEle>(mField);
1807 auto &pipeline = data.update_roller->getOpPtrVector();
1808 if (data.is_ho_geometry) {
1809 pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
1810 pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
1811 }
1812 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1813
1814 pipeline.push_back(
1816 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1817 "SIGMA", commonDataPtr->contactTractionPtr));
1818 pipeline.push_back(
1820 if (data.debug_info)
1821 pipeline.push_back(
1823 pipeline.push_back(new OpGetContactArea("SIGMA", commonDataPtr));
1824 data.update_roller->getRuleHook = integration_rule_bound;
1825 }
1826
1827 auto create_reactions_element = [&]() {
1829 Range &reaction_ents = commonDataPtr->reactionEnts;
1830
1831 auto get_reactions_meshset_ents = [&]() {
1834 if (it->getMeshsetId() == data.reaction_id) {
1835 Range ents;
1836 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
1837 true);
1838 CHKERR mField.get_moab().get_connectivity(ents, reaction_ents, true);
1839 CHKERR mField.get_moab().get_adjacencies(
1840 ents, 1, false, reaction_ents, moab::Interface::UNION);
1841 reaction_ents.merge(ents);
1842 break;
1843 }
1844
1845 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1846 reaction_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));
1852 if (!ents_size /*&& data.reaction_id != -1 */) {
1854 "WORLD", Sev::warning,
1855 "Warning: Provided meshset (-reaction_id: %d) for calculating "
1856 "reactions is EMPTY! Calculating reactions for all FIX_ blocks.",
1858 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
1859 if (bc_mng->checkBlock(bc, "FIX_"))
1860 reaction_ents.merge(*bc.second->getBcEntsPtr());
1861 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1862 reaction_ents);
1863 }
1864
1866 };
1867
1868 CHKERR get_reactions_meshset_ents();
1869
1870 double defs[] = {0, 0, 0};
1871 CHKERR mField.get_moab().tag_get_handle("_REACTION_TAG", 3, MB_TYPE_DOUBLE,
1872 commonDataPtr->reactionTag,
1873 MB_TAG_CREAT | MB_TAG_SPARSE, defs);
1874
1875 // create calculate reactions element
1876 data.calc_reactions = boost::make_shared<DomainEle>(mField);
1877
1878 std::vector<int> ghosts{0, 1, 2};
1879 commonDataPtr->reactionsVec = createSmartGhostVector(
1880 mField.get_comm(), (mField.get_comm_rank() ? 0 : 3), 3,
1881 (mField.get_comm_rank() ? 3 : 0), &*ghosts.begin());
1882
1883 data.calc_reactions->getRuleHook = integration_rule;
1884
1885 auto &pipeline = data.calc_reactions->getOpPtrVector();
1886 pipeline.push_back(
1888 if (data.is_ho_geometry)
1889 CHKERR add_ho_ops(pipeline);
1890
1891 pipeline.push_back(
1893 if (data.is_large_strain) {
1894 pipeline.push_back(new OpLogStrain("U", commonDataPtr));
1895 } else
1896 pipeline.push_back(new OpStrain("U", commonDataPtr));
1897
1898 if (data.with_plasticity) {
1899 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1900 "EP", commonDataPtr->plasticStrainPtr));
1901 pipeline.push_back(
1903 } else
1904 pipeline.push_back(new OpStress("U", commonDataPtr, commonDataPtr->mtD));
1905
1906 // FIXME: include body forces
1907 // pipeline.push_back(new OpBodyForce("U", commonDataPtr, gravity));
1908 // if (data.is_rotating) {
1909 // pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
1910 // pipeline.push_back(new OpCentrifugalForce2("U", commonDataPtr));
1911 // }
1912
1913 if (data.is_large_strain) {
1914 if (data.is_neohooke)
1915 pipeline.push_back(
1917 else
1918 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1919 "U", commonDataPtr, commonDataPtr->mtD));
1920 pipeline.push_back(
1922 } else
1923 pipeline.push_back(
1925
1927 };
1928
1930 CHKERR create_reactions_element();
1931
1933}
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,...)
Definition: LogManager.hpp:311
@ ROW
Definition: definitions.h:123
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
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.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixTensorTimesGradU< 3 > OpMiXLambdaGradURhs
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< 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
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
DEPRECATED auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
struct OpSchurPreconditionMassContact< OpMassPrecHTensorBound > OpSchurPreconditionMassContactBound
struct OpSchurPreconditionMassContact< OpMassPrecHTensorVol > OpSchurPreconditionMassContactVol
static constexpr int approx_order
boost::shared_ptr< DomainEle > calc_reactions
boost::shared_ptr< BoundaryEle > update_roller
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
FieldApproximationBase base
boost::shared_ptr< DomainSideEle > volSideFe
boost::shared_ptr< DomainEle > feLambdaLhs
boost::shared_ptr< DomainEle > feLambdaRhs
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.
transform Hdiv base fluxes from reference element to physical triangle
Set indices on entities on finite element.
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
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)
[Common data]

◆ postProcess()

MoFEMErrorCode ContactPlasticity::postProcess ( )
private

Definition at line 2579 of file multifield_plasticity.cpp.

2579{ return 0; }

◆ runProblem()

MoFEMErrorCode ContactPlasticity::runProblem ( )

Definition at line 591 of file multifield_plasticity.cpp.

591 {
593 CHKERR setUP();
595 CHKERR bC();
596 CHKERR OPs();
597 CHKERR tsSolve();
601}
MoFEMErrorCode postProcess()
MoFEMErrorCode createCommonData()
MoFEMErrorCode checkResults()

◆ setUP()

MoFEMErrorCode ContactPlasticity::setUP ( )
private

Definition at line 834 of file multifield_plasticity.cpp.

834 {
837
838 auto get_base = [this](bool has_hexes = false) {
840
841 // Select base
842 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
843 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
844 enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
845 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz",
846 "bernstein"};
847 PetscInt choice_base_value = has_hexes ? DEMKOWICZ : AINSWORTH;
848 if (has_hexes)
849 choice_base_value = DEMKOWICZ;
850 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
851 LASBASETOPT, &choice_base_value, PETSC_NULL);
852
853 switch (choice_base_value) {
854 case AINSWORTH:
856 MOFEM_LOG("WORLD", Sev::inform)
857 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
858 break;
859 case DEMKOWICZ:
861 MOFEM_LOG("WORLD", Sev::inform)
862 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
863 break;
864 case BERNSTEIN:
866 MOFEM_LOG("WORLD", Sev::inform)
867 << "Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
868 break;
869 default:
870 base = LASTBASE;
871 break;
872 }
874 };
875
876 auto skin_ents = getEntsOnMeshSkin();
877 data.skinEnts = skin_ents;
878
879 Range tets, nodes;
880 int nb_tets, nb_hexes;
881 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
882 CHKERR mField.get_moab().get_connectivity(&*tets.begin(), 1, nodes);
883 CHKERR mField.get_moab().get_number_entities_by_type(0, MBTET, nb_tets, true);
884 CHKERR mField.get_moab().get_number_entities_by_type(0, MBHEX, nb_hexes,
885 true);
886
887 if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {
888 data.is_ho_geometry = PETSC_TRUE;
889 CHKERR simple->addDataField("MESH_NODE_POSITIONS", H1, base, 3);
890 CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2);
891 // auto all_edges_ptr = boost::make_shared<Range>();
892 // CHKERR mField.get_moab().get_entities_by_type(0, MBEDGE, *all_edges_ptr,
893 // true);
894 // CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2,
895 // all_edges_ptr.get());
896 }
897
898 CHKERR get_base(nb_hexes);
899 // Add field
900 CHKERR simple->addDomainField("U", H1, base, 3);
901 CHKERR simple->addBoundaryField("U", H1, base, 3);
902 CHKERR simple->setFieldOrder("U", data.order);
903
904 if (data.with_contact) {
905 CHKERR simple->addDomainField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
906 CHKERR simple->addBoundaryField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
907 CHKERR simple->setFieldOrder("SIGMA", 0);
908 CHKERR simple->setFieldOrder("SIGMA", data.order - 1, &skin_ents);
909 }
910 if (data.with_plasticity) {
912 CHKERR simple->addDomainField("TAU", space_test, USER_BASE, 1);
913 CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
914 CHKERR simple->addDomainField("EP", space_test, USER_BASE, 6);
915 CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
916
917 auto *pipeline_mng = mField.getInterface<PipelineManager>();
918
919 // this creates domain_elements
920 pipeline_mng->setDomainLhsIntegrationRule(
921 [](int, int, int approx_order) { return 1; });
922 pipeline_mng->setDomainRhsIntegrationRule(
923 [](int, int, int approx_order) { return 1; });
924
925 // get access to "TAU" data structure
926 for (string name : {"TAU", "EP"}) {
927
928 auto field_ptr = mField.get_field_structure(name);
929 // get table associating number of dofs to entities depending on
930 // approximation order set on those entities.
931 auto field_order_table =
932 const_cast<Field *>(field_ptr)->getFieldOrderTable();
933
934 // function set zero number of dofs
935 auto get_tau_bubble_order_zero = [](int p) { return 0; };
936 // function set non-zero number of dofs on tetrahedrons
937 auto get_tau_bubble_order_tet = [](int p) -> int {
938 return std::min(p + 1, 2);
939 };
940
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();
946 }
947
948 } else {
949 CHKERR simple->addDomainField("TAU", space_test, base, 1);
950 CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
951 CHKERR simple->addDomainField("EP", space_test, base, 6);
952 CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
953 }
954 }
955
956 if (data.is_rotating && data.with_plasticity && false) {
957 CHKERR simple->addSkeletonField("SKELETON_LAMBDA", L2, base, 7);
958 CHKERR simple->setFieldOrder("SKELETON_LAMBDA",
959 std::max(0, data.order - 1));
960 }
961
962 CHKERR simple->defineFiniteElements();
963
964 // Add Neumann forces elements
968
969 simple->getOtherFiniteElements().push_back("FORCE_FE");
970 simple->getOtherFiniteElements().push_back("PRESSURE_FE");
971 if (data.is_rotating && data.with_plasticity && false)
972 CHKERR simple->setSkeletonAdjacency();
973
974 // PetscBool is_partitioned = PETSC_TRUE;
975 CHKERR simple->defineProblem(PETSC_TRUE);
976 CHKERR simple->buildFields();
977
978 if (nb_hexes) {
979 CHKERR mField.set_field_order(0, MBQUAD, "U", 0);
980 CHKERR mField.set_field_order(0, MBHEX, "U", 0);
982 }
983 if (data.is_ho_geometry) {
984 Projection10NodeCoordsOnField ent_method_material(mField,
985 "MESH_NODE_POSITIONS");
986 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
987 }
988
989 CHKERR simple->buildFiniteElements();
990 if (data.is_rotating && data.with_plasticity && false) {
993 }
994 CHKERR simple->buildProblem();
996 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
997
999}
static Index< 'p', 3 > p
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
@ L2
field with C-1 continuity
Definition: definitions.h:88
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.
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.
FieldSpace space_test
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
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
Provide data structure for (tensor) field approximation.
Projection of edge entities with one mid-node on hierarchical basis.

◆ tsSolve()

MoFEMErrorCode ContactPlasticity::tsSolve ( )
private

Definition at line 1935 of file multifield_plasticity.cpp.

1935 {
1937
1940 ISManager *is_manager = mField.getInterface<ISManager>();
1941
1942 auto preProc =
1943 boost::make_shared<FePrePostProcess>(mField, bcData, commonDataPtr);
1944
1945 auto create_custom_ts = [&]() {
1946 auto set_dm_section = [&](auto dm) {
1948 PetscSection section;
1949 CHKERR mField.getInterface<ISManager>()->sectionCreate(
1950 simple->getProblemName(), &section);
1951 CHKERR DMSetDefaultSection(dm, section);
1952 CHKERR DMSetDefaultGlobalSection(dm, section);
1953 CHKERR PetscSectionDestroy(&section);
1955 };
1956 auto dm = simple->getDM();
1957
1958 CHKERR set_dm_section(dm);
1959 boost::shared_ptr<FEMethod> null;
1960 preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
1961 new TimeForceScale(data.move_history, false));
1962
1963 preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
1964 new TimeForceScale(data.contact_history, false));
1965 // // here we use accelerogram as position data
1966 // for (auto &roller : (*cache).rollerDataVec) {
1967 // if (!roller.positionDataParamName.empty())
1968 // preProc->methodsOpForRollersPosition.push_back(
1969 // boost::shared_ptr<MethodForForceScaling>(
1970 // new TimeAccelerogram(roller.positionDataParamName)));
1971 // if (!roller.directionDataParamName.empty())
1972 // preProc->methodsOpForRollersDirection.push_back(
1973 // boost::shared_ptr<MethodForForceScaling>(
1974 // new TimeAccelerogram(roller.directionDataParamName)));
1975 // }
1976
1977 // array<Range, 3> bc_ents;
1978 // for (auto &bdata : bcData)
1979 // bc_ents[bdata.comp].merge(bdata.ents);
1980
1981 // auto eraseRows = boost::make_shared<EraseRows>(
1982 // mField, simple->getProblemName(), bc_ents[0], bc_ents[1],
1983 // bc_ents[2]);
1984
1985 if (true) {
1986
1987 // Add element to calculate lhs of stiff part
1988 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null,
1989 dirichlet_bc_ptr, null);
1990 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, preProc,
1991 null);
1992
1993 auto postproc_bound_el = [&]() {
1995 if (auto bmc_ptr = block_mat_container_ptr.lock()) {
1996 auto &bmc = *bmc_ptr;
1997 bmc.clear();
1998 // print schur complement after assembling the boundary contributions
1999 // if (commonDataPtr->STauTau)
2000 // CHKERR MatView(commonDataPtr->STauTau,
2001 // PETSC_VIEWER_DRAW_WORLD);
2002 }
2004 };
2005
2007 pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
2008
2009 if (pipeline_mng->getBoundaryLhsFE())
2010 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getBoundaryFEName(),
2011 pipeline_mng->getBoundaryLhsFE(), null,
2012 pipeline_mng->getBoundaryLhsFE());
2013
2014 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(),
2015 pipeline_mng->getDomainLhsFE(), null, null);
2017 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feLambdaLhs,
2018 null, null);
2019
2020 if (pipeline_mng->getSkeletonLhsFE())
2021 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getSkeletonFEName(),
2022 pipeline_mng->getSkeletonLhsFE(), null,
2023 null);
2024
2025 // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
2026 // dirichlet_bc_ptr);
2027 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
2028 preProc);
2029
2030 if (pipeline_mng->getDomainRhsFE()) {
2031
2032 // add dirichlet boundary conditions
2033 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2034 dirichlet_bc_ptr, null);
2035 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2036 preProc, null);
2038 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
2039 feLambdaRhs, null, null);
2040
2041 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
2042 pipeline_mng->getDomainRhsFE(), null,
2043 null);
2044
2045 if (pipeline_mng->getBoundaryRhsFE())
2046 CHKERR DMMoFEMTSSetIFunction(dm, simple->getBoundaryFEName(),
2047 pipeline_mng->getBoundaryRhsFE(), null,
2048 null);
2049 if (pipeline_mng->getSkeletonRhsFE())
2050 CHKERR DMMoFEMTSSetIFunction(dm, simple->getSkeletonFEName(),
2051 pipeline_mng->getSkeletonRhsFE(), null,
2052 null);
2053 // Add surface forces
2054 for (auto fit = neumann_forces.begin(); fit != neumann_forces.end();
2055 fit++) {
2056 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2057 &fit->second->getLoopFe(), NULL, NULL);
2058 }
2059
2060 // Add edge forces
2061 for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
2062 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2063 &fit->second->getLoopFe(), NULL, NULL);
2064 }
2065
2066 // Add nodal forces
2067 for (auto fit = nodal_forces.begin(); fit != nodal_forces.end();
2068 fit++) {
2069 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2070 &fit->second->getLoopFe(), NULL, NULL);
2071 }
2072
2073 // CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2074 // null,
2075 // dirichlet_bc_ptr);
2076 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null, null,
2077 preProc);
2078 }
2079 }
2080
2081 auto print_active_set = [&](std::array<int, 2> &lgauss_pts, string name,
2084 std::vector<int> gauss_pts(2, 0);
2085 auto dm = mField.getInterface<Simple>()->getDM();
2086 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
2087 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2088
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;
2093 };
2094
2095 auto postproc_dom = [&]() {
2097 CHKERR print_active_set(commonDataPtr->stateArrayPlast, "plastic",
2098 mField);
2100 };
2101 auto postproc_bound = [&]() {
2103 CHKERR print_active_set(commonDataPtr->stateArrayCont, "contact", mField);
2105 };
2106
2107 if (data.debug_info) {
2109 pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
2110
2111 if (data.with_contact)
2112 pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
2113 }
2114
2115 auto ts = MoFEM::createTS(mField.get_comm());
2116 CHKERR TSSetDM(ts, dm);
2117 return ts;
2118 };
2119
2120 auto solver = create_custom_ts();
2121
2122 CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
2123
2124 auto dm = simple->getDM();
2125 auto D = smartCreateDMVector(dm);
2126
2127 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
2128 CHKERR TSSetFromOptions(solver);
2129 CHKERR TSSetSolution(solver, D);
2130
2131 CHKERR TSSetUp(solver);
2132
2133 // constexpr bool run_fieldsplit_for_boundary = true;
2134 if (true) {
2135 SNES snes;
2136 CHKERR TSGetSNES(solver, &snes);
2137 KSP ksp;
2138 CHKERR SNESGetKSP(snes, &ksp);
2139 PC pC;
2140 CHKERR KSPGetPC(ksp, &pC);
2141 PetscBool is_pcfs = PETSC_FALSE;
2142 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
2143 SmartPetscObj<AO> sigma_ao; // contact ordering
2144
2145 auto make_is_field_map = [&]() {
2146 PetscSection section;
2147 CHKERR DMGetDefaultSection(dm, &section);
2148
2149 int num_fields;
2150 CHKERR PetscSectionGetNumFields(section, &num_fields);
2151
2152 const MoFEM::Problem *prb_ptr;
2153 CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
2154
2155 map<std::string, SmartPetscObj<IS>> is_map;
2156 for (int ff = 0; ff != num_fields; ff++) {
2157
2158 const char *field_name;
2159 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
2160
2161 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
2162 prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
2163 is_map[field_name]);
2164 };
2165 return is_map;
2166 };
2167 auto is_map = make_is_field_map();
2168
2169 auto set_fieldsplit_on_bc = [&](PC &bc_pc, string name_prb) {
2171
2172 PetscBool is_bc_field_split;
2173 PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
2174 &is_bc_field_split);
2175
2176 auto block_prefix = simple->getProblemName();
2177 auto bc_mng = mField.getInterface<BcManager>();
2178
2179 auto is_all_bc =
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,
2182 1, is_all_bc);
2183 is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_X", "U", name_prb, 0,
2184 0, is_all_bc);
2185 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_X", "U", name_prb, 1,
2186 2, is_all_bc);
2187 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_Z", "U", name_prb, 0,
2188 1, is_all_bc);
2189 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_ALL", "U", name_prb,
2190 0, 2, is_all_bc);
2191 is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_ALL", "U", name_prb, 0,
2192 2, is_all_bc);
2193
2194 int is_all_bc_size;
2195 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
2196
2197 if (is_bc_field_split && is_all_bc_size) {
2198 MOFEM_LOG("WORLD", Sev::inform)
2199 << "Fieldsplit preconditioner for boundary block size: "
2200 << is_all_bc_size;
2201 CHKERR PCFieldSplitSetIS(bc_pc, NULL,
2202 is_all_bc); // boundary block
2203 //FIXME: it could be beneficial to have 'bc' and 'non_bc' in naming
2204
2205 // CHKERR PCFieldSplitSetIS(bc_pc, "bc",
2206 // is_all_bc); // boundary block
2207 // IS is_no_bc;
2208 // CHKERR PCFieldSplitGetISByIndex(bc_pc, 1, &is_no_bc);
2209 // CHKERR PCFieldSplitGetISByIndex(bc_pc, 0, &is_no_bc);
2210 // PetscBool is_equal;
2211 // CHKERR ISEqual(IS is1, IS is2, PetscBool * flg);
2212 // CHKERR PCFieldSplitSetIS(bc_pc, "non_bc",
2213 // is_no_bc); // non-boundary block
2214 // CHKERR ISView(is_no_bc, PETSC_VIEWER_STDOUT_SELF);
2215 // CHKERR ISDestroy(&is_no_bc);
2216 CHKERR PCSetUp(bc_pc);
2217 }
2218 // else {
2219 // CHKERR PCSetFromOptions(bc_pc);
2220 // CHKERR PCSetType(bc_pc, PCLU);
2221 // CHKERR PCSetUp(bc_pc);
2222 // }
2223
2225 };
2226
2227 auto set_global_mat_and_vec = [&]() {
2229 auto fe = mField.getInterface<PipelineManager>()->getDomainLhsFE();
2230 CHKERR DMCreateMatrix_MoFEM(dm, &fe->ts_B);
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);
2235 };
2236
2237 // Setup fieldsplit (block) solver - optional: yes/no
2238 if (is_pcfs == PETSC_TRUE && !data.is_fieldsplit_bc_only) {
2239
2240 CHKERR set_global_mat_and_vec();
2241
2242 PetscBool is_l0_field_split = PETSC_TRUE; // level0 fieldsplit
2243 PetscBool is_l1_field_split = PETSC_FALSE; // level1 fieldsplit
2244
2245 // PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l0_field_split);
2246
2247 auto set_fieldsplit_contact = [&]() {
2249
2250 const MoFEM::Problem *fs_sig_ptr;
2252 if (auto sig_data = fs_sig_ptr->subProblemData) {
2253
2254 is_map["E_IS_SIG"] = sig_data->rowIs;
2255 sigma_ao = sig_data->getSmartRowMap();
2256 // CHKERR AOView(sigma_ao, PETSC_VIEWER_STDOUT_SELF);
2257 CHKERR PCFieldSplitSetIS(pC, NULL, is_map["SIGMA"]);
2258 CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_SIG"]);
2259 // CHKERR PCFieldSplitSetType(pC, PC_COMPOSITE_ADDITIVE);
2260 // CHKERR PCFieldSplitSetType(pC,
2261 // PC_COMPOSITE_MULTIPLICATIVE);
2262 CHKERR PCSetUp(pC);
2263
2264 // SmartPetscObj<Mat> mat_test;
2265 // CHKERR DMCreateMatrix_MoFEM(dmContact, mat_test);
2266
2267 PetscInt n;
2268 KSP *ct_ksp;
2269 CHKERR PCFieldSplitGetSubKSP(pC, &n, &ct_ksp);
2270 PC l1_pc;
2271 CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
2272 PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
2273 &is_l1_field_split);
2274 // important!
2275 // in case of contact analysis we swap the first preconditioner
2276 pC = l1_pc;
2277 CHKERR PetscFree(ct_ksp);
2278 }
2280 };
2281
2282 //FIXME: contact with fieldsplit does not work in parallel, requires further investigation
2283 //TODO: will come back into this in the future
2284 if (data.with_contact && is_l0_field_split && false)
2285 CHKERR set_fieldsplit_contact();
2286
2287 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
2288 if (data.with_plasticity && is_l1_field_split) {
2289
2290 const MoFEM::Problem *schur_ep_ptr;
2291 CHKERR DMMoFEMGetProblemPtr(dmEP, &schur_ep_ptr);
2292 if (auto ep_data = schur_ep_ptr->subProblemData) {
2293
2294 is_map["E_IS_EP"] = ep_data->rowIs;
2295 // CHKERR ISView(is_map["EP"], PETSC_VIEWER_STDOUT_SELF);
2296 if (sigma_ao)
2297 CHKERR AOApplicationToPetscIS(sigma_ao, is_map["EP"]);
2298 // CHKERR ISSort(is_map["EP"]);
2299
2300 CHKERR PCFieldSplitSetIS(pC, "ep", is_map["EP"]);
2301 CHKERR PCFieldSplitSetIS(pC, "etau", is_map["E_IS_EP"]);
2302
2303 PCCompositeType pc_type;
2304 CHKERR PCFieldSplitGetType(pC, &pc_type);
2305
2306 if (pc_type == PC_COMPOSITE_SCHUR) {
2308 commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
2309 if (sigma_ao)
2310 commonDataPtr->aoSEpEp = sigma_ao;
2311
2312 commonDataPtr->blockMatContainerPtr =
2313 boost::make_shared<BlockMatContainer>();
2315 commonDataPtr->blockMatContainerPtr;
2316
2317 CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
2318 commonDataPtr->SEpEp);
2319
2320 CHKERR PCSetUp(pC);
2321 }
2322
2323 PetscInt n;
2324 KSP *tt_ksp;
2325 CHKERR PCFieldSplitGetSubKSP(pC, &n, &tt_ksp);
2326 PC tau_pc;
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);
2331
2332 if (is_tau_field_split) {
2333 const MoFEM::Problem *schur_tau_ptr;
2334 CHKERR DMMoFEMGetProblemPtr(dmTAU, &schur_tau_ptr);
2335 if (auto tau_data = schur_tau_ptr->subProblemData) {
2336
2337 // CHKERR tau_data->getRowIs(is_map["E_IS_TAU"]);
2338 is_map["E_IS_TAU"] = tau_data->rowIs;
2339
2340 AO ep_ao;
2341 CHKERR ep_data->getRowMap(&ep_ao);
2342
2343 if (sigma_ao)
2344 CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
2345 CHKERR AOApplicationToPetscIS(ep_ao, is_map["TAU"]);
2346
2347 CHKERR PCFieldSplitSetIS(tau_pc, "tau", is_map["TAU"]);
2348 CHKERR PCFieldSplitSetIS(tau_pc, "elastic", is_map["E_IS_TAU"]);
2349
2350 CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
2351 if (pc_type == PC_COMPOSITE_SCHUR) {
2353 commonDataPtr->aoSTauTau = tau_data->getSmartRowMap();
2354
2355 CHKERR PCFieldSplitSetSchurPre(tau_pc,
2356 PC_FIELDSPLIT_SCHUR_PRE_USER,
2357 commonDataPtr->STauTau);
2358 }
2359 CHKERR PCSetUp(tau_pc);
2360
2361 PetscInt bc_n;
2362 KSP *bc_ksp;
2363 CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
2364 PC bc_pc;
2365 CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
2366
2367 CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->getName());
2368 CHKERR PetscFree(bc_ksp);
2369 }
2370 }
2371 CHKERR PetscFree(tt_ksp);
2372 }
2373 }
2374 } else if (is_pcfs && data.is_fieldsplit_bc_only) {
2375
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());
2387
2388 }
2389 // else
2390 // {
2391 // SNES snes;
2392 // CHKERR TSGetSNES(solver, &snes);
2393 // KSP ksp;
2394 // CHKERR SNESGetKSP(snes, &ksp);
2395 // PC pCc;
2396 // CHKERR KSPGetPC(ksp, &pCc);
2397 // CHKERR PCSetFromOptions(pCc);
2398 // CHKERR PCSetType(pCc, PCLU);
2399 // }
2400 }
2401
2402 auto set_section_monitor = [&]() {
2404 SNES snes;
2405 CHKERR TSGetSNES(solver, &snes);
2406 PetscViewerAndFormat *vf;
2407 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2408 PETSC_VIEWER_DEFAULT, &vf);
2409 CHKERR SNESMonitorSet(
2410 snes,
2411 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
2412 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
2414 };
2415
2416 auto create_post_process_element = [&]() {
2418 postProcFe = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
2419 postProcFe->generateReferenceElementMesh();
2421 postProcFe->getUserPolynomialBase() =
2422 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
2423 postProcFe->getOpPtrVector().push_back(
2425
2426 postProcFe->getOpPtrVector().push_back(
2428 if (data.is_large_strain) {
2429 postProcFe->getOpPtrVector().push_back(
2430 new OpLogStrain("U", commonDataPtr));
2431 } else
2432 postProcFe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
2433 if (data.with_plasticity) {
2434 postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
2435 "TAU", commonDataPtr->plasticTauPtr));
2436 postProcFe->getOpPtrVector().push_back(
2438 "EP", commonDataPtr->plasticStrainPtr));
2439 postProcFe->getOpPtrVector().push_back(
2441 postProcFe->getOpPtrVector().push_back(
2443 } else
2444 postProcFe->getOpPtrVector().push_back(
2445 new OpStress("U", commonDataPtr, commonDataPtr->mtD));
2446
2447 if (data.is_large_strain) {
2448 if (data.is_neohooke)
2449 postProcFe->getOpPtrVector().push_back(
2451 else
2452 postProcFe->getOpPtrVector().push_back(
2454 commonDataPtr->mtD));
2455 }
2456 if (data.with_contact) {
2457 postProcFe->getOpPtrVector().push_back(
2459 "SIGMA", commonDataPtr->contactStressDivergencePtr));
2460 postProcFe->getOpPtrVector().push_back(
2462 "SIGMA", commonDataPtr->contactStressPtr));
2463
2464 postProcFe->getOpPtrVector().push_back(
2465 new OpPostProcContact("SIGMA", postProcFe->postProcMesh,
2466 postProcFe->mapGaussPts, commonDataPtr));
2467 }
2469 postProcFe->getOpPtrVector().push_back(
2470 new OpPostProcElastic<true>("U", postProcFe->postProcMesh,
2471 postProcFe->mapGaussPts, commonDataPtr));
2472 else
2473 postProcFe->getOpPtrVector().push_back(
2474 new OpPostProcElastic<false>("U", postProcFe->postProcMesh,
2475 postProcFe->mapGaussPts, commonDataPtr));
2476
2478 postProcFe->getOpPtrVector().push_back(
2479 new OpPostProcPlastic("U", postProcFe->postProcMesh,
2480 postProcFe->mapGaussPts, commonDataPtr));
2481
2483 // FIXME: use addFieldValuesPostProc("U", "REACTIONS", vec_rhs); instead
2484 postProcFe->getOpPtrVector().push_back(
2485 new OpSaveReactionForces(mField, "U", postProcFe->postProcMesh,
2486 postProcFe->mapGaussPts, commonDataPtr));
2487 postProcFe->addFieldValuesPostProc("U", "DISPLACEMENT");
2488 if (data.is_ho_geometry)
2489 postProcFe->addFieldValuesPostProc("MESH_NODE_POSITIONS");
2491 };
2492
2493 auto scatter_create = [&](auto coeff) {
2495 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
2496 ROW, "U", coeff, coeff, is);
2497 int loc_size;
2498 CHKERR ISGetLocalSize(is, &loc_size);
2499 Vec v;
2500 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
2501 VecScatter scatter;
2502 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
2503 return std::make_tuple(SmartPetscObj<Vec>(v),
2504 SmartPetscObj<VecScatter>(scatter));
2505 };
2506
2507 auto set_time_monitor = [&]() {
2509 boost::shared_ptr<MMonitor> monitor_ptr(
2510 new MMonitor(dm, mField, commonDataPtr, postProcFe, uXScatter,
2512 boost::shared_ptr<ForcesAndSourcesCore> null;
2513 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
2514 monitor_ptr, null, null);
2516 };
2517
2518 CHKERR set_section_monitor();
2519 CHKERR create_post_process_element();
2520
2521 if (data.is_restart) {
2522
2523 CHKERR PetscPrintf(PETSC_COMM_WORLD,
2524 "Reading vector in binary from %s file...\n",
2525 data.restart_file_str.c_str());
2526 PetscViewer viewer;
2527 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
2528 data.restart_file_str.c_str(), FILE_MODE_READ,
2529 &viewer);
2530 CHKERR VecLoad(D, viewer);
2531
2532 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
2534 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
2535
2536 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
2537 string s = data.restart_file_str;
2538 Tokenizer tok{s, boost::char_separator<char>("_")};
2539 auto it = tok.begin();
2540 const int restart_step = std::stoi((++it)->c_str());
2541 data.restart_step = restart_step;
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);
2547 }
2548
2549 uXScatter = scatter_create(0);
2550 uYScatter = scatter_create(1);
2551 uZScatter = scatter_create(2);
2552 CHKERR set_time_monitor();
2553
2554 // Add iteration post-proc for debugging
2555 // boost::shared_ptr<FEMethod> null_fe;
2556 // SNES snes;
2557 // CHKERR TSGetSNES(solver, &snes);
2558 // auto write_fe = boost::make_shared<FEMethod>();
2559 // write_fe->postProcessHook = [&]() {
2560 // MoFEMFunctionBegin;
2561 // int iter;
2562 // CHKERR SNESGetIterationNumber(snes, &iter);
2563 // CHKERR postProcFe->writeFile(
2564 // "it_out_plastic_" + boost::lexical_cast<std::string>(iter) + ".h5m");
2565 // MoFEMFunctionReturn(0);
2566 // };
2567 // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), postProcFe,
2568 // null_fe, write_fe);
2569
2570 CHKERR TSSolve(solver, D);
2571
2572 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
2573 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
2574 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
2575
2577}
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
FTensor::Index< 'n', SPACE_DIM > n
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
Definition: DMMoFEM.cpp:786
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
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
Definition: DMMoFEM.cpp:839
double D
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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.
Definition: DMMoFEM.cpp:1042
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1011
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
auto createTS(MPI_Comm comm)
constexpr auto field_name
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProcFe
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.
Definition: BcManager.cpp:242
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
intrusive_ptr for managing petsc objects

Member Data Documentation

◆ base

FieldApproximationBase ContactPlasticity::base

Definition at line 157 of file multifield_plasticity.cpp.

◆ bcData

boost::ptr_vector<DataFromMove> ContactPlasticity::bcData
private

Definition at line 141 of file multifield_plasticity.cpp.

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > ContactPlasticity::boundaryMarker
private

Definition at line 139 of file multifield_plasticity.cpp.

◆ commonDataPtr

boost::shared_ptr<OpContactTools::CommonData> ContactPlasticity::commonDataPtr
private

Definition at line 127 of file multifield_plasticity.cpp.

◆ data

ContactPlasticity::ProblemData ContactPlasticity::data
static

Definition at line 244 of file multifield_plasticity.cpp.

◆ dirichlet_bc_ptr

boost::shared_ptr<DirichletDisplacementRemoveDofsBc> ContactPlasticity::dirichlet_bc_ptr
private

Definition at line 137 of file multifield_plasticity.cpp.

◆ dmContact

SmartPetscObj<DM> ContactPlasticity::dmContact
private

Definition at line 154 of file multifield_plasticity.cpp.

◆ dmEP

SmartPetscObj<DM> ContactPlasticity::dmEP
private

Definition at line 152 of file multifield_plasticity.cpp.

◆ dmTAU

SmartPetscObj<DM> ContactPlasticity::dmTAU
private

Definition at line 153 of file multifield_plasticity.cpp.

◆ edge_forces

boost::ptr_map<std::string, EdgeForce> ContactPlasticity::edge_forces
private

Definition at line 136 of file multifield_plasticity.cpp.

◆ feLambdaLhs

boost::shared_ptr<DomainEle> ContactPlasticity::feLambdaLhs
private

Definition at line 147 of file multifield_plasticity.cpp.

◆ feLambdaRhs

boost::shared_ptr<DomainEle> ContactPlasticity::feLambdaRhs
private

Definition at line 146 of file multifield_plasticity.cpp.

◆ invJac

MatrixDouble ContactPlasticity::invJac
private

Definition at line 126 of file multifield_plasticity.cpp.

◆ jAc

MatrixDouble ContactPlasticity::jAc
private

Definition at line 126 of file multifield_plasticity.cpp.

◆ mField

MoFEM::Interface& ContactPlasticity::mField
private

Definition at line 117 of file multifield_plasticity.cpp.

◆ neumann_forces

boost::ptr_map<std::string, NeumannForcesSurface> ContactPlasticity::neumann_forces
private

Definition at line 134 of file multifield_plasticity.cpp.

◆ nodal_forces

boost::ptr_map<std::string, NodalForce> ContactPlasticity::nodal_forces
private

Definition at line 135 of file multifield_plasticity.cpp.

◆ postProcFe

boost::shared_ptr<PostProcVolumeOnRefinedMesh> ContactPlasticity::postProcFe
private

Definition at line 128 of file multifield_plasticity.cpp.

◆ uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uXScatter
private

Definition at line 129 of file multifield_plasticity.cpp.

◆ uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uYScatter
private

Definition at line 130 of file multifield_plasticity.cpp.

◆ uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > ContactPlasticity::uZScatter
private

Definition at line 131 of file multifield_plasticity.cpp.

◆ volSideFe

boost::shared_ptr<DomainSideEle> ContactPlasticity::volSideFe
private

Definition at line 142 of file multifield_plasticity.cpp.


The documentation for this struct was generated from the following file: