Commit a80d59d4 authored by JasonPries's avatar JasonPries
Browse files

Complete implementation of rotational motion, contour selection from Mesh

Implement more rotational motion tests
parent 1cd258a3
......@@ -44,7 +44,7 @@ def plot_scalar(fname):
plt.figure()
plt.set_cmap('viridis')
plt.gca().set_aspect('equal')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,vals,edgecolors='k')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,vals,edgecolors='k',shading='gouraud')
plt.colorbar()
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
......
......@@ -130,28 +130,28 @@ bool Mesh::edges_are_valid() const {
if (orientation(e)) {
Point p0 = base(e);
Point p1 = dc.curve()->point(dc.S0);
if (dist(p0,p1) > tol) {
if (dist(p0, p1) > tol) {
result = false;
break;
}
p0 = tip(e);
p1 = dc.curve()->point(dc.S1);
if (dist(p0,p1) > tol) {
if (dist(p0, p1) > tol) {
result = false;
break;
}
} else {
Point p0 = base(e);
Point p1 = dc.curve()->point(dc.S1);
if (dist(p0,p1) > tol) {
if (dist(p0, p1) > tol) {
result = false;
break;
}
p0 = tip(e);
p1 = dc.curve()->point(dc.S0);
if (dist(p0,p1) > tol) {
if (dist(p0, p1) > tol) {
result = false;
break;
}
......@@ -165,7 +165,7 @@ bool Mesh::edges_are_valid() const {
bool Mesh::find_attached(Point const p, size_t &e_out) {
double tol = length(e_out) * FLT_EPSILON;
if (dist(tip(e_out),p) < tol) {
if (dist(tip(e_out), p) < tol) {
return true;
}
......@@ -174,7 +174,7 @@ bool Mesh::find_attached(Point const p, size_t &e_out) {
if (e_out != twin(e_out)) {
e_out = next(twin(e_out));
while (e_out != e_in) {
if (dist(tip(e_out),p) < tol) {
if (dist(tip(e_out), p) < tol) {
return true;
} else if (e_out != twin(e_out)) {
e_out = next(twin(e_out));
......@@ -187,7 +187,7 @@ bool Mesh::find_attached(Point const p, size_t &e_out) {
if (e_out == twin(e_out)) {
e_out = prev(e_in);
while (e_out != twin(e_out)) {
if (dist(base(e_out),p) < tol) {
if (dist(base(e_out), p) < tol) {
e_out = twin(e_out);
return true;
} else {
......@@ -630,7 +630,7 @@ void Mesh::get_triangles() {
DartConstraint const dc = dart_constraint(bc->dart(0));
edge_queue.push_back(contour->orientation(0) ? dc.forward_dart() : dc.reverse_dart());
while(edge_queue.size() > 0) {
while (edge_queue.size() > 0) {
Edge &et = Edges[edge_queue.back()];
Edge &en = Edges[et.next()];
Edge &ep = Edges[et.prev()];
......@@ -640,12 +640,12 @@ void Mesh::get_triangles() {
if (et.Mark) {
Triangles.emplace_back(et.self(), contour_index);
if(!en.is_constrained() && en.Mark) {
if (!en.is_constrained() && en.Mark) {
edge_queue.push_back(en.twin());
en.Mark = false;
}
if(!ep.is_constrained() && ep.Mark) {
if (!ep.is_constrained() && ep.Mark) {
edge_queue.push_back(ep.twin());
ep.Mark = false;
}
......@@ -676,8 +676,8 @@ void Mesh::insert_internal_boundaries() {
// Find interior curves
std::vector<size_t> interior_index;
for (size_t i = 0; i!= BoundaryConstraints.size(); ++i) {
auto const & bc = BoundaryConstraints[i];
for (size_t i = 0; i != BoundaryConstraints.size(); ++i) {
auto const &bc = BoundaryConstraints[i];
bool on_exterior = false;
for (size_t i = 0; i != Boundary->size(); ++i) {
if (bc->curve() == Boundary->curve(i)) {
......@@ -834,7 +834,7 @@ void Mesh::save_as(std::string path, std::string file_name) const {
void Mesh::sort_constraints() {
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(),
[](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) {
return (size_t)(x->curve().get()) < (size_t)(y->curve().get());
return (size_t) (x->curve().get()) < (size_t) (y->curve().get());
}
);
}
......@@ -986,7 +986,6 @@ void Mesh::triangulate_boundary_polygon() {
}
LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
// TODO: Check for visibility of vertex? That someone contradicts the purpose of this method as a general point location function (as opposed to purely for mesh refinement)
double xp = p.X;
double yp = p.Y;
......@@ -1031,7 +1030,7 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
} else if (dist2 < tol_l) {
ei = prev(ei);
return LocateTriangleResult::Point;
} else if (area01 > tol_a && area12 > tol_a && area20 > tol_a) {
} else if (area01 > tol_a && area12 > tol_a && area20 > tol_a) {
return LocateTriangleResult::Interior;
} else if (area01 < -tol_a && ei != twin(ei)) {
ei = twin(ei);
......@@ -1520,13 +1519,38 @@ Point Mesh::circumcenter(size_t ei) const {
}
std::shared_ptr<BoundaryConstraint> Mesh::boundary_constraint(std::shared_ptr<Curve const> const &c) const {
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, size_t const y) { return (size_t)(x->curve().get()) < y; };
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, size_t const y) { return (size_t) (x->curve().get()) < y; };
std::shared_ptr<BoundaryConstraint> result = *std::lower_bound(BoundaryConstraints.begin(), BoundaryConstraints.end(), (size_t)c.get(), comp);
std::shared_ptr<BoundaryConstraint> result = *std::lower_bound(BoundaryConstraints.begin(), BoundaryConstraints.end(), (size_t) c.get(), comp);
if((size_t)(result->curve().get()) == (size_t)(c.get())) {
if ((size_t) (result->curve().get()) == (size_t) (c.get())) {
return result;
} else {
return nullptr;
}
}
std::shared_ptr<Contour const> Mesh::select_contour(Point const p) {
size_t e{0};
LocateTriangleResult result = locate_triangle(p, e);
if (result != LocateTriangleResult::Interior) {
return std::make_shared<Contour const>();
}
std::sort(Triangles.begin(), Triangles.end(), [](DartTriangle const &t0, DartTriangle const &t1) { return t0.edge() < t1.edge(); });
auto iter = std::lower_bound(Triangles.begin(), Triangles.end(), e, [](DartTriangle const &t0, size_t te) { return t0.edge() < te; });
size_t trip = 0;
while (iter->edge() != e && trip < 2) {
++trip;
e = next(e);
iter = std::lower_bound(Triangles.begin(), Triangles.end(), e, [](DartTriangle const &t0, size_t te) { return t0.edge() < te; });
}
if ((*iter).edge() != e) { // Fail
throw std::exception();
}
return Contours[iter->contour()];
}
\ No newline at end of file
......@@ -101,6 +101,10 @@ public:
// TODO: rename these methods something like, dart_constraint_from_edge_index, curve_from_edge_index
DartConstraint const dart_constraint_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; };
std::shared_ptr<Contour const> select_contour(Point p);
std::shared_ptr<Contour const> contour(size_t i) const { return Contours[i]; };
std::shared_ptr<Curve const> curve_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); };
std::shared_ptr<Curve const> curve(size_t i) const {return BoundaryConstraints[i]->curve(); };
......
......@@ -3,16 +3,16 @@ project(Oersted_Physics)
set(SOURCE_FILES
./include/Physics.hpp
#./src/PhysicsCommon.h
./src/BoundaryCondition.h ./src/BoundaryCondition.cpp
./src/DiscreteBoundary.h ./src/DiscreteBoundary.cpp
./src/BoundaryCondition.h ./src/BoundaryCondition.cpp
./src/Forcing.h ./src/Forcing.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
......@@ -21,10 +21,6 @@ set(SOURCE_FILES
./src/Physics.h ./src/Physics.cpp
#./src/Operator.h ./src/Operator.cpp
./src/Region.h ./src/Region.cpp
./src/Triangle.h ./src/Triangle.cpp
)
......
#ifndef OERSTED_PHYSICS_HPP
#define OERSTED_PHYSICS_HPP
#include "../src/DiscreteBoundary.h"
#include "../src/BoundaryCondition.h"
#include "../src/DiscreteBoundary.h"
#include "../src/DiscreteRegion.h"
#include "../src/FiniteElementMesh.h"
#include "../src/Forcing.h"
#include "../src/MaterialProperties.h"
......@@ -10,7 +11,6 @@
#include "../src/Node.h"
#include "../src/Physics.h"
#include "../src/PhysicsConstants.h"
#include "../src/Region.h"
#include "../src/Triangle.h"
#endif //OERSTED_PHYSICS_HPP
......@@ -19,7 +19,7 @@ class BoundaryCondition {
public:
virtual ~BoundaryCondition() {};
virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;
virtual void apply(Eigen::SparseMatrix<double> &bc_matrix) const = 0;
virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
};
......@@ -36,11 +36,11 @@ public:
ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<2, ElementOrder> const &fem) {
for (auto const &b : boundaries) {
auto fb = fem.boundary(b);
if (fb) { Boundaries.push_back(fb); }
if (fb[0]) { Boundaries.push_back(fb[0]); }
}
}
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
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) {
......@@ -66,43 +66,38 @@ 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} {
double_t sgn = sign();
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) {
auto const &n0 = fem.boundary(pbp.curve0())->nodes();
auto const &n1 = fem.boundary(pbp.curve1())->nodes();
if (pbp.orientation()) {
for (size_t i = 0; i != n0.size(); ++i) {
Map.emplace_back(n0[i], n1[i], sgn);
}
} else {
size_t j = n1.size();
for (size_t i = 0; i != n0.size(); ++i) {
Map.emplace_back(n0[i], n1[--j], sgn);
}
}
b0.push_back(fem.boundary(pbp.curve0())[0]);
b1.push_back(fem.boundary(pbp.curve1())[0]);
orientation.push_back(pbp.orientation());
}
// Remove duplicates
std::sort(Map.begin(), Map.end());
auto last = std::unique(Map.begin(), Map.end());
Map.erase(last, Map.end());
init(b0, b1, orientation);
}
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); };
}
// TODO: needs testing
PeriodicBoundaryCondition(std::vector<std::shared_ptr<DiscreteBoundary<2>>> b0, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b1, std::vector<bool> orientation, bool antiperiodic) : Antiperiodic{antiperiodic} {
init(b0, b1, orientation);
}
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
//double_t sgn = sign();
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.push_back(Eigen::Triplet<double>(v.first(), v.second(), v.value()));
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 {
......@@ -113,10 +108,50 @@ public:
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<2>>> b0, std::vector<std::shared_ptr<DiscreteBoundary<2>>> 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); };
}
}
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
......@@ -126,26 +161,41 @@ class SlidingInterface : public BoundaryCondition {
template<size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
SlidingInterface(std::vector<size_t> first_vars, std::vector<size_t> second_vars, std::vector<double> vals, bool antiperiodic = false)
: First{first_vars}, Second{second_vars}, Value{vals}, Antiperiodic{antiperiodic} {};
SlidingInterface(std::shared_ptr<DiscreteBoundary<2>> b0, std::shared_ptr<DiscreteBoundary<2>> b1, bool antiperiodic = false) : Antiperiodic{antiperiodic} {
std::vector<size_t> const &first = b0->nodes();
First = std::vector<size_t>(first.begin(), first.end() - 1);
SlidingInterface(std::shared_ptr<DiscontinuousBoundary<2>> db, bool antiperiodic = false)
: First{db->nodes()}, Second{db->dnodes()}, Value{db->values()}, Antiperiodic{antiperiodic} {};
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.end()) = -*(Value.end());
*(Value.rbegin()) = -*(Value.rbegin());
}
}
// TODO: SlidingInterface& operator--() {...};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
for (size_t i = 0; i != Value.size(); ++i) {
triplets.push_back(Eigen::Triplet<double>(First[i], Second[i], Value[i]));
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 {
......@@ -154,10 +204,18 @@ public:
}
}
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> Value;
std::vector<double_t> Value;
bool Antiperiodic;
};
......
......@@ -10,10 +10,6 @@ template<size_t Dimension>
class DiscreteBoundary {
};
template<size_t Dimension>
class DiscontinuousBoundary : public DiscreteBoundary<Dimension> {
};
template<>
class DiscreteBoundary<2> {
public:
......@@ -35,61 +31,20 @@ public:
size_t const &node(size_t i) const { return Nodes[i]; };
void node(size_t i, size_t n) { Nodes[i] = n; };
std::shared_ptr<Curve const> curve() const { return CurvePtr; };
bool is_closed() const { return Closed; };
operator size_t() const { return (size_t) (CurvePtr.get()); };
protected:
std::shared_ptr<Curve const> CurvePtr;
std::vector<size_t> Nodes;
bool Closed;
};
template<>
class DiscontinuousBoundary<2> : public DiscreteBoundary<2> {
public:
DiscontinuousBoundary() {};
DiscontinuousBoundary(std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals, bool closed = false) : DiscreteBoundary{nodes, closed}, DiscontinuousNodes{dnodes}, Values{vals} {};
DiscontinuousBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals, bool closed = false) : DiscreteBoundary{cptr, nodes, closed}, DiscontinuousNodes{dnodes}, Values{vals} {};
std::vector<size_t> dnodes() { return DiscontinuousNodes; };
std::vector<size_t> const &dnodes() const { return DiscontinuousNodes; };
std::vector<double_t> values() { return Values; };
std::vector<double_t> const &values() const { return Values; };
size_t dnode(size_t i) { return DiscontinuousNodes[i]; };
size_t const &dnode(size_t i) const { return DiscontinuousNodes[i]; };
protected:
std::vector<size_t> DiscontinuousNodes;
std::vector<double_t> Values;
};
/*
template <size_t Dimension>
class BoundaryPair {};
template <>
class BoundaryPair<2> {
public:
BoundaryPair(std::shared_ptr<Boundary<2>> first, std::shared_ptr<Boundary<2>> second) : First{first}, Second{second} {};
std::shared_ptr<Boundary<2>> const &first() { return First; };
std::shared_ptr<Boundary<2>> const &second() { return Second; };
protected:
std::shared_ptr<Boundary<2>> First;
std::shared_ptr<Boundary<2>> Second;
};
*/
class VariableMap {
public:
VariableMap(size_t first, size_t second, double value) : First{std::min(first, second)}, Second{std::max(first, second)}, Value{value} {};
......
#include "DiscreteRegion.h"
......@@ -3,18 +3,26 @@
#include <vector>
#include "Mesh.hpp"
#include "MaterialProperties.h"
template<size_t Dimension>
class Region {
class DiscreteRegion {
};
template<>
class Region<2> {
class DiscreteRegion<2> {
public:
Region() : Material{Air()} {};
Region(std::vector<size_t> elements) : Elements{elements}, Material{Air()} {};
Region(std::vector<size_t> elements, MaterialProperties material) : Elements{elements}, Material{material} {};
DiscreteRegion() : Region(nullptr), Material{Air()} {};
DiscreteRegion(std::vector<size_t> elements) : Region(nullptr), Elements{elements}, Material{Air()} {};
DiscreteRegion(std::vector<size_t> elements, MaterialProperties material) : Region(nullptr), Elements{elements}, Material{material} {};
DiscreteRegion(std::shared_ptr<Contour const> region) : Region{region}, Material{Air()} {};
DiscreteRegion(std::shared_ptr<Contour const> region, MaterialProperties material) : Region{region}, Material{material} {};
std::vector<size_t> const &elements() const { return Elements; };
......@@ -26,7 +34,11 @@ public:
void push_back(size_t element) { Elements.push_back(element); };
std::shared_ptr<Contour const> region() const { return Region; };
protected:
std::shared_ptr<Contour const> Region;
std::vector<size_t> Elements;
MaterialProperties Material;
......
......@@ -34,7 +34,7 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
Regions.reserve(m.size_contours());
for (size_t i = 0; i != m.size_contours(); ++i) {
Regions.push_back(std::make_shared<Region<2>>()); // TODO: Assign material properties here
Regions.push_back(std::make_shared<DiscreteRegion<2>>(m.contour(i))); // TODO: Assign material properties here
}
Triangles.reserve(m.size_triangles());
......@@ -43,4 +43,7 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
Regions[dt.contour()]->push_back(i);
Triangles.emplace_back(i, m.nodes_from_triangle(dt));
}
sort_boundaries();
sort_regions();
};
\ No newline at end of file
......@@ -10,7 +10,7 @@
#include "MatrixGroup.h"
#include "Node.h"
#include "Triangle.h"
#include "Region.h"
#include "DiscreteRegion.h"