v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
ContactOps::Monitor Struct Reference

#include <users_modules/tutorials/adv-1/src/PostProcContact.hpp>

Inheritance diagram for ContactOps::Monitor:
[legend]
Collaboration diagram for ContactOps::Monitor:
[legend]

Public Member Functions

 Monitor (SmartPetscObj< DM > &dm, bool use_mfront=false, boost::shared_ptr< GenericElementInterface > mfront_interface=nullptr, bool is_axisymmetric=false, int atom_test=0)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode operator() ()
 
MoFEMErrorCode postProcess ()
 
MoFEMErrorCode setScatterVectors (std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter)
 
 Monitor (SmartPetscObj< DM > &dm, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode operator() ()
 
MoFEMErrorCode postProcess ()
 

Private Attributes

SmartPetscObj< DM > dM
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
boost::shared_ptr< moab::Core > postProcMesh = boost::make_shared<moab::Core>()
 
boost::shared_ptr< PostProcEleDomainpostProcDomainFe
 
boost::shared_ptr< PostProcEleBdypostProcBdyFe
 
boost::shared_ptr< BoundaryEleintegrateTraction
 
moab::Core mbVertexPostproc
 
moab::Interface & moabVertex
 
double lastTime
 
double deltaTime
 
int sTEP
 
bool useMFront
 
bool isAxisymmetric
 
boost::shared_ptr< GenericElementInterfacemfrontInterface
 
int atomTest
 

Detailed Description

Definition at line 28 of file PostProcContact.hpp.

Constructor & Destructor Documentation

◆ Monitor() [1/2]

ContactOps::Monitor::Monitor ( SmartPetscObj< DM > &  dm,
bool  use_mfront = false,
boost::shared_ptr< GenericElementInterface mfront_interface = nullptr,
bool  is_axisymmetric = false,
int  atom_test = 0 
)
inline

Definition at line 30 of file PostProcContact.hpp.

36
37 MoFEM::Interface *m_field_ptr;
38 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
39
40 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
41
42 struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
43 OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
44 : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
45 mPtr(m_ptr), scale(s) {}
46 MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
47 *mPtr *= 1./scale;
48 return 0;
49 }
50
51 private:
52 boost::shared_ptr<MatrixDouble> mPtr;
53 double scale;
54 };
55
56 auto push_domain_ops = [&](auto &pip) {
57 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
58 pip, {H1, HDIV}, "GEOMETRY")),
59 "Apply base transform");
60 auto henky_common_data_ptr =
61 commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
62 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
63 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
64 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
65 "SIGMA", contact_stress_ptr));
66 pip.push_back(new OpScale(contact_stress_ptr, scale));
67 return std::make_tuple(henky_common_data_ptr, contact_stress_ptr);
68 };
69
70 auto push_bdy_ops = [&](auto &pip) {
71 // evaluate traction
72 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
73 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
74 "U", common_data_ptr->contactDispPtr()));
75 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
76 "SIGMA", common_data_ptr->contactTractionPtr()));
77 pip.push_back(new OpScale(common_data_ptr->contactTractionPtr(), scale));
78 using C = ContactIntegrators<BoundaryEleOp>;
79 pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
80 common_data_ptr));
81 return common_data_ptr;
82 };
83
84 auto get_domain_pip = [&](auto &pip)
85 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
86 if constexpr (SPACE_DIM == 3) {
87 auto op_loop_side = new OpLoopSide<SideEle>(
88 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
89 boost::make_shared<
90 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
91 pip.push_back(op_loop_side);
92 return op_loop_side->getOpPtrVector();
93 } else {
94 return pip;
95 }
96 };
97
98 auto get_post_proc_domain_fe = [&]() {
99 auto post_proc_fe =
100 boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
101 auto &pip = post_proc_fe->getOpPtrVector();
102
103 auto [henky_common_data_ptr, contact_stress_ptr] =
104 push_domain_ops(get_domain_pip(pip));
105
106 auto u_ptr = boost::make_shared<MatrixDouble>();
107 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
108 auto X_ptr = boost::make_shared<MatrixDouble>();
109 pip.push_back(
110 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
111
112
113
114 pip.push_back(
115
116 new OpPPMap(
117
118 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
119
120 {},
121 {
122
123 {"U", u_ptr}, {"GEOMETRY", X_ptr}
124
125 },
126 {
127
128 {"SIGMA", contact_stress_ptr},
129
130 {"G", henky_common_data_ptr->matGradPtr},
131
132 {"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
133
134 },
135 {}
136
137 )
138
139 );
140
141 if (SPACE_DIM == 3) {
142
143 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
144 pip, {HDIV}, "GEOMETRY")),
145 "Apply transform");
146 auto common_data_ptr = push_bdy_ops(pip);
147
148 pip.push_back(
149
150 new OpPPMap(
151
152 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
153
154 {{"SDF", common_data_ptr->sdfPtr()},
155 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
156
157 {
158
159 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
160 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
161
162 },
163
164 {},
165
166 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
167
168 )
169
170 );
171 }
172
173 return post_proc_fe;
174 };
175
176 auto get_post_proc_bdy_fe = [&]() {
177 auto post_proc_fe =
178 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
179 auto &pip = post_proc_fe->getOpPtrVector();
180
181 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
182 pip, {HDIV}, "GEOMETRY")),
183 "Apply transform");
184 auto common_data_ptr = push_bdy_ops(pip);
185
186 // create OP which run element on side
187 auto op_loop_side = new OpLoopSide<SideEle>(
188 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
189 boost::make_shared<
190 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
191 pip.push_back(op_loop_side);
192
193 auto [henky_common_data_ptr, contact_stress_ptr] =
194 push_domain_ops(op_loop_side->getOpPtrVector());
195
196 auto X_ptr = boost::make_shared<MatrixDouble>();
197 pip.push_back(
198 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
199
200 pip.push_back(
201
202 new OpPPMap(
203
204 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
205
206 {{"SDF", common_data_ptr->sdfPtr()},
207 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
208
209 {{"U", common_data_ptr->contactDispPtr()},
210 {"GEOMETRY", X_ptr},
211 {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
212 {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
213
214 },
215
216 {},
217
218 {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
219
220 )
221
222 );
223
224 return post_proc_fe;
225 };
226
227 auto get_integrate_traction = [&]() {
228 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
229 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
231 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
232 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
233 "Apply transform");
234 // We have to integrate on curved face geometry, thus integration weight
235 // have to adjusted.
236 integrate_traction->getOpPtrVector().push_back(
237 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
238 integrate_traction->getRuleHook = [](int, int, int approx_order) {
239 return 2 * approx_order + geom_order - 1;
240 };
241
243 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
244 integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
245 "push operators to calculate traction");
246
247 return integrate_traction;
248 };
249
250 postProcDomainFe = get_post_proc_domain_fe();
251 if constexpr (SPACE_DIM == 2)
252 postProcBdyFe = get_post_proc_bdy_fe();
253 integrateTraction = get_integrate_traction();
254 }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
int atom_test
Definition: contact.cpp:139
PetscBool use_mfront
Definition: contact.cpp:136
PetscBool is_axisymmetric
Definition: contact.cpp:137
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:400
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition: plastic.cpp:170
int geom_order
Order if fixed.
Definition: plastic.cpp:188
static constexpr int approx_order
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
SmartPetscObj< DM > dM
boost::shared_ptr< GenericElementInterface > mfrontInterface
boost::shared_ptr< BoundaryEle > integrateTraction
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
moab::Interface & moabVertex
boost::shared_ptr< moab::Core > postProcMesh
Deprecated interface functions.
Post post-proc data at points from hash maps.

◆ Monitor() [2/2]

ContactOps::Monitor::Monitor ( SmartPetscObj< DM > &  dm,
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  ux_scatter,
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  uy_scatter,
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  uz_scatter 
)
inline

Definition at line 30 of file PostProcContact.hpp.

34 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
35 uZScatter(uz_scatter), moabVertex(mbVertexPostproc), sTEP(0) {
36
37 MoFEM::Interface *m_field_ptr;
38 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
39
40 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
41
42 auto push_domain_ops = [&](auto &pip) {
43 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
53 };
54
55 auto push_bdy_ops_traction = [&](auto &pip) {
56 // evaluate traction
57 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
58 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
59 "SIGMA", common_data_ptr->contactTractionPtr()));
60 return common_data_ptr;
61 };
62
63 auto push_bdy_ops_sdf = [&](auto &pip) {
64 // evaluate traction
65 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
66 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
67 "U", common_data_ptr->contactDispPtr()));
68 using C = ContactIntegrators<BoundaryEleOp>;
69 pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
70 common_data_ptr));
71 return common_data_ptr;
72 };
73
74 auto get_domain_pip = [&](auto &pip)
75 -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
76 if constexpr (SPACE_DIM == 3) {
77 auto op_loop_side = new OpLoopSide<SideEle>(
78 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
79 boost::make_shared<
80 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
81 pip.push_back(op_loop_side);
82 return op_loop_side->getOpPtrVector();
83 } else {
84 return pip;
85 }
86 };
87
88 auto get_post_proc_domain_fe = [&]() {
89 auto post_proc_fe =
90 boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
91 auto &pip = post_proc_fe->getOpPtrVector();
92
93 auto [henky_common_data_ptr, contact_stress_ptr] =
94 push_domain_ops(get_domain_pip(pip));
95
96 auto u_ptr = boost::make_shared<MatrixDouble>();
97 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
98 auto X_ptr = boost::make_shared<MatrixDouble>();
99 pip.push_back(
100 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
101
102 pip.push_back(
103
104 new OpPPMap(
105
106 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
107
108 {},
109 {
110
111 {"U", u_ptr}, {"GEOMETRY", X_ptr}
112
113 },
114 {
115
116 {"SIGMA", contact_stress_ptr},
117
118 {"G", henky_common_data_ptr->matGradPtr},
119
120 {"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
121
122 },
123 {}
124
125 )
126
127 );
128
129 if (SPACE_DIM == 3) {
130
131 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
132 pip, {HDIV}, "GEOMETRY")),
133 "Apply transform");
134 auto common_data_traction_ptr = push_bdy_ops_traction(pip);
135 auto common_data_sdf_ptr = push_bdy_ops_sdf(pip);
136
137 pip.push_back(
138
139 new OpPPMap(
140
141 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
142
143 {{"SDF", common_data_sdf_ptr->sdfPtr()}},
144
145 {
146
147 {"TRACTION_CONTACT",
148 common_data_traction_ptr->contactTractionPtr()},
149 {"GRAD_SDF", common_data_sdf_ptr->gradSdfPtr()}
150
151 },
152
153 {},
154
155 {{"HESS_SDF", common_data_sdf_ptr->hessSdfPtr()}}
156
157 )
158
159 );
160 }
161
162 return post_proc_fe;
163 };
164
165 auto get_post_proc_bdy_fe = [&]() {
166 auto post_proc_fe =
167 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
168 auto &pip = post_proc_fe->getOpPtrVector();
169
170 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
171 pip, {HDIV}, "GEOMETRY")),
172 "Apply transform");
173 auto common_data_traction_ptr = push_bdy_ops_traction(pip);
174 auto common_data_sdf_ptr = push_bdy_ops_sdf(pip);
175
176 // create OP which run element on side
177 auto op_loop_side = new OpLoopSide<SideEle>(
178 *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
179 boost::make_shared<
180 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
181 pip.push_back(op_loop_side);
182
183 auto [henky_common_data_ptr, contact_stress_ptr] =
184 push_domain_ops(op_loop_side->getOpPtrVector());
185
186 auto X_ptr = boost::make_shared<MatrixDouble>();
187 pip.push_back(
188 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
189
190 pip.push_back(
191
192 new OpPPMap(
193
194 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
195
196 {{"SDF", common_data_sdf_ptr->sdfPtr()}},
197
198 {{"U", common_data_sdf_ptr->contactDispPtr()},
199 {"GEOMETRY", X_ptr},
200 {"TRACTION_CONTACT",
201 common_data_traction_ptr->contactTractionPtr()},
202 {"GRAD_SDF", common_data_sdf_ptr->gradSdfPtr()}
203
204 },
205
206 {},
207
208 {{"HESS_SDF", common_data_sdf_ptr->hessSdfPtr()}}
209
210 )
211
212 );
213
214 return post_proc_fe;
215 };
216
217 auto get_integrate_traction = [&]() {
218 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
219 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
221 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
222 integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
223 "Apply transform");
224 // We have to integrate on curved face geometry, thus integration weight
225 // have to adjusted.
226 integrate_traction->getOpPtrVector().push_back(
227 new OpSetHOWeightsOnSubDim<SPACE_DIM>());
228 integrate_traction->getRuleHook = [](int, int, int approx_order) {
229 return 2 * approx_order + geom_order - 1;
230 };
231
233 (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
234 integrate_traction->getOpPtrVector(), "SIGMA")),
235 "push operators to calculate traction");
236
237 return integrate_traction;
238 };
239
240 postProcDomainFe = get_post_proc_domain_fe();
241 if constexpr (SPACE_DIM == 2)
242 postProcBdyFe = get_post_proc_bdy_fe();
243 integrateTraction = get_integrate_traction();
244 }
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter

Member Function Documentation

◆ operator()() [1/2]

MoFEMErrorCode ContactOps::Monitor::operator() ( )
inline

Definition at line 257 of file PostProcContact.hpp.

257{ return 0; }

◆ operator()() [2/2]

MoFEMErrorCode ContactOps::Monitor::operator() ( )
inline

Definition at line 247 of file PostProcContact.hpp.

247{ return 0; }

◆ postProcess() [1/2]

MoFEMErrorCode ContactOps::Monitor::postProcess ( )
inline

Definition at line 259 of file PostProcContact.hpp.

259 {
261 MoFEM::Interface *m_field_ptr;
262 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
263
264 auto post_proc = [&]() {
266
267 if (!useMFront) {
268 auto post_proc_begin =
269 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
271 auto post_proc_end =
272 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
274
276 post_proc_begin->getFEMethod());
277 if (!postProcBdyFe) {
278 postProcDomainFe->copyTs(*this); // this here is a Monitor
280 } else {
281 postProcDomainFe->copyTs(*this); // this here is a Monitor
282 postProcBdyFe->copyTs(*this);
285 }
287 post_proc_end->getFEMethod());
288
289 CHKERR post_proc_end->writeFile(
290 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
291 } else {
292 mfrontInterface->updateElementVariables();
293 mfrontInterface->postProcessElement(ts_step);
294 }
295
297 };
298
299 auto calculate_traction = [&] {
301 CHKERR VecZeroEntries(CommonData::totalTraction);
303 CHKERR VecAssemblyBegin(CommonData::totalTraction);
304 CHKERR VecAssemblyEnd(CommonData::totalTraction);
306 };
307
308 auto calculate_reactions = [&]() {
310
311 auto res = createDMVector(dM);
312
313 auto assemble_domain = [&]() {
315 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
316 auto &pip = fe_rhs->getOpPtrVector();
317 fe_rhs->f = res;
318
319 auto integration_rule = [](int, int, int approx_order) {
320 return 2 * approx_order + geom_order - 1;
321 };
322 fe_rhs->getRuleHook = integration_rule;
323 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
324 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
325 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
327 };
328
329 auto assemble_boundary = [&]() {
331 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
332 auto &pip = fe_rhs->getOpPtrVector();
333 fe_rhs->f = res;
334
335 auto integration_rule = [](int, int, int approx_order) {
336 return 2 * approx_order + geom_order - 1;
337 };
338 fe_rhs->getRuleHook = integration_rule;
339
340 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
341 "GEOMETRY");
342 // We have to integrate on curved face geometry, thus integration weight
343 // have to adjusted.
344 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
345
346 auto u_disp = boost::make_shared<MatrixDouble>();
347 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
348 pip.push_back(
349 new OpSpringRhs("U", u_disp, [this](double, double, double) {
350 return spring_stiffness;
351 }));
352
353 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
354
356 };
357
358 CHKERR assemble_domain();
359 CHKERR assemble_boundary();
360
361 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
362 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
363 m_field_ptr]() {
365 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
366 *m_field_ptr, fe_post_proc_ptr, res)();
368 };
369 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
370 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
371
373 };
374
375 auto print_max_min = [&](auto &tuple, const std::string msg) {
377 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
378 INSERT_VALUES, SCATTER_FORWARD);
379 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
380 INSERT_VALUES, SCATTER_FORWARD);
381 double max, min;
382 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
383 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
384 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e min %3.4e max %3.4e",
385 msg.c_str(), ts_t, min, max);
387 };
388
389 auto print_traction = [&](const std::string msg) {
391 MoFEM::Interface *m_field_ptr;
392 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
393 if (!m_field_ptr->get_comm_rank()) {
394 const double *t_ptr;
395 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
396 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
397 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
398 if (atomTest == 1 && fabs(ts_t - 1.0) < 1e-3) {
399 double hertz_tract = 158.73;
400 double tol = 4e-3;
401 if (fabs(t_ptr[1] - hertz_tract) / hertz_tract > tol) {
402 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
403 "atom test %d diverged! %3.4e != %3.4e", atom_test,
404 t_ptr[1], hertz_tract);
405 }
406 }
407 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
408 }
410 };
411
412 int se = 1;
413 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every", &se, PETSC_NULL);
414
415 if (!(ts_step % se)) {
416 MOFEM_LOG("CONTACT", Sev::inform)
417 << "Write file at time " << ts_t << " write step " << sTEP;
418 CHKERR post_proc();
419 }
420 CHKERR calculate_traction();
421 CHKERR calculate_reactions();
422
423 CHKERR print_max_min(uXScatter, "Ux");
424 CHKERR print_max_min(uYScatter, "Uy");
425 if (SPACE_DIM == 3)
426 CHKERR print_max_min(uZScatter, "Uz");
427 CHKERR print_traction("Contact force");
428
429 ++sTEP;
430
432 }
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
double spring_stiffness
Definition: contact.cpp:130
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
auto integration_rule
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
double tol
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
virtual int get_comm_rank() const =0

◆ postProcess() [2/2]

MoFEMErrorCode ContactOps::Monitor::postProcess ( )
inline

Definition at line 249 of file PostProcContact.hpp.

249 {
251 MoFEM::Interface *m_field_ptr;
252 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
253
254 auto post_proc = [&]() {
256
257 auto post_proc_begin =
258 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
260 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
261 *m_field_ptr, postProcMesh);
262
263 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
264 if (!postProcBdyFe) {
265 postProcDomainFe->ts_t = this->ts_t; // this here is a Monitor
267 } else {
268 postProcDomainFe->ts_t = this->ts_t; // this here is a Monitor
269 postProcBdyFe->ts_t = this->ts_t;
272 }
273 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
274
275 CHKERR post_proc_end->writeFile(
276 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
278 };
279
280 auto calculate_traction = [&] {
282 CHKERR VecZeroEntries(CommonData::totalTraction);
284 CHKERR VecAssemblyBegin(CommonData::totalTraction);
285 CHKERR VecAssemblyEnd(CommonData::totalTraction);
287 };
288
289 auto calculate_reactions = [&]() {
291
292 auto res = createDMVector(dM);
293
294 auto assemble_domain = [&]() {
296 auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
297 auto &pip = fe_rhs->getOpPtrVector();
298 fe_rhs->f = res;
299
300 auto integration_rule = [](int, int, int approx_order) {
301 return 2 * approx_order + geom_order - 1;
302 };
303 fe_rhs->getRuleHook = integration_rule;
304 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
305 *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform);
306 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
308 };
309
310 auto assemble_boundary = [&]() {
312 auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
313 auto &pip = fe_rhs->getOpPtrVector();
314 fe_rhs->f = res;
315
316 auto integration_rule = [](int, int, int approx_order) {
317 return 2 * approx_order + geom_order - 1;
318 };
319 fe_rhs->getRuleHook = integration_rule;
320
321 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
322 "GEOMETRY");
323 // We have to integrate on curved face geometry, thus integration weight
324 // have to adjusted.
325 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
326
327 auto u_disp = boost::make_shared<MatrixDouble>();
328 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
329 pip.push_back(
330 new OpSpringRhs("U", u_disp, [this](double, double, double) {
331 return spring_stiffness;
332 }));
333
334 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
335
337 };
338
339 CHKERR assemble_domain();
340 CHKERR assemble_boundary();
341
342 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
343 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
344 m_field_ptr]() {
346 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
347 *m_field_ptr, fe_post_proc_ptr, res)();
349 };
350 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
351 CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
352
354 };
355
356 auto print_max_min = [&](auto &tuple, const std::string msg) {
358 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
359 INSERT_VALUES, SCATTER_FORWARD);
360 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
361 INSERT_VALUES, SCATTER_FORWARD);
362 double max, min;
363 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
364 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
365 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e min %3.4e max %3.4e",
366 msg.c_str(), ts_t, min, max);
368 };
369
370 auto print_traction = [&](const std::string msg) {
372 MoFEM::Interface *m_field_ptr;
373 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
374 if (!m_field_ptr->get_comm_rank()) {
375 const double *t_ptr;
376 CHKERR VecGetArrayRead(CommonData::totalTraction, &t_ptr);
377 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
378 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
379 CHKERR VecRestoreArrayRead(CommonData::totalTraction, &t_ptr);
380 }
382 };
383
384 MOFEM_LOG("CONTACT", Sev::inform)
385 << "Write file at time " << ts_t << " write step " << sTEP;
386
387 CHKERR post_proc();
388 CHKERR calculate_traction();
389 CHKERR calculate_reactions();
390
391 CHKERR print_max_min(uXScatter, "Ux");
392 CHKERR print_max_min(uYScatter, "Uy");
393 if (SPACE_DIM == 3)
394 CHKERR print_max_min(uZScatter, "Uz");
395 CHKERR print_traction("Contact force");
396
397 ++sTEP;
398
400 }

◆ preProcess() [1/2]

MoFEMErrorCode ContactOps::Monitor::preProcess ( )
inline

Definition at line 256 of file PostProcContact.hpp.

256{ return 0; }

◆ preProcess() [2/2]

MoFEMErrorCode ContactOps::Monitor::preProcess ( )
inline

Definition at line 246 of file PostProcContact.hpp.

246{ return 0; }

◆ setScatterVectors()

MoFEMErrorCode ContactOps::Monitor::setScatterVectors ( std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  ux_scatter,
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  uy_scatter,
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >  uz_scatter 
)
inline

Definition at line 434 of file PostProcContact.hpp.

437 {
439 uXScatter = ux_scatter;
440 uYScatter = uy_scatter;
441 uZScatter = uz_scatter;
443 }

Member Data Documentation

◆ atomTest

int ContactOps::Monitor::atomTest
private

Definition at line 468 of file PostProcContact.hpp.

◆ deltaTime

double Monitor::deltaTime
private

Definition at line 462 of file PostProcContact.hpp.

◆ dM

SmartPetscObj< DM > Monitor::dM
private

Definition at line 446 of file PostProcContact.hpp.

◆ integrateTraction

boost::shared_ptr< BoundaryEle > Monitor::integrateTraction
private

Definition at line 456 of file PostProcContact.hpp.

◆ isAxisymmetric

bool ContactOps::Monitor::isAxisymmetric
private

Definition at line 466 of file PostProcContact.hpp.

◆ lastTime

double Monitor::lastTime
private

Definition at line 461 of file PostProcContact.hpp.

◆ mbVertexPostproc

moab::Core Monitor::mbVertexPostproc
private

Definition at line 458 of file PostProcContact.hpp.

◆ mfrontInterface

boost::shared_ptr<GenericElementInterface> ContactOps::Monitor::mfrontInterface
private

Definition at line 467 of file PostProcContact.hpp.

◆ moabVertex

moab::Interface & Monitor::moabVertex
private

Definition at line 459 of file PostProcContact.hpp.

◆ postProcBdyFe

boost::shared_ptr< PostProcEleBdy > Monitor::postProcBdyFe
private

Definition at line 454 of file PostProcContact.hpp.

◆ postProcDomainFe

boost::shared_ptr< PostProcEleDomain > Monitor::postProcDomainFe
private

Definition at line 453 of file PostProcContact.hpp.

◆ postProcMesh

boost::shared_ptr< moab::Core > Monitor::postProcMesh = boost::make_shared<moab::Core>()
private

Definition at line 451 of file PostProcContact.hpp.

◆ sTEP

int Monitor::sTEP
private

Definition at line 463 of file PostProcContact.hpp.

◆ useMFront

bool ContactOps::Monitor::useMFront
private

Definition at line 465 of file PostProcContact.hpp.

◆ uXScatter

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

Definition at line 447 of file PostProcContact.hpp.

◆ uYScatter

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

Definition at line 448 of file PostProcContact.hpp.

◆ uZScatter

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

Definition at line 449 of file PostProcContact.hpp.


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