PoissonOperators.hpp

#ifndef __POISSONOPERTATORS_HPP__
#define __POISSONOPERTATORS_HPP__

#include <stdlib.h>
#include <BasicFiniteElements.hpp>

namespace PoissonOps {

using VolEle = FaceElementForcesAndSourcesCore;
using OpVolEle = FaceElementForcesAndSourcesCore::UserDataOperator;

using FaceEle = MoFEM::EdgeElementForcesAndSourcesCore;
using OpFaceEle = FaceEle::UserDataOperator;

using EntData = DataForcesAndSourcesCore::EntData;

// Operators

FTensor::Index<'i', 2> i;

const double body_source = 1;

struct OpLhsUU : public OpVolEle {
public:
  OpLhsUU(std::string domain_field1, std::string domain_field2)
      : OpVolEle(domain_field1, domain_field2, OpVolEle::OPROWCOL) {
    sYmm = true;
  }

  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type, EntData &row_data,
                        EntData &col_data) {
    MoFEMFunctionBegin;

    const int nb_row_dofs = row_data.getIndices().size();
    const int nb_col_dofs = col_data.getIndices().size();

    if (nb_row_dofs && nb_col_dofs) {
      locLhs.resize(nb_row_dofs, nb_col_dofs, false);
      locLhs.clear();
      const int nb_integration_pts = getGaussPts().size2();
      auto t_w = getFTensor0IntegrationWeight();

      auto t_row_diff_base = row_data.getFTensor1DiffN<2>();

      const double vol = getMeasure();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = t_w * vol;
        for (int rr = 0; rr != nb_row_dofs; ++rr) {
          auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
          for (int cc = 0; cc != nb_col_dofs; ++cc) {
            locLhs(rr, cc) += a * t_row_diff_base(i) * t_col_diff_base(i);
            ++t_col_diff_base;
          } // endFor cc
          ++t_row_diff_base;
        } // endFor rr
        ++t_w;
      } // endFor gg
      CHKERR MatSetValues(getFEMethod()->ksp_B, row_data, col_data,
                          &locLhs(0, 0), ADD_VALUES);
      if (row_side != col_side || row_type != col_type) {
        transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
        noalias(transLocLhs) = trans(locLhs);
        CHKERR MatSetValues(getFEMethod()->ksp_B, col_data, row_data,
                            &transLocLhs(0, 0), ADD_VALUES);
      } // (row_side != col_side || row_type != col_type)
    }   // endIf (nb_row_dofs && nb_col_dofs)

    MoFEMFunctionReturn(0);
  }

private:
  MatrixDouble locLhs, transLocLhs;
};

struct OpRhsU : public OpVolEle {
public:
  OpRhsU(std::string domain_field) : OpVolEle(domain_field, OpVolEle::OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;

    const int nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      locRhs.resize(nb_dofs, false);
      locRhs.clear();
      const int nb_integration_pts = getGaussPts().size2();

      auto t_base = data.getFTensor0N();
      auto t_w = getFTensor0IntegrationWeight();
      const double vol = getMeasure();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = vol * t_w;

        for (int rr = 0; rr != nb_dofs; ++rr) {
          locRhs[rr] += a * (t_base * body_source);
          ++t_base;
        } // endFor rr
        ++t_w;
      } // endFor gg
      CHKERR VecSetOption(getFEMethod()->ksp_f, VEC_IGNORE_NEGATIVE_INDICES,
                          PETSC_TRUE);
      CHKERR VecSetValues(getFEMethod()->ksp_f, data, &*locRhs.begin(),
                          ADD_VALUES);
    } // endIf (nb_dofs)

    MoFEMFunctionReturn(0);
  }

private:
  VectorDouble locRhs;
};

}; // namespace PoissonOps

#endif //__POISSONOPERTATORS_HPP__