DataOperators.hpp

Go to the documentation of this file.

/** \file DataOperators.hpp

  Number of structures preforming operations on integration points.

  For example:
  - calculate Jacobian
  - calculate Piola-Transform on shape functions
  - calculate normals and tangent vectors at integration points on triangle

 */

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */

#ifndef __DATAOPERATORS_HPP
#define __DATAOPERATORS_HPP

using namespace boost::numeric;

namespace MoFEM {

/** \brief base operator to do operations at Gauss Pt. level
 * \ingroup mofem_forces_and_sources
 */
struct DataOperator {

  DataOperator(const bool symm = true);

  virtual ~DataOperator() = default;

  /** \brief Operator for bi-linear form, usually to calculate values on left
   * hand side
   */
  virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                                EntityType col_type,
                                DataForcesAndSourcesCore::EntData &row_data,
                                DataForcesAndSourcesCore::EntData &col_data) {
    MoFEMFunctionBeginHot;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "doWork function is not implemented for this operator");
    MoFEMFunctionReturnHot(0);
  }

  virtual MoFEMErrorCode opLhs(DataForcesAndSourcesCore &row_data,
                               DataForcesAndSourcesCore &col_data);

  /** \brief Operator for linear form, usually to calculate values on right hand
   * side
   */
  virtual MoFEMErrorCode doWork(int side, EntityType type,
                                DataForcesAndSourcesCore::EntData &data) {
    MoFEMFunctionBeginHot;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "doWork function is not implemented for this operator");
    MoFEMFunctionReturnHot(0);
  }

  virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data,
                               const bool error_if_no_base = false);

  bool sYmm; ///< If true assume that matrix is symmetric structure

  std::array<bool, MBMAXTYPE>
      doEntities; ///< If true operator is executed for entity.

  // Deprecated variables. Use doEntities instead. I keep them for back
  // compatibility with some older modules. It will be removed in some future.

  bool &doVertices; ///< \deprectaed If false skip vertices
  bool &doEdges;    ///< \deprectaed If false skip edges
  bool &doQuads;    ///< \deprectaed
  bool &doTris;     ///< \deprectaed
  bool &doTets;     ///< \deprectaed
  bool &doPrisms;   ///< \deprectaed

  /**
   * \brief Get if operator uses symmetry of DOFs or not
   *
   * If symmetry is used, only not repeating combinations of entities are
   * looped.  For an example pair of (Vertex, Edge_0) and (Edge_0, Vertex) will
   * calculate the same matrices only transposed. Implementing that this can be
   * exploited by integrating only one pair.
   *
   * @return true if symmetry
   */
  inline bool getSymm() const { return sYmm; }

  /// set if operator is executed taking in account symmetry
  inline void setSymm() { sYmm = true; }

  /// unset if operator is executed for  non symmetric problem
  inline void unSetSymm() { sYmm = false; }

private:
  template <bool Symm>
  inline MoFEMErrorCode opLhs(DataForcesAndSourcesCore &row_data,
                              DataForcesAndSourcesCore &col_data);

  template <bool ErrorIfNoBase>
  inline MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data,
                              const std::array<bool, MBMAXTYPE> &do_entities);
};

/**
 * \brief Transform local reference derivatives of shape function to global
 derivatives

 * \ingroup mofem_forces_and_sources
 */
struct OpSetInvJacH1 : public DataOperator {

  FTensor::Tensor2<double *, 3, 3> tInvJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;

  OpSetInvJacH1(MatrixDouble3by3 &inv_jac)
      : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
                &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
                &inv_jac(2, 2)) {}

  MatrixDouble diffNinvJac;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/**
 * \brief brief Transform local reference derivatives of shape function to
 global derivatives

 * \ingroup mofem_forces_and_sources
 */
struct OpSetInvJacHdivAndHcurl : public DataOperator {

  FTensor::Tensor2<double *, 3, 3> tInvJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetInvJacHdivAndHcurl(MatrixDouble3by3 &inv_jac)
      : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
                &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
                &inv_jac(2, 2)) {}

  MatrixDouble diffHdivInvJac;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/**
 * \brief transform local reference derivatives of shape function to global
 derivatives if higher order geometry is given

 * \ingroup mofem_forces_and_sources
*/
struct OpSetHoInvJacH1 : public DataOperator {

  MatrixDouble &invHoJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  OpSetHoInvJacH1(MatrixDouble &inv_ho_jac) : invHoJac(inv_ho_jac) {}

  MatrixDouble diffNinvJac;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/**
 * \brief transform local reference derivatives of shape function to global
 derivatives if higher order geometry is given
 *

 * \ingroup mofem_forces_and_sources
*/
struct OpSetHoInvJacHdivAndHcurl : public DataOperator {

  MatrixDouble &invHoJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetHoInvJacHdivAndHcurl(MatrixDouble &inv_ho_jac) : invHoJac(inv_ho_jac) {}

  MatrixDouble diffHdivInvJac;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief apply contravariant (Piola) transfer to Hdiv space

Contravariant Piola transformation
\f[
\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
=
\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
\f]

* \ingroup mofem_forces_and_sources

*/
struct OpSetContravariantPiolaTransform : public DataOperator {

  double &vOlume;

  FTensor::Tensor2<double *, 3, 3> tJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetContravariantPiolaTransform(double &volume, MatrixDouble3by3 &jac)
      : vOlume(volume),
        // jAc(jac),
        tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
             &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}

  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief Apply contravariant (Piola) transfer to Hdiv space for HO geometr

* \ingroup mofem_forces_and_sources
*/
struct OpSetHoContravariantPiolaTransform : public DataOperator {

  VectorDouble &detHoJac;
  MatrixDouble &hoJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetHoContravariantPiolaTransform(VectorDouble &det_jac, MatrixDouble &jac)
      : detHoJac(det_jac), hoJac(jac) {}

  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief Apply covariant (Piola) transfer to Hcurl space for HO geometry
 */
struct OpSetHoCovariantPiolaTransform : public DataOperator {

  MatrixDouble &hoInvJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetHoCovariantPiolaTransform(MatrixDouble &inv_jac) : hoInvJac(inv_jac) {}

  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief apply covariant transfer to Hcurl space

Contravariant Piola transformation
\f[
\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
=
\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
\f]


* \ingroup mofem_forces_and_sources
*/
struct OpSetCovariantPiolaTransform : public DataOperator {

  FTensor::Tensor2<double *, 3, 3> tInvJac;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;

  OpSetCovariantPiolaTransform(MatrixDouble3by3 &inv_jac)
      : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
                &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
                &inv_jac(2, 2)) {}

  MatrixDouble piolaN;
  MatrixDouble piolaDiffN;

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief Get field values and gradients at Gauss points
 * \ingroup mofem_forces_and_sources
 */
template <int RANK, int DIM> struct OpGetDataAndGradient : public DataOperator {

  MatrixDouble &dataAtGaussPts;
  MatrixDouble &dataGradAtGaussPts;

  OpGetDataAndGradient(MatrixDouble &data_at_gauss_pt,
                       MatrixDouble &data_grad_at_gauss_pt)
      : dataAtGaussPts(data_at_gauss_pt),
        dataGradAtGaussPts(data_grad_at_gauss_pt) {}

  /**
   * Return tensor associated with matrix storing values
   */
  template <int R>
  FTensor::Tensor1<double *, R> getValAtGaussPtsTensor(MatrixDouble &data) {
    THROW_MESSAGE("Not implemented");
  }

  /**
   * Return tensor associated with matrix storing gradient values
   */
  template <int R, int D>
  FTensor::Tensor2<double *, R, D> getGradAtGaussPtsTensor(MatrixDouble &data) {
    THROW_MESSAGE("Not implemented");
  }

  /**
   * \brief Calculate gradient and values at integration points
   * @param  side side of entity on element
   * @param  type type of entity
   * @param  data data stored on entity (dofs values, dofs indices, etc.)
   * @return      error code
   */
  MoFEMErrorCode calculateValAndGrad(int side, EntityType type,
                                     DataForcesAndSourcesCore::EntData &data) {
    MoFEMFunctionBeginHot;
    const int nb_base_functions = data.getN().size2();
    bool constant_diff = false;
    if (type == MBVERTEX && data.getDiffN().size1() * data.getDiffN().size2() ==
                                DIM * nb_base_functions) {
      constant_diff = true;
    }
    const int nb_dofs = data.getFieldData().size();
    for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
      double *data_ptr, *n_ptr, *diff_n_ptr;
      n_ptr = &data.getN()(gg, 0);
      if (constant_diff) {
        diff_n_ptr = &data.getDiffN()(0, 0);
      } else {
        diff_n_ptr = &data.getDiffN()(gg, 0);
      }
      data_ptr = &*data.getFieldData().data().begin();
      for (int rr = 0; rr < RANK; rr++, data_ptr++) {
        dataAtGaussPts(gg, rr) +=
            cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
        double *diff_n_ptr2 = diff_n_ptr;
        for (unsigned int dd = 0; dd < DIM; dd++, diff_n_ptr2++) {
          dataGradAtGaussPts(gg, DIM * rr + dd) +=
              cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
        }
      }
    }
    MoFEMFunctionReturnHot(0);
  }

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {

    MoFEMFunctionBegin;

    if (data.getFieldData().size() == 0) {
      MoFEMFunctionReturnHot(0);
    }

    unsigned int nb_dofs = data.getFieldData().size();
    if (nb_dofs == 0)
      MoFEMFunctionReturnHot(0);

    if (nb_dofs % RANK != 0) {
      SETERRQ4(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
               type, side, nb_dofs, RANK);
    }
    if (nb_dofs / RANK > data.getN().size2()) {
      std::cerr << side << " " << type << " "
                << ApproximationBaseNames[data.getBase()] << std::endl;
      std::cerr << data.getN() << std::endl;
      SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
               nb_dofs, data.getN().size2());
    }

    if (type == MBVERTEX) {
      dataAtGaussPts.resize(data.getN().size1(), RANK, false);
      dataGradAtGaussPts.resize(data.getN().size1(), RANK * DIM, false);
      dataAtGaussPts.clear();
      dataGradAtGaussPts.clear();
    }

    CHKERR calculateValAndGrad(side, type, data);

    MoFEMFunctionReturn(0);
  }
};

/**
 * \brief Specialization for field with 3 coefficients in 3 dimension
 */
template <>
template <>
FTensor::Tensor1<double *, 3>
OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(MatrixDouble &data);

/**
 * \brief Specialization for field with 3 coefficients in 3 dimension
 */
template <>
template <>
FTensor::Tensor2<double *, 3, 3>
OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data);

/**
 * \brief Specialization for field with 3 coefficients in 3 dimension
 */
template <>
MoFEMErrorCode OpGetDataAndGradient<3, 3>::calculateValAndGrad(
    int side, EntityType type, DataForcesAndSourcesCore::EntData &data);

/**
 * \brief Specialization for field with for scalar field in 3 dimension
 */
template <>
MoFEMErrorCode OpGetDataAndGradient<1, 3>::calculateValAndGrad(
    int side, EntityType type, DataForcesAndSourcesCore::EntData &data);

/** \brief Calculate normals at Gauss points of triangle element
 * \ingroup mofem_forces_and_source
 */
struct OpGetCoordsAndNormalsOnFace : public DataOperator {

  MatrixDouble &cOords_at_GaussPt;
  MatrixDouble &nOrmals_at_GaussPt;
  MatrixDouble &tAngent1_at_GaussPt;
  MatrixDouble &tAngent2_at_GaussPt;

  OpGetCoordsAndNormalsOnFace(MatrixDouble &coords_at_gausspt,
                              MatrixDouble &normals_at_gausspt,
                              MatrixDouble &tangent1_at_gausspt,
                              MatrixDouble &tangent2_at_gausspt)
      : cOords_at_GaussPt(coords_at_gausspt),
        nOrmals_at_GaussPt(normals_at_gausspt),
        tAngent1_at_GaussPt(tangent1_at_gausspt),
        tAngent2_at_GaussPt(tangent2_at_gausspt) {}

  MatrixDouble sPin;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);

  MoFEMErrorCode calculateNormals();
};

/** \brief calculate normals at Gauss points of triangle element
 * \ingroup mofem_forces_and_sources
 */
struct OpGetCoordsAndNormalsOnPrism : public DataOperator {

  MatrixDouble &cOords_at_GaussPtF3;
  MatrixDouble &nOrmals_at_GaussPtF3;
  MatrixDouble &tAngent1_at_GaussPtF3;
  MatrixDouble &tAngent2_at_GaussPtF3;
  MatrixDouble &cOords_at_GaussPtF4;
  MatrixDouble &nOrmals_at_GaussPtF4;
  MatrixDouble &tAngent1_at_GaussPtF4;
  MatrixDouble &tAngent2_at_GaussPtF4;

  OpGetCoordsAndNormalsOnPrism(
      MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3,
      MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3,
      MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4,
      MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
      : cOords_at_GaussPtF3(coords_at_gaussptf3),
        nOrmals_at_GaussPtF3(normals_at_gaussptf3),
        tAngent1_at_GaussPtF3(tangent1_at_gaussptf3),
        tAngent2_at_GaussPtF3(tangent2_at_gaussptf3),
        cOords_at_GaussPtF4(coords_at_gaussptf4),
        nOrmals_at_GaussPtF4(normals_at_gaussptf4),
        tAngent1_at_GaussPtF4(tangent1_at_gaussptf4),
        tAngent2_at_GaussPtF4(tangent2_at_gaussptf4) {}

  MatrixDouble sPin;
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);

  MoFEMErrorCode calculateNormals();
};

/** \brief transform Hdiv base fluxes from reference element to physical
 * triangle \ingroup mofem_forces_and_sources
 */
struct OpSetContravariantPiolaTransformOnFace : public DataOperator {

  const VectorDouble
      *normalRawPtr; //< Normal of the element for linear geometry
  const MatrixDouble *normalsAtGaussPtsRawPtr; //< Normals at integration points
                                               //for nonlinear geometry

  /**
   * @brief Shift in vector for linear geometry
   *
   * Normal can have size larger than three, for example normal for contact
   * prism and flat prims element have six comonents, for top and bottom
   * triangle of the prims.
   *
   */
  int normalShift;

  OpSetContravariantPiolaTransformOnFace(const VectorDouble &normal,
                                         const MatrixDouble &normals_at_pts,
                                         const int normal_shift = 0)
      : normalRawPtr(&normal), normalsAtGaussPtsRawPtr(&normals_at_pts),
        normalShift(normal_shift) {}

  OpSetContravariantPiolaTransformOnFace(
      const VectorDouble *normal_raw_ptr = nullptr,
      const MatrixDouble *normals_at_pts_ptr = nullptr,
      const int normal_shift = 0)
      : normalRawPtr(normal_raw_ptr),
        normalsAtGaussPtsRawPtr(normals_at_pts_ptr), normalShift(normal_shift) {
  }

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief transform Hcurl base fluxes from reference element to physical
 * triangle \ingroup mofem_forces_and_sources
 */
struct OpSetCovariantPiolaTransformOnFace : public DataOperator {

  const VectorDouble &nOrmal;
  const MatrixDouble &normalsAtGaussPts;
  const VectorDouble &tAngent0;
  const MatrixDouble &tangent0AtGaussPt;
  const VectorDouble &tAngent1;
  const MatrixDouble &tangent1AtGaussPt;

  OpSetCovariantPiolaTransformOnFace(const VectorDouble &normal,
                                     const MatrixDouble &normals_at_pts,
                                     const VectorDouble &tangent0,
                                     const MatrixDouble &tangent0_at_pts,
                                     const VectorDouble &tangent1,
                                     const MatrixDouble &tangent1_at_pts)
      : nOrmal(normal), normalsAtGaussPts(normals_at_pts), tAngent0(tangent0),
        tangent0AtGaussPt(tangent0_at_pts), tAngent1(tangent1),
        tangent1AtGaussPt(tangent1_at_pts) {}

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief Calculate tangent vector on edge form HO geometry approximation
 * \ingroup mofem_forces_and_sources
 */
struct OpGetHoTangentOnEdge : public DataOperator {

  MatrixDouble &tAngent;

  OpGetHoTangentOnEdge(MatrixDouble &tangent) : tAngent(tangent) {}

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

/** \brief transform Hcurl base fluxes from reference element to physical edge
 * \ingroup mofem_forces_and_sources
 */
struct OpSetCovariantPiolaTransformOnEdge : public DataOperator {

  const VectorDouble &tAngent;
  const MatrixDouble &tangentAtGaussPt;

  OpSetCovariantPiolaTransformOnEdge(const VectorDouble &tangent,
                                     const MatrixDouble &tangent_at_pts)
      : tAngent(tangent), tangentAtGaussPt(tangent_at_pts) {}

  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
};

} // namespace MoFEM

#endif //__DATAOPERATORS_HPP

/**
 * \defgroup mofem_forces_and_sources Forces and sources
 **/