MinimalSurfaceElement.hpp

Go to the documentation of this file.

/** \file MinimalSurfaceElement.hpp
  \brief Implementation of minimal surface element
  \ingroup minimal_surface_area
 
  Implementation is based on:
  <https://www.dealii.org/developer/doxygen/deal.II/step_15.html>
 
*/
 
/* 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 __MINIMAL_SURFACE_ELEMENT_HPP__
#define __MINIMAL_SURFACE_ELEMENT_HPP__
 
/** \brief Minimal surface equation
\ingroup minimal_surface_area
 
*/
namespace MinimalSurfaceEquation {
 
/** \brief Implementation of minimal area element
\ingroup minimal_surface_area
 
*/
struct MinimalSurfaceElement {
 
  MoFEM::Interface &mField;
 
  struct SurfaceElement : public FaceElementForcesAndSourcesCore {
    SurfaceElement(MoFEM::Interface &m_field)
        : FaceElementForcesAndSourcesCore(m_field) {}
 
    /** \brief Set integrate rule
     */
    int getRule(int order) { return 2 * (order - 1) + 1; };
  };
 
  /** \brief Bondary problem
 
  \f[
  \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L \; \mathbf{u} = \int_L
  \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L \f]
 
  */
  struct EdgeElement : public EdgeElementForcesAndSourcesCore {
    EdgeElement(MoFEM::Interface &m_field)
        : EdgeElementForcesAndSourcesCore(m_field) {}
 
    /** \brief Set integrate rule
     */
    int getRule(int order) { return 2 * order; };
  };
 
  SurfaceElement feRhs; ///< To calculate right hand side
  SurfaceElement feLhs; ///< To calculate left hand side
  EdgeElement feBcEdge; ///< Used to calculate dofs on boundary
 
  MinimalSurfaceElement(MoFEM::Interface &m_field)
      : mField(m_field), feRhs(m_field), feLhs(m_field), feBcEdge(m_field) {}
 
  /** \brief Keep date between operators
   */
  struct CommonData {
    MatrixDouble gradU;      ///< nb_gauss_pts x 2 gradients at integration pts
    VectorDouble normGradU2; ///< size of nb_gauss_pts, norm of gradient
    VectorDouble
        aN; ///< size of nb_gauss_pts, \f$a_n = \frac{1}{\sqrt{1+\|u_i\|^2}}\f$
    MatrixDouble aNbyGradU;     ///< nb_gauss_pts x 2, \f$a_n \nabla u\f$
    MatrixDouble aNpow3byGradU; ///< nb_gauss_pts x 2, \f$a_n^3 \nabla u\f$
  };
 
  /** \brief Integrate vector on rhs,
 
  \f[
  \mathbf{F} = \int_L \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L
  \f]
 
  */
  struct OpAssmebleBcRhs
      : public EdgeElementForcesAndSourcesCore::UserDataOperator {
    Vec vF;
    OpAssmebleBcRhs(const string field_name, Vec v_f)
        : EdgeElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPROW),
          vF(v_f) {}
 
    /** \brief Function on boundary
      Inherit this class and overload this function to change bc.
    */
    virtual double evalFunction(double x, double y) {
      return sin(2 * M_PI * (x + y));
      // return sin(M_PI*x);
    }
 
    VectorDouble nF;
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      nF.resize(nb_dofs, false);
      nF.clear();
      int nb_gauss_pts = data.getN().size1();
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        double val = getLength() * getGaussPts()(1, gg);
        double x = getCoordsAtGaussPts()(gg, 0);
        double y = getCoordsAtGaussPts()(gg, 1);
        noalias(nF) += val * evalFunction(x, y) * data.getN(gg);
      }
      CHKERR VecSetValues(vF, data.getIndices().size(),
                          &*data.getIndices().data().begin(), &nF[0],
                          ADD_VALUES);
      MoFEMFunctionReturn(0);
    }
  };
 
  /** \brief Integrate vector on lhs,
 
  \f[
  \mathbf{A} = \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L
  \f]
 
  */
  struct OpAssmebleBcLhs
      : public EdgeElementForcesAndSourcesCore::UserDataOperator {
    Mat mA;
    OpAssmebleBcLhs(const string field_name, Mat m_a)
        : EdgeElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPROWCOL),
          mA(m_a) {}
    MatrixDouble K, transK;
    PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data) {
      MoFEMFunctionBegin;
 
      int row_nb_dofs = row_data.getIndices().size();
      if (row_nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      int col_nb_dofs = col_data.getIndices().size();
      if (col_nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      K.resize(row_nb_dofs, col_nb_dofs, false);
      K.clear();
      int nb_gauss_pts = row_data.getN().size1();
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        double val = getLength() * getGaussPts()(1, gg);
        VectorAdaptor row_n = row_data.getN(gg);
        VectorAdaptor col_n = col_data.getN(gg);
        noalias(K) += val * outer_prod(row_n, col_n);
      }
      CHKERR MatSetValues(mA, row_nb_dofs,
                          &*row_data.getIndices().data().begin(), col_nb_dofs,
                          &*col_data.getIndices().data().begin(),
                          &*K.data().begin(), ADD_VALUES);
      // Matrix is symmetric, assemble off diagonal part
      if (row_type != col_type || row_side != col_side) {
        transK.resize(col_nb_dofs, row_nb_dofs, false);
        noalias(transK) = trans(K);
        CHKERR MatSetValues(mA, col_nb_dofs,
                            &*col_data.getIndices().data().begin(), row_nb_dofs,
                            &*row_data.getIndices().data().begin(),
                            &*transK.data().begin(), ADD_VALUES);
      }
 
      MoFEMFunctionReturn(0);
    }
  };
 
  /** \brief Evaluate function values and gradients at Gauss Pts
   */
  struct OpGetDataAtGaussPts
      : public FaceElementForcesAndSourcesCore::UserDataOperator {
    CommonData &commonData;
    OpGetDataAtGaussPts(const string field_name, CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPCOL),
          commonData(common_data) {}
 
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      int nb_gauss_pts = data.getN().size1();
 
      // Element loops over entities, it start from vertices
      if (type == MBVERTEX) {
        // Setting size of common data vectors
        commonData.gradU.resize(nb_gauss_pts, 2);
        commonData.gradU.clear();
      }
 
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        MatrixAdaptor grad_shape_fun = data.getDiffN(gg);
        ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU, gg);
        noalias(grad_at_gauss_pt) +=
            prod(trans(grad_shape_fun), data.getFieldData());
      }
 
      MoFEMFunctionReturn(0);
    }
  };
 
  /** \brief Evaluate function values and gradients at Gauss Pts
   */
  struct OpCalculateCoefficientsAtGaussPts
      : public FaceElementForcesAndSourcesCore::UserDataOperator {
    CommonData &commonData;
    bool rhsAndNotLhs;
    OpCalculateCoefficientsAtGaussPts(const string field_name,
                                      CommonData &common_data,
                                      bool rhs_and_not_lhs)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPCOL),
          commonData(common_data), rhsAndNotLhs(rhs_and_not_lhs) {}
 
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      // Element loops over entities, it start from vertices
      if (type == MBVERTEX) {
        int nb_gauss_pts = data.getN().size1();
        commonData.normGradU2.resize(nb_gauss_pts, false);
        for (int gg = 0; gg != nb_gauss_pts; gg++) {
          ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
                                                           gg);
          commonData.normGradU2[gg] =
              inner_prod(grad_at_gauss_pt, grad_at_gauss_pt);
        }
        commonData.aN.resize(nb_gauss_pts, false);
        for (int gg = 0; gg != nb_gauss_pts; gg++) {
          commonData.aN[gg] = 1. / sqrt(1 + commonData.normGradU2[gg]);
        }
        if (rhsAndNotLhs) {
          // needed at right hand side when residual is calculated
          commonData.aNbyGradU.resize(nb_gauss_pts, 2, false);
          for (int gg = 0; gg != nb_gauss_pts; gg++) {
            ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
                                                             gg);
            for (int rr = 0; rr != 2; rr++) {
              commonData.aNbyGradU(gg, rr) =
                  commonData.aN[gg] * grad_at_gauss_pt[rr];
            }
          }
        } else {
          // needed at left hand side when matrix is calculated
          commonData.aNpow3byGradU.resize(nb_gauss_pts, 2, false);
          for (int gg = 0; gg != nb_gauss_pts; gg++) {
            ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
                                                             gg);
            for (int rr = 0; rr != 2; rr++) {
              commonData.aNpow3byGradU(gg, rr) =
                  pow(commonData.aN[gg], 3) * grad_at_gauss_pt[rr];
            }
          }
        }
      }
 
      // doVerticesRow = true;
      // doEdgesRow = false;
      // doQuadsRow = false;
      // doTrisRow = false;
      // doTetsRow = false;
      // doPrismsRow = false;
 
      MoFEMFunctionReturn(0);
    }
  };
 
  /** \brief Assemble residual
   */
  struct OpAssembleResidaul
      : public FaceElementForcesAndSourcesCore::UserDataOperator {
    CommonData &commonData;
    OpAssembleResidaul(const string field_name, CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPROW),
          commonData(common_data) {}
    VectorDouble rEsidual;
    PetscErrorCode doWork(int side, EntityType type,
                          DataForcesAndSourcesCore::EntData &data) {
      MoFEMFunctionBegin;
 
      int nb_dofs = data.getIndices().size();
      if (nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      rEsidual.resize(nb_dofs, false);
      rEsidual.clear();
      int nb_gauss_pts = data.getN().size1();
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        double val = getGaussPts()(2, gg) * getArea();
        ublas::matrix_row<MatrixDouble> an_by_grad_u(commonData.aNbyGradU, gg);
        noalias(rEsidual) += val * prod(data.getDiffN(gg), an_by_grad_u);
      }
      CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
                          &*data.getIndices().data().begin(), &rEsidual[0],
                          ADD_VALUES);
      MoFEMFunctionReturn(0);
    }
  };
 
  /** \brief Assemble tangent matrix
   */
  struct OpAssembleTangent
      : public FaceElementForcesAndSourcesCore::UserDataOperator {
    CommonData &commonData;
    OpAssembleTangent(const string field_name, CommonData &common_data)
        : FaceElementForcesAndSourcesCore::UserDataOperator(
              field_name, UserDataOperator::OPROWCOL),
          commonData(common_data) {
      sYmm = false; // matrix is not symmetric
    }
    MatrixDouble kTangent;
    VectorDouble auxVec;
    MatrixDouble auxMat;
    PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          DataForcesAndSourcesCore::EntData &row_data,
                          DataForcesAndSourcesCore::EntData &col_data) {
      MoFEMFunctionBegin;
 
      int row_nb_dofs = row_data.getIndices().size();
      if (row_nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      int col_nb_dofs = col_data.getIndices().size();
      if (col_nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
      kTangent.resize(row_nb_dofs, col_nb_dofs, false);
      kTangent.clear();
      auxVec.resize(col_nb_dofs, false);
      auxMat.resize(col_nb_dofs, 2, false);
      int nb_gauss_pts = row_data.getN().size1();
      for (int gg = 0; gg != nb_gauss_pts; gg++) {
        // cerr << "gg " << gg << endl;
        double val = getGaussPts()(2, gg) * getArea();
        MatrixAdaptor row_diff_n = row_data.getDiffN(gg);
        MatrixAdaptor col_diff_n = col_data.getDiffN(gg);
        double an = commonData.aN[gg];
        ublas::matrix_row<MatrixDouble> grad_u(commonData.gradU, gg);
        ublas::matrix_row<MatrixDouble> an_pow3_by_grad_u(
            commonData.aNpow3byGradU, gg);
        // cerr << "grad_u " << grad_u << endl;
        // cerr << "an_pow3_by_grad_u " << an_pow3_by_grad_u << endl;
        noalias(auxVec) = prod(col_diff_n, an_pow3_by_grad_u);
        // cerr << "auxVec " << auxVec << endl;
        noalias(auxMat) = outer_prod(auxVec, grad_u);
        // cerr << "auxMat " << auxMat << endl;
        // cerr << row_diff_n << endl;
        // cerr << col_diff_n << endl;
        // cerr << prod(row_diff_n,trans(val*(an*col_diff_n+auxMat))) << endl;
        noalias(kTangent) +=
            val * prod(row_diff_n, trans(an * col_diff_n - auxMat));
        // cerr << "kTangent " << kTangent << endl;
      }
      CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
                          &*row_data.getIndices().data().begin(), col_nb_dofs,
                          &*col_data.getIndices().data().begin(),
                          &*kTangent.data().begin(), ADD_VALUES);
      MoFEMFunctionReturn(0);
    }
  };
};
} // namespace MinimalSurfaceEquation
 
#endif // __MINIMAL_SURFACE_ELEMENT_HPP__
 
/***************************************************************************/ /**
 * \defgroup minimal_surface_area Minimal surface area
 * \ingroup user_modules
 ******************************************************************************/