adJointImpl.hpp

Go to the documentation of this file.

/** \file adJointImpl.hpp
  * \brief Operators managing HO geometry derivatives
 
*/
 
#ifndef __ADJOINT_IMPL_HPP__
#define __ADJOINT_IMPL_HPP__
 
namespace MoFEM {
 
//! [declarations of templates]
 
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpGetCoFactorImpl;
 
//! [declarations of templates]
 
//! [definitions of templates]
 
/**
 * @brief Implementation of cofactor matrix derivative operator for Gauss integration
 * 
 * This operator computes the derivative of the cofactor matrix with respect to 
 * geometry parameters. The cofactor matrix is fundamental in finite element
 * transformations and sensitivity analysis.
 * 
 * Mathematical background:
 * For a Jacobian matrix \f$\mathbf{J}\f$, the cofactor matrix is:
 * \f[
 * \text{cof}(\mathbf{J}) = \det(\mathbf{J}) \cdot \mathbf{J}^{-T}
 * \f]
 * 
 * The derivative computation involves:
 * \f[
 * \frac{d}{d\alpha}\text{cof}(\mathbf{J}) = \text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha}) \cdot \det(\mathbf{J})
 * \f]
 * 
 * This is implemented as: `(inv_jac(j,i) * diff_jac(i,j)) * det`
 * where the Einstein summation convention is used for repeated indices.
 */
template <int SPACE_DIM, typename OpBase>
struct OpGetCoFactorImpl<SPACE_DIM, IntegrationType::GAUSS, OpBase>
    : public OpBase {
 
  OpGetCoFactorImpl(boost::shared_ptr<MatrixDouble> jac_ptr,
                    boost::shared_ptr<MatrixDouble> diff_jac_ptr,
                    boost::shared_ptr<VectorDouble> diff_out_ptr
 
                    )
      : OpBase(NOSPACE, OpBase::OPSPACE), jacPtr(jac_ptr),
        diffJacPtr(diff_jac_ptr), diffOutPtr(diff_out_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
protected:
  boost::shared_ptr<MatrixDouble> jacPtr;      ///< Pointer to Jacobian matrix
  boost::shared_ptr<MatrixDouble> diffJacPtr;  ///< Pointer to Jacobian derivative matrix
  boost::shared_ptr<VectorDouble> diffOutPtr; ///< Pointer to output derivative vector
};
 
//! [definitions of templates]
 
//! [implementations of templates functions]
 
/**
 * @brief Compute cofactor derivative at integration points
 * 
 * This function implements the cofactor derivative computation:
 * \f[
 * \frac{d}{d\alpha}\text{cof}(\mathbf{J}) = \text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha}) \cdot \det(\mathbf{J})
 * \f]
 * 
 * The computation steps are:
 * 1. Compute determinant: \f$\det(\mathbf{J})\f$
 * 2. Compute inverse Jacobian: \f$\mathbf{J}^{-1}\f$
 * 3. Compute trace: \f$\text{trace}(\mathbf{J}^{-T} \cdot \frac{d\mathbf{J}}{d\alpha})\f$
 * 4. Multiply: \f$\text{trace} \times \det(\mathbf{J})\f$
 * 
 * @param side Side of the element
 * @param type Entity type
 * @param data Field data
 * @return MoFEMErrorCode
 */
template <int SPACE_DIM, typename OpBase>
MoFEMErrorCode
OpGetCoFactorImpl<SPACE_DIM, IntegrationType::GAUSS, OpBase>::doWork(
    int side, EntityType type, EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  auto nb_integration_pts = OpBase::getGaussPts().size2();
  diffOutPtr->resize(nb_integration_pts, false);
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*jacPtr);
  auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*diffJacPtr);
  auto t_diff_out = getFTensor0FromVec(*diffOutPtr);
  for (int gg = 0; gg != nb_integration_pts; gg++) {
    // Compute determinant of Jacobian
    auto t_det = determinantTensor(t_jac);
    // Compute inverse Jacobian matrix
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
    CHKERR invertTensor(t_jac, t_det, t_inv_jac);
    // Compute cofactor derivative: trace(J^{-T} * dJ/dalpha) * det(J)
    // Note: t_inv_jac(j,i) gives J^{-T}, t_diff_jac(i,j) gives dJ/dalpha
    t_diff_out = (t_inv_jac(j, i) * t_diff_jac(i, j)) * t_det; 
    ++t_jac;
    ++t_diff_jac;
    ++t_diff_out;
  }
  MoFEMFunctionReturn(0);
}
 
//! [implementations of templates functions]
 
}; // namespace MoFEM
 
#endif