v0.15.0
Loading...
Searching...
No Matches
adJointImpl.hpp
Go to the documentation of this file.
1/** \file adJointImpl.hpp
2 * \brief Operators managing HO geometry derivatives
3
4*/
5
6#ifndef __ADJOINT_IMPL_HPP__
7#define __ADJOINT_IMPL_HPP__
8
9namespace MoFEM {
10
11//! [declarations of templates]
12
13template <int SPACE_DIM, IntegrationType I, typename OpBase>
15
16//! [declarations of templates]
17
18//! [definitions of templates]
19
20/**
21 * @brief Implementation of cofactor matrix derivative operator for Gauss integration
22 *
23 * This operator computes the derivative of the cofactor matrix with respect to
24 * geometry parameters. The cofactor matrix is fundamental in finite element
25 * transformations and sensitivity analysis.
26 *
27 * Mathematical background:
28 * For a Jacobian matrix \f$\mathbf{J}\f$, the cofactor matrix is:
29 * \f[
30 * \text{cof}(\mathbf{J}) = \det(\mathbf{J}) \cdot \mathbf{J}^{-T}
31 * \f]
32 *
33 * The derivative computation involves:
34 * \f[
35 * \frac{d}{d\alpha}\text{cof}(\mathbf{J}) = \text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha}) \cdot \det(\mathbf{J})
36 * \f]
37 *
38 * This is implemented as: `(inv_jac(j,i) * diff_jac(i,j)) * det`
39 * where the Einstein summation convention is used for repeated indices.
40 */
41template <int SPACE_DIM, typename OpBase>
43 : public OpBase {
44
45 OpGetCoFactorImpl(boost::shared_ptr<MatrixDouble> jac_ptr,
46 boost::shared_ptr<MatrixDouble> diff_jac_ptr,
47 boost::shared_ptr<VectorDouble> diff_out_ptr
48
49 )
50 : OpBase(NOSPACE, OpBase::OPSPACE), jacPtr(jac_ptr),
51 diffJacPtr(diff_jac_ptr), diffOutPtr(diff_out_ptr) {}
52
53 MoFEMErrorCode doWork(int side, EntityType type,
55
56protected:
57 boost::shared_ptr<MatrixDouble> jacPtr; ///< Pointer to Jacobian matrix
58 boost::shared_ptr<MatrixDouble> diffJacPtr; ///< Pointer to Jacobian derivative matrix
59 boost::shared_ptr<VectorDouble> diffOutPtr; ///< Pointer to output derivative vector
60};
61
62//! [definitions of templates]
63
64//! [implementations of templates functions]
65
66/**
67 * @brief Compute cofactor derivative at integration points
68 *
69 * This function implements the cofactor derivative computation:
70 * \f[
71 * \frac{d}{d\alpha}\text{cof}(\mathbf{J}) = \text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha}) \cdot \det(\mathbf{J})
72 * \f]
73 *
74 * The computation steps are:
75 * 1. Compute determinant: \f$\det(\mathbf{J})\f$
76 * 2. Compute inverse Jacobian: \f$\mathbf{J}^{-1}\f$
77 * 3. Compute trace: \f$\text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha})\f$
78 * 4. Multiply: \f$\text{trace} \times \det(\mathbf{J})\f$
79 *
80 * @param side Side of the element
81 * @param type Entity type
82 * @param data Field data
83 * @return MoFEMErrorCode
84 */
85template <int SPACE_DIM, typename OpBase>
88 int side, EntityType type, EntitiesFieldData::EntData &data) {
90 auto nb_integration_pts = OpBase::getGaussPts().size2();
91 diffOutPtr->resize(nb_integration_pts, false);
94 auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*jacPtr);
95 auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*diffJacPtr);
96 auto t_diff_out = getFTensor0FromVec(*diffOutPtr);
97 for (int gg = 0; gg != nb_integration_pts; gg++) {
98 // Compute determinant of Jacobian
99 auto t_det = determinantTensor(t_jac);
100 // Compute inverse Jacobian matrix
102 CHKERR invertTensor(t_jac, t_det, t_inv_jac);
103 // Compute cofactor derivative: trace(J^{-T} * dJ/dalpha) * det(J)
104 // Note: t_inv_jac(j,i) gives J^{-T}, t_diff_jac(i,j) gives dJ/dalpha
105 t_diff_out = (t_inv_jac(j, i) * t_diff_jac(i, j)) * t_det;
106 ++t_jac;
107 ++t_diff_jac;
108 ++t_diff_out;
109 }
111}
112
113//! [implementations of templates functions]
114
115}; // namespace MoFEM
116
117#endif
#define FTENSOR_INDEX(DIM, I)
constexpr int SPACE_DIM
@ NOSPACE
Definition definitions.h:83
#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()
#define CHKERR
Inline error check.
IntegrationType
Form integrator integration types.
@ GAUSS
Gaussian quadrature integration.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Data on single entity (This is passed as argument to DataOperator::doWork)
boost::shared_ptr< VectorDouble > diffOutPtr
Pointer to output derivative vector.
boost::shared_ptr< MatrixDouble > jacPtr
Pointer to Jacobian matrix.
boost::shared_ptr< MatrixDouble > diffJacPtr
Pointer to Jacobian derivative matrix.
OpGetCoFactorImpl(boost::shared_ptr< MatrixDouble > jac_ptr, boost::shared_ptr< MatrixDouble > diff_jac_ptr, boost::shared_ptr< VectorDouble > diff_out_ptr)
[declarations of templates]