poisson_2d_homogeneous.hpp

Go to the documentation of this file.

/**
 * \file poisson_2d_homogeneous.hpp
 * \example poisson_2d_homogeneous.hpp
 *
 * Solution of poisson equation. Direct implementation of User Data Operators
 * for teaching proposes.
 *
 * \note In practical application we suggest use form integrators to generalise
 * and simplify code. However, here we like to expose user to ways how to
 * implement data operator from scratch.
 */
 
// Define name if it has not been defined yet
#ifndef __POISSON_2D_HOMOGENEOUS_HPP__
#define __POISSON_2D_HOMOGENEOUS_HPP__
#include <MoFEM.hpp>
using namespace MoFEM;
// Use of alias for some specific functions
// We are solving Poisson's equation in 2D so Face element is used
using EntData = EntitiesFieldData::EntData;
 
template <int DIM>
using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
 
using AssemblyDomainEleOp =
    FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
 
// Namespace that contains necessary UDOs, will be included in the main program
namespace Poisson2DHomogeneousOperators {
 
// Declare FTensor index for 2D problem
FTensor::Index<'i', SPACE_DIM> i;
 
//function type
typedef boost::function<double(const double, const double, const double)>
    ScalarFunc;
 
struct OpDomainLhsMatrixK : public AssemblyDomainEleOp {
public:
  OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
      : AssemblyDomainEleOp(row_field_name, col_field_name,
                            DomainEleOp::OPROWCOL) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const int nb_row_dofs = row_data.getIndices().size();
    const int nb_col_dofs = col_data.getIndices().size();
 
    this->locMat.resize(nb_row_dofs, nb_col_dofs, false);
    this->locMat.clear();
 
    // get element area
    const double area = getMeasure();
 
    // get number of integration points
    const int nb_integration_points = getGaussPts().size2();
    // get integration weights
    auto t_w = getFTensor0IntegrationWeight();
 
    // get derivatives of base functions on row
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
 
    // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
    for (int gg = 0; gg != nb_integration_points; gg++) {
      const double a = t_w * area;
 
      for (int rr = 0; rr != nb_row_dofs; ++rr) {
        // get derivatives of base functions on column
        auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
        for (int cc = 0; cc != nb_col_dofs; cc++) {
          this->locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
 
          // move to the derivatives of the next base functions on column
          ++t_col_diff_base;
        }
 
        // move to the derivatives of the next base functions on row
        ++t_row_diff_base;
      }
 
      // move to the weight of the next integration point
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
struct OpDomainRhsVectorF : public AssemblyDomainEleOp {
public:
  OpDomainRhsVectorF(std::string field_name, ScalarFunc source_term_function)
      : AssemblyDomainEleOp(field_name, field_name, DomainEleOp::OPROW),
        sourceTermFunc(source_term_function) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const int nb_dofs = data.getIndices().size();
 
    this->locF.resize(nb_dofs, false);
    this->locF.clear();
 
    // get element area
    const double area = getMeasure();
 
    // get number of integration points
    const int nb_integration_points = getGaussPts().size2();
    // get integration weights
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    // get base function
    auto t_base = data.getFTensor0N();
 
    // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
    for (int gg = 0; gg != nb_integration_points; gg++) {
      const double a = t_w * area;
 
      // get source term at the integration point coordinates
      double body_source =
          sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
      for (int rr = 0; rr != nb_dofs; rr++) {
        this->locF[rr] += t_base * body_source * a;
 
        // move to the next base function
        ++t_base;
      }
 
      // move to the weight of the next integration point
      ++t_w;
      // move to the coordinates of the next integration point
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  ScalarFunc sourceTermFunc;
};
 
}; // namespace Poisson2DHomogeneousOperators
 
#endif //__POISSON_2D_HOMOGENEOUS_HPP__