Task 5, 5* some fixes and cleanup

sheet5/1/main.cpp

@@ -10,7 +10,7 @@
 #include <cmath>
 using namespace std;
 
-void benchmark(vector<double> &x, vector<double> &y, unsigned int N, unsigned int NLOOPS)
+void benchmark(vector<double> &x, vector<double> &y, unsigned int NLOOPS)
 {
     double sk = 0.0;
     for (int i = 0; i < NLOOPS; ++i)
@@ -92,8 +92,7 @@ int main(int argc, char const *argv[])
     for (int i = 0; i < NLOOPS; ++i)
     {
         sk = scalar(x, y);
-        sk = scalar_trans(x, y);
-        //sk = norm(x);
+        sk = norm(x);
     }
 
     double t1 = omp_get_wtime() - tstart;             // OpenMP
@@ -124,7 +123,7 @@ int main(int argc, char const *argv[])
     cout << "done\n";
     cout << vr << endl;
 
-    N=2;
+    N=200;
     //Data (re-)inizialiion
     for (unsigned int i = 0; i < N; ++i)
     {
@@ -147,9 +146,9 @@ int main(int argc, char const *argv[])
 
          omp_set_schedule(omp_sched_static, 0);
         tstart = omp_get_wtime();
-        benchmark(x, y, N, NLOOPS);
-        t1 = omp_get_wtime()/NLOOPS;
-        cout << "static (chunk 0) "<< (t1-tstart) << endl;
+        benchmark(x, y, NLOOPS);
+        t1 = (omp_get_wtime()-tstart)/NLOOPS;
+        cout << "static (chunk 0) "<< t1 << endl;
         for(int i=0; i<= 5; i++)
         {
             
@@ -159,30 +158,30 @@ int main(int argc, char const *argv[])
             // STATIC
             omp_set_schedule(omp_sched_static, chunk);
             tstart = omp_get_wtime();
-            benchmark(x, y, N, NLOOPS);
-            t1 = omp_get_wtime()/NLOOPS;
-            std::cout << "static:  " << (t1 - tstart) << " s\n";
+            benchmark(x, y, NLOOPS);
+            t1 = (omp_get_wtime()-tstart)/NLOOPS;
+            std::cout << "static:  "<< t1 << " s\n";
 
             // DYNAMIC
             omp_set_schedule(omp_sched_dynamic, chunk);
             tstart = omp_get_wtime();
-            benchmark(x, y, N, NLOOPS);
-            t1 = omp_get_wtime()/NLOOPS;
-            std::cout << "dynamic: " << (t1 - tstart) << " s\n";
+            benchmark(x, y, NLOOPS);
+            t1 = (omp_get_wtime()-tstart)/NLOOPS;
+            std::cout << "dynamic: "<< t1 << " s\n";
 
             // GUIDED
             omp_set_schedule(omp_sched_guided, chunk);
             tstart = omp_get_wtime();
-            benchmark(x, y, N, NLOOPS);
-            t1 = omp_get_wtime()/NLOOPS;
-            std::cout << "guided:  " << (t1 - tstart) << " s\n";
+            benchmark(x, y, NLOOPS);
+            t1 = (omp_get_wtime()-tstart)/NLOOPS;
+            std::cout << "guided:  "<< t1 << " s\n";
 
             // AUTO
             omp_set_schedule(omp_sched_auto, chunk);
             tstart = omp_get_wtime();
-            benchmark(x, y, N, NLOOPS);
-            t1 = omp_get_wtime()/NLOOPS;
-            std::cout << "auto:    " << (t1 - tstart) << " s\n";
+            benchmark(x, y, NLOOPS);
+            t1 = (omp_get_wtime()-tstart)/NLOOPS;
+            std::cout << "auto:    "<< t1 << " s\n";
             cout << endl;
         }
         cout << endl;
@@ -194,7 +193,9 @@ int main(int argc, char const *argv[])
      
     
     vector<int> vec = reduction_vec_append(N);
-    for(int i=0; i< N; i++)
+    unsigned int n = vec.size();
+    cout << "vec.size() = " << n << "\n";
+    for(int i=0; i< n; i++)
     {
         cout << vec[i] << ", ";
     }

sheet5/1/mylib.cpp

@@ -91,7 +91,7 @@ vector<int> reduction_vec(int n)
 
 vector<int> reduction_vec_append(int n)
 {
-    vector<int> vec(n);
+    vector<int> vec;
     #pragma omp parallel default(none) shared(cout,n) reduction(VecAppend:vec)
     {
         int tid = omp_get_thread_num();
@@ -99,10 +99,7 @@ vector<int> reduction_vec_append(int n)
         vector<int> local(n);
         iota(local.begin(), local.end(), tid);
 
-        #pragma omp critical
-        cout << tid << " : " << local.size() << endl;
-
-        vec = local;  
+        vec = local;
     }
 
     return vec;

sheet5/3/mayer_primes.h:Zone.Identifier

@@ -1,3 +0,0 @@
-[ZoneTransfer]
-ZoneId=3
-HostUrl=https://imsc.uni-graz.at/haasegu/Lectures/Math2CPP/Examples/goldbach/mayer_primes.h

sheet5/4/main.cpp

@@ -183,50 +183,63 @@ int main(int argc, char **argv)
 
 
 
+const int NLOOPS = 20;   
+
 for (int k = 3; k <= 8; ++k)
+{
+    unsigned int n = pow(10.0, k);
+
+    vector<double> x(n), y(n);
+    for (unsigned int  i = 0; i < n; ++i)
     {
-        size_t n = (size_t)pow(10, k);
-
-        vector<double> x(n), y(n);
-        for (unsigned int i = 0; i < n; ++i)
-        {
-            double xi= (i % 219) + 1;
-            x[i] = xi;
-            y[i] = 1.0 / xi;
-        }
-
-        // ---- SUM benchmark (sequential) ----
-        double t0 = omp_get_wtime();
-        double s1 = benchmark_A_sum_old(x);
-        double t_sum_seq = omp_get_wtime() - t0;
-
-        // ---- SUM benchmark (parallel) ----
-        t0 = omp_get_wtime();
-        double s2 = benchmark_A_sum(x);
-        double t_sum_omp = omp_get_wtime() - t0;
-
-        double sum_speedup = t_sum_seq / t_sum_omp;
-
-
-        // ---- INNER PRODUCT benchmark (sequential) ----
-        t0 = omp_get_wtime();
-        double ip1 = benchmark_A_old(x, y);
-        double t_inner_seq = omp_get_wtime() - t0;
-
-        // ---- INNER PRODUCT benchmark (parallel) ----
-        t0 = omp_get_wtime();
-        double ip2 = benchmark_A(x, y);
-        double t_inner_omp = omp_get_wtime() - t0;
-
-        double inner_speedup = t_inner_seq / t_inner_omp;
-
-        // ---- Print results ----
-        cout << k << endl;
-            cout << t_sum_seq << ", " << t_sum_omp << ", " << sum_speedup << endl;
-           cout << t_inner_seq << ", " << t_inner_omp << ", " << inner_speedup << endl;
-        cout << endl;
+        double xi = (i % 219) + 1;
+        x[i] = xi;
+        y[i] = 1.0 / xi;
     }
 
+    double s1_guard = 0.0, s2_guard = 0.0;
+    double ip1_guard = 0.0, ip2_guard = 0.0;
+
+    // ---- SUM benchmark (sequential) ----
+    double t0 = omp_get_wtime();
+    for (int r = 0; r < NLOOPS; ++r)
+        s1_guard += benchmark_A_sum_old(x);
+    double t_sum_seq = (omp_get_wtime() - t0) / NLOOPS;
+
+    // ---- SUM benchmark (parallel) ----
+    t0 = omp_get_wtime();
+    for (int r = 0; r < NLOOPS; ++r)
+        s2_guard += benchmark_A_sum(x);
+    double t_sum_omp = (omp_get_wtime() - t0) / NLOOPS;
+
+    double sum_speedup = t_sum_seq / t_sum_omp;
+
+    // ---- INNER PRODUCT benchmark (sequential) ----
+    t0 = omp_get_wtime();
+    for (int r = 0; r < NLOOPS; ++r)
+        ip1_guard += benchmark_A_old(x, y);
+    double t_inner_seq = (omp_get_wtime() - t0) / NLOOPS;
+
+    // ---- INNER PRODUCT benchmark (parallel) ----
+    t0 = omp_get_wtime();
+    for (int r = 0; r < NLOOPS; ++r)
+        ip2_guard += benchmark_A(x, y);
+    double t_inner_omp = (omp_get_wtime() - t0) / NLOOPS;
+
+    double inner_speedup = t_inner_seq / t_inner_omp;
+
+    // ---- Print results ----
+    std::cout << "k = " << k << " (n = 10^" << k << " = " << n << ")\n";
+    std::cout << "SUM    seq: " << t_sum_seq   << " s,  omp: " << t_sum_omp
+              << " s,  speedup = " << sum_speedup   << '\n';
+    std::cout << "INNER  seq: " << t_inner_seq << " s,  omp: " << t_inner_omp
+              << " s,  speedup = " << inner_speedup << '\n';
+    std::cout << "guards: "
+              << s1_guard << ", " << s2_guard << ", "
+              << ip1_guard << ", " << ip2_guard << "\n\n";
+}
+
+
 
 
 return 0;

sheet5/5/getmatrix.cpp

@@ -8,6 +8,7 @@
 #include <iostream>
 #include <list>
 #include <vector>
+#include <omp.h>
 using namespace std;
 
 
@@ -168,33 +169,64 @@ void CRS_Matrix::Debug() const
 
 void CRS_Matrix::CalculateLaplace(vector<double> &f)
 {
+    double t_start = omp_get_wtime();
+
+    double t0 ,t1,t2;
     assert(_mesh.NdofsElement() == 3);               // only for triangular, linear elements
     //cout << _nnz << " vs. " << _id[_nrows] << "  " << _nrows<< endl;
     assert(_nnz == _id[_nrows]);
 
-    for (int k = 0; k < _nrows; ++k)
+    #pragma omp parallel
     {
-        _sk[k] = 0.0;
-    }
-    for (int k = 0; k < _nrows; ++k)
-    {
-        f[k] = 0.0;
-    }
+        #pragma omp for
+        for (int k = 0; k < _nrows; ++k)
+            {
+                _sk[k] = 0.0;
+            }
 
-    double ske[3][3], fe[3];
-    //  Loop over all elements
-    auto const nelem = _mesh.Nelems();
-    auto const &ia   = _mesh.GetConnectivity();
-    auto const &xc   = _mesh.GetCoords();
+        #pragma omp barrier
+        #pragma omp single
+        t0 = omp_get_wtime();
+        
+
+        #pragma omp for   
+        for (int k = 0; k < _nrows; ++k)
+        {
+            f[k] = 0.0;
+        }
+
+        #pragma omp barrier
+        #pragma omp single
+        t1 = omp_get_wtime();
+       
+
+        double ske[3][3], fe[3];
+        //  Loop over all elements
+        auto const nelem = _mesh.Nelems();
+        auto const &ia   = _mesh.GetConnectivity();
+        auto const &xc   = _mesh.GetCoords();
+
+        #pragma omp barrier
+        #pragma omp single
+        t2 = omp_get_wtime();
+        
+        #pragma omp for
+        for (int i = 0; i < nelem; ++i)
+        {
+            CalcElem(ia.data() + 3 * i, xc.data(), ske, fe);
+            AddElem_3(ia.data() + 3 * i, ske, fe, f);
+        }
+        //Debug();
 
-    for (int i = 0; i < nelem; ++i)
-    {
-        CalcElem(ia.data() + 3 * i, xc.data(), ske, fe);
-        AddElem_3(ia.data() + 3 * i, ske, fe, f);
     }
-    //Debug();
+    double t3 = omp_get_wtime();
+    cout << "Zero matrix: " << (t0 - t_start) << " sec\n";
+    cout << "Zero RHS: " << (t1 - t0)     << " sec\n";
+    cout << "Element loop: " << (t3 - t2)    << " sec\n";
+    cout << "Total assembly:" << (t3 - t_start) << " sec\n";
+   
 
-    return;
+    
 }
 
 void CRS_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f)
@@ -286,7 +318,8 @@ void CRS_Matrix::Defect(vector<double> &w,
 {
     assert( _nrows == static_cast<int>(w.size()) );
     assert( w.size() == u.size() && u.size() == f.size() );
-
+    
+    #pragma omp parallel for
     for (int row = 0; row < _nrows; ++row)
     {
         double wi = f[row];
@@ -340,8 +373,10 @@ void CRS_Matrix::AddElem_3(int const ial[3], double const ske[3][3], double cons
                 assert(ip >= 0);
             }
 #endif
+            #pragma omp atomic
             _sk[ip] += ske[i][j];
         }
+        #pragma omp atomic
         f[ii] += fe[i];
     }
 }

sheet5/5/jacsolve.cpp

@@ -23,6 +23,7 @@ void JacobiSolve(CRS_Matrix const &SK, vector<double> const &f, vector<double> &
 
     cout << endl << " Start Jacobi solver for " << nrows << " d.o.f.s"  << endl;
     //  Choose initial guess
+    #pragma omp parallel for
     for (int k = 0; k < nrows; ++k)
     {
         u[k] = 0.0;                          //  u := 0
@@ -44,6 +45,7 @@ void JacobiSolve(CRS_Matrix const &SK, vector<double> const &f, vector<double> &
     // Iteration sweeps
     int iter  = 0;
     double sigma = sigma0;
+    // i dont parallelize this cause its iteration dependent
     while ( sigma > tol2 * sigma0 && maxiter > iter)
     {
         ++iter;

sheet5/5/main.cpp

@@ -11,6 +11,7 @@
 #include <chrono>             // timing
 #include <cmath>
 #include <iostream>
+#include <omp.h>
 using namespace std;
 using namespace std::chrono;  // timing
 
@@ -19,7 +20,6 @@ int main(int, char ** )
 {
     const int numprocs = 1;
     const int myrank   = 0;
-
     if (myrank == 0)
     {
         cout << "\n There are " << numprocs << " processes running.\n \n";
@@ -121,6 +121,7 @@ int main(int, char ** )
 
             //mesh.Write_ascii_matlab("uv.txt", uv);
             //mesh.Visualize(uv);
+            
         }
         return 0;
     }

sheet5/5/vdop.cpp

@@ -11,7 +11,7 @@ void vddiv(vector<double> &x, vector<double> const &y,
 {
     assert( x.size() == y.size() && y.size() == z.size() );
     size_t n = x.size();
-
+    #pragma omp parallel for
     for (size_t k = 0; k < n; ++k)
     {
         x[k] = y[k] / z[k];
@@ -26,7 +26,7 @@ void vdaxpy(std::vector<double> &x, std::vector<double> const &y,
 {
     assert( x.size() == y.size() && y.size() == z.size() );
     size_t n = x.size();
-
+    #pragma omp parallel for
     for (size_t k = 0; k < n; ++k)
     {
         x[k] = y[k] + alpha * z[k];
@@ -41,6 +41,7 @@ double dscapr(std::vector<double> const &x, std::vector<double> const &y)
     size_t n = x.size();
 
     double    s = 0.0;
+     #pragma omp parallel for reduction(+:s)
     for (size_t k = 0; k < n; ++k)
     {
         s += x[k] * y[k];

sheet5/5_star/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\
+	getmatrix.cpp jacsolve.cpp userset.cpp
+# 	dexx.cpp debugd.cpp skalar.cpp  vecaccu.cpp accudiag.cpp
+
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PROGRAM	= ${MAIN}.${COMPILER}
+
+# uncomment the next to lines for debugging and detailed performance analysis
+CXXFLAGS += -g
+# -pg slows down the code on my laptop when using CLANG_
+#LINKFLAGS += -pg
+#CXXFLAGS += -Q --help=optimizers
+#CXXFLAGS += -fopt-info
+
+include ../${COMPILER}default.mk
+
+#############################################################################
+# additional specific cleaning in this directory
+clean_all::
+	@rm -f t.dat*
+
+
+#############################################################################
+# 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

sheet5/5_star/ToDo.txt

@@ -0,0 +1,5 @@
+// Jan 15, 2019
+
+geom.h:75        void SetValues(std::vector<double> &v) const;  // GH: TODO with functor
+Set vector values using a functor ff(x,y).
+ See solution in Progs/cds

sheet5/5_star/ascii_read_meshvector.m

@@ -0,0 +1,43 @@
+function [ xc, ia, v ] = ascii_read_meshvector( fname )
+%
+% Loads the 2D triangular mesh (coordinates, vertex connectivity) 
+%  together with values on its vertices from an ASCII file.
+%   Matlab indexing is stored  (starts with 1).
+% 
+% The input file format is compatible 
+%  with Mesh_2d_3_matlab:Write_ascii_matlab(..) in jacobi_oo_stl/geom.h
+%
+%
+%  IN: fname - filename
+% OUT: xc    - coordinates
+%      ia    - mesh connectivity
+%      v     - solution vector
+
+DELIMETER = ' ';
+
+fprintf('Read file  %s\n',fname)
+
+% Read mesh constants
+nn = dlmread(fname,DELIMETER,[0 0 0 3]);  %% row_1, col_1, row_2, col_2  in C indexing!!!
+nnode = nn(1);
+ndim  = nn(2);
+nelem = nn(3);
+nvert = nn(4);
+
+% Read coordinates
+row_start = 0+1;
+row_end   = 0+nnode;
+xc = dlmread(fname,DELIMETER,[row_start 0 row_end ndim-1]);
+
+% Read connectivity
+row_start = row_end+1;
+row_end   = row_end+nelem;
+ia = dlmread(fname,DELIMETER,[row_start 0 row_end nvert-1]);
+
+% Read solution
+row_start = row_end+1;
+row_end   = row_end+nnode;
+v = dlmread(fname,DELIMETER,[row_start 0 row_end 0]);
+end
+
+

sheet5/5_star/ascii_write_mesh.m

@@ -0,0 +1,49 @@
+function ascii_write_mesh( xc, ia, e, basename)
+%
+% Saves the 2D triangular mesh in the minimal way (only coordinates, vertex connectivity, minimal boundary edge info)
+%  in an ASCII file.
+%  Matlab indexing is stored  (starts with 1).
+% 
+% The output file format is compatible with Mesh_2d_3_matlab:Mesh_2d_3_matlab(std::string const &fname) in jacobi_oo_stl/geom.h
+%
+% IN:
+% coordinates  xc: [2][nnode]
+% connectivity ia: [4][nelem]   with  t(4,:) are the subdomain numbers
+% edges         e: [7][nedges]  boundary edges
+%                              e([1,2],:) - start/end vertex of edge
+%                              e([3,4],:) - start/end values
+%                              e(5,:)     - segment number
+%                              e([6,7],:) - left/right subdomain
+%        basename: file name without extension
+% 
+% Data have been generated via <https://de.mathworks.com/help/pde/ug/initmesh.html initmesh>.
+%
+fname = [basename, '.txt'];
+
+nnode = int32(size(xc,2));
+ndim  = int32(size(xc,1));
+nelem = int32(size(ia,2));
+nvert_e = int32(3);
+
+
+dlmwrite(fname,nnode,'delimiter','\t','precision',16)                 % number of nodes
+dlmwrite(fname,ndim,'-append','delimiter','\t','precision',16)        % space dimension
+dlmwrite(fname,nelem,'-append','delimiter','\t','precision',16)       % number of elements
+dlmwrite(fname,nvert_e,'-append','delimiter','\t','precision',16)     % number of vertices per element
+
+% dlmwrite(fname,xc(:),'-append','delimiter','\t','precision',16)       % coordinates
+dlmwrite(fname,xc([1,2],:).','-append','delimiter','\t','precision',16) % coordinates
+
+% no subdomain info transferred
+tmp=int32(ia(1:3,:));
+% dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16)      % connectivity in Matlab indexing
+dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16)    % connectivity in Matlab indexing
+
+% store only start and end point of boundary edges,
+nbedges = size(e,2);
+dlmwrite(fname,nbedges,'-append','delimiter','\t','precision',16)     % number boundary edges
+tmp=int32(e(1:2,:));
+% dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16)    % boundary edges in Matlab indexing
+dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16)  % boundary edges in Matlab indexing
+
+end

sheet5/5_star/geom.cpp

@@ -0,0 +1,522 @@
+// see:   http://llvm.org/docs/CodingStandards.html#include-style
+#include "geom.h"
+
+#include <algorithm>
+#include <cassert>
+#include <fstream>
+#include <iostream>
+#include <list>
+#include <string>
+#include <vector>
+using namespace std;
+
+Mesh::Mesh(int ndim, int nvert_e, int ndof_e)
+    : _nelem(0), _nvert_e(nvert_e), _ndof_e(ndof_e), _nnode(0), _ndim(ndim), _ia(0), _xc(0)
+{
+}
+
+Mesh::~Mesh()
+{}
+
+void Mesh::SetValues(std::vector<double> &v, const std::function<double(double, double)> &func) const
+{
+    int const nnode = Nnodes();            // number of vertices in mesh
+    assert( nnode == static_cast<int>(v.size()) );
+    for (int k = 0; k < nnode; ++k)
+    {
+        v[k] = func( _xc[2 * k], _xc[2 * k + 1] );
+    }
+}
+
+
+void Mesh::Debug() const
+{
+    cout << "\n ############### Debug  M E S H  ###################\n";
+    cout << "\n ...............    Coordinates       ...................\n";
+    for (int k = 0; k < _nnode; ++k)
+    {
+        cout << k << " : " ;
+        for (int i = 0; i < _ndof_e; ++i )
+        {
+            cout << _xc[2*k+i] << "  ";
+        }
+        cout << endl;
+    }
+    cout << "\n ...............    Elements        ...................\n";
+    for (int k = 0; k < _nelem; ++k)
+    {
+        cout << k << " : ";
+        for (int i = 0; i < _ndof_e; ++i )
+            cout << _ia[_ndof_e * k + i] << "  ";
+        cout << endl;
+    }
+    return;
+}
+
+void Mesh::Write_ascii_matlab(std::string const &fname, std::vector<double> const &v) const
+{
+    assert(Nnodes() ==  static_cast<int>(v.size()));  // fits vector length to mesh information?
+
+    ofstream fout(fname);                             // open file ASCII mode
+    if ( !fout.is_open() )
+    {
+        cout << "\nFile " << fname << " has not been opened.\n\n" ;
+        assert( fout.is_open() && "File not opened."  );
+    }
+
+    string const DELIMETER(" ");    // define the same delimeter as in matlab/ascii_read*.m
+    int const    OFFSET(1);         // convert C-indexing to matlab
+
+    // Write data: #nodes, #space dimensions, #elements, #vertices per element
+    fout << Nnodes() << DELIMETER << Ndims() << DELIMETER << Nelems() << DELIMETER << NverticesElements() << endl;
+
+    // Write cordinates: x_k, y_k   in seperate lines
+    assert( Nnodes()*Ndims() ==  static_cast<int>(_xc.size()));
+    for (int k = 0, kj = 0; k < Nnodes(); ++k)
+    {
+        for (int j = 0; j < Ndims(); ++j, ++kj)
+        {
+            fout << _xc[kj] << DELIMETER;
+        }
+        fout << endl;
+    }
+
+    // Write connectivity: ia_k,0, ia_k,1 etc  in seperate lines
+    assert( Nelems()*NverticesElements() ==  static_cast<int>(_ia.size()));
+    for (int k = 0, kj = 0; k < Nelems(); ++k)
+    {
+        for (int j = 0; j < NverticesElements(); ++j, ++kj)
+        {
+            fout << _ia[kj] + OFFSET << DELIMETER;     // C to matlab
+        }
+        fout << endl;
+    }
+
+    // Write vector
+    for (int k = 0; k < Nnodes(); ++k)
+    {
+        fout << v[k] << endl;
+    }
+
+    fout.close();
+    return;
+}
+
+
+void Mesh::Visualize(std::vector<double> const &v) const
+{
+    // define external command
+    const string exec_m("matlab -nosplash < visualize_results.m");                 // Matlab
+    //const string exec_m("octave --no-window-system --no-gui visualize_results.m"); // Octave
+    //const string exec_m("flatpak run org.octave.Octave visualize_results.m");      // Octave (flatpak): desktop GH
+
+    const string fname("uv.txt");
+    Write_ascii_matlab(fname, v);
+
+    int ierror = system(exec_m.c_str());                                 // call external command
+
+    if (ierror != 0)
+    {
+        cout << endl << "Check path to Matlab/octave on your system" << endl;
+    }
+    cout << endl;
+    return;
+}
+
+
+
+
+// #####################################################################
+Mesh_2d_3_square::Mesh_2d_3_square(int nx, int ny, int myid, int procx, int procy)
+    : Mesh(2, 3, 3), // two dimensions, 3 vertices, 3 dofs
+      _myid(myid), _procx(procx), _procy(procy), _neigh{{-1, -1, -1, -1}}, _color(0),
+_xl(0.0), _xr(1.0), _yb(0.0), _yt(1.0), _nx(nx), _ny(ny)
+{
+    //void IniGeom(int const myid, int const procx, int const procy, int neigh[], int &color)
+    int const ky = _myid / _procx;
+    int const kx = _myid % _procy;  //    MOD(myid,procx)
+    // Determine the neighbors of domain/rank myid
+    _neigh[0] = (ky == 0)       ?  -1 : _myid - _procx;    //   South
+    _neigh[1] = (kx == _procx - 1) ?  -1 : _myid + 1;      //   East
+    _neigh[2] = (ky == _procy - 1) ?  -1 : _myid + _procx;  //   North
+    _neigh[3] = (kx == 0)       ?  -1 : _myid - 1;        //   West
+
+    _color = (kx + ky) & 1 ;
+
+    // subdomain is part of unit square
+    double const hx = 1. / _procx;
+    double const hy = 1. / _procy;
+    _xl = kx * hx;                      //  left
+    _xr = (kx + 1) * hx;                //  right
+    _yb = ky * hy;                      //  bottom
+    _yt = (ky + 1) * hy;                //  top
+
+    // Calculate coordinates
+    int const nnode = (_nx + 1) * (_ny + 1); // number of nodes
+    Resize_Coords(nnode, 2);                 // coordinates in 2D [nnode][ndim]
+    GetCoordsInRectangle(_nx, _ny, _xl, _xr, _yb, _yt, GetCoords().data());
+
+    // Calculate element connectivity (linear triangles)
+    int const nelem = 2 * _nx * _ny;         // number of elements
+    Resize_Connectivity(nelem, 3);           // connectivity matrix [nelem][3]
+    GetConnectivityInRectangle(_nx, _ny, GetConnectivity().data());
+
+    return;
+}
+
+
+void Mesh_2d_3_square::SetU(std::vector<double> &u) const
+{
+    int dx    = _nx + 1;
+    for (int j = 0; j <= _ny; ++j)
+    {
+        int k = j * dx;
+        for (int i = 0; i <= _nx; ++i, ++k)
+        {
+            u[k] = 0.0;
+        }
+    }
+
+}
+
+void Mesh_2d_3_square::SetF(std::vector<double> &f) const
+{
+    int dx    = _nx + 1;
+    for (int j = 0; j <= _ny; ++j)
+    {
+        int k = j * dx;
+        for (int i = 0; i <= _nx; ++i, ++k)
+        {
+            f[k] = 1.0;
+        }
+    }
+
+}
+
+
+std::vector<int> Mesh_2d_3_square::Index_DirichletNodes() const
+{
+    int const dx = 1,
+              dy = _nx + 1,
+              bl  = 0,
+              br  = _nx,
+              tl  = _ny * (_nx + 1),
+              tr  = (_ny + 1) * (_nx + 1) - 1;
+    int const start[4] = { bl, br, tl, bl},
+                         end[4] = { br, tr, tr, tl},
+                                  step[4] = { dx, dy, dx, dy};
+
+    vector<int> idx(0);
+    for (int j = 0; j < 4; j++)
+    {
+        if (_neigh[j] < 0)
+        {
+            for (int i = start[j]; i <= end[j]; i += step[j])
+            {
+                idx.push_back(i);        // node i is Dirichlet node
+            }
+        }
+    }
+    //    remove multiple elements
+    sort(idx.begin(), idx.end());                          // sort
+    idx.erase( unique(idx.begin(), idx.end()), idx.end() ); // remove duplicate data
+
+    return idx;
+}
+
+void Mesh_2d_3_square::SaveVectorP(std::string const &name, vector<double> const &u) const
+{
+//  construct the file name for subdomain myid
+    const string tmp( std::to_string(_myid / 100) + to_string((_myid % 100) / 10) + to_string(_myid % 10) );
+
+    const string namep = name + "." + tmp;
+    ofstream ff(namep.c_str());
+    ff.precision(6);
+    ff.setf(ios::fixed, ios::floatfield);
+
+    // assumes tensor product grid in unit square; rowise numbered (as generated in class constructor)
+    // output is provided for tensor product grid visualization ( similar to Matlab-surf() )
+    auto const &xc = GetCoords();
+    int k = 0;
+    for (int j = 0; j <= _ny; ++j)
+    {
+        for (int i = 0; i <= _nx; ++i, ++k)
+            ff << xc[2 * k + 0] << "   " << xc[2 * k + 1] << "   " << u[k] << endl;
+        ff << endl;
+    }
+
+    ff.close();
+    return;
+}
+
+void Mesh_2d_3_square::GetCoordsInRectangle(int const nx, int const ny,
+        double const xl, double const xr, double const yb, double const yt,
+        double xc[])
+{
+    const double hx = (xr - xl) / nx,
+                 hy = (yt - yb) / ny;
+
+    int k = 0;
+    for (int j = 0; j <= ny; ++j)
+    {
+        const double y0 = yb + j * hy;
+        for (int i = 0; i <= nx; ++i, k += 2)
+        {
+            xc[k  ] = xl + i * hx;
+            xc[k + 1] = y0;
+        }
+    }
+
+    return;
+}
+
+void Mesh_2d_3_square::GetConnectivityInRectangle(int const nx, int const ny, int ia[])
+{
+    const int dx = nx + 1;
+    int k  = 0;
+    int l  = 0;
+    for (int j = 0; j < ny; ++j, ++k)
+    {
+        for (int i = 0; i < nx; ++i, ++k)
+        {
+            ia[l  ] = k;
+            ia[l + 1] = k + 1;
+            ia[l + 2] = k + dx + 1;
+            l += 3;
+            ia[l  ] = k;
+            ia[l + 1] = k + dx;
+            ia[l + 2] = k + dx + 1;
+            l += 3;
+        }
+    }
+    return;
+}
+
+// #################### still some old code (--> MPI) ############################
+
+
+//  Copies the values of w corresponding to the boundary
+//  South (ib==1), East (ib==2), North (ib==3), West (ib==4)
+
+void GetBound(int const ib, int const nx, int const ny, double const w[], double s[])
+{
+    const int //dx = 1,
+    dy = nx + 1,
+    bl  = 0,
+    br  = nx,
+    tl  = ny * (nx + 1),
+    tr  = (ny + 1) * (nx + 1) - 1;
+    switch (ib)
+    {
+        case 1:
+        {
+            for (int i = bl, j = 0; i <= br; ++i, ++j)
+                s[j] = w[i];
+            break;
+        }
+        case 3:
+        {
+            for (int i = tl, j = 0; i <= tr; ++i, ++j)
+                s[j] = w[i];
+            break;
+        }
+        case 4:
+        {
+            for (int i = bl, j = 0; i <= tl; i += dy, ++j)
+                s[j] = w[i];
+            break;
+        }
+        case 2:
+        {
+            for (int i = br, j = 0; i <= tr; i += dy, ++j)
+                s[j] = w[i];
+            break;
+        }
+        default:
+        {
+            cout << endl << "Wrong parameter  ib in " << __FILE__ << ":" << __LINE__ << endl;
+        }
+    }
+    return;
+}
+
+// ----------------------------------------------------------------------------------------------------------
+// Computes w:  = w + s at nodes on  the boundary
+// South (ib == 1), East (ib == 2), North (ib == 3), West (ib == 4)
+
+void AddBound(int const ib, int const nx, int const ny, double w[], double const s[])
+{
+    int const dy = nx + 1,
+              bl  = 0,
+              br  = nx,
+              tl  = ny * (nx + 1),
+              tr  = (ny + 1) * (nx + 1) - 1;
+    switch (ib)
+    {
+        case 1:
+        {
+            for (int i = bl, j = 0; i <= br; ++i, ++j)
+                w[i] += s[j];
+            break;
+        }
+        case 3:
+        {
+            for (int i = tl, j = 0; i <= tr; ++i, ++j)
+                w[i] += s[j];
+            break;
+        }
+        case 4:
+        {
+            for (int i = bl, j = 0; i <= tl; i += dy, ++j)
+                w[i] += s[j];
+            break;
+        }
+        case 2:
+        {
+            for (int i = br, j = 0; i <= tr; i += dy, ++j)
+                w[i] += s[j];
+            break;
+        }
+        default:
+        {
+            cout << endl << "Wrong parameter  ib in " << __FILE__ << ":" << __LINE__ << endl;
+        }
+    }
+    return;
+}
+
+
+// ####################################################################
+
+Mesh_2d_3_matlab::Mesh_2d_3_matlab(string const &fname)
+    : Mesh(2, 3, 3),  // two dimensions, 3 vertices, 3 dofs
+      bedges(0)
+{
+    ifstream ifs(fname);
+    if (!(ifs.is_open() && ifs.good()))
+    {
+        cerr << "Mesh_2d_3_matlab: Error cannot open file " << fname << endl;
+        assert(ifs.is_open());
+    }
+
+    int const OFFSET(1);             // Matlab to C indexing
+    cout << "ASCI file  " << fname << "  opened" << endl;
+
+    // Read some mesh constants
+    int nnode, ndim, nelem, nvert_e;
+    ifs >> nnode >> ndim >> nelem >> nvert_e;
+    cout << nnode << "  " << ndim << "  " << nelem << "  " << nvert_e << endl;
+    assert(ndim == 2 && nvert_e == 3);
+
+    // Allocate memory
+    Resize_Coords(nnode, ndim);                 // coordinates in 2D [nnode][ndim]
+    Resize_Connectivity(nelem, nvert_e);        // connectivity matrix [nelem][nvert]
+
+    // Read ccordinates
+    auto &xc = GetCoords();
+    for (int k = 0; k < nnode * ndim; ++k)
+    {
+        ifs >> xc[k];
+    }
+
+    // Read connectivity
+    auto &ia = GetConnectivity();
+    for (int k = 0; k < nelem * nvert_e; ++k)
+    {
+        ifs >> ia[k];
+        ia[k] -= OFFSET;                // Matlab to C indexing
+    }
+
+    // additional read of boundary information (only start/end point)
+    int nbedges;
+    ifs >> nbedges;
+
+    bedges.resize(nbedges * 2);
+    for (int k = 0; k < nbedges * 2; ++k)
+    {
+        ifs >> bedges[k];
+        bedges[k] -= OFFSET;            // Matlab to C indexing
+    }
+
+    return;
+}
+
+// binary
+//{
+//ifstream ifs(fname, ios_base::in | ios_base::binary);
+//if(!(ifs.is_open() && ifs.good()))
+//{
+//cerr << "ReadBinMatrix: Error cannot open file " << file << endl;
+//assert(ifs.is_open());
+//}
+//cout << "ReadBinMatrix: file opened" << file << endl;
+
+
+//}
+
+// binaryIO.cpp
+//void read_binMatrix(const string& file, vector<int> &cnt, vector<int> &col, vector<double> &ele)
+//{
+
+//ifstream ifs(file, ios_base::in | ios_base::binary);
+
+//if(!(ifs.is_open() && ifs.good()))
+//{
+//cerr << "ReadBinMatrix: Error cannot open file " << file << endl;
+//assert(ifs.is_open());
+//}
+//cout << "ReadBinMatrix: Opened file " << file << endl;
+
+//int _size;
+
+//ifs.read(reinterpret_cast<char*>(&_size), sizeof(int));   // old: ifs.read((char*)&_size, sizeof(int));
+//cnt.resize(_size);
+//cout << "ReadBinMatrix: cnt size: " << _size << endl;
+
+//ifs.read((char*)&_size, sizeof(int));
+//col.resize(_size);
+//cout << "ReadBinMatrix: col size: " << _size << endl;
+
+//ifs.read((char*)&_size, sizeof(int));
+//ele.resize(_size);
+//cout << "ReadBinMatrix: ele size: " << _size << endl;
+
+
+//ifs.read((char*)cnt.data(), cnt.size() * sizeof(int));
+//ifs.read((char*)col.data(), col.size() * sizeof(int));
+//ifs.read((char*)ele.data(), ele.size() * sizeof(double));
+
+//ifs.close();
+//cout << "ReadBinMatrix: Finished reading matrix.." << endl;
+
+//}
+
+
+std::vector<int> Mesh_2d_3_matlab::Index_DirichletNodes() const
+{
+    vector<int> idx(bedges);                              // copy
+
+    sort(idx.begin(), idx.end());                         // sort
+    idx.erase( unique(idx.begin(), idx.end()), idx.end() ); // remove duplicate data
+
+    return idx;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

sheet5/5_star/geom.h

@@ -0,0 +1,381 @@
+#ifndef GEOM_FILE
+#define GEOM_FILE
+#include <array>
+#include <functional>             // function; C++11
+#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)
+     */
+    explicit Mesh(int ndim, int nvert_e = 0, int ndof_e = 0);
+
+    /**
+     * 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();
+
+    /**
+     * 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 convectivity vector [nelems*ndofs].
+     */
+    const std::vector<int>  &GetConnectivity() const
+    {
+        return _ia;
+    }
+
+    /**
+     * Access/Change connectivity information (g1,g2,g3)_i.
+     * @return convectivity 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;
+
+    /**
+     * Prints the information for a finite element mesh
+     */
+    void Debug() const;
+
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions
+     * @return index vector.
+     */
+    virtual std::vector<int> Index_DirichletNodes() const = 0;
+
+    /**
+     * Write vector @p v toghether 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;
+
+    /**
+     * 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;
+
+
+protected:
+    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;
+    }
+
+private:
+    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
+};
+
+/**
+ * 2D finite element mesh of the square consiting 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;
+
+    /**
+      * 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 dicretization 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] 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 coodinates 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
+};
+
+// #################### still some old code (--> MPI) ############################
+/**
+ * Copies the values of @p w corresponding to boundary @p ib
+ * onto vector s.  South (ib==1), East (ib==2), North (ib==3), West (ib==4).
+ * The vector @p s has to be long enough!!
+ * @param[in] ib    my local boundary
+ * @param[in] nx    number of discretization intervals in x-direction
+ * @param[in] ny    number of discretization intervals in y-direction
+ * @param[in] w     vector for all nodes of local discretization
+ * @param[out] s    short vector with values on boundary @p ib
+*/
+// GH_NOTE: Absicherung bei s !!
+void GetBound(int ib, int nx, int ny, double const w[], double s[]);
+
+
+/**
+ * Computes @p w := @p w + @p s  at the interface/boundary nodes on the
+ * boundary @p ib .  South (ib==1), East (ib==2), North (ib==3), West (ib==4)
+ * @param[in] ib    my local boundary
+ * @param[in] nx    number of discretization intervals in x-direction
+ * @param[in] ny    number of discretization intervals in y-direction
+ * @param[in,out] w vector for all nodes of local discretization
+ * @param[in] s     short vector with values on boundary @p ib
+*/
+void AddBound(int ib, int nx, int ny, double w[], double const s[]);
+
+// #################### Mesh from Matlab ############################
+/**
+ * 2D finite element mesh of the square consiting of linear triangular elements.
+ */
+class Mesh_2d_3_matlab: public Mesh
+{
+public:
+    /**
+     * Reads mesh data from a binary file.
+     *
+     * File format, see ascii_write_mesh.m
+     *
+     * @param[in] fname file name
+    */
+    explicit Mesh_2d_3_matlab(std::string const &fname);
+
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions.
+     * @return index vector.
+      *
+      * @warning All boundary nodes are considered as Dirchlet nodes.
+     */
+    std::vector<int> Index_DirichletNodes() const override;
+
+private:
+    /**
+     * Determines the indices of those vertices with Dirichlet boundary conditions
+     * @return index vector.
+     */
+    int Nnbedges() const
+    {
+        return static_cast<int>(bedges.size());
+    }
+
+    std::vector<int> bedges;     //!< boundary edges [nbedges][2] storing start/end vertex
+
+};
+
+#endif