diff --git a/Sheet7/Ex_9to13/accu_template/Makefile b/Sheet7/Ex_9to13/accu_template/Makefile new file mode 100644 index 0000000..d8de97f --- /dev/null +++ b/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 diff --git a/Sheet7/Ex_9to13/accu_template/geom.h b/Sheet7/Ex_9to13/accu_template/geom.h new file mode 100644 index 0000000..22dd912 --- /dev/null +++ b/Sheet7/Ex_9to13/accu_template/geom.h @@ -0,0 +1,712 @@ +#ifndef GEOM_FILE +#define GEOM_FILE +#include +#include // function; C++11 +#include +#include // shared_ptr +#include +#include + +/** + * 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 + * weak-vtables. + */ + 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 &GetConnectivity() const + { + return _ia; + } + + /** + * Access/Change connectivity information (g1,g2,g3)_i. + * @return connectivity vector [nelems*ndofs]. + */ + std::vector &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 &GetCoords() const + { + return _xc; + } + + /** + * Access/Change coordinates of vertices (x,y)_i. + * @return coordinates vector [nnodes*2]. + */ + std::vector &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 &v, const std::function &func) const; + void SetBoundaryValues(std::vector &v, const std::function &func) const; + void SetDirchletValues(std::vector &v, const std::function &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 Index_DirichletNodes() const; + virtual std::vector 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 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 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 &GetEdgeConnectivity() const + { + return _ea; + } + + /** + * Access/Change edge connectivity information (e1,e2,e3)_i. + * @return edge connectivity vector [nelems*_nedge_e]. + */ + std::vector &GetEdgeConnectivity() + { + return _ea; + } + + /** + * Read edge information (v1,v2)_i. + * @return edge connectivity vector [_nedge*2]. + */ + const std::vector &GetEdges() const + { + return _edges; + } + + /** + * Access/Change edge information (v1,v2)_i. + * @return edge connectivity vector [_nedge*2]. + */ + std::vector &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> 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 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(_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 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> 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> 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 _ia; //!< element connectivity + std::vector _xc; //!< coordinates + +protected: + // B.C. + std::vector _bedges; //!< boundary edges [nbedges][2] storing start/end vertex +// 2020-01-08 + std::vector _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 _edges; //!< edges of mesh (vertices ordered ascending) + std::vector _ea; //!< edge based element connectivity + // B.C. + std::vector _ebedges; //!< boundary edges [nbedges] + +private: + const std::vector _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 const &ibref = std::vector(0)); + RefinedMesh(Mesh const &cmesh, std::vector const &ibref); + //RefinedMesh(Mesh const &cmesh, std::vector 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(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 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 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 const &old2new); + + +private: + //Mesh const & _cmesh; //!< coarse mesh + std::vector const _ibref; //!< refinement info + int _nref; //!< number of regular refinements performed + std::vector _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> _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 &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 &f) const; + + /** + * Determines the indices of those vertices with Dirichlet boundary conditions + * @return index vector. + */ + std::vector Index_DirichletNodes() const override; + std::vector 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 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 _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 diff --git a/Sheet7/Ex_9to13/accu_template/main.cpp b/Sheet7/Ex_9to13/accu_template/main.cpp new file mode 100644 index 0000000..fa43227 --- /dev/null +++ b/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 +#include +#include +#include // MPI +//#include // 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 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 xl_int(mesh.Nnodes(), 1); + + // check accumulation (by console output) + mesh.VecAccu(xl_int); + + vector 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 xl_new(mesh.Nnodes(), 1.0); + + mesh.Average(xl_new); + + + + + + + + + + + + + MPI_Finalize(); + return 0; +} + + diff --git a/Sheet7/Ex_9to13/accu_template/par_geom.cpp b/Sheet7/Ex_9to13/accu_template/par_geom.cpp new file mode 100644 index 0000000..786d817 --- /dev/null +++ b/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 +#include +#include +#include +#include // contains clock() +#include +#include +#include +#include // accumulate() +#include +#include + +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 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 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 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 const &t2d) +{ + // number of local elements + const int l_ne = count(t2d.cbegin(), t2d.cend(), myrank); + //cout << myrank << ":: " << lne << endl; + vector 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 l_bedges; + vector 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 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(_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 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 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 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 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 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> 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 sendbuf (MPI_IN_PLACE: used also as recvbuf) + // vector sendcounts (the same as for recv) + // vector sdispls (the same as for recv) + // + // We need to map the interface vertices onto the sendbuffer: + // vector loc_itf local index of interface vertex lk + // vector gloc_itf global index of interface vertex lk + // vector buf2loc local indices of sendbuffer positions (the same as for recv) + + // ---- Determine sendcounts[] and sdipls[] from neigh_itf[] + //vector _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 _sdispls(NumProcs(), 0); + _sdispls.resize(NumProcs(), 0); + partial_sum(_sendcounts.cbegin(), _sendcounts.cend() - 1, _sdispls.begin() + 1); + //vector _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 _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 _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 x(Nnodes(),-1.0); + VecAccu(x); + cout << " VecAccu (InPlace) successful!\n"; + + + return; +} + +bool ParMesh::CheckInterfaceExchange_InPlace() const +{ + vector 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 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) 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 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 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 bool + { return std::abs(a-b) x(Nnodes(),-1.0); + for (size_t i=0; i 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 bool + { return std::abs(a-b) x(Nnodes(),-1.0); + for (size_t i=0; i 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 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 bool + { return std::abs(a-b) &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 &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 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 &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]; + } +} + diff --git a/Sheet7/Ex_9to13/accu_template/par_geom.h b/Sheet7/Ex_9to13/accu_template/par_geom.h new file mode 100644 index 0000000..f134056 --- /dev/null +++ b/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 +#include // function; C++11 +#include +#include +#include // shared_ptr +#include // MPI +#include +#include + +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 + * weak-vtables. + */ + 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 &w) const; + void VecAccu(std::vector& x) const; + + /** Inner product + * @param[in] x vector + * @param[in] y vector + * @return resulting Euclidian inner product + */ + double dscapr(std::vector const &x, std::vector const &y) const + { + return par_scalar(x, y, _icomm); + } + + int GlobalNodes() const; + void Average(std::vector &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 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 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 _v_l2g; //!< vertices: local to global mapping + std::vector _t_l2g; //!< triangles: local to global mapping + std::map _v_g2l; //!< vertices: global to local mapping + std::map _t_g2l; //!< triangles: global to local mapping + + //std::vector e_l2g; //!< edges: local to global mapping + + std::vector _valence; //!< valence of local vertices, i.e. number of subdomains they belong to + // MPI_Alltoallv needs: + mutable std::vector _sendbuf; //!< send buffer a n d receiving buffer (MPI_IN_PLACE) + std::vector _sendcounts; //!< number of data to send to each MPI rank (the same as for recv) + std::vector _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 _loc_itf; //!< local index of interface vertex lk + std::vector _gloc_itf; //!< global index of interface vertex lk + std::vector _buf2loc; //!< local indices of sendbuffer positions (the same as for recv) + + +}; + + +#endif