.

Sheet7/E14/jacob_template/geom.h

@@ -0,0 +1,712 @@
+#ifndef GEOM_FILE
+#define GEOM_FILE
+#include <array>
+#include <functional>             // function; C++11
+#include <iostream>
+#include <memory>                  // shared_ptr
+#include <string>
+#include <vector>
+
+/**
+ * Basis class for finite element meshes.
+ */
+class Mesh
+{
+public:
+    /**
+     * Constructor initializing the members with default values.
+      *
+     * @param[in] ndim  space dimensions (dimension for coordinates)
+      * @param[in] nvert_e  number of vertices per element (dimension for connectivity)
+      * @param[in] ndof_e   degrees of freedom per element (= @p nvert_e for linear elements)
+      * @param[in] nedge_e  number of edges per element (= @p nvert_e for linear elements in 2D)
+     */
+    explicit Mesh(int ndim, int nvert_e = 0, int ndof_e = 0, int nedge_e = 0);
+
+    __attribute__((noinline))
+    Mesh(Mesh const &) = default;
+
+    Mesh &operator=(Mesh const &) = delete;
+
+
+    /**
+     * Destructor.
+     *
+     * See clang warning on
+     * <a href="https://stackoverflow.com/questions/28786473/clang-no-out-of-line-virtual-method-definitions-pure-abstract-c-class/40550578">weak-vtables</a>.
+     */
+    virtual ~Mesh();
+
+    /**
+     * Reads mesh data from a binary file.
+     *
+     * File format, see ascii_write_mesh.m
+     *
+     * @param[in] fname file name
+    */
+    explicit Mesh(std::string const &fname);
+
+    /**
+     * Reads mesh data from a binary file.
+     *
+     * File format, see ascii_write_mesh.m
+     *
+     * @param[in] fname file name
+    */
+    void ReadVertexBasedMesh(std::string const &fname);
+
+    /**
+     * Number of finite elements in (sub)domain.
+     * @return number of elements.
+     */
+    int Nelems() const
+    {
+        return _nelem;
+    }
+
+    /**
+     * Global number of vertices for each finite element.
+     * @return number of vertices per element.
+     */
+    int NverticesElements() const
+    {
+        return _nvert_e;
+    }
+
+    /**
+     * Global number of degrees of freedom (dof) for each finite element.
+     * @return degrees of freedom per element.
+     */
+    int NdofsElement() const
+    {
+        return _ndof_e;
+    }
+
+    /**
+     * Number of vertices in mesh.
+     * @return number of vertices.
+     */
+    int Nnodes() const
+    {
+        return _nnode;
+    }
+
+    /**
+     * Space dimension.
+     * @return number of dimensions.
+     */
+    int Ndims() const
+    {
+        return _ndim;
+    }
+
+    /**
+     * (Re-)Allocates memory for the element connectivity and redefines the appropriate dimensions.
+      *
+     * @param[in] nelem    number of elements
+      * @param[in] nvert_e  number of vertices per element
+     */
+    void Resize_Connectivity(int nelem, int nvert_e)
+    {
+        SetNelem(nelem);               // number of elements
+        SetNverticesElement(nvert_e);  // vertices per element
+        _ia.resize(nelem * nvert_e);
+    }
+
+    /**
+     * Read connectivity information (g1,g2,g3)_i.
+     * @return connectivity vector [nelems*ndofs].
+     */
+    const std::vector<int>  &GetConnectivity() const
+    {
+        return _ia;
+    }
+
+    /**
+     * Access/Change connectivity information (g1,g2,g3)_i.
+     * @return connectivity vector [nelems*ndofs].
+     */
+    std::vector<int>  &GetConnectivity()
+    {
+        return _ia;
+    }
+
+    /**
+     * (Re-)Allocates memory for the element connectivity and redefines the appropriate dimensions.
+      *
+     * @param[in] nnodes    number of nodes
+      * @param[in] ndim      space dimension
+     */
+    void Resize_Coords(int nnodes, int ndim)
+    {
+        SetNnode(nnodes);       // number of nodes
+        SetNdim(ndim);          // space dimension
+        _xc.resize(nnodes * ndim);
+    }
+
+    /**
+     * Read coordinates of vertices (x,y)_i.
+     * @return coordinates vector [nnodes*2].
+     */
+    const std::vector<double> &GetCoords() const
+    {
+        return _xc;
+    }
+
+    /**
+     * Access/Change coordinates of vertices (x,y)_i.
+     * @return coordinates vector [nnodes*2].
+     */
+    std::vector<double> &GetCoords()
+    {
+        return _xc;
+    }
+
+    /**
+     * Calculate values in vector @p v via function @p func(x,y)
+     * @param[in] v     vector
+      * @param[in] func  function of (x,y) returning a double value.
+     */
+    void SetValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
+    void SetBoundaryValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
+    void SetDirchletValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
+
+    /**
+     * Prints the information for a finite element mesh
+     */
+    void Debug() const;
+
+    /**
+     * Prints the edge based information for a finite element mesh
+     */
+    void DebugEdgeBased() const;
+
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions
+     * @return index vector.
+     */
+    virtual std::vector<int> Index_DirichletNodes() const;
+    virtual std::vector<int> Index_BoundaryNodes() const;
+
+    /**
+     * Write vector @p v together with its mesh information to an ASCii file @p fname.
+      *
+      * The data are written in C-style.
+      *
+      * @param[in] fname  file name
+     * @param[in] v      vector
+     */
+    void Write_ascii_matlab(std::string const &fname, std::vector<double> const &v) const;
+
+    /**
+     * Exports the mesh information to ASCii files  @p basename + {_coords|_elements}.txt.
+      *
+      * The data are written in C-style.
+      *
+      * @param[in] basename  first part of file names
+     */
+    void Export_scicomp(std::string const &basename) const;
+
+    /**
+     * Visualize @p v together with its mesh information via matlab or octave.
+      *
+      * Comment/uncomment those code lines in method Mesh:Visualize (geom.cpp)
+      * that are supported on your system.
+      *
+     * @param[in] v      vector
+      *
+      * @warning matlab files ascii_read_meshvector.m  visualize_results.m
+      *          must be in the executing directory.
+     */
+    virtual void Visualize(std::vector<double> const &v) const;
+
+    /**
+     * Global number of edges.
+     * @return number of edges in mesh.
+     */
+    int Nedges() const
+    {
+        return _nedge;
+    }
+
+    /**
+     * Global number of edges for each finite element.
+     * @return number of edges per element.
+     */
+    int NedgesElements() const
+    {
+        return _nedge_e;
+    }
+
+    /**
+     * Read edge connectivity information (e1,e2,e3)_i.
+     * @return edge connectivity vector [nelems*_nedge_e].
+     */
+    const std::vector<int>  &GetEdgeConnectivity() const
+    {
+        return _ea;
+    }
+
+    /**
+     * Access/Change edge connectivity information (e1,e2,e3)_i.
+     * @return edge connectivity vector [nelems*_nedge_e].
+     */
+    std::vector<int>  &GetEdgeConnectivity()
+    {
+        return _ea;
+    }
+
+    /**
+     * Read edge information (v1,v2)_i.
+     * @return edge connectivity vector [_nedge*2].
+     */
+    const std::vector<int>  &GetEdges() const
+    {
+        return _edges;
+    }
+
+    /**
+     * Access/Change edge information (v1,v2)_i.
+     * @return edge connectivity vector [_nedge*2].
+     */
+    std::vector<int>  &GetEdges()
+    {
+        return _edges;
+    }
+
+    /**
+     * Determines all node to node connections from the vertex based mesh.
+      *
+     * @return vector[k][] containing all connections of vertex k, including to itself.
+     */
+    std::vector<std::vector<int>> Node2NodeGraph() const
+    {
+        //// Check version 2 wrt. version 1
+        //auto v1=Node2NodeGraph_1();
+        //auto v2=Node2NodeGraph_2();
+        //if ( equal(v1.cbegin(),v1.cend(),v2.begin()) )
+        //{
+        //std::cout << "\nidentical Versions\n";
+        //}
+        //else
+        //{
+        //std::cout << "\nE R R O R   in Versions\n";
+        //}
+
+        //return Node2NodeGraph_1();
+        return Node2NodeGraph_2();        // 2 times faster than version 1
+    }
+
+    /**
+     * Accesses the father-of-nodes relation.
+      *
+     * @return  vector of length 0 because no relation available.
+      *
+     */
+    virtual std::vector<int> const &GetFathersOfVertices() const
+    {
+        return _dummy;
+    }
+
+    /**
+     * Deletes all edge connectivity information (saves memory).
+     */
+    void Del_EdgeConnectivity();
+
+protected:
+    //public:
+    void SetNelem(int nelem)
+    {
+        _nelem = nelem;
+    }
+
+    void SetNverticesElement(int nvert)
+    {
+        _nvert_e = nvert;
+    }
+
+    void SetNdofsElement(int ndof)
+    {
+        _ndof_e = ndof;
+    }
+
+    void SetNnode(int nnode)
+    {
+        _nnode = nnode;
+    }
+
+    void SetNdim(int ndim)
+    {
+        _ndim = ndim;
+    }
+
+    void SetNedge(int nedge)
+    {
+        _nedge = nedge;
+    }
+
+    /**
+     * Reads vertex based mesh data from a binary file.
+     *
+     * File format, see ascii_write_mesh.m
+     *
+     * @param[in] fname file name
+    */
+    void ReadVectexBasedMesh(std::string const &fname);
+
+    /**
+     * The vertex based mesh data are used to derive the edge based data.
+     *
+     *  @warning Exactly 3 vertices, 3 edges per element are assumed (linear triangle in 2D)
+    */
+    void DeriveEdgeFromVertexBased()
+    {
+        //DeriveEdgeFromVertexBased_slow();
+        //DeriveEdgeFromVertexBased_fast();
+        DeriveEdgeFromVertexBased_fast_2();
+    }
+    void DeriveEdgeFromVertexBased_slow();
+    void DeriveEdgeFromVertexBased_fast();
+    void DeriveEdgeFromVertexBased_fast_2();
+
+
+
+    /**
+     * The edge based mesh data are used to derive the vertex based data.
+      *
+      *  @warning Exactly 3 vertices, 3 edges per element are assumed (linear triangle in 2D)
+    */
+    void DeriveVertexFromEdgeBased();
+
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions
+     * @return index vector.
+     */
+    int Nnbedges() const
+    {
+        return static_cast<int>(_bedges.size());
+    }
+
+    /**
+     * Checks whether the array dimensions fit to their appropriate size parameters
+     * @return index vector.
+     */
+    virtual bool Check_array_dimensions() const;
+
+    /**
+     * Permutes the vertex information in an edge based mesh.
+      *
+     * @param[in] old2new   new indices of original vertices.
+     */
+    void PermuteVertices_EdgeBased(std::vector<int> const &old2new);
+
+private:
+    /**
+     * Determines all node to node connections from the vertex based mesh.
+      *
+     * @return vector[k][] containing all connections of vertex k, including to itself.
+     */
+    std::vector<std::vector<int>> Node2NodeGraph_1() const;  // is correct
+
+    /**
+     * Determines all node to node connections from the vertex based mesh.
+      *
+      * Faster than @p Node2NodeGraph_1().
+      *
+     * @return vector[k][] containing all connections of vertex k, including to itself.
+     */
+    std::vector<std::vector<int>> Node2NodeGraph_2() const;  // is correct
+
+    //private:
+protected:
+    int _nelem;         //!< number elements
+    int _nvert_e;       //!< number of vertices per element
+    int _ndof_e;        //!< degrees of freedom (d.o.f.) per element
+    int _nnode;         //!< number nodes/vertices
+    int _ndim;          //!< space dimension of the problem (1, 2, or 3)
+    std::vector<int> _ia;    //!< element connectivity
+    std::vector<double> _xc; //!< coordinates
+
+protected:
+    // B.C.
+    std::vector<int> _bedges;     //!< boundary edges [nbedges][2] storing start/end vertex
+// 2020-01-08
+    std::vector<int> _sdedges;    //!< boundary edges [nbedges][2] with left/right subdomain number
+
+    //private:
+protected:
+    // edge based connectivity
+    int _nedge;              //!< number of edges in mesh
+    int _nedge_e;            //!< number of edges per element
+    std::vector<int> _edges; //!< edges of mesh (vertices ordered ascending)
+    std::vector<int> _ea;    //!< edge based element connectivity
+    // B.C.
+    std::vector<int> _ebedges; //!< boundary edges [nbedges]
+
+private:
+    const std::vector<int> _dummy; //!< empty dummy vector
+
+};
+
+
+// *********************************************************************
+
+class RefinedMesh: public Mesh
+{
+public:
+    /**
+     * Constructs a refined mesh according to the marked elements in @p ibref.
+     *
+     * If the vector @p ibref has size 0 then all elements will be refined.
+     *
+     * @param[in] cmesh  original mesh for coarsening.
+     * @param[in] ibref  vector containing True/False regarding refinement for each element
+     *
+     */
+    //explicit RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref = std::vector<bool>(0));
+    RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref);
+    //RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref);
+
+    /**
+    * Constructs a refined mesh by regulare refinement of all elements.
+    *
+    * @param[in] cmesh  original mesh for coarsening.
+    *
+    */
+    explicit RefinedMesh(Mesh const &cmesh)
+        : RefinedMesh(cmesh, std::vector<bool>(0))
+    {}
+
+
+    RefinedMesh(RefinedMesh const &) = delete;
+    //RefinedMesh(RefinedMesh const&&) = delete;
+
+    RefinedMesh &operator=(RefinedMesh const &) = delete;
+    //RefinedMesh& operator=(RefinedMesh const&&) = delete;
+
+    /**
+     * Destructor.
+     */
+    virtual ~RefinedMesh() override;
+
+    /**
+     * Refines the mesh according to the marked elements.
+     *
+     * @param[in] ibref  vector containing True/False regarding refinement for each element
+     *
+     * @return the refined mesh
+     *
+     */
+    Mesh RefineElements(std::vector<bool> const &ibref);
+
+    /**
+     * Refines all elements in the actual mesh.
+     *
+     * @param[in] nref  number of regular refinements to perform
+     *
+     */
+    void RefineAllElements(int nref = 1);
+
+    /**
+     * Accesses the father-of-nodes relation.
+      *
+     * @return  father-of-nodes relation [nnodes][2]
+      *
+     */
+    std::vector<int> const &GetFathersOfVertices() const override
+    {
+        return _vfathers;
+    }
+
+protected:
+    /**
+     * Checks whether the array dimensions fit to their appropriate size parameters
+     * @return index vector.
+     */
+    bool Check_array_dimensions() const override;
+
+    /**
+     * Permutes the vertex information in an edge based mesh.
+      *
+     * @param[in] old2new   new indices of original vertices.
+     */
+    void PermuteVertices_EdgeBased(std::vector<int> const &old2new);
+
+
+private:
+    //Mesh const              & _cmesh; //!< coarse mesh
+    std::vector<bool> const   _ibref; //!< refinement info
+    int                       _nref;  //!< number of regular refinements performed
+    std::vector<int>          _vfathers; //!< stores the 2 fathers of each vertex (equal fathers denote original coarse vertex)
+
+};
+
+// *********************************************************************
+
+class gMesh_Hierarchy
+{
+public:
+    /**
+     * Constructs mesh hierarchy of @p nlevel levels starting with coarse mesh @p cmesh.
+      * The coarse mesh @p cmesh will be @p nlevel-1 times geometrically refined.
+      *
+      * @param[in] cmesh   initial coarse mesh
+     * @param[in] nlevel  number levels in mesh hierarchy
+      *
+     */
+    gMesh_Hierarchy(Mesh const &cmesh, int nlevel);
+
+    size_t size() const
+    {
+        return _gmesh.size();
+    }
+
+    /**
+     * Access to mesh @p lev from mesh hierarchy.
+      *
+     * @return mesh @p lev
+      * @warning An out_of_range exception might be thrown.
+      *
+     */
+    Mesh const &operator[](int lev) const
+    {
+        return *_gmesh.at(lev);
+    }
+
+    /**
+     * Access to finest mesh in mesh hierarchy.
+      *
+     * @return finest mesh
+      *
+     */
+    Mesh const &finest() const
+    {
+        return *_gmesh.back();
+    }
+
+    /**
+     * Access to coarest mesh in mesh hierarchy.
+      *
+     * @return coarsest mesh
+      *
+     */
+    Mesh const &coarsest() const
+    {
+        return *_gmesh.front();
+    }
+
+private:
+    std::vector<std::shared_ptr<Mesh>> _gmesh; //!< mesh hierarchy from coarse ([0]) to fine.
+
+};
+
+
+
+// *********************************************************************
+/**
+ * 2D finite element mesh of the square consisting of linear triangular elements.
+ */
+class Mesh_2d_3_square: public Mesh
+{
+public:
+    /**
+     * Generates the f.e. mesh for the unit square.
+     *
+     * @param[in] nx    number of discretization intervals in x-direction
+     * @param[in] ny    number of discretization intervals in y-direction
+     * @param[in] myid  my MPI-rank / subdomain
+     * @param[in] procx number of ranks/subdomains in x-direction
+     * @param[in] procy number of processes in y-direction
+    */
+    Mesh_2d_3_square(int nx, int ny, int myid = 0, int procx = 1, int procy = 1);
+
+    /**
+     * Destructor
+     */
+    ~Mesh_2d_3_square() override;
+
+    /**
+     * Set solution vector based on a tensor product grid in the rectangle.
+     * @param[in] u solution vector
+     */
+    void SetU(std::vector<double> &u) const;
+
+    /**
+     * Set right hand side (rhs) vector on a tensor product grid in the rectangle.
+     * @param[in] f rhs vector
+     */
+    void SetF(std::vector<double> &f) const;
+
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions
+     * @return index vector.
+     */
+    std::vector<int> Index_DirichletNodes() const override;
+    std::vector<int> Index_BoundaryNodes() const override;
+
+    /**
+      * Stores the values of vector @p u of (sub)domain into a file @p name for further processing in gnuplot.
+      * The file stores rowise the x- and y- coordinates together with the value from  @p u .
+      * The domain [@p xl, @p xr] x [@p yb, @p yt] is discretized into @p nx x @p ny intervals.
+      *
+      * @param[in] name  basename of file name (file name will be extended by the rank number)
+      * @param[in] u     local vector
+      *
+      * @warning   Assumes tensor product grid in unit square; rowise numbered
+      *            (as generated in class constructor).
+      *            The output is provided for tensor product grid visualization
+      *            ( similar to Matlab-surf() ).
+      *
+      * @see Mesh_2d_3_square
+      */
+    void SaveVectorP(std::string const &name, std::vector<double> const &u) const;
+
+    // here will still need to implement in the class
+    //  GetBound(), AddBound()
+    //  or better a generalized way with indices and their appropriate ranks for MPI communication
+
+private:
+    /**
+      * Determines the coordinates of the discretization nodes of the domain [@p xl, @p xr] x [@p yb, @p yt]
+      * which is discretized into @p nx x @p ny intervals.
+      * @param[in] nx    number of discretization intervals in x-direction
+      * @param[in] ny    number of discretization intervals in y-direction
+      * @param[in] xl    x-coordinate of left boundary
+      * @param[in] xr    x-coordinate of right boundary
+      * @param[in] yb    y-coordinate of lower boundary
+      * @param[in] yt    y-coordinate of upper boundary
+      * @param[out] xc   coordinate vector of length 2n with x(2*k,2*k+1) as coordinates of node k
+      */
+
+    void GetCoordsInRectangle(int nx, int ny, double xl, double xr, double yb, double yt,
+                              double xc[]);
+    /**
+      * Determines the element connectivity of linear triangular elements of a FEM discretization
+      * of a rectangle using @p nx x @p ny equidistant intervals for discretization.
+      * @param[in] nx    number of discretization intervals in x-direction
+      * @param[in] ny    number of discretization intervals in y-direction
+      * @param[out] ia   element connectivity matrix with ia(3*s,3*s+1,3*s+2) as node numbers od element s
+      */
+    void GetConnectivityInRectangle(int nx, int ny, int ia[]);
+
+private:
+    int _myid;          //!< my MPI rank
+    int _procx;         //!< number of MPI ranks in x-direction
+    int _procy;         //!< number of MPI ranks in y-direction
+    std::array<int, 4> _neigh; //!< MPI ranks of neighbors (negative: no neighbor but b.c.)
+    int _color;         //!< red/black coloring (checker board) of subdomains
+
+    double _xl;         //!< x coordinate of lower left  corner of square
+    double _xr;         //!< x coordinate of lower right corner of square
+    double _yb;         //!< y coordinate or lower left  corner of square
+    double _yt;         //!< y coordinate of upper right corner of square
+    int    _nx;         //!< number of intervals in x-direction
+    int    _ny;         //!< number of intervals in y-direction
+};
+
+// *********************************************************************
+
+
+
+
+#endif

Sheet7/E14/jacob_template/getmatrix.cpp

@@ -0,0 +1,720 @@
+#include "getmatrix.h"
+#include "userset.h"
+
+#include <omp.h>
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <ctime>                  // contains clock()
+#include <iomanip>
+#include <iostream>
+#include <list>
+#include <vector>
+using namespace std;
+
+// ####################################################################
+
+Matrix::Matrix(int const nrows, int const ncols)
+ : _nrows(nrows), _ncols(ncols), _dd(0)
+{}
+
+//Matrix::Matrix()
+ //: Matrix(0,0)
+//{}
+
+
+Matrix::~Matrix()
+{}
+
+//vector<double> const & Matrix::GetDiag() const
+//{
+    //if ( 0==_dd.size() )
+    //{
+        //_dd.resize(Nrows());
+        //this->GetDiag(_dd);
+    //}
+    //assert( Nrows()==static_cast<int>(_dd.size()) );
+    //return _dd;
+//}
+
+// ####################################################################
+
+CRS_Matrix::CRS_Matrix()
+ : Matrix(0,0), _nnz(0), _id(0), _ik(0), _sk(0)
+{}
+
+CRS_Matrix::~CRS_Matrix()
+{}
+
+void CRS_Matrix::Mult(vector<double> &w, vector<double> const &u) const
+{
+    assert( _ncols==static_cast<int>(u.size()) );   // compatibility of inner dimensions
+    assert( _nrows==static_cast<int>(w.size()) );   // compatibility of outer dimensions
+
+    for (int row = 0; row < _nrows; ++row)
+    {
+        double wi = 0.0;
+        for (int ij = _id[row]; ij < _id[row + 1]; ++ij)
+        {
+            wi += _sk[ij] * u[ _ik[ij] ];
+        }
+        w[row] = wi;
+    }
+    return;
+}
+
+void CRS_Matrix::Defect(vector<double> &w,
+                        vector<double> const &f, vector<double> const &u) const
+{
+    assert( _ncols==static_cast<int>(u.size()) );   // compatibility of inner dimensions
+    assert( _nrows==static_cast<int>(w.size()) );   // compatibility of outer dimensions
+    assert( w.size()==f.size() );
+
+#pragma omp parallel for
+    for (int row = 0; row < _nrows; ++row)
+    {
+        double wi = f[row];
+        for (int ij = _id[row]; ij < _id[row + 1]; ++ij)
+        {
+            wi -= _sk[ij] * u[ _ik[ij] ];
+        }
+        w[row] = wi;
+    }
+    return;
+}
+
+void CRS_Matrix::GetDiag(vector<double> &d) const
+{
+    // be carefull when using a rectangular matrix
+    int const nm = min(_nrows, _ncols);
+
+    assert( nm==static_cast<int>(d.size()) );  // instead of stopping we could resize d and warn the user
+
+    for (int row = 0; row < nm; ++row)
+    {
+        const int ia = fetch(row, row); // Find diagonal entry of row
+        assert(ia >= 0);
+        d[row] = _sk[ia];
+    }
+    return;
+}
+
+inline
+int CRS_Matrix::fetch(int const row, int const col) const
+{
+    int const id2 = _id[row + 1];    // end   and
+    int       ip  = _id[row];        // start of recent row (global index)
+
+    while (ip < id2 && _ik[ip] != col)  // find index col (global index)
+    {
+        ++ip;
+    }
+    if (ip >= id2)
+    {
+        ip = -1;
+#ifndef NDEBUG                 // compiler option -DNDEBUG switches off the check
+        cout << "No column  " << col << "  in row  " << row << endl;
+        assert(ip >= id2);
+#endif
+    }
+    return ip;
+}
+
+void CRS_Matrix::Debug() const
+{
+//  ID points to first entry of row
+//  no symmetry assumed
+    cout << "\nMatrix  (" << _nrows << " x " << _ncols << "  with  nnz = " << _id[_nrows] << ")\n";
+
+    for (int row = 0; row < _nrows; ++row)
+    {
+        cout << "Row " << row << " : ";
+        int const id1 = _id[row];
+        int const id2 = _id[row + 1];
+        for (int j = id1; j < id2; ++j)
+        {
+            cout.setf(ios::right, ios::adjustfield);
+            cout << "[" << setw(2) << _ik[j] << "]  " << setw(4) << _sk[j] << "  ";
+        }
+        cout << endl;
+    }
+    return;
+}
+
+
+
+bool CRS_Matrix::Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const
+{
+    bool bn = (nnode==_nrows);       // number of rows
+    if (!bn)
+    {
+        cout << "#########   Error: " << "number of rows" << endl;
+    }
+
+    bool bz = (id[nnode]==_nnz);     // number of non zero elements
+    if (!bz)
+    {
+        cout << "#########   Error: " << "number of non zero elements" << endl;
+    }
+
+    bool bd = equal(id,id+nnode+1,_id.cbegin());  // row starts
+    if (!bd)
+    {
+        cout << "#########   Error: " << "row starts" << endl;
+    }
+
+    bool bk = equal(ik,ik+id[nnode],_ik.cbegin());  // column indices
+    if (!bk)
+    {
+        cout << "#########   Error: " << "column indices" << endl;
+    }
+
+    bool bv = equal(sk,sk+id[nnode],_sk.cbegin());  // values
+    if (!bv)
+    {
+        cout << "#########   Error: " << "values" << endl;
+    }
+
+    return bn && bz && bd && bk && bv;
+}
+
+
+// ####################################################################
+
+FEM_Matrix::FEM_Matrix(Mesh const & mesh)
+ : CRS_Matrix(), _mesh(mesh)
+{
+    Derive_Matrix_Pattern();
+    return;
+}
+
+FEM_Matrix::~FEM_Matrix()
+{}
+
+void FEM_Matrix::Derive_Matrix_Pattern_fast()
+{
+    cout << "\n############   FEM_Matrix::Derive_Matrix_Pattern ";
+    double tstart = clock();
+    int const nelem(_mesh.Nelems());
+    int const ndof_e(_mesh.NdofsElement());
+    auto const &ia(_mesh.GetConnectivity());
+//  Determine the number of matrix rows
+    _nrows = *max_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem);
+    ++_nrows;                                 // node numberng: 0 ... nnode-1
+    assert(*min_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem) == 0); // numbering starts with 0 ?
+
+// CSR data allocation
+    _id.resize(_nrows + 1);                  // Allocate memory for CSR row pointer
+//##########################################################################
+    auto const v2v=_mesh.Node2NodeGraph();
+    _nnz=0;                           // number of connections
+    _id[0] = 0;                       // start of matrix row zero
+    for (size_t v = 0; v<v2v.size(); ++v )
+    {
+        _id[v+1] = _id[v] + v2v[v].size();
+        _nnz += v2v[v].size();
+    }
+    assert(_nnz == _id[_nrows]);
+    _sk.resize(_nnz);                        // Allocate memory for CSR column index vector
+
+// CSR data allocation
+    _ik.resize(_nnz);                        // Allocate memory for CSR column index vector
+// Copy column indices
+    int kk = 0;
+    for (size_t v = 0; v<v2v.size(); ++v )
+    {
+        for (size_t vi=0; vi<v2v[v].size(); ++vi)
+        {
+             _ik[kk] = v2v[v][vi];
+             ++kk;
+        }
+    }
+    _ncols = *max_element(_ik.cbegin(), _ik.cend());  // maximal column number
+    ++_ncols;                                         // node numbering: 0 ... nnode-1
+    //cout << _nrows << "  " << _ncols << endl;
+    assert(_ncols==_nrows);
+
+    double duration = (clock() - tstart) / CLOCKS_PER_SEC;
+    cout << "finished in  " <<  duration  << " sec.    ########\n";
+
+    return;
+}
+
+
+void FEM_Matrix::Derive_Matrix_Pattern_slow()
+{
+    cout << "\n############   FEM_Matrix::Derive_Matrix_Pattern slow ";
+    double tstart = clock();
+    int const nelem(_mesh.Nelems());
+    int const ndof_e(_mesh.NdofsElement());
+    auto const &ia(_mesh.GetConnectivity());
+//  Determine the number of matrix rows
+    _nrows = *max_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem);
+    ++_nrows;                                 // node numberng: 0 ... nnode-1
+    assert(*min_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem) == 0); // numbering starts with 0 ?
+
+//  Collect for each node those nodes it is connected to (multiple entries)
+//  Detect the neighboring nodes
+    vector< list<int> > cc(_nrows);             //  cc[i] is the  list of nodes a node i is connected to
+    for (int i = 0; i < nelem; ++i)
+    {
+        int const idx = ndof_e * i;
+        for (int k = 0; k < ndof_e; ++k)
+        {
+            list<int> &cck = cc[ia[idx + k]];
+            cck.insert( cck.end(), ia.cbegin() + idx, ia.cbegin() + idx + ndof_e );
+        }
+    }
+//  Delete the multiple entries
+    _nnz = 0;
+    for (auto &it : cc)
+    {
+        it.sort();
+        it.unique();
+        _nnz += it.size();
+        // cout << it.size() << " :: "; copy(it->begin(),it->end(), ostream_iterator<int,char>(cout,"  ")); cout << endl;
+    }
+
+// CSR data allocation
+    _id.resize(_nrows + 1);                  // Allocate memory for CSR row pointer
+    _ik.resize(_nnz);                        // Allocate memory for CSR column index vector
+
+//  copy CSR data
+    _id[0] = 0;                              // begin of first row
+    for (size_t i = 0; i < cc.size(); ++i)
+    {
+        //cout << i << "   " << nid.at(i) << endl;;
+        const list<int> &ci = cc[i];
+        const auto nci = static_cast<int>(ci.size());
+        _id[i + 1] = _id[i] + nci; // begin of next line
+        copy(ci.begin(), ci.end(), _ik.begin() + _id[i] );
+    }
+
+    assert(_nnz == _id[_nrows]);
+    _sk.resize(_nnz);                        // Allocate memory for CSR column index vector
+
+    _ncols = *max_element(_ik.cbegin(), _ik.cend());  // maximal column number
+    ++_ncols;                                 // node numbering: 0 ... nnode-1
+    //cout << _nrows << "  " << _ncols << endl;
+    assert(_ncols==_nrows);
+
+    double duration = (clock() - tstart) / CLOCKS_PER_SEC;
+    cout << "finished in  " <<  duration  << " sec.    ########\n";
+
+    return;
+}
+
+
+void FEM_Matrix::CalculateLaplace(vector<double> &f)
+{
+    cout << "\n############   FEM_Matrix::CalculateLaplace ";
+    //double tstart = clock();
+    double tstart = omp_get_wtime();                  // OpenMP    
+    assert(_mesh.NdofsElement() == 3);               // only for triangular, linear elements
+    //cout << _nnz << " vs. " << _id[_nrows] << "  " << _nrows<< endl;
+    assert(_nnz == _id[_nrows]);
+
+    for (int k = 0; k < _nrows; ++k)
+    {
+        _sk[k] = 0.0;
+    }
+    for (int k = 0; k < _nrows; ++k)
+    {
+        f[k] = 0.0;
+    }
+
+    double ske[3][3], fe[3];
+    //  Loop over all elements
+    auto const nelem = _mesh.Nelems();
+    auto const &ia   = _mesh.GetConnectivity();
+    auto const &xc   = _mesh.GetCoords();
+
+#pragma omp parallel for private(ske,fe)
+    for (int i = 0; i < nelem; ++i)
+    {
+        CalcElem(ia.data()+3 * i, xc.data(), ske, fe);
+        //AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
+        AddElem_3(ia.data()+3 * i, ske, fe, f);
+    }
+
+    //double duration = (clock() - tstart) / CLOCKS_PER_SEC;
+    double duration = omp_get_wtime() - tstart;             // OpenMP
+    cout << "finished in  " <<  duration  << " sec.    ########\n";
+    //Debug();
+
+    return;
+}
+
+//void FEM_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f)
+//{
+    //double const PENALTY = 1e6;
+    //auto const idx = _mesh.Index_DirichletNodes();
+    //int const nidx = idx.size();
+
+    //for (int i=0; i<nidx; ++i)
+    //{
+        //int const k = idx[i];
+		//int const id1 = fetch(k, k); // Find diagonal entry of k
+		//assert(id1 >= 0);
+		//_sk[id1] += PENALTY;		// matrix weighted scaling feasible
+        //f[k] += PENALTY * u[k];
+    //}
+
+    //return;
+//}
+
+void FEM_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f)
+{
+    auto const idx = _mesh.Index_DirichletNodes();
+    int const nidx = idx.size();
+
+    for (int i=0; i<nidx; ++i)
+    {
+        int const row = idx[i];
+        for (int ij=_id[row]; ij<_id[row+1]; ++ij)
+        {
+            int const col=_ik[ij];
+            if (col==row)
+            {
+                _sk[ij] = 1.0;
+                f[row]  = u[row];
+            }
+            else
+            {
+                int const id1 = fetch(col, row); // Find entry (col,row)
+		        assert(id1 >= 0);
+                f[col] -= _sk[id1]*u[row];
+                _sk[id1] = 0.0;
+                _sk[ij]  = 0.0;
+            }
+        }
+    }
+
+    return;
+}
+
+
+
+void FEM_Matrix::AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], vector<double> & f)
+{
+    for (int i = 0; i < 3; ++i)
+    {
+        const int ii  = ial[i];           // row ii (global index)
+        for (int j = 0; j < 3; ++j)       // no symmetry assumed
+        {
+            const int jj = ial[j];        // column jj (global index)
+            const int ip = fetch(ii,jj);        // find column entry jj in row ii
+#ifndef NDEBUG                 // compiler option -DNDEBUG switches off the check
+            if (ip<0)          // no entry found !!
+            {
+                cout << "Error in AddElem: (" << ii << "," << jj << ") ["
+                     << ial[0] << "," << ial[1] << "," << ial[2] << "]\n";
+                assert(ip>=0);
+            }
+#endif
+#pragma omp atomic
+            _sk[ip] += ske[i][j];
+        }
+#pragma omp atomic
+        f[ii] += fe[i];
+    }
+}
+
+// ####################################################################
+
+
+//Prolongation::Prolongation(Mesh const & cmesh, Mesh const & fmesh)
+ //: CRS_Matrix(), _cmesh(cmesh), _fmesh(fmesh)
+//{
+    //Derive_Matrix_Pattern();
+    //return;
+//}
+
+
+//void Prolongation::Derive_Matrix_Pattern()
+//{
+    //cout << "\n *****  Subject to ongoing imlementation.  *****\n";
+//}
+
+// ####################################################################
+// *********************************************************************
+
+
+//  general routine for lin. triangular elements
+
+void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3])
+//void CalcElem(const int* __restrict__ ial, const double* __restrict__ xc, double* __restrict__ ske[3], double* __restrict__ fe)
+{
+    const int  i1  = 2 * ial[0],   i2 = 2 * ial[1],   i3 = 2 * ial[2];
+    const double x13 = xc[i3 + 0] - xc[i1 + 0],  y13 = xc[i3 + 1] - xc[i1 + 1],
+                 x21 = xc[i1 + 0] - xc[i2 + 0],  y21 = xc[i1 + 1] - xc[i2 + 1],
+                 x32 = xc[i2 + 0] - xc[i3 + 0],  y32 = xc[i2 + 1] - xc[i3 + 1];
+    const double jac = fabs(x21 * y13 - x13 * y21);
+
+    ske[0][0] = 0.5 / jac * (y32 * y32 + x32 * x32);
+    ske[0][1] = 0.5 / jac * (y13 * y32 + x13 * x32);
+    ske[0][2] = 0.5 / jac * (y21 * y32 + x21 * x32);
+    ske[1][0] = ske[0][1];
+    ske[1][1] = 0.5 / jac * (y13 * y13 + x13 * x13);
+    ske[1][2] = 0.5 / jac * (y21 * y13 + x21 * x13);
+    ske[2][0] = ske[0][2];
+    ske[2][1] = ske[1][2];
+    ske[2][2] = 0.5 / jac * (y21 * y21 + x21 * x21);
+
+    const double xm    = (xc[i1 + 0] + xc[i2 + 0] + xc[i3 + 0]) / 3.0,
+                 ym    = (xc[i1 + 1] + xc[i2 + 1] + xc[i3 + 1]) / 3.0;
+    //fe[0] = fe[1] = fe[2] = 0.5 * jac * FunctF(xm, ym) / 3.0;
+    fe[0] = fe[1] = fe[2] = 0.5 * jac * fNice(xm, ym) / 3.0;
+}
+
+void CalcElem_Masse(int const ial[3], double const xc[], double const cm, double ske[3][3])
+{
+    const int  i1  = 2 * ial[0],   i2 = 2 * ial[1],   i3 = 2 * ial[2];
+    const double x13 = xc[i3 + 0] - xc[i1 + 0],  y13 = xc[i3 + 1] - xc[i1 + 1],
+                 x21 = xc[i1 + 0] - xc[i2 + 0],  y21 = xc[i1 + 1] - xc[i2 + 1];
+                 //x32 = xc[i2 + 0] - xc[i3 + 0],  y32 = xc[i2 + 1] - xc[i3 + 1];
+    const double jac = fabs(x21 * y13 - x13 * y21);
+
+    ske[0][0] += jac/12.0;
+    ske[0][1] += jac/24.0;
+    ske[0][2] += jac/24.0;
+    ske[1][0] += jac/24.0;
+    ske[1][1] += jac/12.0;
+    ske[1][2] += jac/24.0;
+    ske[2][0] += jac/24.0;
+    ske[2][1] += jac/24.0;
+    ske[2][2] += jac/12.0;
+
+    return;
+}
+
+
+// general routine for lin. triangular elements,
+// non-symm. matrix
+// node numbering in element:  a s c e n d i n g   indices !!
+// GH: deprecated
+void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
+             int const id[], int const ik[], double sk[], double f[])
+{
+    for (int i = 0; i < 3; ++i)
+    {
+        const int ii  = ial[i],           // row ii (global index)
+                  id1 = id[ii],           // start and
+                  id2 = id[ii + 1];       // end of row ii in matrix
+        int ip  = id1;
+        for (int j = 0; j < 3; ++j)       // no symmetry assumed
+        {
+            const int jj = ial[j];
+            bool not_found = true;
+            do       // find entry jj (global index) in row ii
+            {
+                not_found = (ik[ip] != jj);
+                ++ip;
+            }
+            while (not_found && ip < id2);
+
+#ifndef NDEBUG                 // compiler option -DNDEBUG switches off the check
+            if (not_found)     // no entry found !!
+            {
+                cout << "Error in AddElem: (" << ii << "," << jj << ") ["
+                     << ial[0] << "," << ial[1] << "," << ial[2] << "]\n";
+                assert(!not_found);
+            }
+#endif
+            sk[ip - 1] += ske[i][j];
+        }
+        f[ii] += fe[i];
+    }
+}
+
+
+// ----------------------------------------------------------------------------
+
+// #####################################################################
+
+BisectInterpolation::BisectInterpolation()
+: Matrix( 0, 0 ), _iv(), _vv()
+{
+}
+
+BisectInterpolation::BisectInterpolation(std::vector<int> const & fathers)
+: Matrix( static_cast<int>(fathers.size())/2, 1+*max_element(fathers.cbegin(),fathers.cend()) ),
+  _iv(fathers), _vv(fathers.size(),0.5)
+{
+}
+
+BisectInterpolation::~BisectInterpolation()
+{}
+
+void BisectInterpolation::GetDiag(vector<double> &d) const
+{
+    assert( Nrows()==static_cast<int>(d.size()) );
+
+    for (int k=0; k<Nrows(); ++k)
+    {
+        if ( _iv[2*k]==_iv[2*k+1] )
+        {
+            d[k] = 1.0;
+        }
+        else
+        {
+            d[k] = 0.0;
+        }
+    }
+    return;
+}
+
+void BisectInterpolation::Mult(vector<double> &wf, vector<double> const &uc) const
+{
+    assert( Nrows()==static_cast<int>(wf.size()) );
+    assert( Ncols()==static_cast<int>(uc.size()) );
+
+#pragma omp parallel for
+    for (int k=0; k<Nrows(); ++k)
+    {
+        wf[k] = _vv[2*k]*uc[_iv[2*k]] + _vv[2*k+1]*uc[_iv[2*k+1]];
+    }
+    return;
+}
+
+//void BisectInterpolation::MultT(vector<double> const &wf, vector<double> &uc) const
+//{
+    //assert( Nrows()==static_cast<int>(wf.size()) );
+    //assert( Ncols()==static_cast<int>(uc.size()) );
+//// GH: atomic slows down the code ==> use different storage for MultT operation (CRS-matrix?)
+////#pragma omp parallel for
+    //for (int k=0; k<Ncols(); ++k)  uc[k] = 0.0;
+//#pragma omp parallel for
+    //for (int k=0; k<Nrows(); ++k)
+    //{
+//#pragma omp atomic
+        //uc[_iv[2*k]  ] += _vv[2*k  ]*wf[k];
+//#pragma omp atomic
+        //uc[_iv[2*k+1]] += _vv[2*k+1]*wf[k];
+    //}
+    //return;
+//}
+
+void BisectInterpolation::MultT(vector<double> const &wf, vector<double> &uc) const
+{
+    assert( Nrows()==static_cast<int>(wf.size()) );
+    assert( Ncols()==static_cast<int>(uc.size()) );
+// GH: atomic slows down the code ==> use different storage for MultT operation (CRS-matrix?)
+//#pragma omp parallel for
+    for (int k=0; k<Ncols(); ++k)  uc[k] = 0.0;
+#pragma omp parallel for
+    for (int k=0; k<Nrows(); ++k)
+    {
+        if (_iv[2*k]!=_iv[2*k+1])
+        {
+#pragma omp atomic
+            uc[_iv[2*k]  ] += _vv[2*k  ]*wf[k];
+#pragma omp atomic
+            uc[_iv[2*k+1]] += _vv[2*k+1]*wf[k];
+        }
+        else
+        {
+#pragma omp atomic
+            uc[_iv[2*k]  ] += 2.0*_vv[2*k  ]*wf[k]; // uses a property of class BisectInterpolation
+        }
+    }
+    return;
+}
+
+void BisectInterpolation::Defect(vector<double> &w,
+                   vector<double> const &f, vector<double> const &u) const
+{
+    assert( Nrows()==static_cast<int>(w.size()) );
+    assert( Ncols()==static_cast<int>(u.size()) );
+    assert( w.size()==f.size() );
+
+    for (int k=0; k<Nrows(); ++k)
+    {
+        w[k] = f[k] - _vv[2*k]*u[_iv[2*k]] + _vv[2*k+1]*u[_iv[2*k+1]];
+    }
+    return;
+}
+
+void BisectInterpolation::Debug() const
+{
+    for (int k=0; k<Nrows(); ++k)
+    {
+        cout << k << " : fathers(" << _iv[2*k] << "," << _iv[2*k+1] << ")    ";
+        cout << "weights(" << _vv[2*k] << "," << _vv[2*k+1] << endl;
+    }
+    cout << endl;
+    return;
+}
+
+int BisectInterpolation::fetch(int row, int col) const
+{
+    int idx(-1);
+    if (_iv[2*row  ] == col) idx = 2*row;
+    if (_iv[2*row+1] == col) idx = 2*row+1;
+    assert(idx>=0);
+    return idx;
+}
+
+// #####################################################################
+
+BisectIntDirichlet::BisectIntDirichlet(std::vector<int> const & fathers, std::vector<int> const & idxc_dir)
+ : BisectInterpolation(fathers)
+{
+    vector<bool> bdir(Ncols(), false);               // Indicator for Dirichlet coarse nodes
+    for (size_t kc=0; kc<idxc_dir.size(); ++kc)
+    {
+        bdir.at(idxc_dir[kc]) = true;                // Mark Dirichlet node from coarse mesh
+    }
+
+    for (size_t j=0; j<_iv.size(); ++j)
+    {
+        if ( bdir.at(_iv[j]) )   _vv[j] = 0.0;       // set weight to zero iff (at least) one father is Dirichlet node
+    }
+    return;
+}
+
+
+BisectIntDirichlet::~BisectIntDirichlet()
+{}
+
+
+// #####################################################################
+
+void DefectRestrict(CRS_Matrix const & SK, BisectInterpolation const& P, 
+       vector<double> &fc, vector<double> &ff, vector<double> &uf)
+{
+    assert( P.Nrows()==static_cast<int>(ff.size()) );
+    assert( P.Ncols()==static_cast<int>(fc.size()) );
+    assert( ff.size()==uf.size() );
+    assert( P.Nrows()==SK.Nrows() );
+
+//#pragma omp parallel for
+    for (int k=0; k<P.Ncols(); ++k)  fc[k] = 0.0;
+   
+// GH: atomic slows down the code ==> use different storage for MultT operation (CRS-matrix?)
+#pragma omp parallel for
+    for (int row = 0; row < SK._nrows; ++row)
+    {        
+        double wi = ff[row];
+        for (int ij = SK._id[row]; ij < SK._id[row + 1]; ++ij)
+        {
+            wi -= SK._sk[ij] * uf[ SK._ik[ij] ];
+        }
+        
+        const int i1=P._iv[2*row];
+        const int i2=P._iv[2*row+1];
+        if (i1!=i2)
+        {
+#pragma omp atomic
+            fc[i1] += P._vv[2*row  ]*wi;
+#pragma omp atomic
+            fc[i2] += P._vv[2*row +1]*wi; 
+        }
+        else
+        {
+#pragma omp atomic
+            fc[i1] += 2.0*P._vv[2*row  ]*wi; // uses a property of class BisectInterpolation
+        }
+    }
+    return;
+}

Sheet7/E14/jacob_template/getmatrix.h

@@ -0,0 +1,545 @@
+#ifndef GETMATRIX_FILE
+#define GETMATRIX_FILE
+
+#include "geom.h"
+#include <cassert>
+#include <vector>
+// #####################################################################
+/**
+ * Abstract matrix class.
+ */
+class Matrix
+{
+public:
+    /**
+     * Constructor for abstract matrix class.
+      *
+      * No memory is allocated.
+     *
+     * @param[in] nrows   number of matrix rows.
+      * @param[in] ncols   number of matrix columns.
+    */
+    Matrix(int nrows, int ncols);
+    //Matrix();
+
+    Matrix(Matrix const &) = default;
+    /**
+     * Destructor.
+      *
+      * No memory is allocated.
+    */
+    virtual ~Matrix();
+
+    /**
+     * Checks whether the matrix is a square matrix.
+     *
+     * @return True iff square matrix.
+    */
+    bool isSquare() const
+    {
+        return _nrows == _ncols;
+    }
+
+    /**
+     * Number of rows in matrix.
+     * @return number of rows.
+     */
+    int Nrows() const
+    {
+        return _nrows;
+    }
+
+    /**
+     * Number of columns in matrix.
+     * @return number of columns.
+     */
+    int Ncols() const
+    {
+        return _ncols;
+    }
+
+    /**
+     * Show the matrix entries.
+     */
+    virtual void Debug() const = 0;
+
+    /**
+     * Extracts the diagonal elements of an inherited matrix.
+     *
+     * @param[in,out]  d  (prellocated) vector of diagonal elements
+     */
+    virtual void GetDiag(std::vector<double> &d) const = 0;
+
+    /**
+     * Extracts the diagonal elements of the matrix.
+     *
+     * @return  d  vector of diagonal elements
+     */
+    std::vector<double> const &GetDiag() const
+    {
+        if ( 0 == _dd.size() )        // GH: better?   Nrows()>static_cast<int>(_dd.size())
+        {
+            _dd.resize(Nrows());
+            this->GetDiag(_dd);
+        }
+        assert( Nrows() == static_cast<int>(_dd.size()) );
+        return _dd;
+    }
+
+    /**
+     * Performs the matrix-vector product  w := K*u.
+     *
+     * @param[in,out] w resulting vector (preallocated)
+     * @param[in]     u vector
+     */
+    virtual void Mult(std::vector<double> &w, std::vector<double> const &u) const = 0;
+
+    /**
+    * Calculates the defect/residuum w := f - K*u.
+    *
+    * @param[in,out] w resulting vector (preallocated)
+    * @param[in]     f load vector
+    * @param[in]     u vector
+    */
+    virtual void Defect(
+        std::vector<double> &w,
+        std::vector<double> const &f, std::vector<double> const &u) const = 0;
+
+    /**
+     * Finds in a CRS matrix the access index for an entry at row @p row and column @p col.
+     *
+     * @param[in] row    row index
+     * @param[in] col    column index
+     * @return index for element (@p row, @p col). If no appropriate entry exists then -1 will be returned.
+     *
+     * @warning assert() stops the function in case that matrix element (@p row, @p col) doesn't exist.
+    */
+    virtual int fetch(int row, int col) const = 0;
+
+protected:
+    int _nrows;              //!< number of rows in matrix
+    int _ncols;              //!< number of columns in matrix
+    mutable std::vector<double> _dd; //!< diagonal matrix elements
+};
+
+// #####################################################################
+class BisectInterpolation;  // class forward declaration
+/**
+ * Matrix in CRS format (compressed row storage; also named CSR),
+ * see an <a href="https://en.wikipedia.org/wiki/Sparse_matrix">introduction</a>.
+ */
+class CRS_Matrix: public Matrix
+{
+public:
+    /**
+     * Constructor
+     *
+     */
+    CRS_Matrix();
+
+    CRS_Matrix(CRS_Matrix const &) = default;
+
+
+    /**
+      * Destructor.
+      */
+    virtual ~CRS_Matrix() override;
+    /**
+     * Extracts the diagonal elements of the sparse matrix.
+     *
+     * @param[in,out]  d  (prellocated) vector of diagonal elements
+     */
+    void GetDiag(std::vector<double> &d) const override;
+    ///**
+    //* Extracts the diagonal elements of the sparse matrix.
+    //*
+    //* @return  d  vector of diagonal elements
+    //*/
+    //std::vector<double> const & GetDiag() const override;
+
+    /**
+     * Performs the matrix-vector product  w := K*u.
+     *
+     * @param[in,out] w resulting vector (preallocated)
+     * @param[in]     u vector
+     */
+    void Mult(std::vector<double> &w, std::vector<double> const &u) const override;
+
+    /**
+    * Calculates the defect/residuum w := f - K*u.
+    *
+    * @param[in,out] w resulting vector (preallocated)
+    * @param[in]     f load vector
+    * @param[in]     u vector
+    */
+    void Defect(std::vector<double> &w,
+                std::vector<double> const &f, std::vector<double> const &u) const override;
+
+    /**
+     * Show the matrix entries.
+     */
+    void Debug() const override;
+
+    /**
+     * Finds in a CRS matrix the access index for an entry at row @p row and column @p col.
+     *
+     * @param[in] row    row index
+     * @param[in] col    column index
+     * @return index for element (@p row, @p col). If no appropriate entry exists then -1 will be returned.
+     *
+     * @warning assert() stops the function in case that matrix element (@p row, @p col) doesn't exist.
+    */
+    int fetch(int row, int col) const override;
+
+    /**
+    * Compare @p this CRS matrix with an external CRS matrix stored in C-Style.
+    *
+    * The method prints statements on differences found.
+    *
+    * @param[in]     nnode  row number of external matrix
+    * @param[in]     id     start indices of matrix rows of external matrix
+    * @param[in]     ik     column indices of external matrix
+    * @param[in]     sk     non-zero values of external matrix
+    *
+    * @return true iff all data are identical.
+    */
+    bool Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const;
+
+    /**
+     * Calculates the defect and projects it to the next coarser level @f$ f_C := P^T \cdot (f_F - SK\cdot u_F) @f$.
+     *
+     * @param[in] SK matrix on fine mesh
+     * @param[in] P      prolongation operator
+     * @param[in,out] fc  resulting coarse mesh vector (preallocated)
+     * @param[in] ff r.h.s. on fine mesh
+     * @param[in] uf status vector on fine mesh
+     *
+    */
+    friend void DefectRestrict(CRS_Matrix const &SK, BisectInterpolation const &P,
+                               std::vector<double> &fc, std::vector<double> &ff, std::vector<double> &uf);
+
+protected:
+    //int _nrows;              //!< number of rows in matrix
+    //int _ncols;              //!< number of columns in matrix
+    int _nnz;                //!< number of non-zero entries
+    std::vector<int> _id;    //!< start indices of matrix rows
+    std::vector<int> _ik;    //!< column indices
+    std::vector<double> _sk; //!< non-zero values
+};
+
+
+/**
+ * FEM Matrix in CRS format (compressed row storage; also named CSR),
+ * see an <a href="https://en.wikipedia.org/wiki/Sparse_matrix">introduction</a>.
+ */
+class FEM_Matrix: public CRS_Matrix
+{
+public:
+    /**
+     * Initializes the CRS matrix structure from the given discretization in @p mesh.
+     *
+     * The sparse matrix pattern is generated but the values are 0.
+     *
+     * @param[in] mesh given discretization
+     *
+     * @warning A reference to the discretization @p mesh is stored inside this class.
+     *          Therefore, changing @p mesh outside requires also
+     *          to call method @p Derive_Matrix_Pattern explicitly.
+     *
+     * @see Derive_Matrix_Pattern
+     */
+    explicit FEM_Matrix(Mesh const &mesh);
+
+    FEM_Matrix(FEM_Matrix const &) = default;
+
+    /**
+      * Destructor.
+      */
+    ~FEM_Matrix() override;
+
+    /**
+     * Generates the sparse matrix pattern and overwrites the existing pattern.
+     *
+     * The sparse matrix pattern is generated but the values are 0.
+    */
+    void Derive_Matrix_Pattern()
+    {
+        //Derive_Matrix_Pattern_slow();
+        Derive_Matrix_Pattern_fast();
+    }
+    void Derive_Matrix_Pattern_fast();
+    void Derive_Matrix_Pattern_slow();
+
+
+    /**
+    * Calculates the entries of f.e. stiffness matrix and load/rhs vector @p f for the Laplace operator in 2D.
+    * No memory is allocated.
+    *
+    * @param[in,out] f (preallocated) rhs/load vector
+    */
+    void CalculateLaplace(std::vector<double> &f);
+
+    /**
+     * Applies Dirichlet boundary conditions to stiffness matrix and to load vector @p f.
+     * The <a href="https://www.jstor.org/stable/2005611?seq=1#metadata_info_tab_contents">penalty method</a>
+     * is used for incorporating the given values @p u.
+     *
+     * @param[in]     u (global) vector with Dirichlet data
+     * @param[in,out] f load vector
+     */
+    void ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f);
+
+    ///**
+    //* Extracts the diagonal elements of the sparse matrix.
+    //*
+    //* @param[in,out]  d  (prellocated) vector of diagonal elements
+    //*/
+    //void GetDiag(std::vector<double> &d) const;   // override in MPI parallel
+
+
+    /**
+      * Adds the element stiffness matrix @p ske and the element load vector @p fe
+      * of one triangular element with linear shape functions to the appropriate positions in
+      * the stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik).
+      *
+      * @param[in]     ial   node indices of the three element vertices
+      * @param[in]     ske   element stiffness matrix
+      * @param[in]     fe    element load vector
+      * @param[in,out] f     distributed local vector storing the right hand side
+      *
+      * @warning Algorithm assumes  linear triangular elements (ndof_e==3).
+     */
+    void AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], std::vector<double> &f);
+
+
+private:
+    Mesh const &_mesh;       //!< reference to discretization
+
+};
+
+
+///**
+//* Prolongation matrix in CRS format (compressed row storage; also named CSR),
+//* see an <a href="https://en.wikipedia.org/wiki/Sparse_matrix">introduction</a>.
+//*
+//* The prolongation is applied for each node from the coarse mesh to the fine mesh and
+//* is derived only geometrically (no operator weighted prolongation).
+//*/
+//class Prolongation: public CRS_Matrix
+//{
+//public:
+///**
+//* Intializes the CRS matrix structure from the given discetization in @p mesh.
+//*
+//* The sparse matrix pattern is generated but the values are 0.
+//*
+//* @param[in] cmesh coarse mesh
+//* @param[in] fmesh fine mesh
+//*
+//* @warning A reference to the discretizations @p fmesh  @p cmesh are stored inside this class.
+//*          Therefore, changing these meshes outside requires also
+//*          to call method @p Derive_Matrix_Pattern explicitely.
+//*
+//* @see Derive_Matrix_Pattern
+//*/
+//Prolongation(Mesh const & cmesh, Mesh const & fmesh);
+
+///**
+//* Destructor.
+//*/
+//~Prolongation() override
+//{}
+
+///**
+//* Generates the sparse matrix pattern and overwrites the existing pattern.
+//*
+//* The sparse matrix pattern is generated but the values are 0.
+//*/
+//void Derive_Matrix_Pattern() override;
+
+//private:
+//Mesh const & _cmesh;      //!< reference to coarse discretization
+//Mesh const & _fmesh;      //!< reference to fine discretization
+//};
+
+// *********************************************************************
+
+
+/**
+ * Interpolation matrix for prolongation coarse mesh (C) to a fine mesh (F)
+ * generated by bisecting edges.
+ *
+ * All interpolation weights are 0.5 (injection points contribute twice).
+ */
+class BisectInterpolation: public Matrix
+{
+public:
+    /**
+     * Generates the interpolation matrix for prolongation coarse mesh to a fine mesh
+     * generated by bisecting edges.
+     * The interpolation weights are all 0.5.
+     *
+     * @param[in] fathers vector[nnodes][2] containing
+     *                    the two coarse grid fathers of a fine grid vertex
+     *
+     */
+    explicit BisectInterpolation(std::vector<int> const &fathers);
+    BisectInterpolation();
+
+    BisectInterpolation(BisectInterpolation const &) = default;
+
+    /**
+     * Destructor.
+     */
+    ~BisectInterpolation() override;
+
+    /**
+     * Extracts the diagonal elements of the matrix.
+     *
+     * @param[in,out]  d  (prellocated) vector of diagonal elements
+     */
+    void GetDiag(std::vector<double> &d) const override;
+    ///**
+    //* Extracts the diagonal elements of the sparse matrix.
+    //*
+    //* @return  d  vector of diagonal elements
+    //*/
+    //std::vector<double> const & GetDiag() const override;
+
+    /**
+     * Performs the prolongation  @f$ w_F := P*u_C @f$.
+     *
+     * @param[in,out] wf resulting fine vector (preallocated)
+     * @param[in]     uc coarse vector
+     */
+    void Mult(std::vector<double> &wf, std::vector<double> const &uc) const override;
+
+    /**
+     * Performs the restriction  @f$ u_C := P^T*w_F @f$.
+     *
+     * @param[in]         wf fine vector
+     * @param[in,out]     uc resulting coarse vector (preallocated)
+     */
+    void MultT(std::vector<double> const &wf, std::vector<double> &uc) const;
+
+    /**
+    * Calculates the defect/residuum w := f - P*u.
+    *
+    * @param[in,out] w resulting vector (preallocated)
+    * @param[in]     f load vector
+    * @param[in]     u coarse vector
+    */
+    void Defect(std::vector<double> &w,
+                std::vector<double> const &f, std::vector<double> const &u) const override;
+
+    /**
+     * Show the matrix entries.
+     */
+    void Debug() const override;
+
+    /**
+     * Finds in this matrix the access index for an entry at row @p row and column @p col.
+     *
+     * @param[in] row    row index
+     * @param[in] col    column index
+     * @return index for element (@p row, @p col). If no appropriate entry exists then -1 will be returned.
+     *
+     * @warning assert() stops the function in case that matrix element (@p row, @p col) doesn't exist.
+    */
+    int fetch(int row, int col) const override;
+
+    /**
+     * Calculates the defect and projects it to the next coarser level @f$ f_C := P^T \cdot (f_F - SK\cdot u_F) @f$.
+     *
+     * @param[in] SK matrix on fine mesh
+     * @param[in] P      prolongation operator
+     * @param[in,out] fc  resulting coarse mesh vector (preallocated)
+     * @param[in] ff r.h.s. on fine mesh
+     * @param[in] uf status vector on fine mesh
+     *
+    */
+    friend void DefectRestrict(CRS_Matrix const &SK, BisectInterpolation const &P,
+                               std::vector<double> &fc, std::vector<double> &ff, std::vector<double> &uf);
+
+protected:
+    std::vector<int> _iv;     //!< fathers[nnode][2] of fine grid nodes, double entries denote injection points
+    std::vector<double> _vv;  //!< weights[nnode][2] of fathers for grid nodes
+};
+
+/**
+ * Interpolation matrix for prolongation from coarse mesh (C)) to a fine mesh (F)
+ * generated by bisecting edges.
+ *
+ * We take into account that values at Dirichlet nodes have to be preserved, i.e.,
+ * @f$ w_F = P \cdot I_D \cdot w_C @f$ and @f$ d_C = I_D  \cdot P^T \cdot  d_F@f$
+ * with @f$ I_D @f$ as @f$ n_C \times n_C @f$ diagonal matrix and entries
+ * @f$ I_{D(j,j)} := \left\{{\begin{array}{l@{\qquad}l} 0 & x_{j}\;\;  \textrm{is Dirichlet node} \\ 1 & \textrm{else} \end{array}}\right. @f$
+ *
+ * Interpolation weights are eighter 0.5 or 0.0 in case of coarse Dirichlet nodes
+ * (injection points contribute twice),
+ * Sets weight to zero iff (at least) one father nodes is a Dirichlet node.
+ */
+class BisectIntDirichlet: public BisectInterpolation
+{
+public:
+    /**
+     * Default constructor.
+    */
+    BisectIntDirichlet()
+        : BisectInterpolation()
+    {}
+
+    /**
+     * Constructs interpolation from father-@p row and column @p col.
+     *
+     * @param[in] fathers    two father nodes from each fine node [nnode_f*2].
+     * @param[in] idxc_dir   vector containing the indices of coarse mesh Dirichlet nodes.
+     *
+    */
+    BisectIntDirichlet(std::vector<int> const &fathers, std::vector<int> const &idxc_dir);
+
+    BisectIntDirichlet(BisectIntDirichlet const &) = default;
+
+
+    /**
+     * Destructor.
+     */
+    ~BisectIntDirichlet() override;
+};
+
+
+
+// *********************************************************************
+
+/**
+ * Calculates the element stiffness matrix @p ske and the element load vector @p fe
+ * of one triangular element with linear shape functions.
+ * @param[in]   ial node indices of the three element vertices
+ * @param[in]   xc  vector of node coordinates with x(2*k,2*k+1) as coordinates of node k
+ * @param[out] ske  element stiffness matrix
+ * @param[out] fe   element load vector
+ */
+void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]);
+
+/**
+ * Adds the element stiffness matrix @p ske and the element load vector @p fe
+ * of one triangular element with linear shape functions to the appropriate positions in
+ * the symmetric stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik)
+ *
+ * @param[in] ial   node indices of the three element vertices
+ * @param[in] ske   element stiffness matrix
+ * @param[in] fe    element load vector
+ * @param[out] sk   vector non-zero entries of CSR matrix
+ * @param[in] id    index vector containing the first entry in a CSR row
+ * @param[in] ik    column index vector of CSR matrix
+ * @param[out] f    distributed local vector storing the right hand side
+ *
+ * @warning Algorithm requires indices in connectivity @p ial in ascending order.
+ *          Currently deprecated.
+*/
+void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
+             int const id[], int const ik[], double sk[], double f[]);
+
+
+
+#endif