Commit 9cdbe1be authored by JasonPries's avatar JasonPries
Browse files

Adds initial support for external plotting of scalar fields using matplotlib

Fix minor inconsistency with antiperiodic boundary conditions
parent e21f3e59
...@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt ...@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt
imgtypes = ['.pdf','.png'] imgtypes = ['.pdf','.png']
def plot_sketch(fname): def plot_sketch(fname):
fname += '.oesk'
print('Sketch: ' +fname) print('Sketch: ' +fname)
data = np.genfromtxt(fname, dtype=float, delimiter=',') data = np.genfromtxt(fname, dtype=float, delimiter=',')
plt.plot(data[:,0], data[:,1]) plt.plot(data[:,0], data[:,1])
...@@ -15,6 +16,7 @@ def plot_sketch(fname): ...@@ -15,6 +16,7 @@ def plot_sketch(fname):
plt.clf() plt.clf()
def plot_mesh(fname): def plot_mesh(fname):
fname += '.oeme'
print('Mesh: '+fname) print('Mesh: '+fname)
data = np.genfromtxt(fname, dtype=float, delimiter=',') data = np.genfromtxt(fname, dtype=float, delimiter=',')
x = data[:,[0,1,2,0]].transpose() x = data[:,[0,1,2,0]].transpose()
...@@ -26,14 +28,40 @@ def plot_mesh(fname): ...@@ -26,14 +28,40 @@ def plot_mesh(fname):
plt.savefig(fname+ext) plt.savefig(fname+ext)
plt.clf() plt.clf()
def plot_scalar(fname):
fname += '.oesc'
print('Scalar: '+fname)
tris = np.genfromtxt(fname, dtype=int, delimiter=',',skip_header=0,max_rows=1)
tris = np.genfromtxt(fname, dtype=int, delimiter=',',skip_header=tris[0],max_rows=tris[1])
nodes = np.genfromtxt(fname, dtype=int, delimiter=',',skip_header=1,max_rows=1)
nodes = np.genfromtxt(fname, dtype=float, delimiter=',',skip_header=nodes[0],max_rows=nodes[1])
vals = np.genfromtxt(fname, dtype=int, delimiter=',',skip_header=2,max_rows=1)
vals = np.genfromtxt(fname, dtype=float, delimiter=',',skip_header=vals[0],max_rows=vals[1])
plt.figure()
plt.set_cmap('viridis')
plt.gca().set_aspect('equal')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,vals,edgecolors='k')
plt.colorbar()
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
for ext in imgtypes:
plt.savefig(fname+ext)
plt.clf()
root = os.getcwd() root = os.getcwd()
root += '/build/test/output/' root += '/build/test/output/'
for path, _, files in os.walk(root): for path, _, files in os.walk(root):
for name in files: for name in files:
_, fext = os.path.splitext(name) name, fext = os.path.splitext(name)
fpath = os.path.join(path, name) fpath = os.path.join(path, name)
if fext == '.oesk': if fext == '.oesk':
plot_sketch(fpath) plot_sketch(fpath)
elif fext == '.oeme': elif fext == '.oeme':
plot_mesh(fpath) plot_mesh(fpath)
elif fext == '.oesc':
plot_scalar(fpath)
...@@ -63,7 +63,9 @@ public: ...@@ -63,7 +63,9 @@ public:
size_t second() const { return Second; }; size_t second() const { return Second; };
double value() const { return Value; }; double_t value() const { return Value; };
void value(double_t value) { Value = value; };
inline bool operator==(VariableMap const &rhs) { inline bool operator==(VariableMap const &rhs) {
return (First == rhs.First && Second == rhs.Second && Value == rhs.Value); return (First == rhs.First && Second == rhs.Second && Value == rhs.Value);
...@@ -76,7 +78,7 @@ public: ...@@ -76,7 +78,7 @@ public:
protected: protected:
size_t First; size_t First;
size_t Second; size_t Second;
double Value; double_t Value;
}; };
class BoundaryMap { class BoundaryMap {
......
...@@ -98,10 +98,10 @@ public: ...@@ -98,10 +98,10 @@ public:
} }
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override { void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double_t sgn = sign(); //double_t sgn = sign();
for (VariableMap const &v : Map) { for (VariableMap const &v : Map) {
triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sgn * v.value())); triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), v.value()));
} }
} }
......
...@@ -112,6 +112,56 @@ public: ...@@ -112,6 +112,56 @@ public:
return qp; return qp;
} }
//void save_as(std::string path, std::string file_name) const;
void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const {
if (!boost::filesystem::exists(path)) {
boost::filesystem::create_directories(path);
}
std::fstream fs;
fs.open(path + file_name + ".oesc", std::fstream::out);
// Write header
size_t loc = 3;
fs << loc << ',' << Triangles.size() << std::endl;
loc += Triangles.size();
fs << loc << ',' << Nodes.size() << std::endl;
loc += Nodes.size();
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) {
fs << t.node(i) << ',';
}
fs << t.node(Triangle<Order>::NumNodes - 1) << std::endl;
}
// Write nodes
for (XY const &n : Nodes) {
fs << n.x() << ',' << n.y() << std::endl;
}
// Write values
for (size_t i = 0; i!=v.size(); ++i){
fs << v(i) << std::endl;
}
/*
for (size_t e = 0; e != Edges.size(); ++e) {
Point const v0 = base(e);
Point const v1 = base(next(e));
Point const v2 = base(next(next(e)));
fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n';
}
*/
fs.close();
}
protected: protected:
std::vector<XY> Nodes; std::vector<XY> Nodes;
std::vector<Triangle<Order>> Triangles; std::vector<Triangle<Order>> Triangles;
......
...@@ -62,7 +62,7 @@ public: ...@@ -62,7 +62,7 @@ public:
std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const; std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const;
protected: protected:
size_t Node[NumNodes]; size_t Node[NumNodes]; // TODO: std::array instead?
}; };
/* /*
......
...@@ -139,6 +139,8 @@ TEST(Library_Integration, Full_Circle_Uniform_Current_Density) { ...@@ -139,6 +139,8 @@ TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
} }
TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) { TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
GTEST_NONFATAL_FAILURE_("TODO: antiperiodic boundary conditions");
// Create Sketch // Create Sketch
Sketch sk; Sketch sk;
...@@ -320,6 +322,8 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) { ...@@ -320,6 +322,8 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
// A = a_h * log(r) + c_h in homogeneous outer region // A = a_h * log(r) + c_h in homogeneous outer region
Eigen::VectorXd vv = msph.recover_boundary(v); Eigen::VectorXd vv = msph.recover_boundary(v);
fem.write_scalar(vv, SAVE_DIR, std::string("quarter_circle_mvp"));
double_t a_f = -0.5; double_t a_f = -0.5;
double_t c_f = -2.0 * a_f * std::log(2.0) - a_f; double_t c_f = -2.0 * a_f * std::log(2.0) - a_f;
double_t a_h = 2 * a_f; double_t a_h = 2 * a_f;
......
...@@ -872,6 +872,165 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) { ...@@ -872,6 +872,165 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
ms.add_magnetic_insulation({outer_boundary}); ms.add_magnetic_insulation({outer_boundary});
variable_map[0].value(-1.0);
variable_map[1].value(-1.0);
variable_map.push_back(VariableMap(0, 0, -1.0));
ms.add_periodic_boundary(variable_map, true);
ms.build_matrices();
ms.apply_conditions();
// Initialize matrices
auto J = ms.init_unknown_matrix();
auto v = ms.init_unknown_vector();
auto r = ms.init_unknown_vector();
auto f = ms.init_unknown_vector();
auto Fx = ms.init_element_array();
auto Fy = ms.init_element_array();
auto dFxdx = ms.init_element_array();
auto dFydy = ms.init_element_array();
auto dFxdy = ms.init_element_array();
v.setZero();
ms.calculate_forcing(f);
EXPECT_NEAR(f(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(f(1), M_SQRT2 / 6.0, FLT_EPSILON);
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
// Verify equation is solved
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
for (size_t i = 0; i != r.size(); ++i) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
std::vector<double> v_expected{0.0, 6.60365e-8};
for (size_t i = 0; i != 2; ++i) {
EXPECT_NEAR(v[i], v_expected[i], 1e-13);
}
// Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
Eigen::ArrayXd Bmag(By.size());
Eigen::ArrayXd Bang(By.size());
for (size_t i = 0; i != By.size(); ++i) {
Bmag(i) = hypot(Bx(i), By(i));
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
}
std::vector<double> b_mag_expected{9.33897e-8, 9.33897e-8, 9.33897e-8, 7.14774e-8, 7.14774e-8, 9.33897e-8};
std::vector<double> b_ang_expected{0.0, -90, 0.0, 112.5, 157.5, -90.0};
for (size_t i = 0; i != 6; ++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-13);
EXPECT_NEAR(Bang[i], b_ang_expected[i], FLT_EPSILON);
}
}
TEST_F(TwoRegionQuarterOctagon, magnetostatic_periodic_rotational_motion) {
GTEST_FATAL_FAILURE_("TODO: Add Motion");
// Initialize physics
Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
ms.time(0.0);
// Add conditions
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
ms.add_periodic_boundary(variable_map, false);
ms.build_matrices();
ms.apply_conditions();
// Initialize matrices
auto J = ms.init_unknown_matrix();
auto v = ms.init_unknown_vector();
auto r = ms.init_unknown_vector();
auto f = ms.init_unknown_vector();
auto Fx = ms.init_element_array();
auto Fy = ms.init_element_array();
auto dFxdx = ms.init_element_array();
auto dFydy = ms.init_element_array();
auto dFxdy = ms.init_element_array();
v.setZero();
ms.calculate_forcing(f);
double nodal_current{M_SQRT2 / 6.0};
for (size_t i = 0; i != 3; ++i) {
EXPECT_NEAR(f(i), nodal_current, FLT_EPSILON);
}
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
// Verify equation is solved
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
for (size_t i = 0; i != r.size(); ++i) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
// Fundamental solution is 0.25*r^2 + C in forced region (0)
EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * 0.1);
for (size_t i = 1; i != 3; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
// Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
Eigen::ArrayXd Bmag(By.size());
Eigen::ArrayXd Bang(By.size());
for (size_t i = 0; i != By.size(); ++i) {
Bmag(i) = hypot(Bx(i), By(i));
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
}
std::vector<double> b_mag_expected(6, 3.86994e-7);
std::vector<double> b_ang_expected{112.5, 157.5, 112.5, 112.5, 157.5, 157.5};
for (size_t i = 0; i != 6; ++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-12);
EXPECT_NEAR(Bang[i], b_ang_expected[i], 1e-1);
}
}
TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic_rotational_motion) {
GTEST_FATAL_FAILURE_("TODO: Add Motion");
// Initialize physics
Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
ms.time(0.0);
// Add conditions
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
variable_map.push_back(VariableMap(0, 0, -1.0)); variable_map.push_back(VariableMap(0, 0, -1.0));
ms.add_periodic_boundary(variable_map, true); ms.add_periodic_boundary(variable_map, true);
...@@ -1165,5 +1324,3 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) { ...@@ -1165,5 +1324,3 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
++*position; ++*position;
} }
} }
\ No newline at end of file
// TODO: Test rotation + periodic boundaries
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment