.

Sheet7/Ex_9to13/accu_template/Makefile

@@ -0,0 +1,54 @@
+#
+# use GNU-Compiler tools
+COMPILER=GCC_
+# COMPILER=GCC_SEQ_
+# alternatively from the shell
+# export COMPILER=GCC_
+# or, alternatively from the shell
+# make COMPILER=GCC_
+
+MAIN = main
+SOURCES = ${MAIN}.cpp vdop.cpp geom.cpp par_geom.cpp
+
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PROGRAM	= ${MAIN}.${COMPILER}
+
+# uncomment the next to lines for debugging and detailed performance analysis
+CXXFLAGS += -g
+# -DNDEBUG
+# -pg slows down the code on my laptop when using CLANG_
+LINKFLAGS += -g 
+#-pg
+#CXXFLAGS += -Q --help=optimizers
+#CXXFLAGS += -fopt-info
+
+include ../${COMPILER}default.mk
+
+#############################################################################
+# additional specific cleaning in this directory
+clean_all::
+	@rm -f uv.txt
+
+
+#############################################################################
+# special testing
+# NPROCS	= 4
+#
+TFILE	= t.dat
+# TTMP	= t.tmp
+#
+graph: $(PROGRAM)
+# 	@rm -f $(TFILE).*
+	# next two lines only sequentially
+	./$(PROGRAM)
+	@mv  $(TFILE).000 $(TFILE)
+# 	$(MPIRUN) $(MPIFLAGS) -np $(NPROCS) $(PROGRAM)
+# 	@echo " "; echo "Manipulate data for graphics."; echo " "
+# 	@cat $(TFILE).* > $(TTMP)
+# 	@sort -b -k 2    $(TTMP)   -o $(TTMP).1
+# 	@sort -b -k 1    $(TTMP).1 -o $(TTMP).2
+# 	@awk  -f nl.awk  $(TTMP).2  > $(TFILE)
+# 	@rm -f $(TTMP).* $(TTMP) $(TFILE).*
+#
+	-gnuplot jac.dem

Sheet7/Ex_9to13/accu_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.
+     */
+    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/Ex_9to13/accu_template/main.cpp

@@ -0,0 +1,108 @@
+//		MPI code in C++.
+//		See [Gropp/Lusk/Skjellum, "Using MPI", p.33/41 etc.]
+//		and  /opt/mpich/include/mpi2c++/comm.h  for details
+
+#include "geom.h"
+#include "par_geom.h"
+#include "vdop.h"
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <mpi.h>            // MPI
+//#include <omp.h>            // OpenMP, for E9 it is not necessary
+using namespace std;
+
+
+int main(int argc, char **argv )
+{
+    MPI_Init(&argc, &argv);
+    MPI_Comm const icomm(MPI_COMM_WORLD);
+    //omp_set_num_threads(1);                   // don't use OMP parallelization for a start
+//
+   {
+       int np;
+       MPI_Comm_size(icomm, &np);
+
+       assert(4 == np);                  // example is only provided for 4 MPI processes
+   }
+// #####################################################################
+// ---- Read the f.e. mesh and the mapping of elements to MPI processes
+    //Mesh const mesh_c("square_4.txt");    //    Files square_4.txt and square_4_sd.txt  are needed
+    ParMesh const mesh("square",icomm);
+    
+    int const numprocs = mesh.NumProcs();
+    int const myrank   = mesh.MyRank();
+    if ( 0 == myrank ) {
+        cout << "\n There are " << numprocs << " processes running.\n \n";
+    }
+        
+    int const check_rank=0;               // choose the MPI process you would like to check the mesh
+    //if ( check_rank == myrank ) mesh.Debug();
+    //if ( check_rank == myrank ) mesh.DebugEdgeBased();
+    
+// ---- allocate local vectors and check skalar product and vector accumulation
+    vector<double> xl(mesh.Nnodes(), 1.0);
+    mesh.SetValues(xl, [](double x, double y) -> double {return x * x * std::sin(2.5 * M_PI * y);} );
+    
+    //E9
+    if (myrank == check_rank){
+        cout << "E9" << endl;
+    }
+
+    double ss = mesh.dscapr(xl,xl);
+    if (myrank == check_rank){
+        cout << myrank << " : scalar : " << ss << endl << endl;
+    }
+
+    mesh.VecAccu(xl);
+
+    //E10
+    if (myrank == check_rank){
+        cout << "E10" << endl;
+    }
+    vector<int> xl_int(mesh.Nnodes(), 1);
+
+    // check accumulation (by console output)
+    mesh.VecAccu(xl_int);
+
+    vector<double> coords = mesh.GetCoords();
+    if (check_rank == myrank){
+        for (size_t i = 0; i < coords.size(); i += 2){
+            cout << "(" << coords[i] << ", " << coords[i + 1] << "):\t" << xl_int[i/2] << endl;
+        }
+    }
+
+    //E11
+    if (myrank == check_rank){
+        cout << "E11" << endl;
+    }
+    int global_nodes = mesh.GlobalNodes();
+    if (check_rank == myrank){
+        cout << "Global nodes: " << global_nodes << endl;
+    }
+
+    //E12
+    if (myrank == check_rank){
+        cout << "E12" << endl;
+    }
+    vector<double> xl_new(mesh.Nnodes(), 1.0);
+
+    mesh.Average(xl_new);
+
+
+
+
+
+
+
+
+
+
+
+
+    MPI_Finalize();
+    return 0;
+}
+
+

Sheet7/Ex_9to13/accu_template/par_geom.cpp

@@ -0,0 +1,590 @@
+// see:   http://llvm.org/docs/CodingStandards.html#include-style
+#include "vdop.h"
+//#include "geom.h"
+#include "par_geom.h"
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <ctime>                  // contains clock()
+#include <fstream>
+#include <iostream>
+#include <list>
+#include <numeric>                // accumulate()
+#include <string>
+#include <vector>
+
+using namespace std;
+
+
+ParMesh::ParMesh(int ndim, int nvert_e, int ndof_e, int nedge_e, MPI_Comm const &icomm)
+    : Mesh(ndim, nvert_e, ndof_e, nedge_e),
+      _icomm(icomm), _numprocs(-1), _myrank(-1),
+      _v_l2g(0), _t_l2g(0), _v_g2l{{}}, _t_g2l{{}}, _valence(0),
+      _sendbuf(0), _sendcounts(0), _sdispls(0),
+      _loc_itf(0), _gloc_itf(0), _buf2loc(0)
+{
+    MPI_Comm_size(icomm, &_numprocs);
+    MPI_Comm_rank(icomm, &_myrank);
+}
+
+ParMesh::~ParMesh()
+{}
+
+
+
+ParMesh::ParMesh(std::string const &sname, MPI_Comm const &icomm)
+    : ParMesh(2, 3, 3, 3, icomm) // two dimensions, 3 vertices, 3 dofs, 3 edges per element
+{
+    //const int numprocs = _icomm.Get_size();
+    const string NS    = "_" + to_string(_numprocs);
+    const string fname = sname + NS + ".txt";
+    //cout << "############  " << fname << endl;
+    ReadVertexBasedMesh(fname);
+    cout << "\n End of sequential File read \n";
+    // ------------------------------------------------------------------------------
+    // Until this point  a l l  processes possess  a l l  mesh info in  g l o b a l  numbering
+    //
+    // Now, we have to select the data belonging to my_rank
+    // and we have to create the mapping local to global (l2g) and vice versa (g2l)
+    // ------------------------------------------------------------------------------
+
+    // save the global node mesh (maybe we need it later)
+    DeriveEdgeFromVertexBased();                     //                    and even more
+    Mesh global_mesh(*this);                         // requires a  l o t  of memory
+    Del_EdgeConnectivity();
+
+    // read the subdomain info
+    const string dname = sname + NS + "_sd" + ".txt";
+    vector<int> t2d = ReadElementSubdomains(dname);  // global mapping triangle to subdomain for all elements
+
+    //const int myrank   = _icomm.Get_rank();
+    Transform_Local2Global_Vertex(_myrank, t2d);      // Vertex based mesh: now in  l o c a l  indexing
+
+    DeriveEdgeFromVertexBased();                     // Generate also the  l o c a l  edge based information
+
+    Generate_VectorAdd();
+
+
+    // Now we have to organize the MPI communication of vertices on the subdomain interfaces
+
+    return;
+}
+
+vector<int> ParMesh::ReadElementSubdomains(string const &dname)
+{
+    ifstream ifs(dname);
+    if (!(ifs.is_open() && ifs.good())) {
+        cerr << "ParMesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
+        assert(ifs.is_open());
+    }
+
+    int const OFFSET{1};             // Matlab to C indexing
+    cout << "ASCI file  " << dname << "  opened" << endl;
+
+    // Read some mesh constants
+    int nelem;
+    ifs >> nelem;
+    cout << nelem << "  " << Nelems() << endl;
+    assert( Nelems() == nelem);
+
+    // Allocate memory
+    vector<int> t2d(nelem, -1);
+    // Read element mapping
+    for (int k = 0; k < nelem; ++k) {
+        int tmp;
+        ifs >> tmp;
+        //t2d[k] = tmp - OFFSET;
+        // 2020-01-08 
+        t2d[k] = min(tmp, NumProcs()) - OFFSET;
+    }
+
+    return t2d;
+}
+
+void ParMesh::Transform_Local2Global_Vertex(int const myrank, vector<int> const &t2d)
+{
+    // number of local elements
+    const int l_ne = count(t2d.cbegin(), t2d.cend(), myrank);
+    //cout << myrank <<  "::  " << lne << endl;
+    vector<int> l_ia(l_ne * NverticesElements(), -1); // local elements still with global vertex numbers
+    _t_l2g.resize(l_ne, -1);
+
+    int lk = 0;
+    for (size_t k = 0; k < t2d.size(); ++k) {
+        if (myrank == t2d[k]) {
+            //if (0==myrank)
+            //{
+            //cout << lk << "   k  " <<  t2d[k] << endl;
+            //}
+            l_ia[3 * lk  ] = _ia[3 * k  ];
+            l_ia[3 * lk + 1] = _ia[3 * k + 1];
+            l_ia[3 * lk + 2] = _ia[3 * k + 2]; // local elements still with global vertex numbers
+            _t_l2g[lk] = k;                       // elements: local to global mapping
+            _t_g2l[k] = lk;                       //           global to local
+            ++lk;
+        }
+    }
+    // Checks:
+    assert( count(l_ia.cbegin(), l_ia.cend(), -1)   == 0 );
+    assert( count(_t_l2g.cbegin(), _t_l2g.cend(), -1) == 0 );
+
+    // Vertices: local to global mapping
+    auto tmp = l_ia;
+    sort(tmp.begin(), tmp.end());
+    auto ip = unique(tmp.begin(), tmp.end());
+    tmp.erase(ip, tmp.end());
+    _v_l2g = tmp;                                 // Vertices: local to global mapping
+    for (size_t lkv = 0; lkv < _v_l2g.size(); ++lkv) {
+        _v_g2l[_v_l2g[lkv]] = lkv;                //           global to local
+    }
+
+    // Boundary edges
+    vector<int> l_bedges;
+    vector<int> l_sdedges;
+    for (size_t b = 0; b < _bedges.size(); b += 2) {
+        int const v1 = _bedges[b   ];             // global vertex numbers
+        int const v2 = _bedges[b + 1];
+        try {
+            int const lv1 = _v_g2l.at(v1);         // map[] would add that element
+            int const lv2 = _v_g2l.at(v2);         //       but at() throws an exeption
+            l_bedges.push_back(lv1);
+            l_bedges.push_back(lv2);               // Boundaries: already in local indexing
+            // 2020-01-08
+            l_sdedges.push_back(_sdedges[b  ]);
+            l_sdedges.push_back(_sdedges[b+1]);
+        }
+        catch (std::out_of_range & err) {
+            //cerr << ".";
+        }
+    }
+
+    // number of local vertices
+    const int l_nn = _v_l2g.size();
+    vector<double> l_xc(Ndims()*l_nn);
+    for (int lkk = 0; lkk < l_nn; ++lkk) {
+        int k = _v_l2g.at(lkk);
+        l_xc[2 * lkk  ]   = _xc[2 * k  ];
+        l_xc[2 * lkk + 1] = _xc[2 * k + 1];
+    }
+
+
+    // Now, we represent the vertex mesh in  l o c a l  numbering
+    // elements
+
+    for (size_t i = 0; i < l_ia.size(); ++i) {
+        l_ia[i] = _v_g2l.at(l_ia[i]);              // element vertices: global to local
+    }
+    SetNelem(l_ne);
+    _ia = l_ia;
+    // boundary
+    _bedges = l_bedges;
+    _sdedges = l_sdedges;
+    // coordinates
+    SetNnode(l_nn);
+    _xc = l_xc;
+
+    return;
+}
+
+
+void ParMesh::Generate_VectorAdd()
+{
+    // Some checks
+    int lnn = Nnodes();                           // local number of vertices
+    assert(static_cast<int>(_v_l2g.size()) == lnn);
+    int ierr{-12345};
+
+    // ---- Determine global largest vertex index
+    int gidx_max{-1};                             // global largest vertex index
+    int lmax = *max_element(_v_l2g.cbegin(), _v_l2g.cend());
+    MPI_Allreduce(&lmax, &gidx_max, 1, MPI_INT, MPI_MAX, _icomm);
+    int gidx_min{-1};                             // global smallest vertex index
+    int lmin = *min_element(_v_l2g.cbegin(), _v_l2g.cend());
+    MPI_Allreduce(&lmin, &gidx_min, 1, MPI_INT, MPI_MIN, _icomm);
+    //cout << gidx_min << "  " << gidx_max << endl;
+    assert(0 == gidx_min);                        // global indices have to start with 0
+
+
+    // ---- Determine for all global vertices the number of subdomains it belongs to
+    vector<int> global(gidx_max+1, 0);              // global scalar array for vertices
+    for (auto const gidx : _v_l2g)  global[gidx] = 1;
+    // https://www.mpi-forum.org/docs/mpi-2.2/mpi22-report/node109.htm
+    ierr = MPI_Allreduce(MPI_IN_PLACE, global.data(), global.size(), MPI_INT, MPI_SUM, _icomm);
+    //if (0 == MyRank())  cout << global << endl;
+    //MPI_Barrier(_icomm);
+    //cout << _xc[2*_v_g2l.at(2)] << " , " << _xc[2*_v_g2l.at(2)+1] << endl;
+    //MPI_Barrier(_icomm);
+
+    //  now, global[] contains the number of subdomains a global vertex belongs to
+    if ( count(global.cbegin(), global.cend(), 0) > 0 )
+        cerr << "\n !!!   Non-continuous global vertex indexing  !!!\n";
+
+    // ---- Determine local interface vertices ( <==> global[] > 1 )
+    //      _loc_itf,  neigh_itf
+    //vector<int> loc_itf;                          // local indices of interface vertices on this MPI process
+    for (size_t lk = 0; lk < _v_l2g.size(); ++lk) {
+        int const gk = _v_l2g[lk];                // global index of local vertex lk
+        if ( global[gk] > 1 ) {
+            _loc_itf.push_back(lk);               // local indices of interface vertices on this MPI process
+        }
+    }
+
+    //MPI_Barrier(_icomm);
+    //if (0 == MyRank())  cout << "\n..._loc_itf...\n" << _loc_itf << "\n......\n";
+    //MPI_Barrier(_icomm);
+    // ---- global indices of local interface vertices
+    //auto gloc_itf(_loc_itf);
+    _gloc_itf=_loc_itf;
+    for_each(_gloc_itf.begin(), _gloc_itf.end(), [this] (auto & v) -> void { v = _v_l2g[v];} );
+    //MPI_Barrier(_icomm);
+    //if (0 == MyRank())  cout << "\n..._gloc_itf...\n" << _gloc_itf << "\n......\n";
+    //DebugVector(_gloc_itf,"_gloc_itf");
+
+    // ---- Determine the global length of interfaces
+    vector<int> vnn(NumProcs(), -1);              // number of interface vertices per MPI rank
+    int l_itf(_loc_itf.size());                    // # local interface vertices
+    ierr = MPI_Allgather(&l_itf, 1, MPI_INT, vnn.data(), 1, MPI_INT, _icomm);
+    assert(0 == ierr);
+    //cout << vnn << endl;
+
+    // ---- Now we consider only the inferface vertices
+    int snn = accumulate(vnn.cbegin(), vnn.cend(), 0); // required length of array for global interface indices
+    //cout << snn << "   " << gnn << endl;
+    vector<int> dispnn(NumProcs(), 0) ;           // displacement of interface vertices per MPI rank
+    partial_sum(vnn.cbegin(), vnn.cend() - 1, dispnn.begin() + 1);
+    //cout << dispnn << endl;
+
+    // ---- Get the global indices for all global interfaces
+    vector<int> g_itf(snn, -1);                    // collects all global indices of the global interfaces
+    // https://www.mpich.org/static//docs/v3.0.x/www3/MPI_Gatherv.html
+    ierr = MPI_Gatherv( _gloc_itf.data(), _gloc_itf.size(), MPI_INT,
+                        g_itf.data(), vnn.data(), dispnn.data(), MPI_INT, 0, _icomm);
+    assert(0 == ierr);
+    // https://www.mpich.org/static/docs/v3.1/www3/MPI_Bcast.html
+    ierr = MPI_Bcast(g_itf.data(), g_itf.size(), MPI_INT, 0, _icomm);
+    assert(0 == ierr);                            // Now, each MPI rank has the all global indices of the global interfaces
+    //MPI_Barrier(_icomm);
+    //if (MyRank() == 0) cout << "\n...g_itf...\n" << g_itf << "\n......\n";
+    //MPI_Barrier(_icomm);
+
+    // ----- Determine all MPI ranks a local interface vertex belongs to
+    vector<vector<int>> neigh_itf(_loc_itf.size());// subdomains a local interface vertex belongs to
+    for (size_t lk = 0; lk < _loc_itf.size(); ++lk) {
+        const int gvert = _gloc_itf[lk];           // global index of local interface node lk
+        for (int rank = 0; rank < NumProcs(); ++rank) {
+            auto const startl = g_itf.cbegin() + dispnn[rank];
+            auto const endl   = startl + vnn[rank];
+            if ( find( startl, endl, gvert) != endl) {
+                neigh_itf[lk].push_back(rank);
+            }
+        }
+    }
+
+    // ---- check the available info in  _loc_itf[lk], _gloc_itf[lk], neigh_itf[lk]
+    //MPI_Barrier(_icomm);
+    ////if (MyRank()==0)  cout << "\n...neigh_itf ...\n"  << neigh_itf << endl;
+    //if (MyRank() == 0) {
+        //for (size_t lk = 0; lk < _loc_itf.size(); ++lk ) {
+            //cout << lk << " : local idx " << _loc_itf[lk] << " , global idx " <<  _gloc_itf[lk];
+            //cout << "  with MPI ranks " << neigh_itf[lk] << endl;
+        //}
+    //}
+    //MPI_Barrier(_icomm);
+    
+    // ---- store the valence (e.g., the number of subdomains it belongs to) of all local vertices
+    _valence.resize(Nnodes(),1);
+    for (size_t lk = 0; lk < _loc_itf.size(); ++lk)
+    {
+        _valence[_loc_itf[lk]] = neigh_itf[lk].size();
+    }
+    //DebugVector(_valence,"_valence",_icomm);
+
+    // ---- We ware going to use MPI_Alltoallv for data exchange on interfaces
+    //      https://www.mpi-forum.org/docs/mpi-3.1/mpi31-report/node109.htm#Node109
+    //      https://www.open-mpi.org/doc/v4.0/man3/MPI_Alltoallv.3.php
+    //int MPI_Alltoallv(const void* sendbuf, const int sendcounts[], const int sdispls[], MPI_Datatype sendtype, void* recvbuf, const int recvcounts[], const int rdispls[], MPI_Datatype recvtype, MPI_Comm comm)
+    //
+    //  MPI_Alltoallv needs:
+    //      vector<double> sendbuf             (MPI_IN_PLACE: used also as recvbuf)
+    //      vector<int> sendcounts             (the same as for recv)
+    //      vector<int> sdispls                (the same as for recv)
+    //
+    //  We need to map the interface vertices onto the sendbuffer:
+    //      vector<int> loc_itf                 local  index of interface vertex lk
+    //      vector<int> gloc_itf                global index of interface vertex lk
+    //      vector<int> buf2loc                 local indices of sendbuffer positions (the same as for recv)
+
+    // ---- Determine sendcounts[] and sdipls[] from neigh_itf[]
+    //vector<int> _sendcounts(NumProcs(), 0);
+    _sendcounts.resize(NumProcs(), 0);
+    for (size_t lk = 0; lk < _loc_itf.size(); ++lk ) {
+        auto const &kneigh = neigh_itf[lk];
+        for (size_t ns = 0; ns < kneigh.size(); ++ns) {
+            ++_sendcounts[kneigh[ns]];
+        }
+    }
+    //if (MyRank() == 0)  cout << "\n..._sendcounts ...\n"  << _sendcounts << endl;
+
+    //vector<int> _sdispls(NumProcs(), 0);
+    _sdispls.resize(NumProcs(), 0);
+    partial_sum(_sendcounts.cbegin(), _sendcounts.cend() - 1, _sdispls.begin() + 1);
+    //vector<int> _sdispls(NumProcs()+1, 0);
+    //partial_sum(_sendcounts.cbegin(), _sendcounts.cend(), _sdispls.begin() + 1);
+    //if (MyRank() == 0)  cout << "\n..._sdispls ...\n"  << _sdispls << endl;
+
+    // ---- Determine size of buffer 'nbuffer' and mapping 'buf2loc'
+    int const nbuffer = accumulate(_sendcounts.cbegin(), _sendcounts.cend(), 0);
+    //vector<int> _buf2loc(nbuffer, -1);
+    _buf2loc.resize(nbuffer, -1);
+    int buf_idx = 0;                              // position in buffer
+    for (int rank = 0; rank < NumProcs(); ++rank) {
+        assert( buf_idx ==  _sdispls[rank]);
+        for (size_t lk = 0; lk < _loc_itf.size(); ++lk ) {
+            auto const &kneigh = neigh_itf[lk];
+            if (find(kneigh.cbegin(),kneigh.cend(),rank)!=kneigh.cend())
+            {
+               _buf2loc[buf_idx] =  _loc_itf[lk];
+               ++buf_idx;
+            }
+        }
+    }
+    //if (MyRank() == 0)  cout << "\n...buf2loc ...\n"  << buf2loc << endl;
+    //DebugVector(buf2loc,"buf2loc",_icomm);
+
+    // ---- Allocate send/recv buffer
+    //vector<double> _sendbuf(nbuffer,-1.0);
+    _sendbuf.resize(nbuffer,-1.0);
+
+    assert(CheckInterfaceExchange_InPlace());
+    cout << " Check of data exchange (InPlace) successful!\n";
+    assert(CheckInterfaceExchange());
+    cout << " Check of data exchange successful!\n";
+    assert(CheckInterfaceAdd_InPlace());
+    cout << " Check of data add successful!\n";
+    assert(CheckInterfaceAdd());
+    cout << " Check of data add (InPlace) successful!\n";
+    
+    vector<double> x(Nnodes(),-1.0);
+    VecAccu(x);
+    cout << " VecAccu (InPlace) successful!\n";
+
+    
+    return;
+}
+
+bool ParMesh::CheckInterfaceExchange_InPlace() const
+{
+    vector<double> x(Nnodes(),-1.0);
+    copy(_v_l2g.cbegin(),_v_l2g.cend(),x.begin());       // init x with global vertex indices
+
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = x[_buf2loc.at(ls)];
+    }
+    int ierr = MPI_Alltoallv(MPI_IN_PLACE, _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                          _sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    assert(ierr==0);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+
+    vector<double> y(x);
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) y[_loc_itf.at(lk)] = -1.0;   // only for interface nodes
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        y[_buf2loc.at(ls)] = _sendbuf[ls];
+    }
+
+    double const  eps=1e-10;
+    bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
+                          [eps](double a, double b) -> bool
+                          { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(b))); }
+                   );
+    return bv;
+}
+
+bool ParMesh::CheckInterfaceExchange() const
+{
+    vector<double> x(Nnodes(),-1.0);
+    copy(_v_l2g.cbegin(),_v_l2g.cend(),x.begin());       // init x with global vertex indices
+
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = x[_buf2loc.at(ls)];
+    }
+    vector<double> recvbuf(_sendbuf.size());
+    int ierr = MPI_Alltoallv(_sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                              recvbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+    //DebugVector(recvbuf,"recvbuf",_icomm);
+    assert(ierr==0);
+
+    vector<double> y(x);
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) y[_loc_itf.at(lk)] = -1.0;   // only for interface nodes
+    for(size_t ls = 0; ls<recvbuf.size(); ++ls)
+    {
+        y[_buf2loc.at(ls)] = recvbuf[ls];
+    }
+    //cout << "WRONG  : " << count(y.cbegin(),y.cend(), -1.0) << endl;
+
+    double const  eps=1e-10;
+    bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
+                          [eps](double a, double b) -> bool
+                          { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(b))); }
+                   );
+    return bv;
+}
+
+bool ParMesh::CheckInterfaceAdd_InPlace() const
+{
+    vector<double> x(Nnodes(),-1.0);
+    for (size_t i=0; i<x.size(); ++i)
+    {
+        x[i] = _xc[2*i]+_xc[2*i+1];                      // init x with coordinate values
+    }
+
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = x[_buf2loc.at(ls)];
+    }
+    int ierr = MPI_Alltoallv(MPI_IN_PLACE, _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                          _sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    assert(ierr==0);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+
+    vector<double> y(x);
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) y[_loc_itf.at(lk)] = 0.0;   // only for interface nodes
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        y[_buf2loc.at(ls)] += _sendbuf[ls];
+    }
+    MPI_Barrier(_icomm);
+    //DebugVector(x,"x",_icomm);
+    //DebugVector(y,"y",_icomm);
+    for (size_t i= 0; i<y.size(); ++i) y[i]/=_valence[i];   // divide by valence
+
+    double const  eps=1e-10;
+    bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
+                          [eps](double a, double b) -> bool
+                          { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(b))); }
+                   );
+    return bv;
+}
+
+bool ParMesh::CheckInterfaceAdd() const
+{
+    vector<double> x(Nnodes(),-1.0);
+    for (size_t i=0; i<x.size(); ++i)
+    {
+        //x[i] = _xc[2*i]+_xc[2*i+1];                      // init x with coordinate values
+        x[i] = _v_l2g[i];
+    }
+    
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = x[_buf2loc.at(ls)];
+    }
+    vector<double> recvbuf(_sendbuf.size());
+    int ierr = MPI_Alltoallv(_sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                              recvbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+    //DebugVector(recvbuf,"recvbuf",_icomm);
+    assert(ierr==0);
+ 
+    vector<double> y(x);
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) y[_loc_itf.at(lk)] = 0.0;   // only for interface nodes
+    for(size_t ls = 0; ls<recvbuf.size(); ++ls)
+    {
+        //if (0==MyRank()) cout << ls << ": " << _buf2loc.at(ls) << "  " << y[_buf2loc.at(ls)] << "("<< x[_buf2loc.at(ls)] << ")" << "  " << recvbuf[ls] << "  (" <<  _sendbuf[ls] << ")" << endl;
+        y[_buf2loc.at(ls)] += recvbuf[ls];
+    }
+    MPI_Barrier(_icomm);
+    //DebugVector(x,"x",_icomm);
+    //DebugVector(y,"y",_icomm);
+    for (size_t i= 0; i<y.size(); ++i) y[i]/=_valence[i];   // divide by valence
+
+    double const  eps=1e-10;
+    bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
+                          [eps](double a, double b) -> bool
+                          { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(b))); }
+                   );
+    return bv;
+}
+
+
+// ----------
+
+void ParMesh::VecAccu(std::vector<double> &w) const
+{
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = w[_buf2loc.at(ls)];
+    }
+    int ierr = MPI_Alltoallv(MPI_IN_PLACE, _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                          _sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    assert(ierr==0);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) w[_loc_itf.at(lk)] = 0.0;   // only for interface nodes
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        w[_buf2loc.at(ls)] += _sendbuf[ls];
+    }
+
+    return;
+}
+
+void ParMesh::VecAccu(std::vector<int> &w) const
+{
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        _sendbuf[ls] = w[_buf2loc.at(ls)];
+    }
+    int ierr = MPI_Alltoallv(MPI_IN_PLACE, _sendcounts.data(), _sdispls.data(), MPI_DOUBLE,
+                          _sendbuf.data(), _sendcounts.data(), _sdispls.data(), MPI_DOUBLE, _icomm);
+    assert(ierr==0);
+    //DebugVector(_sendbuf,"_sendbuf",_icomm);
+
+    for(size_t lk = 0; lk<_loc_itf.size(); ++lk) w[_loc_itf.at(lk)] = 0.0;   // only for interface nodes
+    for(size_t ls = 0; ls<_sendbuf.size(); ++ls)
+    {
+        w[_buf2loc.at(ls)] += _sendbuf[ls];
+    }
+
+    return;
+}
+
+int ParMesh::GlobalNodes() const
+{
+    int local_max = *max_element(_v_l2g.begin(), _v_l2g.end());
+    int global_max = 0;
+
+    MPI_Allreduce(&local_max, &global_max, 1, MPI_INT, MPI_MAX, _icomm);
+
+    std::vector<int> global_flag(global_max + 1, 0);
+    for (int gidx : _v_l2g){
+        global_flag[gidx] = 1;
+    }
+
+    MPI_Allreduce(MPI_IN_PLACE,global_flag.data(),global_flag.size(),MPI_INT,MPI_SUM,_icomm);
+
+    int global_nodes = 0;
+    for (int v : global_flag){
+        if (v > 0) ++global_nodes;
+    }
+
+    return global_nodes;
+}
+
+
+void ParMesh::Average(std::vector<double> &w) const
+{
+    VecAccu(w);
+
+    for (size_t lk = 0; lk < _loc_itf.size(); ++lk)
+    {
+        int const local_idx = _loc_itf[lk];
+        w[local_idx] /= _valence[local_idx];
+    }
+}
+

Sheet7/Ex_9to13/accu_template/par_geom.h

@@ -0,0 +1,151 @@
+#ifndef PAR_GEOM_FILE
+#define PAR_GEOM_FILE
+#include "geom.h"
+#include "vdop.h"
+#include <array>
+#include <functional>             // function; C++11
+#include <iostream>
+#include <map>
+#include <memory>                 // shared_ptr
+#include <mpi.h>                  // MPI
+#include <string>
+#include <vector>
+
+class ParMesh: public 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)
+      * @param[in] icomm    MPI communicator
+     */
+    explicit ParMesh(int ndim, int nvert_e = 0, int ndof_e = 0, int nedge_e = 0, MPI_Comm const &icomm = MPI_COMM_WORLD);
+
+    ParMesh(ParMesh const &) = default;
+
+    ParMesh &operator=(ParMesh 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 ~ParMesh();
+
+    /**
+     * Reads mesh data from a binary file.
+     *
+     * @param[in] sname suffix of file name
+     * @param[in] icomm  MPI communicator
+     * @see  ascii_write_mesh.m for the file format.
+    */
+    explicit ParMesh(std::string const &sname, MPI_Comm const &icomm = MPI_COMM_WORLD);
+
+    void VecAccu(std::vector<double> &w) const;
+    void VecAccu(std::vector<int>& x) const;
+
+    /**  Inner product
+     * @param[in] x  vector
+     * @param[in] y  vector
+     * @return       resulting Euclidian inner product <x,y>
+    */
+    double dscapr(std::vector<double> const &x, std::vector<double> const &y) const
+    {
+        return par_scalar(x, y, _icomm);
+    }
+
+    int GlobalNodes() const;
+    void Average(std::vector<double> &w) const;
+
+private:
+    /**
+     * Reads the global triangle to subdomain mapping.
+     *
+     * @param[in] dname file name
+     *
+     * @see ascii_write_subdomains.m for the file format
+    */
+    std::vector<int> ReadElementSubdomains(std::string const &dname);
+
+
+    /**
+     * Transform
+      *
+     * @param[in] myrank MPI rank of this process
+      * @param[in] t2d       global mapping triangle to subdomain for all elements (vertex based)
+    */
+    void Transform_Local2Global_Vertex(int myrank, std::vector<int> const &t2d);
+
+
+    /**
+     * Transform
+    */
+    void Generate_VectorAdd();
+
+    bool CheckInterfaceExchange_InPlace() const;
+    bool CheckInterfaceExchange() const;
+    bool CheckInterfaceAdd_InPlace() const;
+    bool CheckInterfaceAdd() const;
+
+
+public:
+    /**     MPI rank of the calling process in communication group.
+    *
+    * @return       MPI rank of the calling process
+    */
+    int MyRank() const
+    {
+        return _myrank;
+    }
+
+    /**  Number of MPI processes in communication group.
+     *
+     * @return       Number of MPI processes
+    */
+    int NumProcs() const
+    {
+        return _numprocs;
+    }
+
+    /**  Returns recent
+     * @return       MPI communicator
+    */
+    MPI_Comm GetCommunicator() const
+    {
+        return _icomm;
+    }
+
+private:
+    // Don't use  &_icomm  ==> Error
+    MPI_Comm const _icomm;                     //!< MPI communicator for the group of processes
+    int _numprocs;                             //!< number of MPI processes
+    int _myrank;                               //!< my MPI rank
+    std::vector<int> _v_l2g;                   //!< vertices:  local to global mapping
+    std::vector<int> _t_l2g;                   //!< triangles: local to global mapping
+    std::map<int, int> _v_g2l;                 //!< vertices:  global to local mapping
+    std::map<int, int> _t_g2l;                 //!< triangles: global to local mapping
+
+    //std::vector<int> e_l2g;                    //!< edges:     local to global mapping
+
+    std::vector<int> _valence;                 //!< valence of local vertices, i.e. number of subdomains they belong to
+    //  MPI_Alltoallv needs:
+    mutable std::vector<double> _sendbuf;      //!< send buffer  a n d  receiving buffer (MPI_IN_PLACE)
+    std::vector<int>    _sendcounts;           //!< number of data to send to each MPI rank (the same as for recv)
+    std::vector<int>    _sdispls;              //!< offset of data to send to each MPI rank wrt. _senbuffer (the same as for recv)
+    //
+    //  We need to map the interface vertices onto the sendbuffer:
+    std::vector<int> _loc_itf;                 //!< local  index of interface vertex lk
+    std::vector<int> _gloc_itf;                //!< global index of interface vertex lk
+    std::vector<int> _buf2loc;                 //!< local indices of sendbuffer positions (the same as for recv)
+
+
+};
+
+
+#endif