Commit ce10a42c authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Initial commit mpp_src

parents
build/
.idea/
googletest/
#
# Builds static libraries mpp_src and lib_ps
#
cmake_minimum_required(VERSION 3.5.1)
project(mpp_src)
#---------------------------------------------------------------------------------------#
## Google test (Placed here to not compile it with mpiCC)
#add_subdirectory(googletest)
#include_directories(${PROJECT_BINARY_DIR}/googletest/googletest/include)
#include_directories(${PROJECT_BINARY_DIR}/googlemock/googletest/include)
#set(GTEST_LIB gtest gtest_main)
#---------------------------------------------------------------------------------------#
# MPI
find_package(MPI REQUIRED)
include_directories(${MPI_INCLUDE_PATH}) #/usr/include/mpi/
set(CMAKE_C_COMPILER mpicc) #/usr/bin/mpicc
set(CMAKE_CXX_COMPILER mpiCC) #/usr/bin/mpiCC
set(CMAKE_CXX_COMPILE_FLAGS ${CMAKE_CXX_COMPILE_FLAGS} ${MPI_COMPILE_FLAGS})
set(CMAKE_CXX_LINK_FLAGS ${CMAKE_CXX_LINK_FLAGS} ${MPI_LINK_FLAGS})
set(CMAKE_BUILD_TYPE distribution)
set(CMAKE_CXX_FLAGS_DISTRIBUTION "-Ofast")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O0")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O1")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O2")
set(CMAKE_CXX_FLAGS "-fPIC -g -Wno-deprecated -std=gnu++11")
#---------------------------------------------------------------------------------------#
# Blas Lapack
find_package(BLAS REQUIRED)
find_package(LAPACK REQUIRED)
#---------------------------------------------------------------------------------------#
# Manage folder structure build folder
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/vtk)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/log)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/dual)
#---------------------------------------------------------------------------------------#
# Include and manage libraries
include_directories(${PROJECT_SOURCE_DIR}/src)
include_directories(${PROJECT_SOURCE_DIR}/src/LIB_PS)
include_directories(${PROJECT_SOURCE_DIR}superlu/include)
link_directories(${PROJECT_SOURCE_DIR}/superlu/lib)
#---------------------------------------------------------------------------------------#
# Subdirectories
add_subdirectory(src)
This diff is collapsed.
// $Header: https://svn.math.kit.edu/svn/M++/trunk/src/Algebra.h 214 2014-10-06 08:34:49Z maurer $
#ifndef _ALGEBRA_H_
#define _ALGEBRA_H_
#include "Tensor.h"
#include "IO.h"
#include "MatrixGraph.h"
#include "Interface.h"
class BasicVector {
Scalar* a;
int n;
public:
BasicVector (int N);
BasicVector (Scalar b, int N);
BasicVector (const BasicVector& u);
~BasicVector ();
int size () const { return n; }
Scalar operator [] (int i) const { return a[i]; }
Scalar& operator [] (int i) { return a[i]; }
const Scalar* operator () () const { return a; }
Scalar* operator () () { return a; }
BasicVector& Basis (int k);
BasicVector& operator = (const BasicVector& u);
BasicVector& operator = (Scalar b);
BasicVector& operator *= (Scalar b);
BasicVector& operator /= (Scalar b);
BasicVector& operator += (const BasicVector& u);
BasicVector& operator -= (const BasicVector& u);
void Multiply (Scalar b, const BasicVector& u);
void MultiplyPlus (Scalar b, const BasicVector& u);
void MultiplyMinus (Scalar b, const BasicVector& u);
friend ostream& operator << (ostream& s, const BasicVector& u);
friend Scalar operator * (const BasicVector& u, const BasicVector& v);
double Max () const;
double Min () const;
double norm () const;
void CleanZero ();
BasicVector& ref () { return *this; }
const BasicVector& ref () const { return *this; }
};
class Vectors;
class Operator {
public:
virtual void apply_plus (Vector&, Vector&) {
Exit ("apply_plus not implemented");
}
virtual void apply (Vector&, Vector&);
virtual void apply_minus (Vector&, Vector&);
virtual void apply_transpose_plus (Vector&, Vector&) {
Exit ("apply_transpose_plus not implemented");
}
virtual void apply_transpose (Vector& u, Vector& v);
virtual void apply_transpose_minus (Vector& u, Vector& v);
virtual void multiply_plus (Vector&, const Vector&) const {
Exit ("multiply_plus not implemented");
}
virtual void multiply (Vector& u, const Vector& v) const;
virtual void multiply_minus (Vector& u, const Vector& v) const;
virtual void multiply_transpose_plus (Vector&, const Vector&) const {
Exit ("multiply_transpose_plus not implemented");
}
virtual void multiply_transpose (Vector& u, const Vector& v) const;
virtual void multiply_transpose_minus (Vector& u, const Vector& v) const;
virtual ~Operator () {};
virtual void apply_plus (Vectors& u, Vectors& v);
virtual void apply (Vectors&, Vectors&);
virtual void apply_minus (Vectors&, Vectors&);
virtual void apply_transpose_plus (Vectors& u, Vectors& v);
virtual void apply_transpose (Vectors& u, Vectors& v);
virtual void apply_transpose_minus (Vectors& u, Vectors& v);
virtual void multiply_plus (Vectors& u, const Vectors& v) const;
virtual void multiply (Vectors& u, const Vectors& v) const;
virtual void multiply_minus (Vectors& u, const Vectors& v) const;
virtual void multiply_transpose_plus (Vectors& u, const Vectors& v) const;
virtual void multiply_transpose (Vectors& u, const Vectors& v) const;
virtual void multiply_transpose_minus (Vectors& u, const Vectors& v) const;
};
class DirichletFlags {
bool* f;
int n;
bool Delete;
public:
void ClearDirichletFlags ();
DirichletFlags (int N);
DirichletFlags (const DirichletFlags& dir);
~DirichletFlags ();
const bool* D () const { return f; }
bool* D () { return f; }
};
class Vector : public matrixgraph, public BasicVector, public DirichletFlags {
public:
Vector (const matrixgraph& g);
Vector (const Vector& u);
Vector (Scalar b, const Vector& u);
Vector (const constAB<Operator, Vector>& Ov);
int size () const { return BasicVector::size(); }
Vector& operator = (const Vector& u);
Vector& operator = (Scalar b);
Vector& operator *= (Scalar b);
Vector& operator *= (int b);
Vector& operator /= (Scalar b);
Vector& operator += (const Vector& u);
Vector& operator -= (const Vector& u);
Vector& operator = (const constAB<Scalar, Vector>& au);
Vector& operator += (const constAB<Scalar, Vector>& au);
Vector& operator -= (const constAB<Scalar, Vector>& au);
Vector& operator = (const constAB<int, Vector>& au);
Vector& operator = (const constAB<Vector, Operator>& vO);
Vector& operator += (const constAB<Vector, Operator>& vO);
Vector& operator -= (const constAB<Vector, Operator>& vO);
Vector& operator = (const constAB<Operator, Vector>& Ov);
Vector& operator += (const constAB<Operator, Vector>& Ov);
Vector& operator -= (const constAB<Operator, Vector>& Ov);
const Scalar* operator () () const { return BasicVector::operator()(); }
Scalar* operator () () { return BasicVector::operator()(); }
const Scalar* operator () (int i) const { return (*this)() + Index (i); }
Scalar* operator () (int i) { return (*this)() + Index (i); }
const Scalar* operator () (const row& r) const { return (*this) (r.Id()); }
Scalar* operator () (const row& r) { return (*this) (r.Id()); }
const Scalar* operator () (const Point& z) const {
return (*this) (find_row (z));
}
Scalar* operator () (const Point& z) { return (*this) (find_row (z)); }
Scalar operator () (int i, int k) const { return (*this) (i)[k]; }
Scalar& operator () (int i, int k) { return (*this) (i)[k]; }
Scalar operator () (const row& r, int k) const { return (*this) (r)[k]; }
Scalar& operator () (const row& r, int k) { return (*this) (r)[k]; }
Scalar operator () (const Point& z, int k) const { return (*this) (z)[k];}
Scalar& operator () (const Point& z, int k) { return (*this) (z)[k]; }
const bool* D () const { return DirichletFlags::D(); }
bool* D () { return DirichletFlags::D(); }
const bool* D (int i) const { return DirichletFlags::D() + Index (i); }
bool* D (int i) { return DirichletFlags::D() + Index (i); }
const bool* D (const row& r) const { return D (r.Id()); }
bool* D (const row& r) { return D (r.Id()); }
const bool* D (const Point& z) const { return D (find_row (z)); }
bool* D (const Point& z) { return D (find_row (z)); }
bool D (int i, int k) const { return D (i)[k]; }
bool& D (int i, int k) { return D (i)[k]; }
bool D (const row& r, int k) const { return D (r)[k]; }
bool& D (const row& r, int k) { return D (r)[k]; }
bool D (const Point& z, int k) const { return D (z)[k];}
bool& D (const Point& z, int k) { return D (z)[k]; }
void ClearDirichletValues ();
void print () const;
#ifdef NDOUBLE
Vector& operator *= (double b);
Vector& operator = (const constAB<double, Vector>& au);
Vector& operator += (const constAB<double, Vector>& au);
Vector& operator -= (const constAB<double, Vector>& au);
#endif
};
class Vectors: public matrixgraph, public DirichletFlags {
vector<Vector> V;
public:
Vectors (int n, const matrixgraph& g);
Vectors (int n, const Vector& u);
Vectors (const Vectors& u);
Vectors (const constAB<Operator, Vectors>& Ov);
int size() const {return V.size();}
// int size() const {return V[0]->size();}
const Vector& operator[] (int i) const {return V[i];}
Vector& operator[] (int i) {return V[i];}
void erase (int i) {
V.erase (V.begin() + i);
}
Vectors& operator = (const Vector& u);
Vectors& operator = (const Vectors& u);
Vectors& operator = (Scalar b);
Vectors& operator *= (Scalar b);
Vectors& operator *= (int b);
Vectors& operator /= (Scalar b);
Vectors& operator += (const Vector& u);
Vectors& operator += (const Vectors& u);
Vectors& operator -= (const Vector& u);
Vectors& operator -= (const Vectors& u);
Vectors& operator = (const constAB<Scalar, Vectors>& au);
Vectors& operator += (const constAB<Scalar, Vectors>& au);
Vectors& operator -= (const constAB<Scalar, Vectors>& au);
Vectors& operator = (const constAB<int, Vectors>& au);
Vectors& operator = (const constAB<Vectors, Operator>& vO);
Vectors& operator += (const constAB<Vectors, Operator>& vO);
Vectors& operator -= (const constAB<Vectors, Operator>& vO);
Vectors& operator = (const constAB<Operator, Vectors>& Ov);
Vectors& operator += (const constAB<Operator, Vectors>& Ov);
Vectors& operator -= (const constAB<Operator, Vectors>& Ov);
const Scalar* operator () (int n, int i) const { return (V[n]) (i); }
Scalar* operator () (int n, int i) { return (V[n]) (i); }
const Scalar* operator () (int n, const row& r) const { return (V[n]) (r.Id()); }
Scalar* operator () (int n, const row& r) { return (V[n]) (r.Id()); }
Scalar operator () (int n, int i, int k) const { return (V[n]) (i, k); }
Scalar& operator () (int n, int i, int k) { return (V[n]) (i, k); }
Scalar operator () (int n, const row& r, int k) const { return (V[n]) (r, k); }
Scalar& operator () (int n, const row& r, int k) { return (V[n]) (r, k); }
const bool* D () const { return V[0].D(); }
bool* D () { return V[0].D(); }
const bool* D (int i) const { return V[0].D (i); }
bool* D (int i) { return V[0].D (i); }
const bool* D (const row& r) const { return V[0].D (r); }
bool* D (const row& r) { return V[0].D (r); }
const bool* D (const Point& z) const { return V[0].D (z); }
bool* D (const Point& z) { return V[0].D (z); }
bool D (int i, int k) const { return V[0].D (i, k); }
bool& D (int i, int k) { return V[0].D (i, k); }
bool D (const row& r, int k) const { return V[0].D (r, k); }
bool& D (const row& r, int k) { return V[0].D (r, k); }
bool D (const Point& z, int k) const { return V[0].D (z, k);}
bool& D (const Point& z, int k) { return V[0].D (z, k); }
void ClearDirichletValues ();
void print () const;
#ifdef NDOUBLE
Vector& operator *= (double b);
Vector& operator = (const constAB<double, Vector>& au);
Vector& operator += (const constAB<double, Vector>& au);
Vector& operator -= (const constAB<double, Vector>& au);
#endif // of ifdef NDOUBLE
vector<double> norm() {
vector<double> no (V.size());
for (unsigned int i = 0; i < V.size(); ++i)
no[i] = V[i].norm();
return no;
}
};
Vector& operator + (const Vector& u, const Vector& v);
Vector& operator - (const Vector& u, const Vector& v);
Scalar operator * (const Vector& u, const Vector& v);
ostream& operator << (ostream& s, const Vector& u);
double norm (const Vector& u);
#ifndef DCOMPLEX
double norm (const Vector& u, int i);
#endif
Vectors& operator + (const Vectors& u, const Vectors& v);
Vectors& operator - (const Vectors& u, const Vectors& v);
class RowValues {
Scalar* a[MaxNodalPoints];
void fill (Vector& u, const rows& R);
public:
RowValues (Vector& u, const rows& R);
RowValues (Vector& u, const cell& c);
Scalar& operator () (int i, int k = 0) { return a[i][k]; }
};
class RowBndValues {
BFParts BF;
Scalar* a[MaxNodalPoints];
bool* b[MaxNodalPoints];
void fill (Vector& u, const rows& R);
public:
RowBndValues (Vector& u, const cell& c);
Scalar& operator () (int i, int k = 0) { return a[i][k]; }
bool& D (int i, int k = 0) { return b[i][k]; }
bool onBnd () const { return BF.onBnd(); }
int bc (int i) const { return BF[i]; }
};
template <class F> Vector& operator << (Vector& u, const F& f);
Vector& operator << (Vector& u, int (*F)());
void RandomVector (Vector& u, Operator& B, Operator& IA);
void RandomVectors (Vectors& u, Operator& B, Operator& IA);
class Matrix : public matrixgraph, public BasicVector, public Operator {
const Vector& u;
public:
Matrix (const Vector& U) : matrixgraph (U), BasicVector (U.Size()), u (U) {}
Matrix& operator = (Scalar b) { BasicVector::operator= (b); return *this; }
const Vector& GetVector () const { return u; }
int size () const { return matrixgraph::size(); }
const Scalar* operator () () const { return BasicVector::operator()(); }
Scalar* operator () () { return BasicVector::operator()(); }
const Scalar* operator () (int d) const { return (*this)() + Entry (d); }
Scalar* operator () (int d) { return (*this)() + Entry (d); }
const Scalar* operator () (const row& r0, const row& r1) const {
return (*this)() + GetEntry (r0, r1);
}
Scalar* operator () (const row& r0, const row& r1) {
return (*this)() + GetEntry (r0, r1);
}
Matrix& operator = (const Matrix& A);
Matrix& operator += (const constAB<Scalar, Matrix>& aA);
void multiply (Vectors& u, const Vectors& v) const {
u = 0; multiply_plus (u, v);
}
void multiply_plus (Vector&, const Vector&) const;
void multiply_transpose_plus (Vector&, const Vector&) const;
void multiply_plus (Vectors&, const Vectors&) const;
void multiply_transpose_plus (Vectors&, const Vectors&) const;
void copy (Scalar*, int*, int*) const;
void ClearDirichletValues ();
void Symmetric ();
void print () const;
void EliminateDirichlet();
int rowsize() const {return matrixgraph::rowsize();}
};
ostream& operator << (ostream& s, const Matrix& u);
constAB<Operator, Vector> operator * (const Matrix& A, const Vector& v);
constAB<Operator, Vectors> operator * (const Matrix& A, const Vectors& v);
class RowEntries {
Scalar* a[MaxNodalPoints][MaxNodalPoints];
int n[MaxNodalPoints];
void fill (Matrix& A, const rows& R);
public:
RowEntries (Matrix& A, const rows& R);
RowEntries (Matrix& A, const cell& c);
Scalar& operator () (int i, int j, int k = 0, int l = 0) {
return a[i][j][k * n[j] + l];
}
};
#endif // of #ifndef _ALGEBRA_H_
This diff is collapsed.
set(mpp_src Identify.C
MatrixGraph.C
Algebra.C
Shape.C
Interface.C
Discretization.C
Point.C
Parallel.C
TimeDate.C
ctools.C
IO.C
Mesh.C
Cell.C
Distribution.C
Plot.C
SaveLoad.C
Quadrature.C
DoF.C
Sparse.C
LinearSolver.C
Preconditioner.C
TimeSeries.C
ESolver.C
Transfer.C
Element.C
FaceElement.C
Newton.C
Euler.C
Small.C
LapackMatrix.C
SVD.C
LIB_PS_Solver.C)
set(lib_ps
LIB_PS/LIB_PS.C
LIB_PS/LIB_PS_Matrix.C
LIB_PS/LIB_PS_Parallel.C)
add_library(mpp_src STATIC ${mpp_src})
add_library(lib_ps STATIC ${lib_ps})
This diff is collapsed.
// $Header: https://svn.math.kit.edu/svn/M++/trunk/src/Cell.h 214 2014-10-06 08:34:49Z maurer $
#ifndef _CELL_H_
#define _CELL_H_
#include "Point.h"
#include "Shape.h"
#include <vector>
enum CELLTYPE { NONE = 0,
INTERVAL = 122,
TINTERVAL = 144,
TRIANGLE = 233,
QUADRILATERAL = 244, QUADRILATERAL2 = 248,
TETRAHEDRON = 344,
PYRAMID = 355,
PRISM = 366,
HEXAHEDRON = 388, HEXAHEDRON20 = 3820, HEXAHEDRON27 = 3827
};
class Rule {
CELLTYPE tp;
public:
vector<short> node;
vector<short> face;
Rule (CELLTYPE TP, int n, const short* C) : tp (TP) {
node.resize (n);
for (int j = 0; j < n; ++j) node[j] = C[j];
}
Rule (const Rule& R) : tp (R.tp), node (R.node), face (R.face) {}
Rule () : tp (NONE) {}
CELLTYPE type() const { return tp; }
void operator () (const vector<Point>& z, vector<Point>& x) const {
x.resize (node.size());
for (unsigned int j = 0; j < node.size(); ++j) x[j] = z[node[j]];
}
friend ostream& operator << (ostream& s, const Rule& R) {
s << "z ";
for (unsigned int i = 0; i < R.node.size(); ++i) s << R.node[i] << " ";
s << "f ";
for (unsigned int i = 0; i < R.face.size(); ++i) s << R.face[i] << " ";
return s << "tp " << R.type();
}
int operator [] (int i) const { return node[i]; }
};
class Rules : public vector<Rule> {
public:
Rules (int n, const Rule* R) : vector<Rule> (n) {
for (int i = 0; i < n; ++i) (*this)[i] = R[i];
}
};
class Mesh;
class Cell : public vector<Point> {
short subdomain;
protected:
Cell (const vector<Point>& x, short sd) :
vector<Point> (x), subdomain (sd) {}
const Point& z (int i) const {return (*this).vector<Point>::operator[] (i);}
public:
virtual Point operator () () const = 0;
Point Center () const { return (*this)(); }
virtual CELLTYPE Type () const = 0;
virtual CELLTYPE ReferenceType () const = 0;
virtual const Cell* ReferenceCell () const = 0;
short Subdomain () const { return subdomain; }
virtual int Corners () const = 0;
const Point& operator [] (int i) const {
return (*this).vector<Point>::operator[] (i);
}
const Point& Corner (int i) const { return (*this)[i]; }
virtual const Point& LocalCorner (int) const = 0;
virtual const Point& LocalFaceNormal (int) const = 0;
virtual Point LocalToGlobal (const Point&) const = 0;
virtual Transformation GetTransformation (const Point&) const = 0;
virtual int Edges () const = 0;
virtual Point Edge (int) const = 0;
virtual const Point& LocalEdge (int) const = 0;
virtual const Point& LocalFace (int) const = 0;
virtual const Point& LocalCenter () const = 0;
virtual const Point& EdgeCorner (int, int) const = 0;
virtual short edgecorner (int, int) const = 0;
virtual int Faces () const = 0;
virtual Point Face (int) const = 0;
virtual short FarFace (int) const { return -1; }
virtual short MiddleEdge (int i, int j) const { return -1; }
virtual bool faceedgecircuit (int i) const { return false; }
virtual short FaceCorners (int) const = 0;
virtual const Point& FaceCorner (int, int) const = 0;
virtual short facecorner (int, int) const = 0;
virtual short FaceEdges (int) const = 0;
virtual short faceedge (int, int) const = 0;
Point FaceEdge (int i, int j) const { return Edge (faceedge (i, j)); }
short faceedgecorner (int i, int j, int k) const {
return edgecorner (faceedge (i, j), k);
}
const Point& FaceEdgeCorner (int i, int j, int k) const {
return EdgeCorner (faceedge (i, j), k);
}
virtual double LocalFaceArea (int) const = 0;
virtual double FaceArea (int) const = 0;
virtual const Rules& Refine (vector<Point>&) const = 0;
virtual int dim () const = 0;
virtual bool plane () const = 0;
virtual int Children () const = 0;
virtual Point Child (int i) const = 0;
Point Neighbor (const Mesh&, int) const;
virtual ~Cell () {}
};
inline ostream& operator << (ostream& s, const Cell& C) {
for (unsigned int i = 0; i < C.size(); ++i) s << C[i] << "|";
return s << " tp " << int (C.Type()) << " sd " << C.Subdomain();
}
inline ostream& operator << (ostream& s, const Cell* C) { return s << *C; }
class cell : public hash_map<Point, Cell*, Hash>::const_iterator {
typedef hash_map<Point, Cell*, Hash>::const_iterator Iterator;
public:
cell () {}
cell (Iterator c) : Iterator (c) {}
int size () const { return (*this)->second->size(); }
const Point& operator () () const { return (*this)->first; }
const Point& Center () const { return (*this)->first; }
const Cell& operator * () const { return * ((*this)->second); }
CELLTYPE Type () const { return (*this)->second->Type(); }
CELLTYPE ReferenceType () const { return (*this)->second->ReferenceType();}
const Cell* ReferenceCell () { return (*this)->second->ReferenceCell(); }
short Subdomain () const { return (*this)->second->Subdomain(); }
int Corners () const { return (*this)->second->Corners(); }
const Point& Corner (int i) const { return (*this)->second->Corner (i); }
const Point& LocalCorner (int i) const {
return (*this)->second->LocalCorner (i);
}
const Point& LocalEdge (int i) const {
return (*this)->second->LocalEdge (i);
}
const Point& LocalFace (int i) const {
return (*this)->second->LocalFace (i);
}
const Point& LocalCenter () const { return (*this)->second->LocalCenter();}
const Point& LocalFaceNormal (int