LEPNC_diffusion.hpp

Go to the documentation of this file.

// Implementation of the Locally Enriched Polytopal Non-Conforming scheme in 2D for the diffusion model
//
//   { - div(K \grad(zeta(u))) = f,       inside Omega
//   { K \grad(zeta(u)) . nTF = g,       on GammaN
//   {                              u = g,           on GammaD
// 
//  At the moment, only pure Neumann or pure Dirichlet
//
// Author: Jerome Droniou (jerome.droniou@monash.edu)
//
 
#ifndef _LEPNC_DIFFUSION_HPP
#define _LEPNC_DIFFUSION_HPP
 
#include <functional>
#include <utility>
#include <iostream>
 
#include <boost/timer/timer.hpp>
 
// Matrices and linear solvers
#include <Eigen/Sparse>
#include <Eigen/Dense>
//#include "Eigen/MA41.cpp"
 
#include "mesh.hpp"
#include "lepnccore.hpp"
#include "quad2d.hpp"
#include "TestCase/TestCase.hpp"
#include "BoundaryConditions/BoundaryConditions.hpp"
 
namespace HArDCore2D {
 
// ----------------------------------------------------------------------------
//                            Class definition
// ----------------------------------------------------------------------------
/*
* If using this code in a scientific publication, please cite the reference for the LEPNC:
*
* Non-conforming finite elements on polytopal mesh,
* J. Droniou, R. Eymard, T. Gallouët and R. Herbin
* url:
*
*/
 
class LEPNC_diffusion { 
 
// Types
public:
  using solution_function_type = std::function<double(double,double)>;      
  using source_function_type = std::function<double(double,double,Cell*)>;      
  using grad_function_type = std::function<VectorRd(double,double,Cell*)>;      
  using tensor_function_type = std::function<Eigen::Matrix2d(double,double,Cell*)>;     
 
  LEPNC_diffusion(
        LEPNCCore &nc,                                      
        tensor_function_type kappa,     
        size_t deg_kappa,                           
        source_function_type source,  
        BoundaryConditions BC,                  
        solution_function_type exact_solution,  
        grad_function_type grad_exact_solution,     
        std::string solver_type,        
    std::ostream & output = std::cout    
        );
 
  Eigen::VectorXd solve();
 
  double EnergyNorm(const Eigen::VectorXd Xh) const; 
 
  double get_assembly_time() const; 
  double get_solving_time() const;  
  double get_solving_error() const; 
  double get_itime(size_t idx) const;       
 
private:
  Eigen::MatrixXd diffusion_operator(const size_t iT) const;
 
  Eigen::VectorXd load_operator(const size_t iT) const;
 
    const LEPNCCore& nc;
  const tensor_function_type kappa;
    size_t _deg_kappa;
  const source_function_type source;
    const BoundaryConditions m_BC;
  const solution_function_type exact_solution;
  const grad_function_type grad_exact_solution;
    const std::string solver_type;
  std::ostream & m_output;
 
    // To store local diffusion and source terms
    std::vector<Eigen::MatrixXd> DiffT;
    std::vector<Eigen::VectorXd> SourceT;
 
  // Computation statistics
  size_t _assembly_time;
  size_t _solving_time;
  double _solving_error;
    mutable std::vector<size_t> _itime = std::vector<size_t>(10, 0);
 
};
 
LEPNC_diffusion::LEPNC_diffusion(LEPNCCore &nc, tensor_function_type kappa, size_t deg_kappa, source_function_type source, BoundaryConditions BC, solution_function_type exact_solution, grad_function_type grad_exact_solution, std::string solver_type, std::ostream & output)
  : nc(nc),
        kappa(kappa),
        _deg_kappa(deg_kappa),
    source(source),
        m_BC(BC),
        exact_solution(exact_solution),
        grad_exact_solution(grad_exact_solution),
        solver_type(solver_type),
    m_output(output) {
 
    m_output << "[LEPNC_diffusion] Initializing" << std::endl;
 
        //---- Compute diffusion, mass matrix and source (do not change during Newton) ------//
        const Mesh* mesh = nc.get_mesh();
        DiffT.resize(mesh->n_cells());
        SourceT.resize(mesh->n_cells());
        for (size_t iT = 0; iT < mesh->n_cells(); iT++) {
            DiffT[iT] = diffusion_operator(iT);
            SourceT[iT] = load_operator(iT);
        }
}
 
// ASSEMBLE AND SOLVE THE SCHEME
Eigen::VectorXd LEPNC_diffusion::solve() {
 
    boost::timer::cpu_timer timer; 
    boost::timer::cpu_timer timerint; 
  const auto mesh = nc.get_mesh();
    size_t n_edges_dofs = mesh->n_edges();
    size_t n_cell_dofs = mesh->n_cells() * DimPoly<Cell>(1);
 
    // Vector of Dirichlet BC, either on u (if weight>0) or zeta(u) (if weight=0)
    Eigen::VectorXd DirBC = Eigen::VectorXd::Zero(n_cell_dofs + n_edges_dofs);
    for (size_t ibF = 0; ibF < mesh->n_b_edges(); ibF++){
    Edge* edge = mesh->b_edge(ibF);
    if (m_BC.type(*edge)=="dir"){
          size_t iF = edge->global_index();
          QuadratureRule quadF = generate_quadrature_rule(*edge, 5);
          for (QuadratureNode quadrule : quadF){
              DirBC(n_cell_dofs + iF) += quadrule.w * exact_solution(quadrule.x, quadrule.y);
          }
          // Take average
          DirBC(n_cell_dofs + iF) /= mesh->b_edge(ibF)->measure();
    }
    }
 
 
    // ---- RESOLUTION ------ //
 
  // System matrix and RHS (only on the edge dofs)
  Eigen::SparseMatrix<double> GlobMat(n_edges_dofs, n_edges_dofs);
  std::vector<Eigen::Triplet<double>> triplets_GlobMat;
  Eigen::VectorXd GlobRHS = Eigen::VectorXd::Zero(n_edges_dofs);
 
  // Static condensation: matrix and source to recover cell unknowns
  Eigen::SparseMatrix<double> ScMat(n_cell_dofs, n_edges_dofs);
  std::vector<Eigen::Triplet<double>> triplets_ScMat;
  Eigen::VectorXd ScRHS = Eigen::VectorXd::Zero(n_cell_dofs);
 
  // Assemble local contributions
  for (size_t iT = 0; iT < mesh->n_cells(); iT++) {
    Cell* cell = mesh->cell(iT);
    size_t n_edgesT = cell->n_edges();
  
    // Local static condensation of element unknowns
    Eigen::MatrixXd MatT = DiffT[iT];
    Eigen::MatrixXd ATT = MatT.topLeftCorner(DimPoly<Cell>(1), DimPoly<Cell>(1));
    Eigen::MatrixXd ATF = MatT.topRightCorner(DimPoly<Cell>(1), n_edgesT);
    Eigen::MatrixXd AFT = MatT.bottomLeftCorner(n_edgesT, DimPoly<Cell>(1));
    Eigen::MatrixXd AFF = MatT.bottomRightCorner(n_edgesT, n_edgesT);
 
    Eigen::PartialPivLU<Eigen::MatrixXd> invATT;
    invATT.compute(ATT);
  
    Eigen::MatrixXd invATT_ATF = invATT.solve(ATF);
    Eigen::VectorXd invATT_source_cell = invATT.solve(SourceT[iT].head(DimPoly<Cell>(1)));
 
    // Local matrix and right-hand side on the face unknowns
    Eigen::MatrixXd MatF = Eigen::MatrixXd::Zero(n_edgesT, n_edgesT);
    Eigen::VectorXd bF = Eigen::VectorXd::Zero(n_edgesT);
    MatF = AFF - AFT * invATT_ATF;
    bF = SourceT[iT].tail(n_edgesT) - AFT * invATT_source_cell;
              
    // Assemble static condensation operator
    ScRHS.segment(iT * DimPoly<Cell>(1), DimPoly<Cell>(1)) = invATT_source_cell;
    for (size_t i = 0; i < DimPoly<Cell>(1); i++){
      size_t iGlobal = iT * DimPoly<Cell>(1) + i;
      for (size_t j = 0; j < n_edgesT; j++){
        size_t jGlobal = cell->edge(j)->global_index();
        triplets_ScMat.emplace_back(iGlobal, jGlobal, invATT_ATF(i, j));
      }
    }
 
    // GLOBAL MATRIX and RHS on the edge unknowns, using the matrices MatF obtained after static condensation
    for (size_t i = 0; i < n_edgesT; i++){
      size_t iGlobal = cell->edge(i)->global_index();
      for (size_t j = 0; j < n_edgesT; j++){
        size_t jGlobal = cell->edge(j)->global_index();
        triplets_GlobMat.emplace_back(iGlobal, jGlobal, MatF(i, j));
      }
        GlobRHS(iGlobal) += bF(i);
    }
  }
 
  // Assemble the global linear system (without BC), and matrix to recover statically-condensed cell dofs
  GlobMat.setFromTriplets(std::begin(triplets_GlobMat), std::end(triplets_GlobMat));
  ScMat.setFromTriplets(std::begin(triplets_ScMat), std::end(triplets_ScMat));
 
  //    Dirichlet boundary conditions
  size_t n_edge_unknowns = mesh->n_edges() - m_BC.n_dir_edges();
  Eigen::SparseMatrix<double> A = GlobMat.topLeftCorner(n_edge_unknowns, n_edge_unknowns);
  Eigen::VectorXd B = GlobRHS.head(n_edge_unknowns);
 
  // Solve condensed system and recover cell unknowns.
  Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
  solver.compute(A);
  Eigen::VectorXd Xh_edges = solver.solve(B);
 
  Eigen::VectorXd Xh = DirBC; 
  Xh.segment(n_cell_dofs, n_edge_unknowns) = Xh_edges;
  Xh.head(n_cell_dofs) = ScRHS - ScMat * Xh.tail(n_edges_dofs);
 
 
  return Xh;
}
 
 
//******************************** 
//      local diffusion matrix 
//********************************
 
Eigen::MatrixXd LEPNC_diffusion::diffusion_operator(const size_t iT) const {
 
  const auto mesh = nc.get_mesh();
    Cell* cell = mesh->cell(iT);
    const size_t n_edgesT = cell->n_edges();
 
  // Total number of degrees of freedom local to this cell
  size_t local_dofs = n_edgesT + DimPoly<Cell>(1);
    
    // Diffusion in the cell.
    std::function<Eigen::Matrix2d(double,double)> kappaT = [&](double x, double y){
            return kappa(x, y, cell);
    };
 
    // Compute cell quadrature nodes and values of cell basis functions and gradients at these nodes
    size_t doe = 2 + _deg_kappa;
    QuadratureRule quadT = generate_quadrature_rule(*cell, doe, true);
    std::vector<Eigen::ArrayXXd> Dnc_phiT_quadT = nc.grad_nc_basis_quad(iT, quadT);
    // Diffusion tensor at the quadrature nodes
    std::vector<Eigen::Matrix2d> kappaT_quadT(quadT.size());
    std::transform(quadT.begin(), quadT.end(), kappaT_quadT.begin(),
            [&kappaT](QuadratureNode qr) -> Eigen::MatrixXd { return kappaT(qr.x, qr.y); });
 
    return nc.gram_matrix(Dnc_phiT_quadT, Dnc_phiT_quadT, local_dofs, local_dofs, quadT, true, kappaT_quadT);  
}
 
 
//******************************** 
//      local load term 
//********************************
 
Eigen::VectorXd LEPNC_diffusion::load_operator(const size_t iT) const {
    // Load for the cell DOFs (first indices) and face DOFs (last indices)
  const auto mesh = nc.get_mesh();
    Cell* cell = mesh->cell(iT);
    size_t n_edgesT = cell->n_edges();
    size_t local_dofs = DimPoly<Cell>(1) + n_edgesT;
  Eigen::VectorXd b = Eigen::VectorXd::Zero(local_dofs);
 
    // Quadrature nodes and values of cell basis functions at these nodes
    size_t doe = 5;
    QuadratureRule quadT = generate_quadrature_rule(*cell, doe, true);
    size_t nbq = quadT.size();
    std::vector<Eigen::ArrayXd> nc_phiT_quadT = nc.nc_basis_quad(iT, quadT);
 
    // Value of source times quadrature weights at the quadrature nodes
    Eigen::ArrayXd weight_source_quad = Eigen::ArrayXd::Zero(nbq);
    for (size_t iqn = 0; iqn < nbq; iqn++){
        weight_source_quad(iqn) = quadT[iqn].w * source(quadT[iqn].x, quadT[iqn].y, cell);
    }
 
    for (size_t i=0; i < local_dofs; i++){
        b(i) = (weight_source_quad * nc_phiT_quadT[i]).sum();
    }
 
    // Boundary values, if we have a boundary cell and Neumann edges
    if (cell->is_boundary()){
    // Boundary values only on boundary Neumann edges
    for (size_t ilF = 0; ilF < cell->n_edges(); ilF++) {
      Edge* edge = cell->edge(ilF);
      if (m_BC.type(*edge)=="neu"){
                // BC on boundary faces
                if (edge->is_boundary()){
                  const auto& nTF = cell->edge_normal(ilF);
                const auto& phi_i = nc.nc_basis(iT, DimPoly<Cell>(1) + ilF);
                    QuadratureRule quadF = generate_quadrature_rule(*edge, 5);
                    std::function<double(VectorRd)> Kgrad_n = [&](VectorRd p){
              return nTF.dot(kappa(p.x(),p.y(),cell) * grad_exact_solution(p.x(),p.y(),cell)) * phi_i(p.x(),p.y());
                        };
                    for (QuadratureNode qr : quadF){
                        b(DimPoly<Cell>(1) + ilF) += qr.w * Kgrad_n(qr.vector());
                    }
                }
            }
        }
    }
    
  return b;
}
 
 
 
//*******************************************************
// Norms: L2 mass lumped, Energy 
//*******************************************************
 
double LEPNC_diffusion::EnergyNorm(const Eigen::VectorXd Xh) const {
  const Mesh* mesh = nc.get_mesh();
    double value = 0.0;
 
  for (size_t iT = 0; iT < mesh-> n_cells(); iT++) {
        Eigen::VectorXd XTF = nc.nc_restr(Xh, iT);
        value += XTF.transpose() * DiffT[iT] * XTF;
  }
 
  return sqrt(value);
}
 
 
double LEPNC_diffusion::get_assembly_time() const {
  return double(_assembly_time) * pow(10, -9);
}
 
double LEPNC_diffusion::get_solving_time() const {
  return double(_solving_time) * pow(10, -9);
}
 
double LEPNC_diffusion::get_itime(size_t idx) const {
  return double(_itime[idx]) * pow(10, -9);
}
 
double LEPNC_diffusion::get_solving_error() const {
  return _solving_error;
}
 
} // end of namespace HArDCore2D
 
#endif //_LEPNC_DIFFUSION_HPP