Commit 07c274f3 authored by Pries, Jason's avatar Pries, Jason
Browse files

Start of de-templating Physics implementation

    Part of this involves using more shared pointers and vectors (as opposed to arrays) as certain member variables
    Some of that has been done, some work may (probably) still needs to be done
    It's been a few months since this has been looked at so it will require some review

This is a bit of a mess right now
parent b270e42c
......@@ -25,7 +25,9 @@ set(SOURCE_FILES
./src/Triangle.h ./src/Triangle.cpp
src/Solution.cpp src/Solution.h)
./src/Solution.cpp ./src/Solution.h
./src/PostProcessor.cpp ./src/PostProcessor.h)
add_library(physics SHARED ${SOURCE_FILES})
......
......@@ -17,7 +17,7 @@
class BoundaryCondition {
public:
virtual ~BoundaryCondition() {};
virtual ~BoundaryCondition() { std::cout << "Get rid of unnecessary templates" << std::endl; };
virtual void apply(Eigen::SparseMatrix<double> &bc_matrix) const = 0;
......@@ -33,9 +33,10 @@ class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition
public:
ZeroDirichlet(std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundaries) : Boundaries{boundaries} {};
ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<2, ElementOrder> const &fem) {
//ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<ElementOrder, QuadratureOrder> const &fem) {
ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, std::shared_ptr<FiniteElementMeshInterface> fem) {
for (auto const &b : boundaries) {
auto fb = fem.boundary(b);
auto fb = fem->boundary(b);
if (fb[0]) { Boundaries.push_back(fb[0]); }
}
}
......@@ -65,14 +66,15 @@ class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public Bound
public:
PeriodicBoundaryCondition(std::vector<VariableMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2, ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
//PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<ElementOrder, QuadratureOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, std::shared_ptr<FiniteElementMeshInterface> fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
std::vector<std::shared_ptr<DiscreteBoundary<2>>> b0;
std::vector<std::shared_ptr<DiscreteBoundary<2>>> b1;
std::vector<bool> orientation;
for (auto const &pbp : periodic_boundary_pairs) {
b0.push_back(fem.boundary(pbp.curve0())[0]);
b1.push_back(fem.boundary(pbp.curve1())[0]);
b0.push_back(fem->boundary(pbp.curve0())[0]);
b1.push_back(fem->boundary(pbp.curve1())[0]);
orientation.push_back(pbp.orientation());
}
......
#include "FiniteElementMesh.h"
template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
/*
template<size_t Order, size_t QuadratureOrder> FiniteElementMesh<Order, QuadratureOrder>::FiniteElementMesh(Mesh &m) {
// std::vector<XY> nodes
// std::vector<Triangle<Order>> tris
// std::vector<std::shared_ptr<Region<2>>> r
......@@ -11,7 +12,7 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
Nodes.emplace_back(m.point(i));
}
/* TODO: Changes to rotational motion implementation
// TODO: Changes to rotational motion implementation
Initially, perform a straightforward conversion as currently implemented
Define a new Boundary type called (called something like SlidingBoundary)
SlidingBoundary holds a reference to the two boundaries and ensures the correct mapping
......@@ -25,7 +26,7 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
TODO: To make these easier to implement from the user side, Boundary should hold a shared_ptr to the defining curve (from the Mesh object)
Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
*/
//
Boundaries.reserve(m.size_boundary_constraints());
for (size_t i = 0; i!= m.size_boundary_constraints(); ++i) {
......@@ -47,3 +48,4 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
sort_boundaries();
sort_regions();
};
*/
\ No newline at end of file
......@@ -12,35 +12,26 @@
#include "Triangle.h"
#include "DiscreteRegion.h"
template<size_t Dimension, size_t Order>
class FiniteElementMesh {
};
template<size_t Order>
class FiniteElementMesh<2, Order> {
class FiniteElementMeshInterface {
public:
FiniteElementMesh() {};
FiniteElementMeshInterface() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {
FiniteElementMeshInterface(std::vector<XY> nodes, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Regions(r), Boundaries(b) {
sort_boundaries();
sort_regions();
};
FiniteElementMesh(Mesh &m);
virtual ~FiniteElementMeshInterface() {};
size_t size_nodes() const { return Nodes.size(); };
size_t size_elements() const { return Triangles.size(); };
std::vector<Triangle<Order>> const &triangles() const { return Triangles; };
Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
auto &boundaries() { return Boundaries; };
auto &boundary(size_t i) { return Boundaries[i]; };
std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundary(std::shared_ptr<Curve const> const &c) const {
// Boundary may be discontinuous (e.g. multiple overlapping curves). Therefore, multiple DiscreteBoundaries may be returned (in general, upper != lower)
auto lower_comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
auto upper_comp = [](size_t const &y, std::shared_ptr<DiscreteBoundary<2>> const &x) { return y < (size_t) (x->curve().get()); };
......@@ -54,7 +45,174 @@ public:
}
};
std::vector<std::shared_ptr<DiscreteBoundary<2>>> make_discontinuous(std::shared_ptr<Curve const> const &c) {
// TODO: virtual std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Region const> const &r) const = 0;
virtual std::vector<std::shared_ptr<DiscreteBoundary<2>>> make_discontinuous(std::shared_ptr<Curve const> const &c) = 0;
auto const &boundaries() const { return Boundaries; };
auto const &nodes() const { return Nodes; };
auto const &node(size_t i) const { return Nodes[i]; };
auto &regions() { return Regions; };
auto &region(size_t i) { return Regions[i]; };
auto const &regions() const { return Regions; };
auto const &region(size_t i) const { return Regions[i]; };
std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Contour const> const &c) const {
auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, size_t const &y) { return (size_t) (x->region().get()) < y; };
auto iter = std::lower_bound(Regions.begin(), Regions.end(), (size_t) c.get(), comp);
if (iter != Regions.end()) {
return *iter;
} else {
return std::make_shared<DiscreteRegion<2>>();
}
}
size_t nodes_size() const { return Nodes.size(); };
virtual size_t elements_size() const = 0;
virtual void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const = 0;
virtual void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const = 0;
virtual DiagonalMatrixGroup determinant() const = 0;
virtual SparseMatrixGroup basis() const = 0;
virtual DerivativeMatrixGroup derivative() const = 0;
virtual std::vector<std::vector<XY>> quadrature_points() const = 0;
protected:
std::vector<XY> Nodes;
std::vector<std::shared_ptr<DiscreteRegion<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
void sort_boundaries() {
auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, std::shared_ptr<DiscreteBoundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
std::sort(Boundaries.begin(), Boundaries.end(), comp);
}
void sort_regions() {
auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, std::shared_ptr<DiscreteRegion<2>> const &y) { return (size_t) (x->region().get()) < (size_t) (y->region().get()); };
std::sort(Regions.begin(), Regions.end(), comp);
}
};
template<size_t ElementOrder, size_t QuadratureOrder>
class FiniteElementMesh : public FiniteElementMeshInterface {
public:
FiniteElementMesh() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<ElementOrder>> tris, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : FiniteElementMeshInterface{nodes,r,b}, Triangles(tris) {};
FiniteElementMesh(Mesh &m) {
// std::vector<XY> nodes
// std::vector<Triangle<Order>> tris
// std::vector<std::shared_ptr<Region<2>>> r
// std::vector<std::shared_ptr<Boundary<2>>> b
Nodes.reserve(m.size_points()); // TODO: This will be larger when Order > 1
for (size_t i = 0; i!= m.size_points(); ++i) {
Nodes.emplace_back(m.point(i));
}
/* TODO: Changes to rotational motion implementation
Initially, perform a straightforward conversion as currently implemented
Define a new Boundary type called (called something like SlidingBoundary)
SlidingBoundary holds a reference to the two boundaries and ensures the correct mapping
In FiniteElementMesh, add a method called make_sliding(...) or something similar
This method will implement a procedure which duplicates the nodes on the boundary and renumbers the nodes in the appropriate elements
Then, the SlidingInterface boundary condition holds a pointer to the SlidingBoundary instead of two Boundary objects
Perhaps a similar approach can be taken for implementing periodic boundary conditions (e.g. MappedBoundary holds two (or several arrays of) boundaries and enforces the proper invariant)
e.g. make_mapped(...)
Then PeriodicBoundaryCondition holds a pointer to the MappedBoundary instead of two Boundary objects
TODO: To make these easier to implement from the user side, Boundary should hold a shared_ptr to the defining curve (from the Mesh object)
Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
*/
Boundaries.reserve(m.size_boundary_constraints());
for (size_t i = 0; i!= m.size_boundary_constraints(); ++i) {
Boundaries.push_back(std::make_shared<DiscreteBoundary<2>>(m.curve(i), m.boundary_nodes(i)));
}
Regions.reserve(m.size_contours());
for (size_t i = 0; i != m.size_contours(); ++i) {
Regions.push_back(std::make_shared<DiscreteRegion<2>>(m.contour(i))); // TODO: Assign material properties here
}
Triangles.reserve(m.size_triangles());
for (size_t i = 0; i != m.size_triangles(); ++i) {
auto &dt = m.triangle(i);
Regions[dt.contour()]->push_back(i);
Triangles.emplace_back(i, m.nodes_from_triangle(dt));
}
sort_boundaries();
sort_regions();
};
std::vector<Triangle<ElementOrder>> const &triangles() const { return Triangles; };
Triangle<ElementOrder> const &triangle(size_t i) const { return Triangles[i]; };
size_t elements_size() const override { return Triangles.size(); };
DiagonalMatrixGroup determinant() const override {
DiagonalMatrixGroup mat(TriangleQuadrature<QuadratureOrder>::size, Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i]. template determinant<QuadratureOrder>(mat, Nodes);
}
return mat;
}
SparseMatrixGroup basis() const override {
SparseMatrixGroup mat(TriangleQuadrature<QuadratureOrder>::size, Nodes.size(), Triangles.size(), Triangle<ElementOrder>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i]. template basis<QuadratureOrder>(mat, Nodes);
}
return mat;
}
DerivativeMatrixGroup derivative() const override {
DerivativeMatrixGroup df(TriangleQuadrature<QuadratureOrder>::size, Nodes.size(), Triangles.size(), Triangle<ElementOrder>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i]. template derivative<QuadratureOrder>(df, Nodes);
}
return df;
}
std::vector<std::vector<XY>> quadrature_points() const override {
std::vector<std::vector<XY>> qp;
qp.reserve(Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
qp.push_back(Triangles[i]. template quadrature_points<QuadratureOrder>(Nodes));
}
return qp;
}
// TODO: std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Region const> const &r) const override {...}
std::vector<std::shared_ptr<DiscreteBoundary<2>>> make_discontinuous(std::shared_ptr<Curve const> const &c) override {
// TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
std::vector<std::shared_ptr<DiscreteBoundary<2>>> b_vec = boundary(c);
......@@ -90,7 +248,7 @@ public:
double_t dx0 = Nodes[boundary_nodes[ii]].x() - Nodes[boundary_nodes[i]].x();
double_t dy0 = Nodes[boundary_nodes[ii]].y() - Nodes[boundary_nodes[i]].y();
for (Triangle<Order> &t : Triangles) {
for (Triangle<ElementOrder> &t : Triangles) {
for (size_t k = 0; k != 3; ++k) {
if (t.node(k) == boundary_nodes[i]) {
size_t k0 = t.node(0);
......@@ -204,78 +362,7 @@ public:
return b_vec;
}
auto const &boundaries() const { return Boundaries; };
auto const &nodes() const { return Nodes; };
auto const &node(size_t i) const { return Nodes[i]; };
auto &regions() { return Regions; };
auto &region(size_t i) { return Regions[i]; };
auto const &regions() const { return Regions; };
auto const &region(size_t i) const { return Regions[i]; };
std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Contour const> const &c) const {
auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, size_t const &y) { return (size_t) (x->region().get()) < y; };
auto iter = std::lower_bound(Regions.begin(), Regions.end(), (size_t) c.get(), comp);
if (iter != Regions.end()) {
return *iter;
} else {
return std::make_shared<DiscreteRegion<2>>();
}
}
template<size_t Q>
DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
DiagonalMatrixGroup<TriangleQuadrature<Q>::size> mat(Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].determinant<Q>(mat, Nodes);
}
return mat;
}
template<size_t Q>
SparseMatrixGroup<TriangleQuadrature<Q>::size> basis() const {
SparseMatrixGroup<TriangleQuadrature<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].basis<Q>(mat, Nodes);
}
return mat;
}
template<size_t Q>
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> derivative() const {
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].derivative<Q>(df, Nodes);
}
return df;
}
template<size_t Q>
std::vector<std::vector<XY>> quadrature_points() const {
std::vector<std::vector<XY>> qp;
qp.reserve(Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
}
return qp;
}
void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const {
void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const override {
if (!boost::filesystem::exists(path)) {
boost::filesystem::create_directories(path);
}
......@@ -294,11 +381,11 @@ public:
fs << loc << ',' << v.size() << std::endl;
// Write triangles
for (Triangle<Order> const &t : Triangles) {
for (size_t i = 0; i != Triangle<Order>::NumNodes - 1; ++i) {
for (Triangle<ElementOrder> const &t : Triangles) {
for (size_t i = 0; i != Triangle<ElementOrder>::NumNodes - 1; ++i) {
fs << t.node(i) << ',';
}
fs << t.node(Triangle<Order>::NumNodes - 1) << std::endl;
fs << t.node(Triangle<ElementOrder>::NumNodes - 1) << std::endl;
}
// Write nodes
......@@ -315,7 +402,7 @@ public:
fs.close();
}
void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const {
void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const override {
if (!boost::filesystem::exists(path)) {
boost::filesystem::create_directories(path);
}
......@@ -334,11 +421,11 @@ public:
fs << loc << ',' << Fx.size() << std::endl;
// Write triangles
for (Triangle<Order> const &t : Triangles) {
for (size_t i = 0; i != Triangle<Order>::NumNodes - 1; ++i) {
for (Triangle<ElementOrder> const &t : Triangles) {
for (size_t i = 0; i != Triangle<ElementOrder>::NumNodes - 1; ++i) {
fs << t.node(i) << ',';
}
fs << t.node(Triangle<Order>::NumNodes - 1) << std::endl;
fs << t.node(Triangle<ElementOrder>::NumNodes - 1) << std::endl;
}
// Write nodes
......@@ -355,21 +442,7 @@ public:
}
protected:
std::vector<XY> Nodes;
std::vector<Triangle<Order>> Triangles;
std::vector<std::shared_ptr<DiscreteRegion<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
void sort_boundaries() {
auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, std::shared_ptr<DiscreteBoundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
std::sort(Boundaries.begin(), Boundaries.end(), comp);
}
void sort_regions() {
auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, std::shared_ptr<DiscreteRegion<2>> const &y) { return (size_t) (x->region().get()) < (size_t) (y->region().get()); };
std::sort(Regions.begin(), Regions.end(), comp);
}
std::vector<Triangle<ElementOrder>> Triangles;
};
#endif //OERSTED_FINITEELEMENTMESH_H
......@@ -28,8 +28,8 @@ class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Forcing {
public:
HomogeneousForcing(std::function<double(double)> f,
std::vector<std::shared_ptr<DiscreteRegion<2>>> regions,
SparseMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &basis,
DiagonalMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &weights)
SparseMatrixGroup const &basis,
DiagonalMatrixGroup const &weights)
: Function{f}, Regions{regions}, Basis(basis), Weights(weights) {
// TODO: May need to adjust implementation once quadrilaterals are added
......@@ -59,8 +59,8 @@ protected:
std::vector<std::shared_ptr<DiscreteRegion<2>>> Regions;
SparseMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &Basis;
DiagonalMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &Weights;
SparseMatrixGroup const &Basis;
DiagonalMatrixGroup const &Weights;
Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
};
......
......@@ -4,15 +4,11 @@
#include "Eigen"
#include "Eigen/Sparse"
// TODO: template<size_t Q, size_t D> MatrixGroup {};
// TODO: using template<size_t Q> BasisMatrixGroup = MatrixGroup<Q,1>;
// TODO: using template<size_t Q> DerivativeMatrixGroup = MatrixGroup<Q,2>;
template<size_t Q>
class SparseMatrixGroup {
public:
SparseMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) {
for (size_t i = 0; i != Q; ++i) {
SparseMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) {
Matrices.resize(Q);
for (size_t i = 0; i != Matrices.size(); ++i) {
Matrices[i].resize(rows, cols);
Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
}
......@@ -25,20 +21,20 @@ public:
auto &operator[](size_t const q) const { return Matrices[q]; };
void transform(Eigen::SparseMatrix<double> &A) {
for (size_t i = 0; i != Q; ++i) {
for (size_t i = 0; i != Matrices.size(); ++i) {
Matrices[i] = A * Matrices[i];
}
}
protected:
std::array<Eigen::SparseMatrix<double>, Q> Matrices;
std::vector<Eigen::SparseMatrix<double_t>> Matrices;
};
template<size_t Q>
class DiagonalMatrixGroup {
public:
DiagonalMatrixGroup(size_t const dim) {
for (size_t i = 0; i != Q; ++i) {
DiagonalMatrixGroup(size_t Q, size_t const dim) {
Matrices.resize(Q);
for (size_t i = 0; i != Matrices.size(); ++i) {
Matrices[i].resize(dim);
}
}
......@@ -47,14 +43,12 @@ public:
auto &operator()(size_t const q) { return Matrices[q]; };
protected:
//std::array<Eigen::DiagonalMatrix<double, Eigen::Dynamic>, Q> Matrices;
std::array<Eigen::ArrayXd, Q> Matrices;
std::vector<Eigen::ArrayXd> Matrices;
};
template<size_t Q>
class DerivativeMatrixGroup {
public:
DerivativeMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) : Dx{rows, cols, nnz}, Dy{rows, cols, nnz} {};
DerivativeMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) : Dx{Q, rows, cols, nnz}, Dy{Q, rows, cols, nnz} {};
auto const &dx(size_t const q) const { return Dx[q]; };
......@@ -70,8 +64,8 @@ public:
}
protected:
SparseMatrixGroup<Q> Dx;
SparseMatrixGroup<Q> Dy;
SparseMatrixGroup Dx;
SparseMatrixGroup Dy;
};
#endif //OERSTED_MATRIXGROUP_H
......@@ -38,23 +38,8 @@ enum class FieldVariable {
Psi = 26 + 23
};
class PhysicsInterface {
public:
double time() const { return Time; };
double &time() { return Time; };
void time(double t) { Time = t; };
protected:
double Time{0.0};
std::vector<std::shared_ptr<Forcing>> ForcingCondtions;
std::vector<std::shared_ptr<BoundaryCondition>> BoundaryConditions;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class PhysicsData {
class PhysicsInterface {
public:
static constexpr size_t QuadraturePoints = TriangleQuadrature<QuadratureOrder>::size;
......@@ -66,50 +51,59 @@ public:
virtual void apply_conditions() = 0;
*/
PhysicsData(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
ElementWeights{d.triangles().size()},
Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
//PhysicsInterface(FiniteElementMesh<ElementOrder, QuadratureOrder> &d) : Domain{d},
PhysicsInterface(std::shared_ptr<FiniteElementMeshInterface> d) : Domain{d},
ElementWeights{QuadraturePoints, d->elements_size()},
Derivatives{QuadraturePoints, d->nodes_size(), d->elements_size(), Triangle<ElementOrder>::NumNodes},
Basis{QuadraturePoints, d->nodes_size(), d->elements_size(), Triangle<ElementOrder>::NumNodes} {}
double time() const { return Time; };
double &time() { return Time; };
void time(double t) { Time = t; };