v0.13.2
Loading...
Searching...
No Matches
PlasticOpsGeneric.hpp
Go to the documentation of this file.
1
2
3/** \file PlasticOpsGeneric.hpp
4 * \example PlasticOpsGeneric.hpp
5 */
6
7namespace PlasticOps {
8
10 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
12 commonDataPtr(common_data_ptr) {
13 // Opetor is only executed for vertices
14 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
15}
16
17MoFEMErrorCode OpCalculatePlasticSurface::doWork(int side, EntityType type,
18 EntData &data) {
20
21 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
22 auto t_stress =
23 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commonDataPtr->mStressPtr));
24
25 commonDataPtr->plasticSurface.resize(nb_gauss_pts, false);
26 commonDataPtr->plasticFlow.resize(size_symm, nb_gauss_pts, false);
27 auto t_flow =
28 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticFlow);
29
30 for (auto &f : commonDataPtr->plasticSurface) {
31 f = platsic_surface(deviator(t_stress, trace(t_stress)));
32 auto t_flow_tmp = plastic_flow(f,
33
34 deviator(t_stress, trace(t_stress)),
35
37 t_flow(i, j) = t_flow_tmp(i, j);
38 ++t_flow;
39 ++t_stress;
40 }
41
43}
44
46 const std::string field_name, MoFEM::Interface &m_field,
47 boost::shared_ptr<CommonData> common_data_ptr,
48 boost::shared_ptr<MatrixDouble> m_D_ptr)
49 : DomainEleOp(field_name, DomainEleOp::OPROW), mField(m_field),
50 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
51 // Opetor is only executed for vertices
52 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
53}
54
55MoFEMErrorCode OpCalculatePlasticity::doWork(int side, EntityType type,
56 EntData &data) {
58
59 const size_t nb_gauss_pts = getGaussPts().size2();
61 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
62 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
63 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
64 auto t_flow =
65 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticFlow);
66 auto t_plastic_strain_dot =
67 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
68 auto t_stress =
69 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commonDataPtr->mStressPtr));
70
71 auto t_D =
72 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
73 auto t_D_Op = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
74
75 auto t_diff_plastic_strain = diff_tensor();
76 auto t_diff_deviator = diff_deviator(diff_tensor());
77
80 t_flow_dir_dstress(i, j, k, l) =
81 1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l));
82 t_flow_dir_dstrain(i, j, k, l) =
83 t_flow_dir_dstress(i, j, m, n) * t_D_Op(m, n, k, l);
84
85 commonDataPtr->resC.resize(nb_gauss_pts, false);
86 commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
87 commonDataPtr->resCdStrain.resize(size_symm, nb_gauss_pts, false);
88 commonDataPtr->resCdStrainDot.resize(size_symm, nb_gauss_pts, false);
89 commonDataPtr->resFlow.resize(size_symm, nb_gauss_pts, false);
90 commonDataPtr->resFlowDtau.resize(size_symm, nb_gauss_pts, false);
91 commonDataPtr->resFlowDstrain.resize(size_symm * size_symm, nb_gauss_pts,
92 false);
93 commonDataPtr->resFlowDstrainDot.resize(size_symm * size_symm, nb_gauss_pts,
94 false);
95
96 commonDataPtr->resC.clear();
97 commonDataPtr->resCdTau.clear();
98 commonDataPtr->resCdStrain.clear();
99 commonDataPtr->resCdStrainDot.clear();
100 commonDataPtr->resFlow.clear();
101 commonDataPtr->resFlowDtau.clear();
102 commonDataPtr->resFlowDstrain.clear();
103 commonDataPtr->resFlowDstrainDot.clear();
104
105 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
106 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
107 auto t_res_c_dstrain =
108 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
109 auto t_res_c_dstrain_dot =
110 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrainDot);
111 auto t_res_flow =
112 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resFlow);
113 auto t_res_flow_dtau =
114 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resFlowDtau);
115 auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
116 commonDataPtr->resFlowDstrain);
117 auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
118 commonDataPtr->resFlowDstrainDot);
119
120 auto next = [&]() {
121 ++t_tau;
122 ++t_tau_dot;
123 ++t_f;
124 ++t_flow;
125 ++t_plastic_strain_dot;
126 ++t_stress;
127 ++t_res_c;
128 ++t_res_c_dtau;
129 ++t_res_c_dstrain;
130 ++t_res_c_dstrain_dot;
131 ++t_res_flow;
132 ++t_res_flow_dtau;
133 ++t_res_flow_dstrain;
134 ++t_res_flow_dstrain_dot;
135 ++t_w;
136 };
137
138 auto get_avtive_pts = [&]() {
139 int nb_points_avtive_on_elem = 0;
140 int nb_points_on_elem = 0;
141
142 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
143 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
144 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
145 auto t_plastic_strain_dot =
146 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
147
148 for (auto &f : commonDataPtr->plasticSurface) {
149 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
150 const auto ww = w(eqiv, t_tau_dot, t_f, hardening(t_tau));
151 const auto sign_ww = constrian_sign(ww);
152
153 ++nb_points_on_elem;
154 if (sign_ww > 0) {
155 ++nb_points_avtive_on_elem;
156 }
157
158 ++t_tau;
159 ++t_tau_dot;
160 ++t_f;
161 ++t_plastic_strain_dot;
162 }
163
164 int &active_points = commonDataPtr->activityData[0];
165 int &avtive_full_elems = commonDataPtr->activityData[1];
166 int &avtive_elems = commonDataPtr->activityData[2];
167 int &nb_points = commonDataPtr->activityData[3];
168 int &nb_elements = commonDataPtr->activityData[4];
169
170 ++nb_elements;
171 nb_points += nb_points_on_elem;
172 if (nb_points_avtive_on_elem > 0) {
173 ++avtive_elems;
174 active_points += nb_points_avtive_on_elem;
175 if (nb_points_avtive_on_elem == nb_points_on_elem) {
176 ++avtive_full_elems;
177 }
178 }
179
180 if (nb_points_avtive_on_elem != nb_points_on_elem)
181 return 1;
182 else
183 return 0;
184 };
185
186 if (getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
187 get_avtive_pts();
188 }
189
190 for (auto &f : commonDataPtr->plasticSurface) {
191
192 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
193 auto t_diff_eqiv = diff_equivalent_strain_dot(eqiv, t_plastic_strain_dot,
194 t_diff_plastic_strain);
195
196 const auto sigma_y = hardening(t_tau);
197 const auto d_sigma_y = hardening_dtau(t_tau);
198
199 auto ww = w(eqiv, t_tau_dot, t_f, sigma_y);
200 auto abs_ww = constrain_abs(ww);
201 auto sign_ww = constrian_sign(ww);
202
203 auto c = constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww);
204 auto c_dot_tau = diff_constrain_ddot_tau(sign_ww, eqiv);
205 auto c_equiv = diff_constrain_eqiv(sign_ww, eqiv, t_tau_dot);
206 auto c_sigma_y = diff_constrain_dsigma_y(sign_ww);
207 auto c_f = diff_constrain_df(sign_ww);
208
209 auto t_dev_stress = deviator(t_stress, trace(t_stress));
211 t_flow_dir(k, l) = 1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
213 t_flow_dstrain(i, j) = t_flow(k, l) * t_D_Op(k, l, i, j);
214
215 auto get_res_c = [&]() { return c; };
216
217 auto get_res_c_dstrain = [&](auto &t_diff_res) {
218 t_diff_res(i, j) = c_f * t_flow_dstrain(i, j);
219 };
220
221 auto get_res_c_dstrain_dot = [&](auto &t_diff_res) {
222 t_diff_res(i, j) = (getTSa() * c_equiv) * t_diff_eqiv(i, j);
223 };
224
225 auto get_res_c_dtau = [&]() {
226 return getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
227 };
228
229 auto get_res_flow = [&](auto &t_res_flow) {
230 const auto a = sigma_y;
231 const auto b = t_tau_dot;
232 t_res_flow(k, l) = a * t_plastic_strain_dot(k, l) - b * t_flow_dir(k, l);
233 };
234
235 auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
236 const auto da = d_sigma_y;
237 const auto db = getTSa();
238 t_res_flow_dtau(k, l) =
239 da * t_plastic_strain_dot(k, l) - db * t_flow_dir(k, l);
240 };
241
242 auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
243 const auto b = t_tau_dot;
244 t_res_flow_dstrain(m, n, k, l) = -t_flow_dir_dstrain(m, n, k, l) * b;
245 };
246
247 auto get_res_flow_dstrain_dot = [&](auto &t_res_flow_dstrain_dot) {
248 const auto a = sigma_y;
249 t_res_flow_dstrain_dot(m, n, k, l) =
250 (a * getTSa()) * t_diff_plastic_strain(m, n, k, l);
251 };
252
253 t_res_c = get_res_c();
254 get_res_flow(t_res_flow);
255
256 if (getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
257 t_res_c_dtau = get_res_c_dtau();
258 get_res_c_dstrain(t_res_c_dstrain);
259 get_res_c_dstrain_dot(t_res_c_dstrain_dot);
260 get_res_flow_dtau(t_res_flow_dtau);
261 get_res_flow_dstrain(t_res_flow_dstrain);
262 get_res_flow_dstrain_dot(t_res_flow_dstrain_dot);
263 }
264
265 next();
266 }
267
269}
270
272 boost::shared_ptr<CommonData> common_data_ptr,
273 boost::shared_ptr<MatrixDouble> m_D_ptr,
274 const double scale)
276 commonDataPtr(common_data_ptr), scaleStress(scale), mDPtr(m_D_ptr) {
277 // Operator is only executed for vertices
278 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
279}
280
281//! [Calculate stress]
282MoFEMErrorCode OpPlasticStress::doWork(int side, EntityType type,
283 EntData &data) {
285 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
286 commonDataPtr->mStressPtr->resize((SPACE_DIM * (SPACE_DIM + 1)) / 2,
287 nb_gauss_pts);
288 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
289 auto t_strain =
290 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commonDataPtr->mStrainPtr));
291 auto t_plastic_strain =
292 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrain);
293 auto t_stress =
294 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commonDataPtr->mStressPtr));
295
296 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
297 t_stress(i, j) =
298 t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
299 t_stress(i, j) /= scaleStress;
300 ++t_strain;
301 ++t_plastic_strain;
302 ++t_stress;
303 }
304
306}
307//! [Calculate stress]
308
310 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
311 boost::shared_ptr<MatrixDouble> m_D_ptr)
313 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
314
315MoFEMErrorCode
316OpCalculatePlasticFlowRhs::iNtegrate(EntitiesFieldData::EntData &data) {
318
319 const auto nb_integration_pts = getGaussPts().size2();
320 const auto nb_base_functions = data.getN().size2();
321
322 auto t_res_flow =
323 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resFlow);
324
325 auto t_L = symm_L_tensor();
326
327 auto next = [&]() { ++t_res_flow; };
328
329 auto t_w = getFTensor0IntegrationWeight();
330 auto t_base = data.getFTensor0N();
331 auto &nf = AssemblyDomainEleOp::locF;
332 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
333 double alpha = getMeasure() * t_w;
334 ++t_w;
335
337 t_rhs(L) = alpha * (t_res_flow(i, j) * t_L(i, j, L));
338 next();
339
340 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
341 size_t bb = 0;
342 for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
343 t_nf(L) += t_base * t_rhs(L);
344 ++t_base;
345 ++t_nf;
346 }
347 for (; bb < nb_base_functions; ++bb)
348 ++t_base;
349 }
350
352}
353
355 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
356 boost::shared_ptr<MatrixDouble> m_D_ptr)
358 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
359
360MoFEMErrorCode
361OpCalculateConstraintsRhs::iNtegrate(EntitiesFieldData::EntData &data) {
363
364 const size_t nb_integration_pts = getGaussPts().size2();
365 const size_t nb_base_functions = data.getN().size2();
366
367 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
368
369 auto next = [&]() { ++t_res_c; };
370
371 auto t_w = getFTensor0IntegrationWeight();
372 auto &nf = AssemblyDomainEleOp::locF;
373 auto t_base = data.getFTensor0N();
374 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
375 const double alpha = getMeasure() * t_w;
376 ++t_w;
377 const auto res = alpha * t_res_c;
378 next();
379
380 size_t bb = 0;
381 for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
382 nf[bb] += t_base * res;
383 ++t_base;
384 }
385 for (; bb < nb_base_functions; ++bb)
386 ++t_base;
387 }
388
390}
391
393 const std::string row_field_name, const std::string col_field_name,
394 boost::shared_ptr<CommonData> common_data_ptr,
395 boost::shared_ptr<MatrixDouble> m_D_ptr)
396 : AssemblyDomainEleOp(row_field_name, col_field_name,
397 DomainEleOp::OPROWCOL),
398 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
399 sYmm = false;
400}
401
402static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
405 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
406 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
407 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
408}
409
410static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
413 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
414 &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
415 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
416 &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
417 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
418 &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
419 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
420 &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
421 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
422 &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
423 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
424 &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
425}
426
427MoFEMErrorCode
428OpCalculatePlasticFlowLhs_dEP::iNtegrate(EntitiesFieldData::EntData &row_data,
429 EntitiesFieldData::EntData &col_data) {
431
433
434 const auto nb_integration_pts = getGaussPts().size2();
435 const auto nb_row_base_functions = row_data.getN().size2();
436
437 auto t_res_flow_dstrain = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
438 commonDataPtr->resFlowDstrain);
439 auto t_res_flow_dstrain_dot = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(
440 commonDataPtr->resFlowDstrainDot);
441 auto t_L = symm_L_tensor();
442
443 auto next = [&]() {
444 ++t_res_flow_dstrain;
445 ++t_res_flow_dstrain_dot;
446 };
447
448 auto t_w = getFTensor0IntegrationWeight();
449 auto t_row_base = row_data.getFTensor0N();
450 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
451 double alpha = getMeasure() * t_w;
452 ++t_w;
453
455 t_res_mat(O, L) =
456 alpha * (t_L(i, j, O) * ((t_res_flow_dstrain_dot(i, j, k, l) -
457 t_res_flow_dstrain(i, j, k, l)) *
458 t_L(k, l, L)));
459 next();
460
461 size_t rr = 0;
462 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
463 auto t_mat = get_mat_tensor_sym_dtensor_sym(rr, locMat,
465 auto t_col_base = col_data.getFTensor0N(gg, 0);
466 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
467 t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
468 ++t_mat;
469 ++t_col_base;
470 }
471
472 ++t_row_base;
473 }
474
475 for (; rr < nb_row_base_functions; ++rr)
476 ++t_row_base;
477 }
478
480}
481
483 const std::string row_field_name, const std::string col_field_name,
484 boost::shared_ptr<CommonData> common_data_ptr,
485 boost::shared_ptr<MatrixDouble> m_D_ptr)
486 : AssemblyDomainEleOp(row_field_name, col_field_name,
487 DomainEleOp::OPROWCOL),
488 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
489 sYmm = false;
490}
491
492static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
495 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
496}
497
498static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
501 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
502 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
503}
504
506 EntitiesFieldData::EntData &row_data,
507 EntitiesFieldData::EntData &col_data) {
509
510 const auto nb_integration_pts = getGaussPts().size2();
511 const size_t nb_row_base_functions = row_data.getN().size2();
513
514 const auto type = getFEType();
515 const auto nb_nodes = moab::CN::VerticesPerEntity(type);
516
517 auto t_res_flow_dtau =
518 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resFlowDtau);
519
520 auto t_L = symm_L_tensor();
521
522 auto next = [&]() { ++t_res_flow_dtau; };
523
524 auto t_w = getFTensor0IntegrationWeight();
525 auto t_row_base = row_data.getFTensor0N();
526 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
527 double alpha = getMeasure() * t_w;
528 ++t_w;
530 t_res_vec(L) = alpha * (t_res_flow_dtau(i, j) * t_L(i, j, L));
531 next();
532
533 size_t rr = 0;
534 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
535 auto t_mat =
537 auto t_col_base = col_data.getFTensor0N(gg, 0);
538 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
539 t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
540 ++t_mat;
541 ++t_col_base;
542 }
543 ++t_row_base;
544 }
545 for (; rr != nb_row_base_functions; ++rr)
546 ++t_row_base;
547 }
548
550}
551
553 const std::string row_field_name, const std::string col_field_name,
554 boost::shared_ptr<CommonData> common_data_ptr,
555 boost::shared_ptr<MatrixDouble> m_D_ptr)
556 : AssemblyDomainEleOp(row_field_name, col_field_name,
557 DomainEleOp::OPROWCOL),
558 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
559 sYmm = false;
560}
561
564 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
565}
566
569 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
570}
571
572MoFEMErrorCode
573OpCalculateConstraintsLhs_dEP::iNtegrate(EntitiesFieldData::EntData &row_data,
574 EntitiesFieldData::EntData &col_data) {
576
577 const auto nb_integration_pts = getGaussPts().size2();
578 const auto nb_row_base_functions = row_data.getN().size2();
579
580 auto t_c_dstrain =
581 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
582 auto t_c_dstrain_dot =
583 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrainDot);
584
585 auto next = [&]() {
586 ++t_c_dstrain;
587 ++t_c_dstrain_dot;
588 };
589
590 auto t_L = symm_L_tensor();
591
592 auto t_w = getFTensor0IntegrationWeight();
593 auto t_row_base = row_data.getFTensor0N();
594 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
595 const double alpha = getMeasure() * t_w;
596 ++t_w;
597
599 t_res_vec(L) = t_L(i, j, L) * (t_c_dstrain_dot(i, j) - t_c_dstrain(i, j));
600 next();
601
602 auto t_mat =
604 size_t rr = 0;
605 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
606 const auto row_base = alpha * t_row_base;
607 auto t_col_base = col_data.getFTensor0N(gg, 0);
608 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
609 t_mat(L) += (row_base * t_col_base) * t_res_vec(L);
610 ++t_mat;
611 ++t_col_base;
612 }
613 ++t_row_base;
614 }
615 for (; rr != nb_row_base_functions; ++rr)
616 ++t_row_base;
617 }
618
620}
621
623 const std::string row_field_name, const std::string col_field_name,
624 boost::shared_ptr<CommonData> common_data_ptr)
625 : AssemblyDomainEleOp(row_field_name, col_field_name,
626 DomainEleOp::OPROWCOL),
627 commonDataPtr(common_data_ptr) {
628 sYmm = false;
629}
630
632 EntitiesFieldData::EntData &row_data,
633 EntitiesFieldData::EntData &col_data) {
635
636 const auto nb_integration_pts = getGaussPts().size2();
637 const auto nb_row_base_functions = row_data.getN().size2();
638
639 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
640 auto next = [&]() { ++t_res_c_dtau; };
641
642 auto t_w = getFTensor0IntegrationWeight();
643 auto t_row_base = row_data.getFTensor0N();
644 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
645 const double alpha = getMeasure() * t_w;
646 ++t_w;
647
648 const auto res = alpha * (t_res_c_dtau);
649 next();
650
651 auto mat_ptr = locMat.data().begin();
652 size_t rr = 0;
653 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
654 auto t_col_base = col_data.getFTensor0N(gg, 0);
655 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
656 *mat_ptr += t_row_base * t_col_base * res;
657 ++t_col_base;
658 ++mat_ptr;
659 }
660 ++t_row_base;
661 }
662 for (; rr < nb_row_base_functions; ++rr)
663 ++t_row_base;
664 }
665
667}
668}; // namespace PlasticOps
constexpr double a
constexpr int SPACE_DIM
#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
const double c
speed of light (cm/ns)
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
auto symm_L_tensor()
Definition: PlasticOps.hpp:315
auto deviator(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress, double trace)
Definition: PlasticOps.hpp:373
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
double constrian_sign(double x)
Definition: PlasticOps.hpp:488
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
Definition: PlasticOps.hpp:440
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
FTensor::Index< 'j', SPACE_DIM > j
Definition: PlasticOps.hpp:93
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:103
double diff_constrain_eqiv(double sign, double eqiv, double dot_tau)
Definition: PlasticOps.hpp:540
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, SPACE_DIM > &&t_diff_deviator)
Definition: PlasticOps.hpp:445
FTensor::Index< 'O', size_symm > O
Definition: PlasticOps.hpp:107
FTensor::Index< 'L', size_symm > L
Definition: PlasticOps.hpp:106
FTensor::Index< 'l', SPACE_DIM > l
Definition: PlasticOps.hpp:95
FTensor::Index< 'k', SPACE_DIM > k
Definition: PlasticOps.hpp:94
auto diff_tensor()
[Operators definitions]
Definition: PlasticOps.hpp:308
auto diff_constrain_df(double sign)
Definition: PlasticOps.hpp:545
double w(double eqiv, double dot_tau, double f, double sigma_y)
Definition: PlasticOps.hpp:513
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:104
auto diff_constrain_dsigma_y(double sign)
Definition: PlasticOps.hpp:547
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:101
auto diff_deviator(FTensor::Ddg< double, SPACE_DIM, SPACE_DIM > &&t_diff_stress)
Definition: PlasticOps.hpp:387
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: PlasticOps.hpp:92
FTensor::Index< 'm', SPACE_DIM > m
Definition: PlasticOps.hpp:96
double constrain_abs(double x)
Definition: PlasticOps.hpp:502
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:102
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain)
Definition: PlasticOps.hpp:576
double trace(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress)
Definition: PlasticOps.hpp:364
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w)
Definition: PlasticOps.hpp:528
double diff_constrain_ddot_tau(double sign, double eqiv)
Definition: PlasticOps.hpp:536
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_plastic_strain_dot)
Definition: PlasticOps.hpp:567
FTensor::Index< 'n', SPACE_DIM > n
Definition: PlasticOps.hpp:97
double scale
Definition: plastic.cpp:97
constexpr auto size_symm
Definition: plastic.cpp:33
double hardening(double tau)
Definition: plastic.cpp:122
double hardening_dtau(double tau)
Definition: plastic.cpp:126
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:290
OpCalculateConstraintsLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mat_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:302
OpCalculateConstraintsLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpCalculateConstraintsRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:162
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:235
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:248
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
[Calculate stress]
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:151
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:116
OpCalculatePlasticity(const std::string field_name, MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:127
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:128
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
Definition: PlasticOps.hpp:141
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > mDPtr, const double scale=1)
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOps.hpp:139