Commits (12)
......@@ -37,6 +37,10 @@ public:
int Degree() const {
return degree;
}
int Size() const {
return size;
}
};
using LagrangeDiscretization = LagrangeDiscretizationT<>;
......
......@@ -40,6 +40,8 @@ public:
degree(degree), size(size) {}
int Degree() const { return degree; }
int Size() const { return size; }
};
using MultiPartDiscretization = MultiPartDiscretizationT<>;
......
......@@ -7,6 +7,41 @@ ArgyrisBasicElementT<TT, sDim, tDim>::ArgyrisBasicElementT(const VectorMatrixBas
init();
}
template<typename TT, int sDim, int tDim>
std::vector<TT>
ArgyrisBasicElementT<TT, sDim, tDim>::values(const PointT<TT, sDim, tDim> &z) const {
std::vector<TT> v(21);
for (int j = 0; j < 21; ++j) v[j] = this->shape(z, j);
std::vector<TT> V(21);
for (int j = 0; j < 21; ++j)
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, V[j], v, j);
return V;
}
template<typename TT, int sDim, int tDim>
std::vector<VectorFieldT<TT, sDim>>
ArgyrisBasicElementT<TT, sDim, tDim>::derivatives(const PointT<TT, sDim, tDim> &z) const {
vector<VectorFieldT<TT, sDim>> G(21);
for (int j = 0; j < 21; ++j)
G[j] = this->GetTransformation(0) * this->shape.LocalGradient(z, j);
std::vector<VectorFieldT<TT, sDim>> D(21);
for (int j = 0; j < 21; ++j)
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, D[j], G, j);
return D;
}
template<typename TT, int sDim, int tDim>
std::vector<SymTensorT<TT, sDim>>
ArgyrisBasicElementT<TT, sDim, tDim>::hessians(const PointT<TT, sDim, tDim> &z) const {
vector<SymTensorT<TT, sDim>> H(21);
for (int j = 0; j < 21; ++j)
H[j] = this->multiplyTheta(this->shape.LocalHessian(z, j));
vector<SymTensorT<TT, sDim>> D(21);
for (int j = 0; j < 21; ++j)
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, D[j], H, j);
return D;
}
template<typename TT, int sDim, int tDim>
int ArgyrisBasicElementT<TT, sDim, tDim>::get_maxk(int i) const {
if (i < 3) return 6;
......@@ -222,16 +257,11 @@ TT ArgyrisElementT<TT, sDim, tDim>::Value(int q, const Vector &u, int m) const {
template<typename TT, int sDim, int tDim>
TT ArgyrisElementT<TT, sDim, tDim>::Value(const PointT<TT, sDim, tDim> &z, const Vector &u,
int m) const {
vector<TT> V(21);
for (int j = 0; j < 21; ++j) V[j] = this->shape(z, j);
TT U = TT(0.0);
vector<TT> v_z = this->values(z);
TT U{};
for (int i = 0; i < this->size(); ++i)
for (int k = 0; k < this->get_maxk(i); ++k) {
TT s = 0.0;
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, s, V,
this->indexingShape(i, k));
U += u(this->r(i), this->indexing_k(i, k, m)) * s;
}
for (int k = 0; k < this->get_maxk(i); ++k)
U += u(this->r(i), this->indexing_k(i, k, m)) * v_z[this->indexingShape(i, k)];
return U;
}
......@@ -263,17 +293,12 @@ template<typename TT, int sDim, int tDim>
VectorFieldT<TT, sDim> ArgyrisElementT<TT, sDim, tDim>::Derivative(const PointT<TT, sDim, tDim> &z,
const Vector &u,
int m) const {
vector<VectorFieldT<TT, sDim>> G(21);
for (int j = 0; j < 21; ++j)
G[j] = this->GetTransformation(0) * this->shape.LocalGradient(z, j);
vector<VectorFieldT<TT, sDim>> d_z = this->derivatives(z);
VectorFieldT<TT, sDim> D;
for (int i = 0; i < this->size(); ++i)
for (int k = 0; k < this->get_maxk(i); ++k) {
VectorFieldT<TT, sDim> v;
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, v, G,
this->indexingShape(i, k));
D += u(this->r(i), this->indexing_k(i, k, m)) * v;
}
for (int k = 0; k < this->get_maxk(i); ++k)
D += u(this->r(i), this->indexing_k(i, k, m)) * d_z[this->indexingShape(i, k)];
return D;
}
......@@ -302,18 +327,12 @@ template<typename TT, int sDim, int tDim>
SymTensorT<TT, sDim>
ArgyrisElementT<TT, sDim, tDim>::Hessian(const PointT<TT, sDim, tDim> &z, const Vector &u,
int m) const {
vector<SymTensorT<TT, sDim>> H(21);
for (int j = 0; j < 21; ++j)
H[j] = this->multiplyTheta(this->shape.LocalHessian(z, j));
SymTensorT<TT, sDim> D;
SymTensorT<TT, sDim> H;
vector<SymTensorT<TT, sDim>> h_z = this->hessians(z);
for (int i = 0; i < this->size(); ++i)
for (int k = 0; k < this->get_maxk(i); ++k) {
SymTensorT<TT, sDim> v;
applyC(this->C0, this->C1, this->C2, this->C3, this->C4, this->C5, v, H,
this->indexingShape(i, k));
D += u(this->r(i), this->indexing_k(i, k, m)) * v;
}
return D;
for (int k = 0; k < this->get_maxk(i); ++k)
H += u(this->r(i), this->indexing_k(i, k, m)) * h_z[this->indexingShape(i, k)];
return H;
}
template<typename TT, int sDim, int tDim>
......@@ -347,7 +366,6 @@ TT ArgyrisElementT<TT, sDim, tDim>::Laplace(const PointT<TT, sDim, tDim> &z, con
}
template
class ArgyrisBasicElementT<>;
......
......@@ -24,6 +24,12 @@ protected:
ArgyrisBasicElementT(const VectorMatrixBase &base, const cell &c);
std::vector<TT> values(const PointT<TT, sDim, tDim> &z) const;
std::vector<VectorFieldT<TT, sDim>> derivatives(const PointT<TT, sDim, tDim> &z) const;
std::vector<SymTensorT<TT, sDim>> hessians(const PointT<TT, sDim, tDim> &z) const;
public:
PointT<TT, sDim, tDim> Normal(int i) const { return sign[i] * normal[i]; }
......
#include "DivergenceFreeElement.hpp"
template<typename TT, int sDim, int tDim>
DivergenceFreeElementT<TT, sDim, tDim>::DivergenceFreeElementT(
const VectorMatrixBase &base,
......@@ -59,10 +58,13 @@ VelocityT<TT, sDim>
DivergenceFreeElementT<TT, sDim, tDim>::VelocityField(const PointT<TT, sDim, tDim> &z,
const Vector &u,
int m) const {
std::vector<VectorFieldT<TT, sDim>> d_z = this->derivatives(z);
VectorFieldT<TT, sDim> D;
for (int i = 0; i < this->size(); ++i)
for (int k = 0; k < this->get_maxk(i); ++k)
D += u(this->r(i), this->indexing_k(i, k, m)) * VelocityField(z, i, k);
D += u(this->r(i), this->indexing_k(i, k, m)) *
VelocityT<TT, sDim>(-d_z[this->indexingShape(i, k)][1],
d_z[this->indexingShape(i, k)][0]);
return D;
}
......@@ -107,10 +109,15 @@ VelocityGradientT<TT, sDim>
DivergenceFreeElementT<TT, sDim, tDim>::VelocityFieldGradient(const PointT<TT, sDim, tDim> &z,
const Vector &u,
int m) const {
vector<SymTensorT<TT, sDim>> h_z = this->hessians(z);
VelocityGradientT<TT, sDim> H;
for (int i = 0; i < this->size(); ++i)
for (int k = 0; k < this->get_maxk(i); ++k)
H += u(this->r(i), this->indexing_k(i, k, m)) * VelocityFieldGradient(z, i, k);
for (int k = 0; k < this->get_maxk(i); ++k) {
const SymTensorT<TT, sDim> &h_z_ik = h_z[this->indexingShape(i, k)];
H += u(this->r(i), this->indexing_k(i, k, m)) *
VelocityGradientT<TT, sDim>(-h_z_ik(1, 0), -h_z_ik(1, 1),
h_z_ik(0, 0), h_z_ik(0, 1));
}
return H;
}
......
......@@ -39,7 +39,7 @@ public:
int m = 0) const;
};
using DivergenceFreeElement = DivergenceFreeElementT<> ;
using DivergenceFreeElement = DivergenceFreeElementT<>;
#ifdef BUILD_IA
......
#include "VtuPlot.hpp"
#include "Parallel.hpp"
#include "DivergenceFreeElement.hpp"
using namespace std;
SerialVtuPlot::SerialVtuPlot() : VtuPlot() {
......@@ -88,34 +90,37 @@ void SerialVtuPlot::initCells(const Mesh &init) {
}
}
void SerialVtuPlot::addPointData(const string& name, const Vector& data, const std::vector<int> &indices){
void SerialVtuPlot::addPointData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
int dataDim = indices.size()>1? 3 : 1;
ExchangeBuffer exBuffer(data.CommSplit());
for(row r = data.rows(); r!= data.rows_end(); ++r){
if(plotMesh().find_vertex(r())!=plotMesh().vertices_end()){
for (row r = data.rows(); r != data.rows_end(); ++r) {
if (plotMesh().find_vertex(r()) != plotMesh().vertices_end()) {
exBuffer.Send(0) << r();
for (int j : indices)
exBuffer.Send(0) << data(r, j);
}
}
exBuffer.Communicate();
communicatePointData(name, exBuffer, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::communicatePointData(const string &name, ExchangeBuffer &exBuffer,
int dataDim, int numIndices) {
exBuffer.Communicate();
auto newPointData = std::make_unique<std::unordered_map<Point, PointData>>();
for (short q = 0; q < PPM->Size(data.CommSplit()); ++q)
for (short q = 0; q < PPM->Size(commSplit()); ++q)
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
PointData newData(dataDim);
for (int i = 0; i < indices.size(); i++) {
for (int i = 0; i < numIndices; i++) {
double a{};
exBuffer.Receive(q) >> a;
newData.AddData(i, a);
}
for (int i = indices.size(); i < dataDim; i++) {
for (int i = numIndices; i < dataDim; i++) {
newData.AddData(i, 0.0);
}
......@@ -127,36 +132,37 @@ void SerialVtuPlot::addPointData(const string& name, const Vector& data, const s
if (activeArrays.first.empty()) activeArrays.first = name;
}
//int SerialVtuPlot::dataDim(int dim) const { return (dim > 1 ? 3 : 1); }
void SerialVtuPlot::addCellData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
int dataDim = indices.size()>1? 3 : 1;
ExchangeBuffer exBuffer(data.CommSplit());
for(row r = data.rows(); r!= data.rows_end(); ++r){
if(plotMesh().find_cell(r())!=plotMesh().cells_end()){
for (row r = data.rows(); r != data.rows_end(); ++r) {
if (plotMesh().find_cell(r()) != plotMesh().cells_end()) {
exBuffer.Send(0) << r();
for (int j : indices)
exBuffer.Send(0) << data(r, j);
}
}
exBuffer.Communicate();
communicateCellData(name, exBuffer, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::communicateCellData(const string &name, ExchangeBuffer &exBuffer,
int dataDim, int numIndices) {
exBuffer.Communicate();
auto newCellData = std::make_unique<std::unordered_map<Point, PointData>>();
for (short q = 0; q < PPM->Size(data.CommSplit()); ++q) {
for (short q = 0; q < PPM->Size(commSplit()); ++q) {
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
PointData newData(dataDim);
for (int i = 0; i < indices.size(); i++) {
for (int i = 0; i < numIndices; i++) {
double a;
exBuffer.Receive(q) >> a;
newData.AddData(i, a);
}
for(int i = indices.size(); i < dataDim; i++){
for (int i = numIndices; i < dataDim; i++) {
newData.AddData(i, 0.0);
}
newCellData->insert({z, newData});
......@@ -169,13 +175,62 @@ void SerialVtuPlot::addCellData(const string &name, const Vector &data,
if (activeArrays.second.empty()) activeArrays.second = name;
}
void SerialVtuPlot::addArgyrisData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
std::unordered_map<Point, std::vector<double>> pData{};
for (cell c = data.cells(); c != data.cells_end(); ++c) {
ArgyrisElement E(data, c);
for (int i = 0; i < c.Corners(); ++i) {
const Point &corner = c.Corner(i);
if (pData.find(corner) == pData.end()) {
std::vector<double> cornerData{};
for (int m : indices)
cornerData.push_back(E.Value(c.LocalCorner(i), data, m));
pData[corner] = cornerData;
}
}
}
ExchangeBuffer E;
for (const auto &[point, value] : pData) {
E.Send(0) << point;
for (int j = 0; j < indices.size(); ++j)
E.Send(0) << value[j];
}
communicatePointData(name, E, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::addDivergenceFreeData(const string &name, const Vector &data) {
checkMesh(data.GetMesh());
std::unordered_map<Point, Velocity> pData{};
for (cell c = data.cells(); c != data.cells_end(); ++c) {
DivergenceFreeElement E(data, c);
for (int i = 0; i < c.Corners(); ++i) {
const Point &corner = c.Corner(i);
if (pData.find(corner) == pData.end())
pData[corner] = E.VelocityField(c.LocalCorner(i), data);
}
}
ExchangeBuffer E;
for (const auto &[point, value] : pData) {
E.Send(0) << point << value[0] << value[1];
}
communicatePointData(name, E, 3, 2);
}
void SerialVtuPlot::DeformGrid(const Vector &deformation) {
checkMesh(deformation.GetMesh());
ExchangeBuffer exBuffer(deformation.CommSplit());
int dim = deformation.dim();
for(row r = deformation.rows(); r!= deformation.rows_end(); ++r){
if(plotMesh().find_vertex(r())!=plotMesh().vertices_end()){
for (row r = deformation.rows(); r != deformation.rows_end(); ++r) {
if (plotMesh().find_vertex(r()) != plotMesh().vertices_end()) {
exBuffer.Send(0) << r();
for (int j = 0; j < dim; ++j)
exBuffer.Send(0) << deformation(r, j);
......@@ -275,8 +330,7 @@ void SerialVtuPlot::checkMesh(const Mesh &other) {
if(!plotMesh.IsInstantiated()){
plotMesh(other);
initialize(other);
}
else if (&other != &plotMesh()) {
} else if (&other != &plotMesh()) {
THROW("Plotting: Meshes don't match!")
}
}
......
#include "VtuPlot.hpp"
#include "Config.hpp"
#include "LagrangeDiscretization.hpp"
#include "MultiPartDiscretization.hpp"
#include "DoFDiscretization.hpp"
#include "DGDiscretization.hpp"
#include "ArgyrisDiscretization.hpp"
#include "BasicDoFs.hpp"
#include "LagrangeDoF.hpp"
#include "MultiPartDoF.hpp"
#include "DGDoF.hpp"
#include <fstream>
......@@ -215,36 +217,39 @@ void VtuPlot::Flush(PlotStatus status) {
}
}
void VtuPlot::addData(const string &name, const Vector &data, std::vector<int> indices) {
const DoF &dofToPlot = data.GetDoF(0);
const IDiscretization &discToPlot = data.GetDisc();
if(indices.empty()){
int dim = min(3, dofToPlot.size());
int dim = min(3, ((const LagrangeDiscretization &) discToPlot).Size());
for(int i = 0;i<dim;++i){
indices.emplace_back(i);
}
}
if (typeid(dofToPlot) == typeid(LagrangeDoF)
|| typeid(dofToPlot) == typeid(MultiPartDoF)) {
int d = ((const LagrangeDoF &) dofToPlot).Degree();
if (typeid(discToPlot) == typeid(LagrangeDiscretization) ||
typeid(discToPlot) == typeid(MultiPartDiscretization)) {
int d = ((const LagrangeDiscretization &) discToPlot).Degree();
if (d > 0)
addPointData(name, data, indices);
else
addCellData(name, data, indices);
} else if (typeid(dofToPlot) == typeid(CellDoF)) {
} else if (typeid(discToPlot) == typeid(DoFDiscretization<CellDoF>)) {
addCellData(name, data, indices);
} else if (typeid(dofToPlot) == typeid(VertexDoF)) {
} else if (typeid(discToPlot) == typeid(DoFDiscretization<VertexDoF>)) {
addPointData(name, data, indices);
} else if (typeid(dofToPlot) == typeid(DGDoF)) {
} else if (typeid(discToPlot) == typeid(DGDiscretization) ||
typeid(discToPlot) == typeid(DoFDiscretization<DGDoF>)) {
addCellData(name, data, indices);
} else if (typeid(discToPlot) == typeid(ArgyrisDiscretization)) {
addArgyrisData(name, data, indices);
} else if (typeid(discToPlot) == typeid(DivergenceFreeDiscretization)) {
addDivergenceFreeData(name, data);
} else {
mout << "Warning: Plot type for DoF " << typeid(dofToPlot).name() << " not defined."
<< endl;
mout << "Warning: Plot type for Discretization " << typeid(discToPlot).name() << " not defined." << endl;
}
}
int VtuPlot::commSplit() {
return plotMesh().CommSplit();
}
int VtuPlot::commSplit() { return plotMesh().CommSplit(); }
VtuPlot &operator<<(VtuPlot &plot, const Vector &v) {
plot.AddData("value" + std::to_string(plot.DataCount()), v);
......
......@@ -53,7 +53,7 @@ public:
return mesh != nullptr;
}
void Reset(){
void Reset() {
mesh = nullptr;
}
};
......@@ -106,17 +106,23 @@ protected:
void stream(std::ostream &out);
void addData(const string& name, const Vector& data, std::vector<int> indices);
void addData(const string &name, const Vector &data, std::vector<int> indices);
virtual void
addCellData(const string &name, const Vector &data, const std::vector<int> &indices) = 0;
virtual void
addPointData(const string &name, const Vector &data, const std::vector<int> &indices) = 0;
int commSplit();
virtual void addArgyrisData(const string &name, const Vector &data, const std::vector<int> &indices) = 0;
virtual void addDivergenceFreeData(const string &name, const Vector &data) = 0;
int commSplit();;
public:
VtuPlot() : plotMesh() {}
explicit VtuPlot(const Mesh &init) : plotMesh(init){};
explicit VtuPlot(const Mesh &init) : plotMesh(init) {};
/// Creates a Plotter where the Vector init defines the Grid on which the VtuPlotter operates.
explicit VtuPlot(const Vector &init) : VtuPlot(init.GetMesh()) {};
......@@ -126,7 +132,8 @@ public:
}
void FileName(const string &name) { plotName = name; }
const std::string& FileName() const { return plotName; }
const std::string &FileName() const { return plotName; }
void ChangeDataPath(const string &dataPath) {
this->dataPath = dataPath;
......@@ -141,14 +148,16 @@ public:
*/
virtual void DeformGrid(const Vector &diplacement) = 0;
virtual void AddData(const string &name, const Vector& data) {
virtual void AddData(const string &name, const Vector &data) {
addData(name, data, std::vector<int>{});
};
virtual void AddData(const string& name, const Vector& data, int index) {
virtual void AddData(const string &name, const Vector &data, int index) {
addData(name, data, std::vector<int>{index});
}
virtual void AddData(const string& name, const Vector& data, const std::vector<int> &indices){
for (int i : indices){
virtual void AddData(const string &name, const Vector &data, const std::vector<int> &indices) {
for (int i : indices) {
std::string new_name = name + "_" + std::to_string(i);
AddData(new_name, data, i);
}
......@@ -160,7 +169,7 @@ public:
* @param data Contains the data to be plotted.
* @param dim
*/
//virtual void AddPointData(const string &name, const Vector &data) = 0;
//virtual void AddPointData(const string &name, const Vector &data) = 0;
/**
* Adds a CellData Container to the plotted file. As a Vector typically stores values in the vertices of the mesh,
......@@ -190,12 +199,17 @@ public:
void AddProcLoad(const Vector &gridVector) { AddSubdomain(gridVector.GetMesh()); }
virtual std::string FileString();
virtual void PlotFile(const string &filename) = 0;
virtual void PlotFile() { PlotFile(plotName); };
void ClearData();
void ClearMesh();
void Erase();
void Flush(PlotStatus s);
~VtuPlot() {
......@@ -206,19 +220,35 @@ public:
};
class SerialVtuPlot : public VtuPlot {
void initialize(const Mesh& init);
void initialize(const Mesh &init);
void initPoints(const Mesh &init);
void initCells(const Mesh &init);
void checkMesh(const Mesh& init);
void checkMesh(const Mesh &init);
void
addCellData(const string &name, const Vector &data, const std::vector<int> &indices) override;
void addCellData(const string &name, const Vector &data, const std::vector<int> &indices) override;
void
communicateCellData(const string& name, ExchangeBuffer &buffer, int dataDim, int numIndices);
void
addPointData(const string &name, const Vector &data, const std::vector<int> &indices) override;
void
communicatePointData(const string& name, ExchangeBuffer &buffer, int dataDim, int numIndices);
void
addArgyrisData(const string &name, const Vector &data, const std::vector<int> &indices) override;
void addDivergenceFreeData(const string &name, const Vector &data) override;
public:
SerialVtuPlot();
explicit SerialVtuPlot(const Mesh &init);
/// Creates a Plotter where the Vector init defines the Grid on which the VtuPlotter operates.
......@@ -264,10 +294,12 @@ public:
int vtkCellType(CELLTYPE cType);
VtuPlot &operator<<(VtuPlot &plot, const Vector &v);
VtuPlot &operator<<(VtuPlot &plot, PlotStatus status);
VtuPlot &operator<<(VtuPlot &plot, const string &name);
// Use as global macro to quickly save plot
VtuPlot &operator<<(VtuPlot &plot, const std::pair<std::string, PlotStatus>& saveoptions);
VtuPlot &operator<<(VtuPlot &plot, const std::pair<std::string, PlotStatus> &saveoptions);
#endif //VTUPLOT_H
......@@ -28,6 +28,7 @@ void Newton::operator()(const IAssemble &assemble, Vector &u) {
int JU_cnt = 0; // counts iteration steps without Jacobian update
Vector c(r);
line_iter = 0;
for (iter = 0; iter < max_iter; ++iter) {
vout(3) << "r(" << iter << ")= " << endl << r << endl;
if (d < eps) { break; }
......@@ -61,6 +62,7 @@ void Newton::operator()(const IAssemble &assemble, Vector &u) {
u += c;
d = assemble.Residual(u, r);
energy = assemble.Energy(u);
line_iter++;
if (d < d_previous) break;
}
}
......
......@@ -21,7 +21,8 @@ protected:
std::string prefix = "NonLinear";
int iter; // todo remove
int iter{0}; // todo remove
int line_iter{0};
double d, d_0; // todo remove
......@@ -57,6 +58,12 @@ public:
);
}
void PrintIterationInfo() const {
mout.PrintInfo("Newton iterations", verbose,
PrintInfoEntry("Steps", iter, 0),
PrintInfoEntry("line search steps",line_iter, 0));
}
virtual void operator()(const IAssemble &, Vector &) = 0;
virtual std::string Name() const = 0;
......