Operator for linear form, usually to calculate values on right hand side.
139 {
142 constexpr int coordinate_size =
144 constexpr int stretch_size = dim * (dim + 1) / 2;
145 const int nb_integration_pts =
getGaussPts().size2();
149 "Plastic cell %llu has no control-DOF layout",
151
152#ifndef NDEBUG
153 const auto check_matrix =
154 [nb_integration_pts](const boost::shared_ptr<MatrixDouble> &matrix_ptr,
155 const int columns, const char *name) {
157 if (matrix_ptr->size1() != nb_integration_pts ||
158 matrix_ptr->size2() != columns)
160 "Inconsistent %s data in reduced Helmholtz gradient",
161 name);
163 };
168 "rotation-adjoint gradient");
169#endif
170
172 auto t_adjoint_piola =
175 auto t_adjoint_stretch =
178 auto t_adjoint_rotation =
181 auto t_adjoint_rotation_gradient =
184 auto t_reference_p =
dataAtPts->getFTensorApproxP(nb_integration_pts);
185 auto t_elastic_f =
dataAtPts->getFTensorSmallH(nb_integration_pts);
186 auto t_plastic_h =
dataAtPts->getFTensorPlasticH(nb_integration_pts);
187 auto t_plastic_f =
dataAtPts->getFTensorPlasticF(nb_integration_pts);
188 auto t_inv_plastic_f =
189 dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
190 auto t_adjoint_stress_d_p =
191 dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
192 auto t_levi_d_p =
193 dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
194 auto t_rotation_gradient =
195 dataAtPts->getFTensorRotAxisGrad(nb_integration_pts);
197
198 std::array<PetscScalar, coordinate_size> local_gradient{};
203 return std::exp(value);
204 };
205
206 for (int gg = 0; gg != nb_integration_pts; ++gg) {
209 t_eigen_vectors(
i,
j) = t_plastic_h(
i,
j);
211 MB_SUCCESS)
213 "Failed to diagonalise plastic logarithmic stretch while "
214 "assembling its reduced gradient");
215 const int nb_unique = getUniqNb<dim>(t_eigen_values);
216 if (nb_unique < dim)
217 CHKERR sortEigenVals<dim>(t_eigen_values, t_eigen_vectors);
219 t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, nb_unique);
221
224 t_adjoint_rotation_gradient_pull(
k,
j) =
225 t_adjoint_rotation_gradient(
k,
i) * t_inv_plastic_f(
i,
j);
226 t_rotation_gradient_pull(
k,
j) =
227 t_rotation_gradient(
k,
i) * t_inv_plastic_f(
i,
j);
228
229 for (int coordinate = 0; coordinate != coordinate_size; ++coordinate) {
231 t_coordinate_direction(C) = 0.;
232 t_coordinate_direction(coordinate) = 1.;
233 const auto t_direction =
235 t_coordinate_direction);
237 t_d_plastic_f(
i,
j) =
238 t_diff_exp(
i,
j,
k,
l) * t_direction(
k,
l);
240 t_d_inv_plastic_f(
i,
j) =
241 -t_inv_plastic_f(
i,
k) * t_d_plastic_f(
k,
l) *
242 t_inv_plastic_f(
l,
j);
243
244
245
246
247 double residual_derivative_contraction =
248 -t_adjoint_piola(
i,
k) * t_elastic_f(
i,
j) *
250
252 t_d_intermediate_p(
i,
j) =
253 t_reference_p(
i,
k) * t_d_plastic_f(
j,
k) / det_plastic_f;
255 t_d_adjoint_stress(S) =
256 t_adjoint_stress_d_p(
i,
j, S) * t_d_intermediate_p(
i,
j);
257
258
259 residual_derivative_contraction -=
260 det_plastic_f *
261 (t_adjoint_stretch(S) * t_d_adjoint_stress(S));
262
264 t_d_levi(
k) = t_levi_d_p(
k,
i,
j) * t_d_intermediate_p(
i,
j);
265 residual_derivative_contraction -=
266 det_plastic_f * t_adjoint_rotation(
k) * t_d_levi(
k);
267
270 t_d_adjoint_rotation_gradient_pull;
272 t_d_adjoint_rotation_gradient_pull(
k,
j) =
273 t_adjoint_rotation_gradient(
k,
i) *
274 t_d_inv_plastic_f(
i,
j);
275 t_d_rotation_gradient_pull(
k,
j) =
276 t_rotation_gradient(
k,
i) * t_d_inv_plastic_f(
i,
j);
277 residual_derivative_contraction +=
279 (t_d_adjoint_rotation_gradient_pull(
k,
j) *
280 t_rotation_gradient_pull(
k,
j) +
281 t_adjoint_rotation_gradient_pull(
k,
j) *
282 t_d_rotation_gradient_pull(
k,
j));
283 }
284
285
286
287
288 local_gradient[coordinate] -=
289 getMeasure() * t_w * residual_derivative_contraction;
290 }
291
292 ++t_adjoint_piola;
293 ++t_adjoint_stretch;
294 ++t_adjoint_rotation;
295 ++t_adjoint_rotation_gradient;
296 ++t_reference_p;
297 ++t_elastic_f;
298 ++t_plastic_h;
299 ++t_plastic_f;
300 ++t_inv_plastic_f;
301 ++t_adjoint_stress_d_p;
302 ++t_levi_d_p;
303 ++t_rotation_gradient;
304 ++t_w;
305 }
306
308 local_gradient.data(), ADD_VALUES);
310 }
#define FTENSOR_INDEXES(DIM,...)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
boost::function< T(const T)> Fun
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, SPACE_DIM > plasticLogarithmicStretchTensorFromCoordinates(const FTensor::Tensor1< T, plasticLogarithmicStretchCoordinateSize > &t_coordinates)
constexpr int plasticLogarithmicStretchCoordinateSize
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'm', 3 > m
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Public Member Functions inherited from MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator