v0.15.0
Loading...
Searching...
No Matches
PlasticOpsLargeStrains.hpp
Go to the documentation of this file.
1/** \file PlasticOpsLargeStrains.hpp
2 * \example PlasticOpsLargeStrains.hpp
3 */
4
5namespace PlasticOps {
6
7template <int DIM, typename AssemblyDomainEleOp>
10 : public AssemblyDomainEleOp {
12 const std::string row_field_name, const std::string col_field_name,
13 boost::shared_ptr<CommonData> common_data_ptr,
14 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
15 boost::shared_ptr<MatrixDouble> m_D_ptr);
16 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
17 EntitiesFieldData::EntData &col_data);
18
19private:
20 boost::shared_ptr<CommonData> commonDataPtr;
21 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
22 boost::shared_ptr<MatrixDouble> mDPtr;
23 MatrixDouble locMat;
24 MatrixDouble resDiff;
25};
26
27template <int DIM, typename AssemblyDomainEleOp>
30 OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl(
31 const std::string row_field_name, const std::string col_field_name,
32 boost::shared_ptr<CommonData> common_data_ptr,
33 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
34 boost::shared_ptr<MatrixDouble> m_D_ptr)
35 : AssemblyDomainEleOp(row_field_name, col_field_name,
36 AssemblyDomainEleOp::OPROWCOL),
37 commonDataPtr(common_data_ptr),
38 commonHenckyDataPtr(common_hencky_data_ptr), mDPtr(m_D_ptr) {
39 AssemblyDomainEleOp::sYmm = false;
40}
41
42template <int DIM, typename AssemblyDomainEleOp>
44 DIM, GAUSS,
45 AssemblyDomainEleOp>::iNtegrate(EntitiesFieldData::EntData &row_data,
46 EntitiesFieldData::EntData &col_data) {
48
49 FTensor::Index<'i', DIM> i;
50 FTensor::Index<'j', DIM> j;
51 FTensor::Index<'k', DIM> k;
52 FTensor::Index<'l', DIM> l;
53 FTensor::Index<'m', DIM> m;
54 FTensor::Index<'n', DIM> n;
55
56 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
58
59 auto &locMat = AssemblyDomainEleOp::locMat;
60
61 const size_t nb_integration_pts = row_data.getN().size1();
62 const size_t nb_row_base_functions = row_data.getN().size2();
63
64 if (AssemblyDomainEleOp::rowType == MBVERTEX &&
65 AssemblyDomainEleOp::rowSide == 0) {
66 resDiff.resize(DIM * DIM * size_symm, nb_integration_pts, false);
67 auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(resDiff);
68 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
69 auto t_logC_dC =
70 getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
71 auto t_grad =
72 getFTensor2FromMat<DIM, DIM>(*(commonHenckyDataPtr->matGradPtr));
74 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
75 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
77 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
79 t_DLogC_dC(i, j, k, l) = t_D(m, n, k, l) * t_logC_dC(m, n, i, j);
81 t_FDLogC_dC(i, j, k, l) = t_F(i, m) * t_DLogC_dC(m, j, k, l);
83 t_res_diff(i, j, L) = t_FDLogC_dC(i, j, k, l) * t_L(k, l, L);
84 ++t_logC_dC;
85 ++t_grad;
86 ++t_res_diff;
87 }
88 }
89
90 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
91 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
92 auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(resDiff);
93
94 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
95 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
96
97 size_t rr = 0;
98 for (; rr != AssemblyDomainEleOp::nbRows / DIM; ++rr) {
99
100 auto t_mat =
102
104 t_tmp(i, L) = (t_res_diff(i, j, L) * (alpha * t_row_diff_base(j)));
105
106 auto t_col_base = col_data.getFTensor0N(gg, 0);
107 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
108
109 t_mat(i, L) -= (t_col_base * t_tmp(i, L));
110
111 ++t_mat;
112 ++t_col_base;
113 }
114
115 ++t_row_diff_base;
116 }
117
118 for (; rr < nb_row_base_functions; ++rr)
119 ++t_row_diff_base;
120
121 ++t_w;
122 ++t_res_diff;
123 }
124
126}
127
128template <int DIM, typename AssemblyDomainEleOp>
131 : public AssemblyDomainEleOp {
133 const std::string row_field_name, const std::string col_field_name,
134 boost::shared_ptr<CommonData> common_data_ptr,
135 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
136 boost::shared_ptr<MatrixDouble> m_D_ptr);
137 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
138 EntitiesFieldData::EntData &col_data);
139
140private:
141 boost::shared_ptr<CommonData> commonDataPtr;
142 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
143 boost::shared_ptr<MatrixDouble> mDPtr;
144 MatrixDouble resDiff;
145};
146
147template <int DIM, typename AssemblyDomainEleOp>
150 const std::string row_field_name, const std::string col_field_name,
151 boost::shared_ptr<CommonData> common_data_ptr,
152 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
153 boost::shared_ptr<MatrixDouble> m_D_ptr)
154 : AssemblyDomainEleOp(row_field_name, col_field_name,
155 AssemblyDomainEleOp::OPROWCOL),
156 commonDataPtr(common_data_ptr),
157 commonHenckyDataPtr(common_hencky_data_ptr), mDPtr(m_D_ptr) {
158 AssemblyDomainEleOp::sYmm = false;
159}
160
161template <int DIM, typename AssemblyDomainEleOp>
162MoFEMErrorCode
164 iNtegrate(EntitiesFieldData::EntData &row_data,
165 EntitiesFieldData::EntData &col_data) {
167
168 FTensor::Index<'i', DIM> i;
169 FTensor::Index<'j', DIM> j;
170 FTensor::Index<'k', DIM> k;
171 FTensor::Index<'l', DIM> l;
172 FTensor::Index<'m', DIM> m;
173 FTensor::Index<'n', DIM> n;
174 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
176
177 auto &locMat = AssemblyDomainEleOp::locMat;
178
179 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
180 const auto nb_row_base_functions = row_data.getN().size2();
181
182 if (AssemblyDomainEleOp::colType == MBVERTEX &&
183 AssemblyDomainEleOp::colSide == 0) {
184
185 resDiff.resize(size_symm * DIM * DIM, nb_integration_pts, false);
186 auto t_res_diff = getFTensor3FromMat<size_symm, DIM, DIM>(resDiff);
187
188 auto t_res_flow_dstrain =
189 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
190 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
191 auto t_logC_dC =
192 getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
193
194 auto next = [&]() {
195 ++t_res_flow_dstrain;
196 ++t_grad;
197 ++t_logC_dC;
198 ++t_res_diff;
199 };
200
202 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
204 t_diff_grad(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
205
206 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
207 FTensor::Ddg<double, DIM, DIM> t_diff_ls_dlogC_dC;
210
211 t_diff_ls_dlogC_dC(i, j, k, l) =
212 (t_res_flow_dstrain(i, j, m, n)) * (t_logC_dC(m, n, k, l) / 2);
213
214 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
215 t_dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
216
218 t_res_diff(L, i, j) =
219 (t_L(m, n, L) * t_diff_ls_dlogC_dC(m, n, k, l)) * t_dC_dF(k, l, i, j);
220 next();
221 }
222 }
223
224 auto t_res_diff = getFTensor3FromMat<size_symm, DIM, DIM>(resDiff);
225
226 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
227 auto t_row_base = row_data.getFTensor0N();
228 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
229 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
230 ++t_w;
231
232 size_t rr = 0;
233 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
234 const auto row_base = alpha * t_row_base;
235 auto t_mat =
237 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
238 for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; ++cc) {
239 t_mat(L, l) += row_base * (t_res_diff(L, l, k) * t_col_diff_base(k));
240 ++t_mat;
241 ++t_col_diff_base;
242 }
243 ++t_row_base;
244 }
245
246 for (; rr < nb_row_base_functions; ++rr)
247 ++t_row_base;
248
249 ++t_res_diff;
250 }
251
253}
254
255template <int DIM, typename AssemblyDomainEleOp>
257 public AssemblyDomainEleOp {
259 const std::string row_field_name, const std::string col_field_name,
260 boost::shared_ptr<CommonData> common_data_ptr,
261 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
262 boost::shared_ptr<MatrixDouble> m_D_ptr);
263 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
264 EntitiesFieldData::EntData &col_data);
265
266private:
267 boost::shared_ptr<CommonData> commonDataPtr;
268 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
269 boost::shared_ptr<MatrixDouble> mDPtr;
270 MatrixDouble resDiff;
271};
272
273template <int DIM, typename AssemblyDomainEleOp>
276 const std::string row_field_name, const std::string col_field_name,
277 boost::shared_ptr<CommonData> common_data_ptr,
278 boost::shared_ptr<HenckyOps::CommonData> common_hencky_data_ptr,
279 boost::shared_ptr<MatrixDouble> m_D_ptr)
280 : AssemblyDomainEleOp(row_field_name, col_field_name,
281 DomainEleOp::OPROWCOL),
282 commonDataPtr(common_data_ptr),
283 commonHenckyDataPtr(common_hencky_data_ptr), mDPtr(m_D_ptr) {
284 AssemblyDomainEleOp::sYmm = false;
285}
286
287template <int DIM, typename AssemblyDomainEleOp>
288MoFEMErrorCode
290 iNtegrate(EntitiesFieldData::EntData &row_data,
291 EntitiesFieldData::EntData &col_data) {
293
294 FTensor::Index<'i', DIM> i;
295 FTensor::Index<'j', DIM> j;
296 FTensor::Index<'k', DIM> k;
297 FTensor::Index<'l', DIM> l;
298
299 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
300
301 AssemblyDomainEleOp::locMat.resize(AssemblyDomainEleOp::nbRows, AssemblyDomainEleOp::nbCols,
302 false);
303 AssemblyDomainEleOp::locMat.clear();
304
305 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
306 const auto nb_row_base_functions = row_data.getN().size2();
307
308 if (AssemblyDomainEleOp::colType == MBVERTEX &&
309 AssemblyDomainEleOp::colSide == 0) {
310
311 resDiff.resize(DIM * DIM, nb_integration_pts, false);
312 auto t_res_diff = getFTensor2FromMat<DIM, DIM>(resDiff);
313 auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
314 auto t_logC_dC =
315 getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
316 auto t_c_dstrain =
317 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
318
319 auto next = [&]() {
320 ++t_grad;
321 ++t_logC_dC;
322 ++t_c_dstrain;
323 ++t_res_diff;
324 };
325
326 auto t_diff_grad_symmetrise = diff_symmetrize(FTensor::Number<DIM>());
327
328 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
332
333 t_diff_ls_dlog_c(k, l) =
334 (t_c_dstrain(i, j)) * (t_logC_dC(i, j, k, l) / 2);
335 t_F(i, j) = t_grad(i, j) + t_kd(i, j);
336 t_dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
337
339 t_res_diff(i, j) = (t_diff_ls_dlog_c(k, l) * t_dC_dF(k, l, i, j));
340 next();
341 }
342 }
343
344 auto t_res_diff = getFTensor2FromMat<DIM, DIM>(resDiff);
345
346 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
347 auto t_row_base = row_data.getFTensor0N();
348 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
349 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
350 ++t_w;
351
352 auto t_mat =
353 getFTensor1FromPtr<DIM>(AssemblyDomainEleOp::locMat.data().data());
354 size_t rr = 0;
355 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
356 const auto row_base = alpha * t_row_base;
357 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
358 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; cc++) {
359 t_mat(i) += row_base * (t_res_diff(i, j) * t_col_diff_base(j));
360 ++t_mat;
361 ++t_col_diff_base;
362 }
363 ++t_row_base;
364 }
365 for (; rr != nb_row_base_functions; ++rr)
366 ++t_row_base;
367
368 ++t_res_diff;
369 }
370
372}
373
374}; // namespace PlasticOps
Kronecker Delta class.
#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()
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto symm_L_tensor(FTensor::Number< DIM >)
auto diff_symmetrize(FTensor::Number< DIM >)
static FTensor::Tensor2< FTensor::PackPtr< double *, 3 >, 2, 3 > get_mat_vector_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
constexpr auto size_symm
Definition plastic.cpp:42
FTensor::Index< 'm', 3 > m