Operator the left hand side matrix. More...

#include "tutorials/scl-11/src/PoissonDiscontinousGalerkin.hpp"

Inheritance diagram for Poisson2DiscontGalerkinOperators::OpL2LhsPenalty:

Collaboration diagram for Poisson2DiscontGalerkinOperators::OpL2LhsPenalty:

Public Member Functions
	OpL2LhsPenalty (boost::shared_ptr< FaceSideEle > side_fe_ptr)
	Construct a new OpL2LhsPenalty.

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

Private Attributes
boost::shared_ptr< FaceSideEle >	sideFEPtr
	pointer to element to get data on edge/face sides

MatrixDouble	locMat
	local operator matrix

Detailed Description

Operator the left hand side matrix.

Examples: poisson_2d_dis_galerkin.cpp.

Definition at line 115 of file PoissonDiscontinousGalerkin.hpp.

Constructor & Destructor Documentation

◆ OpL2LhsPenalty()

Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::OpL2LhsPenalty ( boost::shared_ptr< FaceSideEle > side_fe_ptr )

inline

Construct a new OpL2LhsPenalty.

Parameters

side_fe_ptr pointer to FE to evaluate side elements

Examples: PoissonDiscontinousGalerkin.hpp.

Definition at line 123 of file PoissonDiscontinousGalerkin.hpp.

124 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr) {}

BoundaryEleOp

BoundaryEle::UserDataOperator BoundaryEleOp

Definition child_and_parent.cpp:40

NOSPACE

@ NOSPACE

Definition definitions.h:83

Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::sideFEPtr

boost::shared_ptr< FaceSideEle > sideFEPtr

pointer to element to get data on edge/face sides

Definition PoissonDiscontinousGalerkin.hpp:266

Member Function Documentation

◆ doWork()

MoFEMErrorCode Poisson2DiscontGalerkinOperators::OpL2LhsPenalty::doWork	(	int	side,
		EntityType	type,
		EntitiesFieldData::EntData &	data )

inline

Examples: PoissonDiscontinousGalerkin.hpp.

Definition at line 126 of file PoissonDiscontinousGalerkin.hpp.

                                                        {
    MoFEMFunctionBegin;
 
    // Collect data from side domian elements
    CHKERR loopSideFaces("dFE", *sideFEPtr);
    const auto in_the_loop =
        sideFEPtr->nInTheLoop; // return number of elements on the side
 
#ifndef NDEBUG
    const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpL2LhsPenalty inTheLoop " << ele_type_name[in_the_loop];
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpL2LhsPenalty sense " << senseMap[0] << " " << senseMap[1];
#endif
 
    // calculate  penalty
    const double s = getMeasure() / (areaMap[0] + areaMap[1]);
    const double p = penalty * s;
 
    // get normal of the face or edge
    auto t_normal = getFTensor1Normal();
    t_normal.normalize();
 
    // get number of integration points on face
    const size_t nb_integration_pts = getGaussPts().size2();
 
    // beta paramter is zero, when penalty method is used, also, takes value 1,
    // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
    const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
 
    // iterate the sides rows
    for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
 
      // gent number of DOFs on the right side. 
      const auto nb_rows = indicesRowSideMap[s0].size();
 
      if (nb_rows) {
 
        // get orientation of the local element edge
        const auto sense_row = senseMap[s0];
 
        // iterate the side cols
        const auto nb_row_base_functions = rowBaseSideMap[s0].size2();
        for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
 
          // get orientation of the edge
          const auto sense_col = senseMap[s1];
 
          // number of dofs, for homogenous approximation this value is
          // constant.
          const auto nb_cols = indicesColSideMap[s1].size();
 
          // resize local element matrix
          locMat.resize(nb_rows, nb_cols, false);
          locMat.clear();
 
          // get base functions, and integration weights
          auto t_row_base = get_ntensor(rowBaseSideMap[s0]);
          auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[s0]);
          auto t_w = getFTensor0IntegrationWeight();
 
          // iterate integration points on face/edge
          for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
            // t_w is integration weight, and measure is element area, or
            // volume, depending if problem is in 2d/3d.
            const double alpha = getMeasure() * t_w;
            auto t_mat = locMat.data().begin();
            
            // iterate rows
            size_t rr = 0;
            for (; rr != nb_rows; ++rr) {
 
              // calculate tetting function 
              FTensor::Tensor1<double, SPACE_DIM> t_vn_plus;
              t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
              FTensor::Tensor1<double, SPACE_DIM> t_vn;
              t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
 
              // get base functions on columns
              auto t_col_base = get_ntensor(colBaseSideMap[s1], gg, 0);
              auto t_diff_col_base =
                  get_diff_ntensor(colDiffBaseSideMap[s1], gg, 0);
 
              // iterate columns
              for (size_t cc = 0; cc != nb_cols; ++cc) {
 
                // calculate variance of tested function 
                FTensor::Tensor1<double, SPACE_DIM> t_un;
                t_un(i) = -p * (t_col_base * t_normal(i) * sense_col -
                                beta * t_diff_col_base(i) / p);
 
                // assemble matrix
                *t_mat -= alpha * (t_vn(i) * t_un(i));
                *t_mat -= alpha * (t_vn_plus(i) * (beta * t_diff_col_base(i)));
 
                // move to next column base and element of matrix
                ++t_col_base;
                ++t_diff_col_base;
                ++t_mat;
              }
 
              // move to next row base
              ++t_row_base;
              ++t_diff_row_base;
            }
 
            // this is obsolete for this particular example, we keep it for
            // generality. in case of multi-physics problems different fields can
            // chare hierarchical base but use different approx. order, so is
            // possible to have more base functions than DOFs on element.
            for (; rr < nb_row_base_functions; ++rr) {
              ++t_row_base;
              ++t_diff_row_base;
            }
 
            ++t_w;
          }
 
          // assemble system
          CHKERR ::MatSetValues(getKSPB(), indicesRowSideMap[s0].size(),
                                &*indicesRowSideMap[s0].begin(),
                                indicesColSideMap[s1].size(),
                                &*indicesColSideMap[s1].begin(),
                                &*locMat.data().begin(), ADD_VALUES);
 
          // stop of boundary element
          if (!in_the_loop)
            MoFEMFunctionReturnHot(0);
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }