15 using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
21 using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
31 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
32 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
33 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
38 CHKERR DMoFEMGetInterfacePtr(
dM, &m_field_ptr);
40 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
42 auto push_domain_ops = [&](
auto &pip) {
44 pip, {
H1,
HDIV},
"GEOMETRY")),
45 "Apply base transform");
46 auto henky_common_data_ptr =
47 commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
48 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform);
49 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
50 pip.push_back(
new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
51 "SIGMA", contact_stress_ptr));
52 return std::make_tuple(henky_common_data_ptr, contact_stress_ptr);
55 auto push_bdy_ops = [&](
auto &pip,
int space_dim) {
58 pip, {
HDIV},
"GEOMETRY")),
61 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
62 pip.push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
63 "SIGMA", common_data_ptr->contactTractionPtr()));
64 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
65 "U", common_data_ptr->contactDispPtr()));
66 return common_data_ptr;
68 return boost::shared_ptr<ContactOps::CommonData>();
71 auto get_domain_pip = [&](
auto &pip)
72 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
74 auto op_loop_side =
new OpLoopSide<SideEle>(
75 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
77 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
78 pip.push_back(op_loop_side);
79 return op_loop_side->getOpPtrVector();
85 auto get_post_proc_domain_fe = [&]() {
87 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
88 auto &pip = post_proc_fe->getOpPtrVector();
90 auto common_data_ptr = push_bdy_ops(pip,
SPACE_DIM - 1);
91 auto [henky_common_data_ptr, contact_stress_ptr] =
92 push_domain_ops(get_domain_pip(pip));
94 auto u_ptr = boost::make_shared<MatrixDouble>();
95 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_ptr));
96 auto X_ptr = boost::make_shared<MatrixDouble>();
98 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
104 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
116 (common_data_ptr) ? common_data_ptr->contactTractionPtr()
123 {
"SIGMA", contact_stress_ptr},
125 {
"G", henky_common_data_ptr->matGradPtr},
127 {
"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
140 auto get_post_proc_bdy_fe = [&]() {
142 boost::make_shared<PostProcEleBdy>(*m_field_ptr,
postProcMesh);
143 auto &pip = post_proc_fe->getOpPtrVector();
145 auto common_data_ptr = push_bdy_ops(pip,
SPACE_DIM);
147 auto op_loop_side =
new OpLoopSide<SideEle>(
148 *m_field_ptr,
"dFE",
SPACE_DIM, Sev::noisy,
150 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
151 pip.push_back(op_loop_side);
153 auto [henky_common_data_ptr, contact_stress_ptr] =
154 push_domain_ops(op_loop_side->getOpPtrVector());
156 auto X_ptr = boost::make_shared<MatrixDouble>();
158 new OpCalculateVectorFieldValues<SPACE_DIM>(
"GEOMETRY", X_ptr));
164 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
168 {{
"U", common_data_ptr->contactDispPtr()},
170 {
"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()}},
183 auto get_integrate_traction = [&]() {
184 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
185 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
187 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
188 integrate_traction->getOpPtrVector(), {HDIV},
"GEOMETRY")),
192 integrate_traction->getOpPtrVector().push_back(
193 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
199 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
200 integrate_traction->getOpPtrVector(),
"SIGMA")),
201 "push operators to calculate traction");
203 return integrate_traction;
206 postProcDomainFe = get_post_proc_domain_fe();
208 postProcBdyFe = get_post_proc_bdy_fe();
209 integrateTraction = get_integrate_traction();
218 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
220 auto post_proc = [&]() {
223 auto post_proc_begin =
224 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
226 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
227 *m_field_ptr, postProcMesh);
229 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
230 if (!postProcBdyFe) {
231 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcDomainFe);
233 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", postProcDomainFe);
234 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcBdyFe);
236 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
238 CHKERR post_proc_end->writeFile(
239 "out_contact_" + boost::lexical_cast<std::string>(sTEP) +
".h5m");
243 auto calculate_traction = [&] {
245 CHKERR VecZeroEntries(CommonData::totalTraction);
246 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", integrateTraction);
247 CHKERR VecAssemblyBegin(CommonData::totalTraction);
248 CHKERR VecAssemblyEnd(CommonData::totalTraction);
252 auto calculate_reactions = [&]() {
255 auto res = createDMVector(dM);
257 auto assemble_domain = [&]() {
259 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
260 auto &pip = fe_rhs->getOpPtrVector();
270 auto mat_acceleration = boost::make_shared<MatrixDouble>();
271 pip.push_back(
new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
272 "U", mat_acceleration));
275 [](
double,
double,
double) {
return rho; }));
278 auto mat_velocity = boost::make_shared<MatrixDouble>();
279 pip.push_back(
new OpCalculateVectorFieldValuesDot<SPACE_DIM>(
282 new OpInertiaForce(
"U", mat_velocity, [](
double,
double,
double) {
287 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
288 *m_field_ptr, pip,
"U",
"MAT_ELASTIC", Sev::inform);
289 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", fe_rhs);
293 auto assemble_boundary = [&]() {
295 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
296 auto &pip = fe_rhs->getOpPtrVector();
304 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {
HDIV},
308 pip.push_back(
new OpSetHOWeightsOnSubDim<SPACE_DIM>());
310 auto u_disp = boost::make_shared<MatrixDouble>();
311 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_disp));
313 new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
322 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
323 auto get_post_proc_hook_rhs = [
this, fe_post_proc_ptr, res,
326 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
327 *m_field_ptr, fe_post_proc_ptr, res)();
330 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
331 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
336 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
338 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
339 INSERT_VALUES, SCATTER_FORWARD);
340 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
341 INSERT_VALUES, SCATTER_FORWARD);
343 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
344 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
345 MOFEM_LOG_C(
"CONTACT", Sev::inform,
"%s time %3.4e min %3.4e max %3.4e",
346 msg.c_str(), ts_t, min, max);
350 auto print_traction = [&](
const std::string msg) {
353 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
356 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
357 MOFEM_LOG_C(
"CONTACT", Sev::inform,
"%s time %3.4e %3.4e %3.4e %3.4e",
358 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
359 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
365 <<
"Write file at time " << ts_t <<
" write step " << sTEP;
368 CHKERR calculate_traction();
369 CHKERR calculate_reactions();
371 CHKERR print_max_min(uXScatter,
"Ux");
372 CHKERR print_max_min(uYScatter,
"Uy");
374 CHKERR print_max_min(uZScatter,
"Uz");
375 CHKERR print_traction(
"Contact force");
383 SmartPetscObj<DM>
dM;
384 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uXScatter;
385 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uYScatter;
386 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uZScatter;
388 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
#define MOFEM_LOG_C(channel, severity, format,...)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int geom_order
Order if fixed.
static constexpr int approx_order
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomain
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleDomain
virtual int get_comm_rank() const =0
Deprecated interface functions.
Post post-proc data at points from hash maps.