v0.15.0
Loading...
Searching...
No Matches
PlasticOpsSmallStrains.hpp
Go to the documentation of this file.
1
2
3/** \file PlasticOpsSmallStrains.hpp
4 * \example PlasticOpsSmallStrains.hpp
5 */
6
7namespace PlasticOps {
8
9template <int DIM, typename AssemblyDomainEleOp>
12 : public AssemblyDomainEleOp {
14 const std::string row_field_name, const std::string col_field_name,
15 boost::shared_ptr<CommonData> common_data_ptr,
16 boost::shared_ptr<MatrixDouble> m_D_ptr);
17 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
18 EntitiesFieldData::EntData &col_data);
19
20private:
21 boost::shared_ptr<CommonData> commonDataPtr;
22 boost::shared_ptr<MatrixDouble> mDPtr;
23};
24
25template <int DIM, typename AssemblyDomainEleOp>
28 const std::string row_field_name, const std::string col_field_name,
29 boost::shared_ptr<CommonData> common_data_ptr,
30 boost::shared_ptr<MatrixDouble> m_D_ptr)
31 : AssemblyDomainEleOp(row_field_name, col_field_name,
32 AssemblyDomainEleOp::OPROWCOL),
33 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
34 AssemblyDomainEleOp::sYmm = false;
35}
36
38get_mat_vector_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number<2>) {
40 &mat(2 * rr + 0, 0), &mat(2 * rr + 0, 1), &mat(2 * rr + 0, 2),
41 &mat(2 * rr + 1, 0), &mat(2 * rr + 1, 1), &mat(2 * rr + 1, 2)};
42}
43
45get_mat_vector_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number<3>) {
47 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
48 &mat(3 * rr + 0, 3), &mat(3 * rr + 0, 4), &mat(3 * rr + 0, 5),
49
50 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
51 &mat(3 * rr + 1, 3), &mat(3 * rr + 1, 4), &mat(3 * rr + 1, 5),
52
53 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2),
54 &mat(3 * rr + 2, 3), &mat(3 * rr + 2, 4), &mat(3 * rr + 2, 5)};
55}
56
57template <int DIM, typename AssemblyDomainEleOp>
58MoFEMErrorCode
60 iNtegrate(EntitiesFieldData::EntData &row_data,
61 EntitiesFieldData::EntData &col_data) {
63
64 FTensor::Index<'i', DIM> i;
65 FTensor::Index<'j', DIM> j;
66 FTensor::Index<'k', DIM> k;
67 FTensor::Index<'l', DIM> l;
68 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
70
71 auto &locMat = AssemblyDomainEleOp::locMat;
72
73 const size_t nb_integration_pts = row_data.getN().size1();
74 const size_t nb_row_base_functions = row_data.getN().size2();
75
76 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
78
80 t_DL(i, j, L) = t_D(i, j, k, l) * t_L(k, l, L);
81
82 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
83 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
84 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
85 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
86 ++t_w;
87
88 size_t rr = 0;
89 for (; rr != AssemblyDomainEleOp::nbRows / DIM; ++rr) {
90
91 auto t_mat =
93
95 t_tmp(i, L) = (t_DL(i, j, L)) * (alpha * t_row_diff_base(j));
96
97 auto t_col_base = col_data.getFTensor0N(gg, 0);
98 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
99
100 t_mat(i, L) -= (t_col_base * t_tmp(i, L));
101
102 ++t_mat;
103 ++t_col_base;
104 }
105
106 ++t_row_diff_base;
107 }
108
109 for (; rr < nb_row_base_functions; ++rr)
110 ++t_row_diff_base;
111 }
112
114}
115
116template <int DIM, typename AssemblyDomainEleOp>
118 : public AssemblyDomainEleOp {
120 const std::string row_field_name, const std::string col_field_name,
121 boost::shared_ptr<CommonData> common_data_ptr,
122 boost::shared_ptr<MatrixDouble> m_D_ptr);
123 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
124 EntitiesFieldData::EntData &col_data);
125
126private:
127 boost::shared_ptr<CommonData> commonDataPtr;
128 boost::shared_ptr<MatrixDouble> mDPtr;
129};
130
131template <int DIM, typename AssemblyDomainEleOp>
134 const std::string row_field_name, const std::string col_field_name,
135 boost::shared_ptr<CommonData> common_data_ptr,
136 boost::shared_ptr<MatrixDouble> m_D_ptr)
137 : AssemblyDomainEleOp(row_field_name, col_field_name,
138 AssemblyDomainEleOp::OPROWCOL),
139 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
140 AssemblyDomainEleOp::sYmm = false;
141}
142
143template <int DIM, typename AssemblyDomainEleOp>
144MoFEMErrorCode
146 EntitiesFieldData::EntData &row_data,
147 EntitiesFieldData::EntData &col_data) {
149
150 FTensor::Index<'i', DIM> i;
151 FTensor::Index<'j', DIM> j;
152 FTensor::Index<'k', DIM> k;
153 FTensor::Index<'l', DIM> l;
154 FTensor::Index<'m', DIM> m;
155 FTensor::Index<'n', DIM> n;
156
157 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
159
160 auto &locMat = AssemblyDomainEleOp::locMat;
161
162 const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
163 const size_t nb_row_base_functions = row_data.getN().size2();
164
165 auto t_res_flow_dstrain =
166 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
167
168 auto next = [&]() { ++t_res_flow_dstrain; };
169
171 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
173 t_diff_grad(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
174
175 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
176 auto t_row_base = row_data.getFTensor0N();
177 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
178
179 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
180 ++t_w;
182 t_res_tens(L, i, j) =
183 alpha * ((t_L(m, n, L) * (t_res_flow_dstrain(m, n, k, l))) *
184 t_diff_grad(k, l, i, j));
185 next();
186
187 size_t rr = 0;
188 for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
189 auto t_mat =
191 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
192 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; ++cc) {
193 t_mat(L, l) += t_row_base * (t_res_tens(L, l, k) * t_col_diff_base(k));
194 ++t_mat;
195 ++t_col_diff_base;
196 }
197 ++t_row_base;
198 }
199
200 for (; rr < nb_row_base_functions; ++rr)
201 ++t_row_base;
202 }
203
205}
206
207template <int DIM, typename AssemblyDomainEleOp>
210 const std::string row_field_name, const std::string col_field_name,
211 boost::shared_ptr<CommonData> common_data_ptr,
212 boost::shared_ptr<MatrixDouble> m_D_ptr)
213 : AssemblyDomainEleOp(row_field_name, col_field_name,
214 DomainEleOp::OPROWCOL),
215 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
216 AssemblyDomainEleOp::sYmm = false;
217}
218
219template <int DIM, typename AssemblyDomainEleOp>
220MoFEMErrorCode
222 EntitiesFieldData::EntData &row_data,
223 EntitiesFieldData::EntData &col_data) {
225
226 FTensor::Index<'i', DIM> i;
227 FTensor::Index<'j', DIM> j;
228 FTensor::Index<'k', DIM> k;
229 FTensor::Index<'l', DIM> l;
230
231 auto &locMat = AssemblyDomainEleOp::locMat;
232
233 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
234 const auto nb_row_base_functions = row_data.getN().size2();
235
236 auto t_c_dstrain =
237 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
238 auto t_diff_grad_symmetrise = diff_symmetrize(FTensor::Number<DIM>());
239
240 auto next = [&]() { ++t_c_dstrain; };
241
242 auto get_mat_scalar_dvector = [&]() {
243 if constexpr (DIM == 2)
244 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&locMat(0, 0),
245 &locMat(0, 1)};
246 else
248 &locMat(0, 0), &locMat(0, 1), &locMat(0, 2)};
249 };
250
251 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
252 auto t_row_base = row_data.getFTensor0N();
253 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
254 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
255 ++t_w;
256
258 t_res_mat(i, j) =
259 ((t_c_dstrain(k, l)) * t_diff_grad_symmetrise(k, l, i, j));
260 next();
261
262 auto t_mat = get_mat_scalar_dvector();
263 size_t rr = 0;
264 for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
265 const double row_base = alpha * t_row_base;
266 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
267 for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; cc++) {
268 t_mat(i) += row_base * (t_res_mat(i, j) * t_col_diff_base(j));
269 ++t_mat;
270 ++t_col_diff_base;
271 }
272 ++t_row_base;
273 }
274 for (; rr != nb_row_base_functions; ++rr)
275 ++t_row_base;
276 }
277
279}
280
281}; // 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