v0.13.2
Loading...
Searching...
No Matches
PostProcContact.hpp
Go to the documentation of this file.
1/**
2 * \file PostProcContact.hpp
3 *
4 *
5 * @copyright Copyright (c) 2023
6 */
7
8namespace ContactOps {
9
10template <int DIM> struct PostProcEleByDim;
11
12template <> struct PostProcEleByDim<2> {
13 using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
14 using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
15 using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
16};
17
18template <> struct PostProcEleByDim<3> {
19 using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
20 using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
21 using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
22};
23
27
29 OpAssembleTraction(boost::shared_ptr<CommonData> common_data_ptr,
30 SmartPetscObj<Vec> total_traction);
31 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
32
33private:
34 boost::shared_ptr<CommonData> commonDataPtr;
35 SmartPetscObj<Vec> totalTraction;
36};
37
38struct Monitor : public FEMethod {
39
40 Monitor(SmartPetscObj<DM> &dm,
41 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
42 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
43 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
44 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
45 uZScatter(uz_scatter), moabVertex(mbVertexPostproc), sTEP(0) {
46
47 MoFEM::Interface *m_field_ptr;
48 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
50 createSmartVectorMPI(m_field_ptr->get_comm(),
51 (m_field_ptr->get_comm_rank() == 0) ? 3 : 0, 3);
52
53 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
54
55 auto get_post_proc_domain_fe = [&]() {
56 auto post_proc_fe =
57 boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
58 auto &pip = post_proc_fe->getOpPtrVector();
59
60 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
61 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
62 henky_common_data_ptr->matGradPtr = common_data_ptr->mGradPtr();
63 henky_common_data_ptr->matDPtr = common_data_ptr->mDPtr();
64
65 auto push_domain_ops = [&](auto &pip) {
66 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
67 pip, {H1, HDIV} /*, "GEOMETRY"*/)),
68 "Apply base transform");
70 ContactOps::addMatBlockOps(*m_field_ptr, pip, "U", "MAT_ELASTIC",
71 common_data_ptr->mDPtr(), Sev::inform),
72 "Set block data");
74 "U", common_data_ptr->mGradPtr()));
75 pip.push_back(
76 new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
77 pip.push_back(
78 new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
79 pip.push_back(
80 new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
81 pip.push_back(
82 new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
83 pip.push_back(
84 new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
85 pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
86 "SIGMA", common_data_ptr->contactStressPtr()));
87 };
88
89 // Evaluate domain on side element
90 if constexpr (SPACE_DIM == 3) {
91 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
92 pip, {HDIV} /*, "GEOMETRY"*/)),
93 "Apply transform");
94 // create OP which run element on side
95 auto op_loop_side =
96 new OpLoopSide<SideEle>(*m_field_ptr, "dFE", SPACE_DIM);
97 // push ops to side element, through op_loop_side operator
98 push_domain_ops(op_loop_side->getOpPtrVector());
99 pip.push_back(op_loop_side);
100 // evaluate traction
101 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
102 "SIGMA", common_data_ptr->contactTractionPtr()));
103 } else {
104 push_domain_ops(pip);
105 }
106
107 auto u_ptr = boost::make_shared<MatrixDouble>();
108 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
109
110 post_proc_fe->getOpPtrVector().push_back(
111
112 new OpPPMap(
113
114 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
115
116 {},
117
118 {{"U", u_ptr},
119
120 // Note: post-process tractions in 3d, i.e. when mesh is
121 // post-process on skin
122 {"t", (SPACE_DIM == 3) ? common_data_ptr->contactTractionPtr()
123 : nullptr}},
124
125 {
126
127 {"SIGMA", common_data_ptr->contactStressPtr()},
128
129 {"G", common_data_ptr->mGradPtr()},
130
131 {"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
132
133 },
134
135 {}
136
137 )
138
139 );
140
141 return post_proc_fe;
142 };
143
144 auto get_post_proc_bdy_fe = [&]() {
145 auto post_proc_fe =
146 boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
147 auto &pip = post_proc_fe->getOpPtrVector();
148
149 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
150
151 CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
152 pip, {HDIV} /*, "GEOMETRY"*/)),
153 "Apply transform");
154 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
155 "U", common_data_ptr->contactDispPtr()));
156 pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
157 "SIGMA", common_data_ptr->contactTractionPtr()));
158
159 pip.push_back(
160
161 new OpPPMap(
162
163 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
164
165 {},
166
167 {{"U", common_data_ptr->contactDispPtr()},
168 {"t", common_data_ptr->contactTractionPtr()}},
169
170 {},
171
172 {}
173
174 )
175
176 );
177
178 return post_proc_fe;
179 };
180
181 auto get_integrate_traction = [&]() {
182 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
183 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
185 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
186 integrate_traction->getOpPtrVector(), {HDIV} /*, "GEOMETRY"*/)),
187 "Apply transfrom");
188 integrate_traction->getOpPtrVector().push_back(
189 new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
190 "SIGMA", common_data_ptr->contactTractionPtr()));
191 integrate_traction->getOpPtrVector().push_back(
192 new OpAssembleTraction(common_data_ptr, totalTraction));
193 integrate_traction->getRuleHook = [](int, int, int approx_order) {
194 return 2 * approx_order;
195 };
196 return integrate_traction;
197 };
198
199 postProcDomainFe = get_post_proc_domain_fe();
200 if constexpr (SPACE_DIM == 2)
201 postProcBdyFe = get_post_proc_bdy_fe();
202 integrateTraction = get_integrate_traction();
203 }
204
205 MoFEMErrorCode preProcess() { return 0; }
206 MoFEMErrorCode operator()() { return 0; }
207
208 MoFEMErrorCode postProcess() {
210 MoFEM::Interface *m_field_ptr;
211 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
212
213 auto post_proc = [&]() {
215
216 auto post_proc_begin =
217 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
218 postProcMesh);
219 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
220 *m_field_ptr, postProcMesh);
221
222 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
223 if (!postProcBdyFe) {
224 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
225 } else {
226 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
227 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
228 }
229 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
230
231 CHKERR post_proc_end->writeFile(
232 "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
234 };
235
236 auto calculate_traction = [&] {
238 CHKERR VecZeroEntries(totalTraction);
239 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
240 CHKERR VecAssemblyBegin(totalTraction);
241 CHKERR VecAssemblyEnd(totalTraction);
243 };
244
245 auto print_max_min = [&](auto &tuple, const std::string msg) {
247 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
248 INSERT_VALUES, SCATTER_FORWARD);
249 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
250 INSERT_VALUES, SCATTER_FORWARD);
251 double max, min;
252 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
253 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
254 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e min %3.4e max %3.4e",
255 msg.c_str(), ts_t, min, max);
257 };
258
259 auto print_traction = [&](const std::string msg) {
261 MoFEM::Interface *m_field_ptr;
262 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
263 if (!m_field_ptr->get_comm_rank()) {
264 const double *t_ptr;
265 CHKERR VecGetArrayRead(totalTraction, &t_ptr);
266 MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
267 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
268 CHKERR VecRestoreArrayRead(totalTraction, &t_ptr);
269 }
271 };
272
273 MOFEM_LOG("CONTACT", Sev::inform)
274 << "Write file at time " << ts_t << " write step " << sTEP;
275
276 CHKERR post_proc();
277 CHKERR calculate_traction();
278
279 CHKERR print_max_min(uXScatter, "Ux");
280 CHKERR print_max_min(uYScatter, "Uy");
281 if (SPACE_DIM == 3)
282 CHKERR print_max_min(uZScatter, "Uz");
283 CHKERR print_traction("Force");
284
285 ++sTEP;
286
288 }
289
290private:
291 SmartPetscObj<DM> dM;
292 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
293 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
294 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
295
296 SmartPetscObj<Vec> totalTraction;
297
298 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
299
300 boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
301 boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
302
303 boost::shared_ptr<BoundaryEle> integrateTraction;
304
306 moab::Interface &moabVertex;
307
308 double lastTime;
309 double deltaTime;
310 int sTEP;
311};
312
313OpAssembleTraction::OpAssembleTraction(
314 boost::shared_ptr<CommonData> common_data_ptr,
315 SmartPetscObj<Vec> total_traction)
317 commonDataPtr(common_data_ptr), totalTraction(total_traction) {}
318
319MoFEMErrorCode OpAssembleTraction::doWork(int side, EntityType type,
320 EntData &data) {
322
323 FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
324
325 auto t_w = getFTensor0IntegrationWeight();
326 auto t_traction =
327 getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
328
329 const auto nb_gauss_pts = getGaussPts().size2();
330 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
331 const double alpha = t_w * getMeasure();
332 t_sum_t(i) += alpha * t_traction(i);
333 ++t_w;
334 ++t_traction;
335 }
336
337 constexpr int ind[] = {0, 1, 2};
338 CHKERR VecSetValues(totalTraction, 3, ind, &t_sum_t(0), ADD_VALUES);
339
341}
342
343} // namespace ContactOps
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
constexpr int SPACE_DIM
#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 MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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 MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:166
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: contact.cpp:755
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static constexpr int approx_order
moab::Interface & moabVertex
MoFEMErrorCode postProcess()
boost::shared_ptr< BoundaryEle > integrateTraction
SmartPetscObj< Vec > totalTraction
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
boost::shared_ptr< PostProcEleDomain > postProcDomainFe
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)
SmartPetscObj< DM > dM
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< moab::Core > postProcMesh
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode preProcess()
boost::shared_ptr< PostProcEleBdy > postProcBdyFe
MoFEMErrorCode operator()()
SmartPetscObj< Vec > totalTraction
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< CommonData > commonDataPtr
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomain
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleDomain
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Post post-proc data at points from hash maps.