Upload files to "ex3_benchmarks"
This commit is contained in:
parent
8d998191a9
commit
50845d49be
5 changed files with 985 additions and 0 deletions
381
ex3_benchmarks/geom.h
Normal file
381
ex3_benchmarks/geom.h
Normal file
|
|
@ -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
|
||||
347
ex3_benchmarks/getmatrix.cpp
Normal file
347
ex3_benchmarks/getmatrix.cpp
Normal file
|
|
@ -0,0 +1,347 @@
|
|||
#include "getmatrix.h"
|
||||
#include "userset.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <list>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
|
||||
// general routine for lin. triangular elements
|
||||
|
||||
void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3])
|
||||
//void CalcElem(const int* __restrict__ ial, const double* __restrict__ xc, double* __restrict__ ske[3], double* __restrict__ fe)
|
||||
{
|
||||
const int i1 = 2 * ial[0], i2 = 2 * ial[1], i3 = 2 * ial[2];
|
||||
const double x13 = xc[i3 + 0] - xc[i1 + 0], y13 = xc[i3 + 1] - xc[i1 + 1],
|
||||
x21 = xc[i1 + 0] - xc[i2 + 0], y21 = xc[i1 + 1] - xc[i2 + 1],
|
||||
x32 = xc[i2 + 0] - xc[i3 + 0], y32 = xc[i2 + 1] - xc[i3 + 1];
|
||||
const double jac = fabs(x21 * y13 - x13 * y21);
|
||||
|
||||
ske[0][0] = 0.5 / jac * (y32 * y32 + x32 * x32);
|
||||
ske[0][1] = 0.5 / jac * (y13 * y32 + x13 * x32);
|
||||
ske[0][2] = 0.5 / jac * (y21 * y32 + x21 * x32);
|
||||
ske[1][0] = ske[0][1];
|
||||
ske[1][1] = 0.5 / jac * (y13 * y13 + x13 * x13);
|
||||
ske[1][2] = 0.5 / jac * (y21 * y13 + x21 * x13);
|
||||
ske[2][0] = ske[0][2];
|
||||
ske[2][1] = ske[1][2];
|
||||
ske[2][2] = 0.5 / jac * (y21 * y21 + x21 * x21);
|
||||
|
||||
const double xm = (xc[i1 + 0] + xc[i2 + 0] + xc[i3 + 0]) / 3.0,
|
||||
ym = (xc[i1 + 1] + xc[i2 + 1] + xc[i3 + 1]) / 3.0;
|
||||
//fe[0] = fe[1] = fe[2] = 0.5 * jac * FunctF(xm, ym) / 3.0;
|
||||
fe[0] = fe[1] = fe[2] = 0.5 * jac * fNice(xm, ym) / 3.0;
|
||||
}
|
||||
|
||||
|
||||
// general routine for lin. triangular elements,
|
||||
// non-symm. matrix
|
||||
// node numbering in element: a s c e n d i n g indices !!
|
||||
// GH: deprecated
|
||||
void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
|
||||
int const id[], int const ik[], double sk[], double f[])
|
||||
{
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
const int ii = ial[i], // row ii (global index)
|
||||
id1 = id[ii], // start and
|
||||
id2 = id[ii + 1]; // end of row ii in matrix
|
||||
int ip = id1;
|
||||
for (int j = 0; j < 3; ++j) // no symmetry assumed
|
||||
{
|
||||
const int jj = ial[j];
|
||||
bool not_found = true;
|
||||
do // find entry jj (global index) in row ii
|
||||
{
|
||||
not_found = (ik[ip] != jj);
|
||||
++ip;
|
||||
}
|
||||
while (not_found && ip < id2);
|
||||
|
||||
#ifndef NDEBUG // compiler option -DNDEBUG switches off the check
|
||||
if (not_found) // no entry found !!
|
||||
{
|
||||
cout << "Error in AddElem: (" << ii << "," << jj << ") ["
|
||||
<< ial[0] << "," << ial[1] << "," << ial[2] << "]\n";
|
||||
assert(!not_found);
|
||||
}
|
||||
#endif
|
||||
sk[ip - 1] += ske[i][j];
|
||||
}
|
||||
f[ii] += fe[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
|
||||
|
||||
|
||||
// ####################################################################
|
||||
|
||||
CRS_Matrix::CRS_Matrix(Mesh const &mesh)
|
||||
: _mesh(mesh), _nrows(0), _nnz(0), _id(0), _ik(0), _sk(0)
|
||||
{
|
||||
Derive_Matrix_Pattern();
|
||||
return;
|
||||
}
|
||||
|
||||
void CRS_Matrix::Derive_Matrix_Pattern()
|
||||
{
|
||||
int const nelem(_mesh.Nelems());
|
||||
int const ndof_e(_mesh.NdofsElement());
|
||||
auto const &ia(_mesh.GetConnectivity());
|
||||
// Determine the number of matrix rows
|
||||
_nrows = *max_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem);
|
||||
++_nrows; // node numberng: 0 ... nnode-1
|
||||
assert(*min_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem) == 0); // numbering starts with 0 ?
|
||||
|
||||
// Collect for each node those nodes it is connected to (multiple entries)
|
||||
// Detect the neighboring nodes
|
||||
vector< list<int> > cc(_nrows); // cc[i] is the list of nodes a node i is connected to
|
||||
for (int i = 0; i < nelem; ++i)
|
||||
{
|
||||
int const idx = ndof_e * i;
|
||||
for (int k = 0; k < ndof_e; ++k)
|
||||
{
|
||||
list<int> &cck = cc.at(ia[idx + k]);
|
||||
cck.insert( cck.end(), ia.cbegin() + idx, ia.cbegin() + idx + ndof_e );
|
||||
}
|
||||
}
|
||||
// Delete the multiple entries
|
||||
_nnz = 0;
|
||||
for (auto &it : cc)
|
||||
{
|
||||
it.sort();
|
||||
it.unique();
|
||||
_nnz += static_cast<int>(it.size());
|
||||
// cout << it.size() << " :: "; copy(it->begin(),it->end(), ostream_iterator<int,char>(cout," ")); cout << endl;
|
||||
}
|
||||
|
||||
// CSR data allocation
|
||||
_id.resize(_nrows + 1); // Allocate memory for CSR row pointer
|
||||
_ik.resize(_nnz); // Allocate memory for CSR column index vector
|
||||
|
||||
// copy CSR data
|
||||
_id[0] = 0; // begin of first row
|
||||
for (size_t i = 0; i < cc.size(); ++i)
|
||||
{
|
||||
//cout << i << " " << nid.at(i) << endl;;
|
||||
const list<int> &ci = cc.at(i);
|
||||
const auto nci = static_cast<int>(ci.size());
|
||||
_id[i + 1] = _id[i] + nci; // begin of next line
|
||||
copy(ci.begin(), ci.end(), _ik.begin() + _id[i] );
|
||||
}
|
||||
|
||||
assert(_nnz == _id[_nrows]);
|
||||
_sk.resize(_nnz); // Allocate memory for CSR column index vector
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
void CRS_Matrix::Debug() const
|
||||
{
|
||||
// ID points to first entry of row
|
||||
// no symmetry assumed
|
||||
cout << "\nMatrix (nnz = " << _id[_nrows] << ")\n";
|
||||
|
||||
for (int row = 0; row < _nrows; ++row)
|
||||
{
|
||||
cout << "Row " << row << " : ";
|
||||
int const id1 = _id[row];
|
||||
int const id2 = _id[row + 1];
|
||||
for (int j = id1; j < id2; ++j)
|
||||
{
|
||||
cout.setf(ios::right, ios::adjustfield);
|
||||
cout << "[" << setw(2) << _ik[j] << "] " << setw(4) << _sk[j] << " ";
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void CRS_Matrix::CalculateLaplace(vector<double> &f)
|
||||
{
|
||||
assert(_mesh.NdofsElement() == 3); // only for triangular, linear elements
|
||||
//cout << _nnz << " vs. " << _id[_nrows] << " " << _nrows<< endl;
|
||||
assert(_nnz == _id[_nrows]);
|
||||
|
||||
for (int k = 0; k < _nrows; ++k)
|
||||
{
|
||||
_sk[k] = 0.0;
|
||||
}
|
||||
for (int k = 0; k < _nrows; ++k)
|
||||
{
|
||||
f[k] = 0.0;
|
||||
}
|
||||
|
||||
double ske[3][3], fe[3];
|
||||
// Loop over all elements
|
||||
auto const nelem = _mesh.Nelems();
|
||||
auto const &ia = _mesh.GetConnectivity();
|
||||
auto const &xc = _mesh.GetCoords();
|
||||
|
||||
for (int i = 0; i < nelem; ++i)
|
||||
{
|
||||
CalcElem(ia.data() + 3 * i, xc.data(), ske, fe);
|
||||
//AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
|
||||
AddElem_3(ia.data() + 3 * i, ske, fe, f);
|
||||
}
|
||||
|
||||
//Debug();
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void CRS_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f)
|
||||
{
|
||||
double const PENALTY = 1e6;
|
||||
auto const idx = _mesh.Index_DirichletNodes();
|
||||
int const nidx = static_cast<int>(idx.size());
|
||||
|
||||
for (int row = 0; row < nidx; ++row)
|
||||
{
|
||||
int const k = idx[row];
|
||||
int const id1 = fetch(k, k); // Find diagonal entry of row
|
||||
assert(id1 >= 0);
|
||||
_sk[id1] += PENALTY; // matrix weighted scaling feasible
|
||||
f[k] += PENALTY * u[k];
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void CRS_Matrix::GetDiag(vector<double> &d) const
|
||||
{
|
||||
assert( _nrows == static_cast<int>(d.size()) );
|
||||
|
||||
for (int row = 0; row < _nrows; ++row)
|
||||
{
|
||||
const int ia = fetch(row, row); // Find diagonal entry of row
|
||||
assert(ia >= 0);
|
||||
d[row] = _sk[ia];
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
bool CRS_Matrix::Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const
|
||||
{
|
||||
bool bn = (nnode == _nrows); // number of rows
|
||||
if (!bn)
|
||||
{
|
||||
cout << "######### Error: " << "number of rows" << endl;
|
||||
}
|
||||
|
||||
bool bz = (id[nnode] == _nnz); // number of non zero elements
|
||||
if (!bz)
|
||||
{
|
||||
cout << "######### Error: " << "number of non zero elements" << endl;
|
||||
}
|
||||
|
||||
bool bd = equal(id, id + nnode + 1, _id.cbegin()); // row starts
|
||||
if (!bd)
|
||||
{
|
||||
cout << "######### Error: " << "row starts" << endl;
|
||||
}
|
||||
|
||||
bool bk = equal(ik, ik + id[nnode], _ik.cbegin()); // column indices
|
||||
if (!bk)
|
||||
{
|
||||
cout << "######### Error: " << "column indices" << endl;
|
||||
}
|
||||
|
||||
bool bv = equal(sk, sk + id[nnode], _sk.cbegin()); // values
|
||||
if (!bv)
|
||||
{
|
||||
cout << "######### Error: " << "values" << endl;
|
||||
}
|
||||
|
||||
return bn && bz && bd && bk && bv;
|
||||
}
|
||||
|
||||
|
||||
void CRS_Matrix::Mult(vector<double> &w, vector<double> const &u) const
|
||||
{
|
||||
assert( _nrows == static_cast<int>(w.size()) );
|
||||
assert( w.size() == u.size() );
|
||||
|
||||
for (int row = 0; row < _nrows; ++row)
|
||||
{
|
||||
double wi = 0.0;
|
||||
for (int ij = _id[row]; ij < _id[row + 1]; ++ij)
|
||||
{
|
||||
wi += _sk[ij] * u[ _ik[ij] ];
|
||||
}
|
||||
w[row] = wi;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void CRS_Matrix::Defect(vector<double> &w,
|
||||
vector<double> const &f, vector<double> const &u) const
|
||||
{
|
||||
assert( _nrows == static_cast<int>(w.size()) );
|
||||
assert( w.size() == u.size() && u.size() == f.size() );
|
||||
|
||||
for (int row = 0; row < _nrows; ++row)
|
||||
{
|
||||
double wi = f[row];
|
||||
for (int ij = _id[row]; ij < _id[row + 1]; ++ij)
|
||||
{
|
||||
wi -= _sk[ij] * u[ _ik[ij] ];
|
||||
}
|
||||
w[row] = wi;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
int CRS_Matrix::fetch(int const row, int const col) const
|
||||
{
|
||||
int const id2 = _id[row + 1]; // end and
|
||||
int ip = _id[row]; // start of recent row (global index)
|
||||
|
||||
while (ip < id2 && _ik[ip] != col) // find index col (global index)
|
||||
{
|
||||
++ip;
|
||||
}
|
||||
if (ip >= id2)
|
||||
{
|
||||
ip = -1;
|
||||
#ifndef NDEBUG // compiler option -DNDEBUG switches off the check
|
||||
cout << "No column " << col << " in row " << row << endl;
|
||||
assert(ip >= id2);
|
||||
#endif
|
||||
}
|
||||
return ip;
|
||||
}
|
||||
|
||||
|
||||
void CRS_Matrix::AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], vector<double> &f)
|
||||
{
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
const int ii = ial[i]; // row ii (global index)
|
||||
for (int j = 0; j < 3; ++j) // no symmetry assumed
|
||||
{
|
||||
const int jj = ial[j]; // column jj (global index)
|
||||
int ip = fetch(ii, jj); // find column entry jj in row ii
|
||||
#ifndef NDEBUG // compiler option -DNDEBUG switches off the check
|
||||
if (ip < 0) // no entry found !!
|
||||
{
|
||||
cout << "Error in AddElem: (" << ii << "," << jj << ") ["
|
||||
<< ial[0] << "," << ial[1] << "," << ial[2] << "]\n";
|
||||
assert(ip >= 0);
|
||||
}
|
||||
#endif
|
||||
_sk[ip] += ske[i][j];
|
||||
}
|
||||
f[ii] += fe[i];
|
||||
}
|
||||
}
|
||||
|
||||
178
ex3_benchmarks/getmatrix.h
Normal file
178
ex3_benchmarks/getmatrix.h
Normal file
|
|
@ -0,0 +1,178 @@
|
|||
#ifndef GETMATRIX_FILE
|
||||
#define GETMATRIX_FILE
|
||||
|
||||
#include "geom.h"
|
||||
#include <vector>
|
||||
|
||||
/**
|
||||
* Calculates the element stiffness matrix @p ske and the element load vector @p fe
|
||||
* of one triangular element with linear shape functions.
|
||||
* @param[in] ial node indices of the three element vertices
|
||||
* @param[in] xc vector of node coordinates with x(2*k,2*k+1) as coodinates of node k
|
||||
* @param[out] ske element stiffness matrix
|
||||
* @param[out] fe element load vector
|
||||
*/
|
||||
void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]);
|
||||
|
||||
/**
|
||||
* Adds the element stiffness matrix @p ske and the element load vector @p fe
|
||||
* of one triangular element with linear shape functions to the appropriate positions in
|
||||
* the symmetric stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik)
|
||||
*
|
||||
* @param[in] ial node indices of the three element vertices
|
||||
* @param[in] ske element stiffness matrix
|
||||
* @param[in] fe element load vector
|
||||
* @param[out] sk vector non-zero entries of CSR matrix
|
||||
* @param[in] id index vector containing the first entry in a CSR row
|
||||
* @param[in] ik column index vector of CSR matrix
|
||||
* @param[out] f distributed local vector storing the right hand side
|
||||
*
|
||||
* @warning Algorithm requires indices in connectivity @p ial in ascending order.
|
||||
* Currently deprecated.
|
||||
*/
|
||||
void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
|
||||
int const id[], int const ik[], double sk[], double f[]);
|
||||
|
||||
|
||||
// #####################################################################
|
||||
/**
|
||||
* Square matrix in CRS format (compressed row storage; also named CSR),
|
||||
* see an <a href="https://en.wikipedia.org/wiki/Sparse_matrix">introduction</a>.
|
||||
*/
|
||||
class CRS_Matrix
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Intializes the CRS matrix structure from the given discetization in @p mesh.
|
||||
*
|
||||
* The sparse matrix pattern is generated but the values are 0.
|
||||
*
|
||||
* @param[in] mesh given discretization
|
||||
*
|
||||
* @warning A reference to the discretization @p mesh is stored inside this class.
|
||||
* Therefore, changing @p mesh outside requires also
|
||||
* to call method @p Derive_Matrix_Pattern explicitely.
|
||||
*
|
||||
* @see Derive_Matrix_Pattern
|
||||
*/
|
||||
explicit CRS_Matrix(Mesh const & mesh);
|
||||
|
||||
/**
|
||||
* Destructor.
|
||||
*/
|
||||
~CRS_Matrix()
|
||||
{}
|
||||
|
||||
/**
|
||||
* Generates the sparse matrix pattern and overwrites the existing pattern.
|
||||
*
|
||||
* The sparse matrix pattern is generated but the values are 0.
|
||||
*/
|
||||
void Derive_Matrix_Pattern();
|
||||
|
||||
/**
|
||||
* Calculates the entries of f.e. stiffness matrix and load/rhs vector @p f for the Laplace operator in 2D.
|
||||
* No memory is allocated.
|
||||
*
|
||||
* @param[in,out] f (preallocated) rhs/load vector
|
||||
*/
|
||||
void CalculateLaplace(std::vector<double> &f);
|
||||
|
||||
/**
|
||||
* Applies Dirichlet boundary conditions to stiffness matrix and to load vector @p f.
|
||||
* The <a href="https://www.jstor.org/stable/2005611?seq=1#metadata_info_tab_contents">penalty method</a>
|
||||
* is used for incorporating the given values @p u.
|
||||
*
|
||||
* @param[in] u (global) vector with Dirichlet data
|
||||
* @param[in,out] f load vector
|
||||
*/
|
||||
void ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f);
|
||||
|
||||
/**
|
||||
* Extracts the diagonal elemenst of the sparse matrix.
|
||||
*
|
||||
* @param[in,out] d (prellocated) vector of diagonal elements
|
||||
*/
|
||||
void GetDiag(std::vector<double> &d) const;
|
||||
|
||||
/**
|
||||
* Performs the matrix-vector product w := K*u.
|
||||
*
|
||||
* @param[in,out] w resulting vector (preallocated)
|
||||
* @param[in] u vector
|
||||
*/
|
||||
void Mult(std::vector<double> &w, std::vector<double> const &u) const;
|
||||
|
||||
/**
|
||||
* Calculates the defect/residuum w := f - K*u.
|
||||
*
|
||||
* @param[in,out] w resulting vector (preallocated)
|
||||
* @param[in] f load vector
|
||||
* @param[in] u vector
|
||||
*/
|
||||
void Defect(std::vector<double> &w,
|
||||
std::vector<double> const &f, std::vector<double> const &u) const;
|
||||
|
||||
/**
|
||||
* Number rows in matrix.
|
||||
* @return number of rows.
|
||||
*/
|
||||
int Nrows() const
|
||||
{return _nrows;}
|
||||
|
||||
/**
|
||||
* Show the matrix entries.
|
||||
*/
|
||||
void Debug() const;
|
||||
|
||||
/**
|
||||
* Finds in a CRS matrix the access index for an entry at row @p row and column @p col.
|
||||
*
|
||||
* @param[in] row row index
|
||||
* @param[in] col column index
|
||||
* @return index for element (@p row, @p col). If no appropriate entry exists then -1 will be returned.
|
||||
*
|
||||
* @warning assert() stops the function in case that matrix element (@p row, @p col) doesn't exist.
|
||||
*/
|
||||
int fetch(int row, int col) const;
|
||||
|
||||
/**
|
||||
* Adds the element stiffness matrix @p ske and the element load vector @p fe
|
||||
* of one triangular element with linear shape functions to the appropriate positions in
|
||||
* the stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik).
|
||||
*
|
||||
* @param[in] ial node indices of the three element vertices
|
||||
* @param[in] ske element stiffness matrix
|
||||
* @param[in] fe element load vector
|
||||
* @param[in,out] f distributed local vector storing the right hand side
|
||||
*
|
||||
* @warning Algorithm assumes linear triangular elements (ndof_e==3).
|
||||
*/
|
||||
void AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], std::vector<double> &f);
|
||||
|
||||
/**
|
||||
* Compare @p this CRS matrix with an external CRS matrix stored in C-Style.
|
||||
*
|
||||
* The method prints statements on differences found.
|
||||
*
|
||||
* @param[in] nnode row number of external matrix
|
||||
* @param[in] id start indices of matrix rows of external matrix
|
||||
* @param[in] ik column indices of external matrix
|
||||
* @param[in] sk non-zero values of external matrix
|
||||
*
|
||||
* @return true iff all data are identical.
|
||||
*/
|
||||
bool Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const;
|
||||
|
||||
private:
|
||||
Mesh const & _mesh; //!< reference to discretization
|
||||
int _nrows; //!< number of rows in matrix
|
||||
int _nnz; //!< number of non-zero entries
|
||||
std::vector<int> _id; //!< start indices of matrix rows
|
||||
std::vector<int> _ik; //!< column indices
|
||||
std::vector<double> _sk; //!< non-zero values
|
||||
|
||||
};
|
||||
|
||||
|
||||
#endif
|
||||
61
ex3_benchmarks/jacsolve.cpp
Normal file
61
ex3_benchmarks/jacsolve.cpp
Normal file
|
|
@ -0,0 +1,61 @@
|
|||
#include "vdop.h"
|
||||
#include "getmatrix.h"
|
||||
#include "jacsolve.h"
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
// #####################################################################
|
||||
// const int neigh[], const int color,
|
||||
// const MPI::Intracomm& icomm,
|
||||
void JacobiSolve(CRS_Matrix const &SK, vector<double> const &f, vector<double> &u)
|
||||
{
|
||||
const double omega = 1.0;
|
||||
const int maxiter = 1000;
|
||||
const double tol = 1e-5, // tolerance
|
||||
tol2 = tol * tol; // tolerance^2
|
||||
|
||||
int nrows = SK.Nrows(); // number of rows == number of columns
|
||||
assert( nrows == static_cast<int>(f.size()) && f.size() == u.size() );
|
||||
|
||||
cout << endl << " Start Jacobi solver for " << nrows << " d.o.f.s" << endl;
|
||||
// Choose initial guess
|
||||
for (int k = 0; k < nrows; ++k)
|
||||
{
|
||||
u[k] = 0.0; // u := 0
|
||||
}
|
||||
|
||||
vector<double> dd(nrows); // matrix diagonal
|
||||
vector<double> r(nrows); // residual
|
||||
vector<double> w(nrows); // correction
|
||||
|
||||
SK.GetDiag(dd); // dd := diag(K)
|
||||
////DebugVector(dd);{int ijk; cin >> ijk;}
|
||||
|
||||
// Initial sweep
|
||||
SK.Defect(r, f, u); // r := f - K*u
|
||||
|
||||
vddiv(w, r, dd); // w := D^{-1}*r
|
||||
double sigma0 = dscapr(w, r); // s0 := <w,r>
|
||||
|
||||
// Iteration sweeps
|
||||
int iter = 0;
|
||||
double sigma = sigma0;
|
||||
while ( sigma > tol2 * sigma0 && maxiter > iter)
|
||||
{
|
||||
++iter;
|
||||
vdaxpy(u, u, omega, w ); // u := u + om*w
|
||||
SK.Defect(r, f, u); // r := f - K*u
|
||||
vddiv(w, r, dd); // w := D^{-1}*r
|
||||
sigma = dscapr(w, r); // s0 := <w,r>
|
||||
// cout << "Iteration " << iter << " : " << sqrt(sigma/sigma0) << endl;
|
||||
}
|
||||
cout << "aver. Jacobi rate : " << exp(log(sqrt(sigma / sigma0)) / iter) << " (" << iter << " iter)" << endl;
|
||||
cout << "final error: " << sqrt(sigma / sigma0) << " (rel) " << sqrt(sigma) << " (abs)\n";
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
18
ex3_benchmarks/jacsolve.h
Normal file
18
ex3_benchmarks/jacsolve.h
Normal file
|
|
@ -0,0 +1,18 @@
|
|||
#ifndef JACSOLVE_FILE
|
||||
#define JACSOLVE_FILE
|
||||
#include "getmatrix.h"
|
||||
#include <vector>
|
||||
|
||||
|
||||
/**
|
||||
* Solves linear system of equations K @p u = @p f via the Jacobi iteration.
|
||||
* We use a distributed symmetric CSR matrix @p SK and initial guess of the
|
||||
* solution is set to 0.
|
||||
* @param[in] SK CSR matrix
|
||||
* @param[in] f distributed local vector storing the right hand side
|
||||
* @param[out] u accumulated local vector storing the solution.
|
||||
*/
|
||||
void JacobiSolve(CRS_Matrix const &SK, std::vector<double> const &f, std::vector<double> &u);
|
||||
|
||||
|
||||
#endif
|
||||
Loading…
Add table
Add a link
Reference in a new issue