Operator to evaluate Dirichlet boundary conditions using DG. More...

#include <users_modules/tutorials/scl-11/src/PoissonDiscontinousGalerkin.hpp>

Inheritance diagram for Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs:

Collaboration diagram for Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs:

Public Member Functions
	OpL2BoundaryRhs (boost::shared_ptr< FaceSideEle > side_fe_ptr, ScalarFun fun)

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

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

ScalarFun	sourceFun
	pointer to function to evaluate value of function on boundary More...

VectorDouble	rhsF
	vector to store local operator right hand side More...

Detailed Description

Operator to evaluate Dirichlet boundary conditions using DG.

Examples: poisson_2d_dis_galerkin.cpp.

Definition at line 274 of file PoissonDiscontinousGalerkin.hpp.

Constructor & Destructor Documentation

◆ OpL2BoundaryRhs()

Poisson2DiscontGalerkinOperators::OpL2BoundaryRhs::OpL2BoundaryRhs	(	boost::shared_ptr< FaceSideEle >	side_fe_ptr,
		ScalarFun	fun
	)

inline

Examples: PoissonDiscontinousGalerkin.hpp.

Definition at line 276 of file PoissonDiscontinousGalerkin.hpp.

277 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr),

278 sourceFun(fun) {}

Member Function Documentation

◆ doWork()

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

inline

Examples: PoissonDiscontinousGalerkin.hpp.

Definition at line 280 of file PoissonDiscontinousGalerkin.hpp.

                                                                   {
     MoFEMFunctionBegin;
  
     // get normal of the face or edge
     CHKERR loopSideFaces("dFE", *sideFEPtr);
     const auto in_the_loop =
         sideFEPtr->nInTheLoop; // return number of elements on the side
  
     // calculate  penalty
     const double s = getMeasure() / (areaMap[0]);
     const double p = penalty * s;
  
     // get normal of the face or edge
     auto t_normal = getFTensor1Normal();
     t_normal.normalize();
  
     auto t_w = getFTensor0IntegrationWeight();
  
     const size_t nb_rows = indicesRowSideMap[0].size();
  
     if (!nb_rows)
       MoFEMFunctionReturnHot(0);
  
     // resize local operator vector
     rhsF.resize(nb_rows, false);
     rhsF.clear();
  
     const size_t nb_integration_pts = getGaussPts().size2();
     const size_t nb_row_base_functions = rowBaseSideMap[0].size2();
  
     // shape funcs
     auto t_row_base = get_ntensor(rowBaseSideMap[0]);
     auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[0]);
     auto t_coords = getFTensor1CoordsAtGaussPts();
  
     const auto sense_row = senseMap[0];
     const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
  
     for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
       const double alpha = getMeasure() * t_w;
  
       const double source_val =
           -p * sourceFun(t_coords(0), t_coords(1), t_coords(2));
  
       auto t_f = rhsF.data().begin();
  
       size_t rr = 0;
       for (; rr != nb_rows; ++rr) {
  
         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);
  
         // assemble operator local vactor
         *t_f -= alpha * t_vn(i) * (source_val * t_normal(i));
  
         ++t_row_base;
         ++t_diff_row_base;
         ++t_f;
       }
  
       for (; rr < nb_row_base_functions; ++rr) {
         ++t_row_base;
         ++t_diff_row_base;
       }
  
       ++t_w;
       ++t_coords;
     }
  
     // assemble local operator vector to global vector
     CHKERR ::VecSetValues(getKSPf(), indicesRowSideMap[0].size(),
                           &*indicesRowSideMap[0].begin(), &*rhsF.data().begin(),
                           ADD_VALUES);
  
     MoFEMFunctionReturn(0);
   }