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. More...

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 More...

MatrixDouble	locMat
	local operator matrix More...

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) {}

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);
   }