Commit 838ff8f9 authored by Pries, Jason's avatar Pries, Jason

Reorganize Physics code to one class per file

parent b4a74b59
......@@ -3,9 +3,9 @@ cmake_minimum_required(VERSION 3.2)
project(Oersted)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage -fopenmp")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 -fopenmp")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
include_directories(./lib/)
include_directories(./lib/Eigen/)
......
......@@ -4,30 +4,29 @@ set(SOURCE_FILES
./include/Physics.hpp
./src/BoundaryCondition.h ./src/BoundaryCondition.cpp
./src/DiscreteBoundary.h ./src/DiscreteBoundary.cpp
./src/DiscreteRegion.h ./src/DiscreteRegion.cpp
./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
./src/Forcing.h ./src/Forcing.cpp
./src/MaterialProperties.h ./src/MaterialProperties.cpp
./src/MatrixGroup.h ./src/MatrixGroup.cpp
./src/PostProcessorInterface.h ./src/PostProcessorInterface.cpp
./src/Node.h ./src/Node.cpp
./src/Physics.h ./src/Physics.cpp
./src/PhysicsConstants.h
./src/Triangle.h ./src/Triangle.cpp
./src/Solution.cpp ./src/Solution.h
./src/PostProcessor.cpp ./src/PostProcessor.h)
./src/ZeroDirichlet.cpp ./src/ZeroDirichlet.h
./src/PeriodicBoundaryCondition.cpp ./src/PeriodicBoundaryCondition.h
./src/SlidingInterface.cpp ./src/SlidingInterface.h
./src/VariableMap.cpp ./src/VariableMap.h
./src/BoundaryMap.cpp ./src/BoundaryMap.h
./src/MappedBoundaryCondition.cpp ./src/MappedBoundaryCondition.h
./src/HomogeneousForcing.cpp ./src/HomogeneousForcing.h
./src/FunctionArguments.cpp ./src/FunctionArguments.h
./src/FiniteElementMeshInterface.cpp ./src/FiniteElementMeshInterface.h
)
add_library(physics SHARED ${SOURCE_FILES})
......
......@@ -11,6 +11,7 @@
#include "../src/Node.h"
#include "../src/Physics.h"
#include "../src/PhysicsConstants.h"
#include "../src/PostProcessorInterface.h"
#include "../src/Triangle.h"
#endif //OERSTED_PHYSICS_HPP
......@@ -17,203 +17,9 @@
class BoundaryCondition {
public:
virtual void apply(Eigen::SparseMatrix<double> &bc_matrix) const = 0;
virtual void apply(Eigen::SparseMatrix<double_t> &bc_matrix) const = 0;
virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
};
class ZeroDirichlet : public BoundaryCondition {
public:
ZeroDirichlet(std::vector<std::shared_ptr<DiscreteBoundary>> boundaries) : Boundaries{boundaries} {};
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);
if (fb[0]) { Boundaries.push_back(fb[0]); }
}
}
void apply(Eigen::SparseMatrix<double> &bc_matrix) const override {};
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (auto const &b : Boundaries) {
for (auto const &i : b->nodes()) {
index.insert(i);
}
}
};
protected:
std::vector<std::shared_ptr<DiscreteBoundary>> Boundaries;
};
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
class PeriodicBoundaryCondition : public BoundaryCondition {
public:
PeriodicBoundaryCondition(std::vector<VariableMap> map, bool antiperiodic) : Map{map}, 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>> b0;
std::vector<std::shared_ptr<DiscreteBoundary>> 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]);
orientation.push_back(pbp.orientation());
}
init(b0, b1, orientation);
}
PeriodicBoundaryCondition(std::vector<std::shared_ptr<DiscreteBoundary>> b0, std::vector<std::shared_ptr<DiscreteBoundary>> b1, std::vector<bool> orientation, bool antiperiodic) : Antiperiodic{antiperiodic} {
init(b0, b1, orientation);
}
void apply(Eigen::SparseMatrix<double> &bc_matrix) const override {
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(bc_matrix.rows());
for (size_t i = 0; i != bc_matrix.cols(); ++i) {
triplets.push_back(Eigen::Triplet<double>(i, i, 1.0));
}
for (VariableMap const &v : Map) {
triplets[v.second()] = Eigen::Triplet<double>(v.first(), v.second(), v.value());
}
Eigen::SparseMatrix<double> this_bc(bc_matrix.rows(), bc_matrix.cols());
this_bc.setFromTriplets(triplets.begin(), triplets.end());
bc_matrix = this_bc * bc_matrix * this_bc;
}
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (VariableMap const &v : Map) {
index.insert(v.second());
}
}
double_t sign() const { return (Antiperiodic ? -1.0 : 1.0); };
std::vector<VariableMap> const &map() const { return Map; };
std::vector<VariableMap> &map() { return Map; };
VariableMap const &map(size_t i) const { return Map[i]; };
VariableMap &map(size_t i) { return Map[i]; };
protected:
std::vector<VariableMap> Map;
bool Antiperiodic;
private:
void init(std::vector<std::shared_ptr<DiscreteBoundary>> b0, std::vector<std::shared_ptr<DiscreteBoundary>> b1, std::vector<bool> orientation) {
double_t sgn = sign();
for (size_t i = 0; i != b0.size(); ++i) {
auto const &n0 = b0[i]->nodes();
auto const &n1 = b1[i]->nodes();
if (orientation[i]) {
for (size_t j = 0; j != n0.size(); ++j) {
Map.emplace_back(n0[j], n1[j], sgn);
}
} else {
size_t k = n1.size();
for (size_t j = 0; j != n0.size(); ++j) {
Map.emplace_back(n0[j], n1[--k], sgn);
}
}
}
// Remove duplicates
std::sort(Map.begin(), Map.end());
auto last = std::unique(Map.begin(), Map.end());
Map.erase(last, Map.end());
if (!Antiperiodic) { // Remove center node if it exists
auto eq = [](VariableMap const &x) { return (x.first() == x.second()); };
auto center = std::find_if(Map.begin(), Map.end(), eq);
if (center != Map.end()) { Map.erase(center); };
}
}
};
class SlidingInterface : public BoundaryCondition {
public:
SlidingInterface(std::shared_ptr<DiscreteBoundary> b0, std::shared_ptr<DiscreteBoundary> b1, bool antiperiodic = false) : Antiperiodic{antiperiodic} {
std::vector<size_t> const &first = b0->nodes();
First = std::vector<size_t>(first.begin(), first.end() - 1);
std::vector<size_t> const &second = b1->nodes();
Second = std::vector<size_t>(second.begin(), second.end() - 1);
Value = std::vector<double_t>(First.size(), 1.0);
}
SlidingInterface &operator++() {
std::rotate(Second.begin(), Second.begin() + 1, Second.end());
if (Antiperiodic) {
std::rotate(Value.begin(), Value.begin() + 1, Value.end());
*(Value.rbegin()) = -*(Value.rbegin());
}
}
SlidingInterface &operator+=(const size_t &m) {
std::rotate(Second.begin(), Second.begin() + m, Second.end());
if (Antiperiodic) {
std::rotate(Value.begin(), Value.begin() + m, Value.end());
auto v_iter = Value.rbegin();
for (size_t i = 0; i != m; ++i) {
*(v_iter) = -*(v_iter);
++v_iter;
}
}
}
// TODO: SlidingInterface& operator--() {...};
// TODO: SlidingInterface& operator-=() {...};
void apply(Eigen::SparseMatrix<double> &bc_matrix) const override {
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(bc_matrix.rows());
for (size_t i = 0; i != bc_matrix.cols(); ++i) {
triplets.push_back(Eigen::Triplet<double>(i, i, 1.0));
}
for (size_t i = 0; i != First.size(); ++i) {
triplets[Second[i]] = Eigen::Triplet<double>(First[i], Second[i], Value[i]);
}
Eigen::SparseMatrix<double> this_bc(bc_matrix.rows(), bc_matrix.cols());
this_bc.setFromTriplets(triplets.begin(), triplets.end());
bc_matrix = this_bc * bc_matrix * this_bc;
}
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (size_t i = 0; i != Second.size(); ++i) {
index.insert(Second[i]);
}
}
std::vector<size_t> const first() const { return First; };
std::vector<size_t> const second() const { return Second; };
std::vector<double_t> const value() const { return Value; };
size_t size() const { return Value.size(); };
protected: // TODO: Use VariableMap instead?
std::vector<size_t> First;
std::vector<size_t> Second;
std::vector<double_t> Value;
bool Antiperiodic;
};
#endif //OERSTED_BOUNDARYCONDITION_H
#include "BoundaryMap.h"
#include <algorithm>
BoundaryMap::BoundaryMap(std::vector<size_t> first, std::vector<size_t> second, std::vector<double_t> value, bool ap) : Antiperiodic{ap} {
Map.reserve(first.size());
for (size_t i = 0; i != first.size(); ++i) {
Map.emplace_back(first[i], second[i], value[i]);
}
}
void BoundaryMap::move(long long delta) {
std::vector<size_t> first;
std::vector<size_t> second;
std::vector<double_t> value;
first.reserve(Map.size());
second.reserve(Map.size());
value.reserve(Map.size());
for (auto const &v : Map) {
first.push_back(v.first());
second.push_back(v.second());
value.push_back(v.value());
}
std::rotate(second.begin(), second.begin() + delta, second.end());
if (Antiperiodic) {
std::rotate(value.begin(), value.begin() + delta, value.end());
auto v_iter = value.rbegin();
for (size_t i = 0; i != delta; ++i) {
*(v_iter) = -*(v_iter);
++v_iter;
}
}
for(size_t i = 0; i != Map.size(); ++i) {
Map[i].second(second[i]);
Map[i].value(value[i]);
}
}
#ifndef OERSTED_BOUNDARYMAP_H
#define OERSTED_BOUNDARYMAP_H
#include <vector>
#include "VariableMap.h"
class BoundaryMap {
public:
BoundaryMap() {};
BoundaryMap(std::vector<VariableMap> map) : Map{map} {}
BoundaryMap(std::vector<size_t> first, std::vector<size_t> second, std::vector<double_t> value, bool ap);
void move(long long delta);
size_t size() const { return Map.size(); }
auto begin() const { return Map.begin(); }
auto end() const { return Map.end(); }
protected:
bool Antiperiodic{false};
std::vector<VariableMap> Map;
};
#endif //OERSTED_BOUNDARYMAP_H
#include "DiscreteBoundary.h"
#include "DiscreteBoundary.h"
\ No newline at end of file
......@@ -4,7 +4,7 @@
#include <vector>
#include <cstddef>
#include <Sketch.hpp>
#include "Sketch.hpp"
class DiscreteBoundary {
public:
......@@ -36,54 +36,10 @@ public:
protected:
std::shared_ptr<Curve const> CurvePtr;
std::vector<size_t> Nodes;
bool Closed;
};
class VariableMap {
public:
VariableMap(size_t first, size_t second, double value) : First{std::min(first, second)}, Second{std::max(first, second)}, Value{value} {};
size_t first() const { return First; };
size_t second() const { return Second; };
double_t value() const { return Value; };
void value(double_t value) { Value = value; };
inline bool operator==(VariableMap const &rhs) {
return (First == rhs.First && Second == rhs.Second && Value == rhs.Value);
}
inline bool operator<(VariableMap const &rhs) {
return (First < rhs.First) || (First == rhs.First && Second < rhs.Second);
}
protected:
size_t First;
size_t Second;
double_t Value;
};
class BoundaryMap {
public:
BoundaryMap(std::vector<VariableMap> map) : Map{map} {};
size_t size() const { return Map.size(); };
VariableMap const &map(size_t i) const { return Map[i]; };
std::vector<VariableMap> const &map() const { return Map; };
size_t first(size_t i) const { return Map[i].first(); };
size_t second(size_t i) const { return Map[i].second(); };
double value(size_t i) const { return Map[i].value(); };
std::vector<size_t> Nodes;
protected:
std::vector<VariableMap> Map;
bool Closed;
};
#endif //OERSTED_BOUNDARY_H
......@@ -7,7 +7,6 @@
#include "MaterialProperties.h"
class DiscreteRegion {
public:
DiscreteRegion() : Region(nullptr), Material{Air()} {};
......
#include "FiniteElementMesh.h"
/*
template<size_t Order, size_t QuadratureOrder> FiniteElementMesh<Order, QuadratureOrder>::FiniteElementMesh(Mesh &m) {
template<size_t P, size_t Q>
FiniteElementMesh<P, Q>::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.reserve(m.size_points()); // TODO: This will be larger when ElementOrder > 1
for (size_t i = 0; i != m.size_points(); ++i) {
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
......@@ -26,16 +26,16 @@ template<size_t Order, size_t QuadratureOrder> FiniteElementMesh<Order, Quadratu
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)));
for (size_t i = 0; i != m.size_boundary_constraints(); ++i) {
Boundaries.push_back(std::make_shared<DiscreteBoundary>(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
Regions.push_back(std::make_shared<DiscreteRegion>(m.contour(i))); // TODO: Assign material properties here
}
Triangles.reserve(m.size_triangles());
......@@ -48,4 +48,286 @@ template<size_t Order, size_t QuadratureOrder> FiniteElementMesh<Order, Quadratu
sort_boundaries();
sort_regions();
};
*/
\ No newline at end of file
template<size_t P, size_t Q>
DiagonalMatrixGroup FiniteElementMesh<P, Q>::determinant() const {
DiagonalMatrixGroup mat(TriangleQuadrature<Q>::size, Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].template determinant<Q>(mat, Nodes);
}
return mat;
}
template<size_t P, size_t Q>
SparseMatrixGroup FiniteElementMesh<P, Q>::basis() const {
SparseMatrixGroup mat(TriangleQuadrature<Q>::size, Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].template basis<Q>(mat, Nodes);
}
return mat;
}
template<size_t P, size_t Q>
DerivativeMatrixGroup FiniteElementMesh<P, Q>::derivative() const {
DerivativeMatrixGroup df(TriangleQuadrature<Q>::size, Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].template derivative<Q>(df, Nodes);
}
return df;
}
template<size_t P, size_t Q>
std::vector<std::vector<XY>> FiniteElementMesh<P, Q>::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].template quadrature_points<Q>(Nodes));
}
return qp;
}
// TODO: std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Region const> const &r) const override {...}
template<size_t P, size_t Q>
std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_discontinuous(std::shared_ptr<Curve const> const &c) {
// TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
std::vector<std::shared_ptr<DiscreteBoundary>> b_vec = boundary(c);
if (b_vec.size() == 0 || b_vec.size() == 2) { // Boundary curve not found || boundary pair found
return b_vec;
} else if (b_vec.size() > 2) { // should not happen
throw std::runtime_error("FiniteElementMesh<2>::make_discontinuous expects to find at most 2 DiscreteBoundary objects");
} // else b.size() == 1 and we need to construct an new boundary
// Copy nodes
std::vector<size_t> boundary_nodes = b_vec[0]->nodes();
std::vector<size_t> discontinuous_nodes;
discontinuous_nodes.reserve(boundary_nodes.size());
for (size_t i = 0; i != boundary_nodes.size(); ++i) {
discontinuous_nodes.push_back(Nodes.size());
Nodes.emplace_back(Nodes[boundary_nodes[i]]);
}
// Construct new DiscontinuousBoundary
bool is_closed = b_vec[0]->is_closed();
std::shared_ptr < DiscreteBoundary > db = std::make_shared<DiscreteBoundary>(b_vec[0]->curve(), discontinuous_nodes, is_closed);
Boundaries.push_back(db);
sort_boundaries();
// Renumber nodes in elements when the element is on the "right" side of the plane defined by the discontinuous boundary edge
// (cross product of edge vector with edge midpoint to incenter vector is negative)
// TODO: This a brute renumbering algorithm with O(t*n) complexity where t = Triangles.size() and n = bnodes.size()
// TODO: Sort + renumber should result in O(n*log(t)) complexity, but all present assumptions about this class assume the the ordering of the elements never changes
// TODO: However, triangles do store an ID which could reverse the sorting
for (size_t i = 0; i != (discontinuous_nodes.size() - 1 + is_closed); ++i) {
size_t ii = (i + 1) % discontinuous_nodes.size();
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<P> &t : Triangles) {
for (size_t k = 0; k != 3; ++k) { // TODO: This essentially assumes Element order = 1
if (t.node(k) == boundary_nodes[i]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
double_t cross = dx0 * dy1 - dy0 * dx1;
if (cross < 0) {
t.node(k, discontinuous_nodes[i]);
}
}
if (t.node(k) == boundary_nodes[ii]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
double_t cross = dx0 * dy1 - dy0 * dx1;