Skip to content

GitLab

Commits (52)
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ a5aba5d3
......@@ -105,10 +105,9 @@ trigger-cardmech:
- $CI_PIPELINE_SOURCE == 'web'
trigger:
project: mpp/cardmech
branch: master
branch: debug
strategy: depend
trigger-navierstokes:
stage: downstream
variables:
......
src/Algebra.C
View file @ a5aba5d3
#include "Algebra.h"
#include "parallel/Parallel.h"
#include "Sparse.h"
#include "dof/MultiPartDoF.hpp"
/*********************/
/* class BasicVector */
......@@ -44,9 +45,8 @@ BasicVector &BasicVector::operator=(Scalar b) {
return *this;
}
BasicVector &BasicVector::operator*=(Scalar b) {
if (b == 1){
if (b == 1) {
return *this;
}
for (int i = 0; i < n; ++i) a[i] *= b;
......@@ -54,7 +54,7 @@ BasicVector &BasicVector::operator*=(Scalar b) {
}
BasicVector &BasicVector::operator/=(Scalar b) {
if (b == 1){
if (b == 1) {
return *this;
}
for (int i = 0; i < n; ++i) a[i] /= b;
......@@ -129,11 +129,10 @@ DirichletFlags::DirichletFlags(int N) : n(N), Delete(true) {
}
DirichletFlags::DirichletFlags(const DirichletFlags &dir) :
f(dir.f), n(dir.n), Delete(false) {}
f(dir.f), n(dir.n), Delete(false) {}
DirichletFlags::~DirichletFlags() { if (Delete) delete[] f; }
constAB<Vector, Operator> operator*(const Vector &v, const Operator &A) {
return constAB<Vector, Operator>(v, A);
}
......@@ -155,19 +154,20 @@ constAB<int, Vector> operator*(const int &a, const Vector &v) {
/* class Vector */
/****************/
Vector::Vector(const matrixgraph &g) :
matrixgraph(g), BasicVector(g.size()), DirichletFlags(g.size()) {}
matrixgraph(g), BasicVector(g.size()), DirichletFlags(g.size()) {}
Vector::Vector(const Vector &u) : matrixgraph(u),
BasicVector(u), DirichletFlags(u) {}
Vector::Vector(Scalar b, const Vector &u) : matrixgraph(u),
BasicVector(b, u.size()), DirichletFlags(u) {}
Vector::Vector(Scalar b, const matrixgraph &g) : matrixgraph(g),
BasicVector(b, g.size()), DirichletFlags(g.size()) {}
BasicVector(b, g.size()), DirichletFlags(g.size()) {}
Vector::Vector(const constAB<Operator, Vector> &Ov) :
matrixgraph(Ov.second()), BasicVector(Ov.second()),
DirichletFlags(Ov.second()) {
matrixgraph(Ov.second()), BasicVector(Ov.second()),
DirichletFlags(Ov.second()) {
Ov.first().multiply(*this, Ov.second());
}
......@@ -293,7 +293,7 @@ void Vector::ClearDirichletValues() {
}
void Vector::print() const {
vector <row> R(Size());
vector<row> R(Size());
for (row r = rows(); r != rows_end(); ++r)
R[r.Id()] = r;
for (int k = 0; k < R.size(); ++k) {
......@@ -369,7 +369,6 @@ constAB<Operator, Vectors> operator*(const Operator &A, const Vectors &v) {
return constAB<Operator, Vectors>(A, v);
}
constAB<Scalar, Vectors> operator*(const Scalar &a, const Vectors &v) {
return constAB<Scalar, Vectors>(a, v);
}
......@@ -378,43 +377,46 @@ constAB<int, Vectors> operator*(const int &a, const Vectors &v) {
return constAB<int, Vectors>(a, v);
}
Vectors::Vectors(int n, const matrixgraph& g): matrixgraph(g), DirichletFlags(g.size()), V(n) {
for (int i =0; i < n; ++i){
Vectors::Vectors(int n, const matrixgraph &g) : matrixgraph(g), DirichletFlags(g.size()), V(n) {
for (int i = 0; i < n; ++i) {
V[i] = new Vector(g);
}
}
Vectors::Vectors(int n, const Vector& u): matrixgraph(u), DirichletFlags(u), V(n) {
for (int i =0; i < n; ++i){
Vectors::Vectors(int n, const Vector &u) : matrixgraph(u), DirichletFlags(u), V(n) {
for (int i = 0; i < n; ++i) {
V[i] = new Vector(u);
}
}
Vectors::Vectors(const Vectors& u): matrixgraph(u[0]), DirichletFlags(u[0]), V(u.size()) {
Vectors::Vectors(const Vectors &u) : matrixgraph(u[0]), DirichletFlags(u[0]), V(u.size()) {
// !!!! Here: matrixgraph(u[0]), DirichletFlags(u[0]) --
// there are no DirichletFlags set for "Vectors"...
for (int i=0; i < V.size(); ++i){
for (int i = 0; i < V.size(); ++i) {
V[i] = new Vector(*(u.V[i]));
*(V[i]) = *(u.V[i]);
}
}
Vectors::Vectors(const constAB<Operator, Vectors> &Ov) :
matrixgraph(Ov.second()), V(Ov.second().size()),
DirichletFlags(Ov.second()) {
for (int i =0; i < V.size(); ++i){
matrixgraph(Ov.second()), V(Ov.second().size()),
DirichletFlags(Ov.second()) {
for (int i = 0; i < V.size(); ++i) {
V[i] = new Vector(Ov.second());
}
Ov.first().multiply(*this, Ov.second());
}
Vectors::~Vectors() {
for (int i=0; i< V.size();++i){
for (int i = 0; i < V.size(); ++i) {
delete V[i];
}
}
void Vectors::resize(int N) {
for (int i=0; i<size(); ++i) delete V[i];
for (int i = 0; i < size(); ++i) delete V[i];
V.resize(N);
for (int i=0; i<size(); ++i) V[i] = new Vector(*this);
for (int i = 0; i < size(); ++i) V[i] = new Vector(*this);
}
void Vectors::resizeData(int N) {
......@@ -488,25 +490,25 @@ Vectors &Vectors::operator-=(const Vectors &u) {
Vectors &Vectors::operator=(const constAB<Scalar, Vectors> &au) {
for (int i = 0; i < au.second().size(); ++i)
(*V[i]).Multiply(au.first(),au.second()[i]);
(*V[i]).Multiply(au.first(), au.second()[i]);
return *this;
}
Vectors &Vectors::operator+=(const constAB<Scalar, Vectors> &au) {
for (int i = 0; i < au.second().size(); ++i)
(*V[i]).MultiplyPlus(au.first(),au.second()[i]);
(*V[i]).MultiplyPlus(au.first(), au.second()[i]);
return *this;
}
Vectors &Vectors::operator-=(const constAB<Scalar, Vectors> &au) {
for (int i = 0; i < au.second().size(); ++i)
(*V[i]).MultiplyMinus(au.first(),au.second()[i]);
(*V[i]).MultiplyMinus(au.first(), au.second()[i]);
return *this;
}
Vectors &Vectors::operator=(const constAB<int, Vectors> &au) {
for (int i = 0; i < au.second().size(); ++i)
(*V[i]).Multiply(au.first(),au.second()[i]);
(*V[i]).Multiply(au.first(), au.second()[i]);
return *this;
}
......@@ -546,14 +548,13 @@ void Vectors::ClearDirichletValues() {
for (int n = 0; n < V.size(); ++n)
a[n] = (*V[n])();
bool *b = D();
for (int i=0; i<V[0]->size(); ++i, ++b) {
for (int i = 0; i < V[0]->size(); ++i, ++b) {
if (*b)
for (int n = 0; n < V.size(); ++n) *(a[n]) = 0;
for (int n = 0; n < V.size(); ++n) ++(a[n]);
}
}
Saver &Vectors::save(Saver &saver) const {
saver << size();
for (int i = 0; i < size(); ++i)
......@@ -801,7 +802,7 @@ Matrix::~Matrix() {
}
}
void Matrix::CreateSparse () {
void Matrix::CreateSparse() {
mout << "mem used: " << PPM->MemoryUsage() << ", before Matrix::CreateSparse" << endl;
if (sm) delete sm;
sm = new SparseMatrix(*this);
......@@ -1182,11 +1183,8 @@ constAB<Operator, Vectors> operator*(const Matrix &A, const Vectors &v) {
return constAB<Operator, Vectors>(A, v);
}
RowBndValues::RowBndValues(Vector &u, const cell &c) : BF(u.GetMesh(), c),
a(0), b(0) {
a(0), b(0) {
rows R(u, c);
a.resize(R.size());
b.resize(R.size());
......@@ -1230,7 +1228,7 @@ MixedRowBndValues::MixedRowBndValues(Vector &u, const cell &c, const rows &R)
}
MixedRowFaceValues::MixedRowFaceValues(int face, Vector &u, const cell &c, const rows &R) :
a(u.NumberOfDoFs()) {
a(u.NumberOfDoFs()) {
for (int n = 0; n < u.NumberOfDoFs(); ++n) {
const vector<NPData> &data_face_n = u.dof::NodalPointFaceData(n, c, face);
a[n].resize(data_face_n.size());
......@@ -1240,7 +1238,7 @@ a(u.NumberOfDoFs()) {
}
MixedRowBndFaceValues::MixedRowBndFaceValues(int face, Vector &u, const cell &c, const rows &R) :
a(u.NumberOfDoFs()), b(u.NumberOfDoFs()) {
a(u.NumberOfDoFs()), b(u.NumberOfDoFs()) {
for (int n = 0; n < u.NumberOfDoFs(); ++n) {
const vector<NPData> &data_face_n = u.dof::NodalPointFaceData(n, c, face);
a[n].resize(data_face_n.size());
......@@ -1327,4 +1325,25 @@ void EGRowEntries::fill(Matrix &A, const rows &R, const rows &R2) {
Scalar &EGRowEntries::operator()(int i, int j, int k, int l) {
return a[i][j][k * n[j] + l];
}
\ No newline at end of file
}
PartRowEntries::PartRowEntries(Matrix &A, const rows &R, int p) : a{} {
int q[8] = {};
for (unsigned int i = 0; i < R.size(); ++i) {
q[i] = 0;
n[i] = R[i].n();
if (n[i] > 1) q[i] = p;
}
for (unsigned int i = 0; i < R.size(); ++i)
for (unsigned int j = 0; j < R.size(); ++j)
a[i][j] = A(R[i], R[j]) + q[i] * n[j] + q[j];
}
PartRowBndValues::PartRowBndValues(Vector &u, const rows &R, const cell &c) : BF(u.GetMesh(), c), a{}, b{} {
for (unsigned int i = 0; i < R.size(); ++i) {
int q = 0;
if (R[i].n() > 1) q = DomainPart(c.Subdomain());
a[i] = u(R[i]) + q;
b[i] = u.D(R[i]) + q;
}
}
src/Algebra.h
View file @ a5aba5d3
......@@ -8,6 +8,7 @@
#include "Interface.h"
#include "utility/SaveLoad.h"
class BasicVector {
Scalar *a;
int n;
......@@ -132,7 +133,6 @@ public:
virtual void multiply_transpose_minus(Vectors &u, const Vectors &v) const;
};
class DirichletFlags {
bool *f;
int n;
......@@ -161,7 +161,6 @@ constAB<Scalar, Vector> operator*(const Scalar &, const Vector &);
constAB<int, Vector> operator*(const int &, const Vector &);
class Vector : public matrixgraph, public BasicVector, public DirichletFlags {
public:
Vector(const matrixgraph &g);
......@@ -284,9 +283,9 @@ public:
#endif
};
inline Saver &operator<<(Saver &saver, const Vector & u) { return u.save(saver); }
inline Saver &operator<<(Saver &saver, const Vector &u) { return u.save(saver); }
inline Loader &operator>>(Loader &loader, Vector & u) { return u.load(loader); }
inline Loader &operator>>(Loader &loader, Vector &u) { return u.load(loader); }
class Vectors;
......@@ -300,7 +299,7 @@ constAB<int, Vectors> operator*(const int &, const Vectors &);
class Vectors : public matrixgraph, public DirichletFlags {
protected:
vector<Vector*> V;
vector<Vector *> V;
public:
Vectors(int n, const matrixgraph &g);
......@@ -309,7 +308,9 @@ public:
Vectors(const Vectors &u);
Vectors(const constAB<Operator, Vectors> &Ov);
~Vectors ();
~Vectors();
int size() const { return V.size(); }
void resize(int N);
......@@ -325,16 +326,19 @@ public:
(V[i]) = 0;
}
const Vector& operator[] (int i) const {return *(V[i]);}
Vector& operator[] (int i) {return *(V[i]);}
const Vector &operator[](int i) const { return *(V[i]); }
Vector &operator[](int i) { return *(V[i]); }
void erase(int i) {
mout << OUT("TODO") << endl;
V.erase(V.begin() + i);
}
void reverse() {
std::reverse(V.begin(),V.end());
std::reverse(V.begin(), V.end());
}
Vectors &operator=(const Vector &u);
Vectors &operator=(const Vectors &u);
......@@ -375,29 +379,50 @@ public:
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); }
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;
......@@ -421,21 +446,23 @@ public:
Loader &load(Loader &loader);
};
inline Saver &operator<<(Saver &saver, const Vectors & u) { return u.save(saver); }
inline Saver &operator<<(Saver &saver, const Vectors &u) { return u.save(saver); }
inline Loader &operator>>(Loader &loader, Vectors & u) { return u.load(loader); }
inline Loader &operator>>(Loader &loader, Vectors &u) { return u.load(loader); }
Vector operator + (const Vector& u, const Vector& v);
Vector operator - (const Vector& u, const Vector& v);
Vector operator+(const Vector &u, const Vector &v);
Vector operator-(const Vector &u, const Vector &v);
Scalar operator*(const Vector &u, const Vector &v);
Logging &operator<<(Logging &s, const Vector &u);
double norm(const Vector &u);
double norm(const double &u);
double norm(double &u);
double norm(double &u);
#ifndef DCOMPLEX
......@@ -443,8 +470,9 @@ double norm(const Vector &u, int i);
#endif
Vectors operator + (const Vectors& u, const Vectors& v);
Vectors operator - (const Vectors& u, const Vectors& v);
Vectors operator+(const Vectors &u, const Vectors &v);
Vectors operator-(const Vectors &u, const Vectors &v);
class RowValues {
vector<Scalar *> a;
......@@ -481,7 +509,6 @@ public:
int bc(int i) const { return BF[i]; }
};
/**
* Sets the RowValues corresponding to the nodal points
* n: index of the shapes (ordered as in the discretization/dof)
......@@ -522,7 +549,6 @@ public:
int bc(int i) const { return BF[i]; }
};
class MixedRowFaceValues {
vector<vector<Scalar *>> a;
public:
......@@ -577,18 +603,25 @@ public:
};
class Matrix;
constAB<Scalar, Matrix> operator*(const Scalar &, const Matrix &);
class SparseMatrix;
class Matrix : public matrixgraph, public BasicVector, public Operator {
const Vector &u;
SparseMatrix* sm;
SparseMatrix *sm;
public:
Matrix(const Vector& U) : matrixgraph(U), BasicVector(U.Size()), u(U), sm(nullptr) {}
Matrix(const Matrix& A) : matrixgraph(A), BasicVector(A), u(A.u), sm(nullptr) {}
virtual ~Matrix ();
const SparseMatrix* GetSparse() const { return sm; }
void CreateSparse ();
Matrix(const Vector &U) : matrixgraph(U), BasicVector(U.Size()), u(U), sm(nullptr) {}
Matrix(const Matrix &A) : matrixgraph(A), BasicVector(A), u(A.u), sm(nullptr) {}
virtual ~Matrix();
const SparseMatrix *GetSparse() const { return sm; }
void CreateSparse();
Matrix &operator=(Scalar b) {
BasicVector::operator=(b);
return *this;
......@@ -660,7 +693,7 @@ public:
RowEntries(Matrix &A, const cell &c) : RowEntries(A, rows(A, c)) {}
RowEntries(Matrix &A, const tcell &tc) : RowEntries(A, rows(A, tc)) {}
RowEntries(Matrix &A, const tcell &tc) : RowEntries(A, rows(A, tc)) {}
Scalar &operator()(int i, int j, int k = 0, int l = 0) {
return a[i][j][k * n[j] + l];
......@@ -680,4 +713,34 @@ public:
}
};
class PartRowBndValues {
BFParts BF;
Scalar *a[MaxNodalPoints];
bool *b[MaxNodalPoints];
public:
PartRowBndValues(Vector &u, const rows &R, 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]; }
};
class PartRowEntries {
Scalar *a[MaxNodalPoints][MaxNodalPoints];
int n[MaxNodalPoints];
public:
PartRowEntries(Matrix &A, const rows &R, int p);
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_
src/CMakeLists.txt
View file @ a5aba5d3
......@@ -42,7 +42,7 @@ set(mpp_src ${src} ${utility} ${parallel} ${cell}
if (USE_CXSC)
include(intervalarithmetic/CMakeLists.txt)
add_library(SRC STATIC ${mpp_src} ${intervalarithmetic})
elseif(USE_SPACETIME)
elseif (USE_SPACETIME)
add_library(SRC STATIC ${mpp_src})
else ()
add_library(SRC STATIC ${mpp_src})
......
src/Distribution.C
View file @ a5aba5d3
......@@ -6,6 +6,7 @@
#include <list>
#include <map>
extern int TimeLevel;
using std::list;
......@@ -210,6 +211,7 @@ void FullOverlap(Mesh &M) {
E.ClearBuffers();
pout << "Faces Size " << M.FaceCount() << endl;
}
/*
class CellBoundaryFaces {
bnd_face bf[6];
......@@ -837,7 +839,7 @@ void OverlapLayer(Mesh &M, Vertices &bndVertices) {
check_c = 1;
for (int fe = 0; fe < c.FaceEdges(fi); ++fe) {
auto v = bndVertices.Find(c.FaceEdge(fi, fe));
if(v == bndVertices.End())
if (v == bndVertices.End())
bndVertices.Insert(c.FaceEdge(fi, fe), new Vertex(0));
else
v->second->inc();
......@@ -1053,7 +1055,7 @@ void SlipFacesOverlap_cdd2d(Mesh &M, vector<Point> &N, int level) {
// }
for (int j = 0; j < c.FaceEdges(i); ++j) {
//if(M.FE[ns].find_faceedge(c.FaceEdge(i,j)) == M.FE[ns].faceedges_end()) continue;
InsertFaceEdgeIn(M.OverlapFE[ns],c.FaceEdge(i, j), c.Face(i));
InsertFaceEdgeIn(M.OverlapFE[ns], c.FaceEdge(i, j), c.Face(i));
//M.ProcSets::Copy(p,c.FaceEdge(i,j));
}
}
......@@ -1142,8 +1144,7 @@ void Distribute(Mesh &M, const string &name) {
++k;
D.MarkCell(C[i], q);
}
}
else if (name == "RCB" || name == "RCBx" || name == "RCBy" || name == "RCBz" ||
} else if (name == "RCB" || name == "RCBx" || name == "RCBy" || name == "RCBz" ||
name == "RCB2D" ||
name == "RCB2Dx" || name == "RCB2Dy") {
int N = M.CellCount();
......@@ -1196,8 +1197,7 @@ void Distribute(Mesh &M, const string &name) {
for (int i = 0; i < C.size(); ++i) {
D.MarkCell(C[i], Dest[i]);
}
}
else Exit(name + " not implemented");
} else Exit(name + " not implemented");
D.Communicate();
int cellsize = M.CellCount();
int mincells = PPM->Min(cellsize);
......
src/Distribution.h
View file @ a5aba5d3
......@@ -59,7 +59,8 @@ inline Buffer &operator<<(Buffer &b, const DataContainer &d) {
}
inline Buffer &operator>>(Buffer &b, DataContainer &d) {
short m; b >> m;
short m;
b >> m;
d = DataContainer(m);
for (int i = 0; i < d.size(); ++i)
b >> d.data[i];
......@@ -76,7 +77,7 @@ inline Buffer &operator>>(Buffer &b, Mesh &M) {
for (int i = 0; i < n; ++i) b >> x[i];
if (M.StoresData()) {
vector<DataContainer> D(n);
for (int i = 0; i < n; ++i){
for (int i = 0; i < n; ++i) {
b >> D[i];
}
M.InsertCell(CELLTYPE(tp), sd, x, D);
......@@ -115,7 +116,6 @@ inline Buffer &operator<<(Buffer &b, const CellBoundaryFaces &bf) {
return b;
}
void Distribute(Mesh &, const string &);
void DistributeOverlap(Mesh &, const string &);
......@@ -128,7 +128,6 @@ void SlipFacesOverlap_cdd(Mesh &M);
void SlipFacesOverlap_cdd2d(Mesh &M, vector<Point> &N, int level);
inline bool Less_x(const cell &c0, const cell &c1) {
const Point &P = c0();
const Point &Q = c1();
......@@ -322,5 +321,4 @@ rcb_geom(int level, vector<int> &D, int begin, int end, vector<cell> &C, char d)
}
}
#endif // of #ifndef _DISTRIBUTION_H_
src/Euler.C
View file @ a5aba5d3
#include "Euler.h"
Euler::Euler(const char *conf, NonlinearSolver &n) :
N(n), verbose(1), extrapolate(1), max_estep(100000),
ReferenceCalculation(0), check_conv(true), SSum(0) {
N(n), verbose(1), extrapolate(1), max_estep(100000),
ReferenceCalculation(0), check_conv(true), SSum(0) {
config.get("EulerVerbose", verbose);
config.get("EulerExtrapolate", extrapolate);
config.get("ReferenceCalculation", ReferenceCalculation);
......@@ -12,8 +13,8 @@ Euler::Euler(const char *conf, NonlinearSolver &n) :
}
Euler::Euler(NonlinearSolver &n) :
N(n), verbose(1), extrapolate(1), max_estep(100000),
ReferenceCalculation(0), check_conv(true), SSum(0) {
N(n), verbose(1), extrapolate(1), max_estep(100000),
ReferenceCalculation(0), check_conv(true), SSum(0) {
config.get("EulerVerbose", verbose);
config.get("EulerExtrapolate", extrapolate);
config.get("ReferenceCalculation", ReferenceCalculation);
......@@ -181,7 +182,7 @@ void Euler::operator()(Mesh &M, TAssemble &A, Vector &u) {
}
LinearEuler::LinearEuler(const char *conf, Solver &s) :
S(s), verbose(1), ReferenceCalculation(0), max_estep(100000), check_conv(true) {
S(s), verbose(1), ReferenceCalculation(0), max_estep(100000), check_conv(true) {
config.get("EulerVerbose", verbose);
config.get("ReferenceCalculation", ReferenceCalculation);
config.get("EulerSteps", max_estep);
......
src/Euler.h
View file @ a5aba5d3
......@@ -4,6 +4,7 @@
#include "solver/Newton.h"
#include "TimeSeries.h"
class Euler {
NonlinearSolver &N;
int verbose;
......
src/GlobalDefinitions.hpp
View file @ a5aba5d3
......@@ -4,6 +4,7 @@
#include <string>
#include <vector>
using std::string;
using std::vector;
......@@ -12,6 +13,8 @@ using std::max;
using std::to_string;
#include <ext/numeric>
using __gnu_cxx::power;
#ifdef GCC_4_4
......@@ -19,9 +22,12 @@ using __gnu_cxx::power;
#include <cstdio>
#include <cstring>
#endif
#include <unordered_map>
#define hash_map std::unordered_map
typedef double Scalar;
......@@ -42,12 +48,14 @@ inline double eval(const Scalar &z) { return z; }
const Scalar iUnit(-1);
constexpr double Eps = 1e-10;
constexpr double Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062;
constexpr int MaxPointTypes = 10;
constexpr int MaxQuadraturePoints = 216;
constexpr int MaxNodalPoints = 216;
#endif // GLOBALDEFINITIONS
src/Identify.C
View file @ a5aba5d3
#include "Identify.h"
#include "parallel/Parallel.h"
void IdentifySet::Add(const Point &x) {
long unsigned int m = size();
if (m == 0) {
......@@ -22,7 +23,6 @@ Logging &operator<<(Logging &s, const IdentifySet &I) {
return s;
}
Logging &operator<<(Logging &s, const identifyset &I) {
return s << I() << " : " << *I << endl;
}
......
src/Identify.h
View file @ a5aba5d3
......@@ -5,6 +5,7 @@
#include "geometry/Point.h"
class IdentifySet : public std::vector<Point> {
public:
IdentifySet() {};
......
src/Interface.C
View file @ a5aba5d3
......@@ -32,7 +32,7 @@ void DirichletConsistent(Vector &u) {
<< endl;
E.Receive(q) >> z;
int i = u.Idx(z);
if (i == -1) pout << OUT(z) << " not found \n";
if (i == -1) pout << OUT(z) << " not found \n";
assert(i != -1);
for (int k = 0; k < u.Dof(i); ++k) {
Scalar a;
......
src/Lapdef.h
View file @ a5aba5d3
......@@ -84,7 +84,6 @@
#ifdef LAPACK
extern "C" void dsytri_(char *UPLO, int *N, void *A, int *LDA, int *IPIV, void *WORK, int *INFO);
extern "C" void dsytrf_(char *UPLO, int *N, void *A, int *LDA, int *IPIV, void *WORK, void *LWORK, int *INFO);
......@@ -224,7 +223,6 @@ extern "C" void
dsyev_(char *JOBZ, char *UPLO, int *N, void *A, int *LDA, void *W, void *WORK, int *LWORK, int *INFO,
int ch1, int ch2);
#else // of #ifdef LAPACK
void dgetrf_(int *M, int *N, void *A, int *LDA, int *IPIV, int *INFO) {
......
src/Plot.C
View file @ a5aba5d3
......@@ -6,6 +6,7 @@
#include "parallel/Parallel.h"
#include "utility/M_IOFiles.hpp"
using namespace std;
string DataPathName("data");
......@@ -23,9 +24,9 @@ ostream &operator<<(ostream &s, const PlotData &p) {
ostream &operator<<(ostream &s, const plotdata &p) {
s << p[0];
for (int i = 1;i<SpaceDimension; ++i)
for (int i = 1; i < SpaceDimension; ++i)
s << " " << p[i];
for(int i=SpaceDimension; i<3;++i)
for (int i = SpaceDimension; i < 3; ++i)
s << " 0.0";
for (int j = 0; j < p->second.size(); ++j) s << " " << p[j];
return s << " 0 0 0 " << endl;
......@@ -221,7 +222,6 @@ void PlotCellDatas::singledata(const Mesh &M, const D &d, int shift) {
}
}
ostream &operator<<(ostream &s, const Plot &plot) {
for (plotdata p = plot.VD.plotdatas(); p != plot.VD.plotdatas_end(); ++p)
s << p;
......@@ -346,7 +346,6 @@ inline int Less0(const void *z0, const void *z1) {
return ((*z_0)[0] > (*z_1)[0]);
}
void Plot::gnu_celldata(const char *name, int i) {
if (!PPM->master()) return;
if (M.dim() > 1) return;
......@@ -366,7 +365,6 @@ void Plot::gnu_celldata(const char *name, int i) {
delete[] zd;
}
void Plot::gnu_deformation(const char *name) {
if (!PPM->master()) return;
string filename = DataPathName + string("/gp/") + name + string(".gp");
......@@ -453,12 +451,10 @@ void Plot::vtk_mesh(ostream &out, int deform) {
out << "CELLS " << C.size() << " " << (C[0].size() + 1) * C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
out << C[i].size() << " ";
switch (C[i].size()) {
case 2:
out << C[i][0] << " " << C[i][1];
switch(C[i].size()) {
case 2:out << C[i][0] << " " << C[i][1];
break;
case 3:
out << C[i][0] << " " << C[i][1] << " " << C[i][2];
case 3:out << C[i][0] << " " << C[i][1] << " " << C[i][2];
break;
case 4:
if (M.dim() == 2)
......@@ -475,19 +471,16 @@ void Plot::vtk_mesh(ostream &out, int deform) {
}
out << "CELL_TYPES " << C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
switch (C[i].size()) {
case 2:
out << "3";
switch(C[i].size()) {
case 2:out << "3";
break;
case 3:
out << "5";
case 3:out << "5";
break;
case 4:
if (M.dim() == 2) out << "8";
else out << "10";
break;
case 8:
out << "11";
case 8:out << "11";
break;
}
}
......@@ -510,30 +503,25 @@ void Plot::vtk_2d_graph(ostream &out, int k) {
out << "CELLS " << C.size() << " " << (C[0].size() + 1) * C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
out << C[i].size() << " ";
switch (C[i].size()) {
case 3:
out << C[i][0] << " " << C[i][1] << " " << C[i][2];
switch(C[i].size()) {
case 3:out << C[i][0] << " " << C[i][1] << " " << C[i][2];
break;
case 4:
out << C[i][0] << " " << C[i][1] << " " << C[i][3] << " " << C[i][2];
case 4:out << C[i][0] << " " << C[i][1] << " " << C[i][3] << " " << C[i][2];
break;
}
}
out << "CELL_TYPES " << C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
switch (C[i].size()) {
case 2:
out << "3";
switch(C[i].size()) {
case 2:out << "3";
break;
case 3:
out << "5";
case 3:out << "5";
break;
case 4:
if (M.dim() == 2) out << "8";
else out << "10";
break;
case 8:
out << "11";
case 8:out << "11";
break;
}
}
......@@ -568,9 +556,8 @@ void Plot::dx_mesh(ostream &out, bool deformed) {
<< C[0].size() << " items "
<< C.size() << " data follows" << endl;
for (int i = 0; i < C.size(); ++i, out << endl)
switch (C[i].size()) {
case 3:
out << C[i][0] << " " << C[i][1] << " " << C[i][2];
switch(C[i].size()) {
case 3:out << C[i][0] << " " << C[i][1] << " " << C[i][2];
break;
case 4:
if (M.dim() == 2)
......@@ -646,7 +633,6 @@ void Plot::vtk_celltensor(ostream &out, int k) {
}
}
void Plot::vtk_scalar(ostream &out, int k) {
out << "SCALARS scalar_value float 1" << endl
<< "LOOKUP_TABLE default" << endl;
......@@ -664,9 +650,7 @@ void Plot::vtk_point_data(ostream &out) {
out << "POINT_DATA " << VD.size() << endl;
}
void Plot::dualgraph(ostream &out) {
}
void Plot::vtk_vertex_vector(const char *name, int k, int deformed) {
......@@ -678,7 +662,7 @@ void Plot::vtk_vertex_vector(const char *name, int k, int deformed) {
vtk_vector(out, k);
}
void Plot::vtk_vertexvector(const string& name, int k, int deformed) {
void Plot::vtk_vertexvector(const string &name, int k, int deformed) {
if (!PPM->master()) return;
string filename = DataPathName + string("/vtk/") + name + string(".vtk");
M_ofstream out(filename.c_str());
......@@ -696,7 +680,7 @@ void Plot::vtk_vertexdata(const char *name, int k, bool deformed) {
vtk_scalar(out, k);
}
void Plot::vtk_vertexdata(const string& name, int k, bool deformed) {
void Plot::vtk_vertexdata(const string &name, int k, bool deformed) {
if (!PPM->master()) return;
string filename = DataPathName + string("/vtk/") + name + string(".vtk");
M_ofstream out(filename.c_str());
......@@ -733,7 +717,7 @@ void Plot::vtk_celldata(const char *name, int k, bool deformed, const char *plot
vtk_celldata(out, k, deformed, plotname);
}
void Plot::vtk_celldata(const string& name, int k, bool deformed, const char *plotname) {
void Plot::vtk_celldata(const string &name, int k, bool deformed, const char *plotname) {
if (!PPM->master()) return;
string filename = DataPathName + string("/vtk/") + name + string(".vtk");
M_ofstream out(filename.c_str());
......@@ -777,7 +761,7 @@ void Plot::vtk_cellvector(const char *name, int k, bool deformed) {
vtk_cellvector(out, k);
}
void Plot::vtk_cellvector(const string& name, int k, bool deformed) {
void Plot::vtk_cellvector(const string &name, int k, bool deformed) {
if (!PPM->master()) return;
string filename = DataPathName + string("/vtk/") + name + string(".vtk");
M_ofstream out(filename.c_str());
......@@ -1285,9 +1269,8 @@ void TPlot::vtk_mesh(ostream &out, int deform) {
out << "CELLS " << C.size() << " " << (C[0].size() + 1) * C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
out << C[i].size() << " ";
switch (C[i].size()) {
case 4:
out << C[i][0] << " " << C[i][1] << " " << C[i][2] << " " << C[i][3];
switch(C[i].size()) {
case 4:out << C[i][0] << " " << C[i][1] << " " << C[i][2] << " " << C[i][3];
break;
case 6:
out << C[i][0] << " " << C[i][1] << " " << C[i][2] << " " << C[i][3] << " " << C[i][4] << " "
......@@ -1301,15 +1284,12 @@ void TPlot::vtk_mesh(ostream &out, int deform) {
}
out << "CELL_TYPES " << C.size() << endl;
for (int i = 0; i < C.size(); ++i, out << endl) {
switch (C[i].size()) {
case 4:
out << "8";
switch(C[i].size()) {
case 4:out << "8";
break;
case 6:
out << "13";
case 6:out << "13";
break;
case 8:
out << "11";
case 8:out << "11";
break;
}
}
......
src/Plot.h
View file @ a5aba5d3
......@@ -8,6 +8,7 @@
#include "mesh/Mesh.h"
#include "mesh/SpaceTimeMesh.hpp"
class PlotData : public vector<double> {
int i;
int c;
......@@ -82,7 +83,6 @@ public:
void data(const Mesh &M, const D &d, int shift = 0);
PlotVertexDatas(const Mesh &M, int N);
};
class PlotCell : public vector<int> {
......@@ -119,14 +119,14 @@ class Plot {
int vtkplot;
double PlotTolerance;
void plotPoint(std::ostream &os, const Point& p){
void plotPoint(std::ostream &os, const Point &p) {
os << p[0];
for (int i = 1;i<SpaceDimension; ++i)
for (int i = 1; i < SpaceDimension; ++i)
if (abs(p[i]) < GeometricTolerance)
os << " 0.0";
else
os << " " << p[i];
for(int i=SpaceDimension; i<3;++i)
for (int i = SpaceDimension; i < 3; ++i)
os << " 0.0";
}
......@@ -241,9 +241,9 @@ public:
const char *plotname = "scalar_value");
void vtk_celldata(const string &name,
int k = 0,
bool deformed = false,
const char *plotname = "scalar_value");
int k = 0,
bool deformed = false,
const char *plotname = "scalar_value");
void vtk_celldata(std::ostream &out,
int k = 0,
......@@ -360,16 +360,17 @@ class TPlot {
int vtkplot;
double PlotTolerance;
protected:
void plotPoint(std::ostream &os, const Point& p){
void plotPoint(std::ostream &os, const Point &p) {
os << p[0];
for (int i = 1;i<SpaceDimension; ++i)
for (int i = 1; i < SpaceDimension; ++i)
if (abs(p[i]) < GeometricTolerance)
os << " 0.0";
else
os << " " << p[i];
for(int i=SpaceDimension; i<3;++i)
for (int i = SpaceDimension; i < 3; ++i)
os << " 0.0";
}
public:
const SpaceTimeMesh &STM;
TPlotVertexDatas VD;
......
src/Quadrature.C
View file @ a5aba5d3
#include "Quadrature.h"
#include <unordered_map>
template<>
QpointT<double, SpaceDimension, TimeDimension>::QpointT() : Quadrature("Qpoint", 1) {
z[0] = Point(0.0, 0.0, 0.0);
......@@ -125,8 +126,7 @@ QintT<double, SpaceDimension, TimeDimension>::QintT(int exactUpTo)
z[8] = Point(0.984080119753813, 0.0, 0.0);
break;
case 18:
case 19:
w[0] = 0.033335672154344;
case 19:w[0] = 0.033335672154344;
z[0] = Point(0.0130467357414141, 0.0, 0.0);
w[1] = 0.0747256745752903;
z[1] = Point(0.0674683166555077, 0.0, 0.0);
......@@ -1117,8 +1117,6 @@ const QuadratureT<> &GetQuadratureT(const CELLTYPE cellType,
Exit("Quadrature not implemented!");
}
const Quadrature &GetQuadrature(const CELLTYPE cellType,
int exactUpTo,
int exactUpTo1,
......@@ -1225,7 +1223,6 @@ const QuadratureT<T, sDim, tDim> &GetQuadrature(const std::string &name) {
const Quadrature &GetQuadrature(const std::string &name) {
return GetQuadrature<double, SpaceDimension, TimeDimension>(name);
}
void clearQuadrature() {
......
src/Quadrature.h
View file @ a5aba5d3
......@@ -128,8 +128,8 @@ const Quadrature &GetQuadrature(const CELLTYPE,
int exactUpTo1 = -1,
int exactUpTo2 = -1);
template <class T, int sDim,int tDim>
const QuadratureT<T,sDim, tDim> &GetQuadrature(const std::string &name);
template<class T, int sDim, int tDim>
const QuadratureT<T, sDim, tDim> &GetQuadrature(const std::string &name);
const Quadrature &GetQuadrature(const std::string &);
......
src/Small.C
View file @ a5aba5d3
<
......@@ -3,6 +3,7 @@
#include "Lapdef.h"
#include "spectrum/Spectrum.hpp"
void PrintScalar(Logging &s, const Scalar &a) {
#ifdef NDOUBLE
s << " " << a;
......@@ -46,7 +47,7 @@ SmallMatrix::SmallMatrix(int N, int M, const SmallMatrix &A) : n(N), m(M) {
}
SmallMatrix::SmallMatrix(const SmallMatrix &A, const SmallMatrix &B) :
n(A.n), m(B.m) {
n(A.n), m(B.m) {
a = new Scalar[n * m];
int rows = min(A.m, B.n);
for (int i = 0; i < n; ++i)
......@@ -59,7 +60,7 @@ SmallMatrix::SmallMatrix(const SmallMatrix &A, const SmallMatrix &B) :
SmallMatrix::SmallMatrix(int _n, int _m,
const SmallMatrix &A, const SmallMatrix &B) :
n(_n), m(_m) {
n(_n), m(_m) {
a = new Scalar[n * m];
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
......@@ -70,7 +71,7 @@ SmallMatrix::SmallMatrix(int _n, int _m,
}
SmallMatrix::SmallMatrix(const SparseMatrix &M, const vector<int> &ind) :
n(int(ind.size())), m(int(ind.size())) {
n(int(ind.size())), m(int(ind.size())) {
a = new Scalar[n * m];
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
......@@ -136,7 +137,6 @@ SmallMatrix::SmallMatrix(const Matrix &A) : n(A.size()), m(A.size()) {
}
}
SmallMatrix::SmallMatrix(const SmallMatrices &A) : n(A.rows()), m(A.cols()) {
a = new Scalar[n * m];
int j = 0;
......@@ -148,7 +148,7 @@ SmallMatrix::SmallMatrix(const SmallMatrices &A) : n(A.rows()), m(A.cols()) {
}
AdjointSmallMatrix::AdjointSmallMatrix(const SmallMatrices &A)
: SmallMatrix(A.cols(), A.rows()) {
: SmallMatrix(A.cols(), A.rows()) {
int j = 0;
for (int s = 0; s < A.size(); ++s) {
for (int k = 0; k < A.cols(s); ++j, ++k)
......@@ -160,7 +160,7 @@ AdjointSmallMatrix::AdjointSmallMatrix(const SmallMatrices &A)
AdjointSmallMatrix::AdjointSmallMatrix(int n, int m,
const SmallMatrix &A,
const SmallMatrix &B)
: SmallMatrix(n, m) {
: SmallMatrix(n, m) {
int K = min(A.rows(), B.rows());
for (int i = 0; i < rows(); ++i)
for (int j = 0; j < cols(); ++j) {
......@@ -704,7 +704,7 @@ SmallVector::SmallVector(const Vector &u) : std::valarray<Scalar>(u.size()) {
}
SmallVector::SmallVector(const SmallMatrix &A, const SmallVector &u)
: std::valarray<Scalar>(A.rows()) {
: std::valarray<Scalar>(A.rows()) {
for (int j = 0; j < size(); ++j) {
Scalar s = 0;
for (int k = 0; k < u.size(); ++k)
......@@ -714,7 +714,7 @@ SmallVector::SmallVector(const SmallMatrix &A, const SmallVector &u)
}
SmallVector::SmallVector(const SmallVector &u, const SmallMatrix &A)
: std::valarray<Scalar>(std::min(int(u.size()), A.cols())) {
: std::valarray<Scalar>(std::min(int(u.size()), A.cols())) {