L2Rotation.hpp

Go to the documentation of this file.

#ifndef ESHELBIAN_L2_ROTATION_HPP
#define ESHELBIAN_L2_ROTATION_HPP
 
namespace EshelbianPlasticity {
 
/** Decompose a tetrahedral L2 vector into its constant and higher modes.
 * The core evaluator owns total-field interpolation and vector/DM lookup.
 * Used before rotation kinematics; the same linear split applies to directions.
 */
struct OpCalculateL2RotationValues
    : MoFEM::OpCalculateVectorFieldValues_General<3, MatrixDouble> {
  using OpBase = MoFEM::OpCalculateVectorFieldValues_General<3, MatrixDouble>;
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
 
  OpCalculateL2RotationValues(const std::string &field, MatrixPtr total,
                              MatrixPtr constant, MatrixPtr higher,
                              SmartPetscObj<Vec> vector = {})
      : OpBase(field, std::move(total), vector, MBTET),
        constantPtr(std::move(constant)), higherPtr(std::move(higher)) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override {
    MoFEMFunctionBegin;
    CHKERR OpBase::doWork(side, type, data);
    if (type != MBTET)
      MoFEMFunctionReturnHot(0);
    const int points = getGaussPts().size2();
    auto a = MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
        *constantPtr, points)();
    auto q = MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
        *higherPtr, points)();
    constantPtr->clear();
    higherPtr->clear();
    const int dofs = data.getFieldData().size();
    if (!dofs || !points)
      MoFEMFunctionReturnHot(0);
    if (dofs < 3 || dofs % 3)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "L2 rotation requires a complete constant vector block");
    FTensor::Tensor1<double, 3> c0;
    FTENSOR_INDEX(3, i);
    if (dataVec) {
      // The core base retains the evaluated vector coefficients after restoring
      // field data. Reuse them so constant and total values share one lookup.
      auto coefficients = getFTensor1FromPtr<3>(&dotVector[0]);
      c0(i) = coefficients(i);
    } else {
      auto coefficients = data.getFTensor1FieldData<3>();
      c0(i) = coefficients(i);
    }
    const double phi0 = data.getN()(0, 0);
#ifndef NDEBUG
    if (!std::isfinite(phi0) || phi0 == 0.)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Invalid constant L2 rotation basis");
    for (int gg = 0; gg != points; ++gg)
      if (std::abs(data.getN()(gg, 0) - phi0) > 1.e-13)
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "First L2 rotation basis must be constant");
#endif
    auto total = MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
        *dataPtr, points)();
    for (int gg = 0; gg != points; ++gg) {
      a(i) = phi0 * c0(i);
      q(i) = total(i) - a(i);
      ++a;
      ++q;
      ++total;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  MatrixPtr constantPtr, higherPtr;
};
 
/** Assemble a rotation covector with different constant/higher test modes.
 * Used for tangent-state and incremental-adjoint loads; the form base owns
 * local sizing and global assembly, including ghost contribution handling.
 */
template <AssemblyType Assembly>
struct OpAssembleL2RotationLoad
    : FormsIntegrators<VolUserDataOperator>::template Assembly<
          Assembly>::OpBase {
  using OpBase = typename FormsIntegrators<
      VolUserDataOperator>::template Assembly<Assembly>::OpBase;
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  OpAssembleL2RotationLoad(const std::string &field, MatrixPtr constant,
                           MatrixPtr higher)
      : OpBase(field, VolUserDataOperator::OPROW),
        constantPtr(std::move(constant)), higherPtr(std::move(higher)) {}
  MoFEMErrorCode iNtegrate(EntData &row_data) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEX(3, i);
    using View = GetFTensor1FromMatType<3, -1, DL>;
    auto values = std::array{MatrixSizeHelper<View, DL>::get(
                                 *constantPtr, OpBase::nbIntegrationPts)(),
                             MatrixSizeHelper<View, DL>::get(
                                 *higherPtr, OpBase::nbIntegrationPts)()};
    auto basis = row_data.getFTensor0N();
    auto weight = OpBase::getFTensor0IntegrationWeight();
    for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
      const double alpha = weight * OpBase::getMeasure();
      auto f = OpBase::template getNf<3>();
      int rr = 0;
      for (; rr != OpBase::nbRows / 3; ++rr) {
        f(i) += (alpha * basis) * values[rr ? 1 : 0](i);
        ++f;
        ++basis;
      }
      for (; rr != OpBase::nbRowBaseFunctions; ++rr)
        ++basis;
      ++values[0];
      ++values[1];
      ++weight;
    }
    MoFEMFunctionReturn(0);
  }
  MatrixPtr constantPtr, higherPtr;
};
 
} // namespace EshelbianPlasticity
#endif