Commit 0e76bba5 authored by JasonPries's avatar JasonPries
Browse files

Partial commit because I'm nervous about losing this work

Work on rotational motion and related boundary refinement
Output for vector plots
Starts on nonlinear material properties
Begin encapsulating solver routines
parent 27c32224
...@@ -2,6 +2,8 @@ import os ...@@ -2,6 +2,8 @@ import os
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from math import *
imgtypes = ['.pdf','.png'] imgtypes = ['.pdf','.png']
def plot_sketch(fname): def plot_sketch(fname):
...@@ -52,6 +54,39 @@ def plot_scalar(fname): ...@@ -52,6 +54,39 @@ def plot_scalar(fname):
plt.savefig(fname+ext) plt.savefig(fname+ext)
plt.clf() plt.clf()
def plot_vector(fname):
fname += '.oeve'
print('Vector: '+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])
mag = vals[:,0].copy()
for i in range(0,len(vals)):
mag[i] = hypot(vals[i,0], vals[i,1])
xc = (nodes[tris[:,0],0]+nodes[tris[:,1],0]+nodes[tris[:,2],0])/3.0
yc = (nodes[tris[:,0],1]+nodes[tris[:,1],1]+nodes[tris[:,2],1])/3.0
plt.figure()
plt.set_cmap('viridis')
plt.gca().set_aspect('equal')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,facecolors = mag,edgecolors='k')
plt.colorbar()
#plt.quiver(xc,yc,vals[:,0],vals[:,1],pivot='mid')
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/'
...@@ -65,3 +100,5 @@ for path, _, files in os.walk(root): ...@@ -65,3 +100,5 @@ for path, _, files in os.walk(root):
plot_mesh(fpath) plot_mesh(fpath)
elif fext == '.oesc': elif fext == '.oesc':
plot_scalar(fpath) plot_scalar(fpath)
elif fext == '.oeve':
plot_vector(fpath)
...@@ -29,14 +29,15 @@ double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const { ...@@ -29,14 +29,15 @@ double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const { double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
std::cout << "//TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued" << std::endl; std::cout << "//TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued" << std::endl;
std::cout << "//TODO: Constrained dart enqueuing needs to be revisited, and need to make sure enqueued edges are unique" << std::endl;
double delta_max{delta_min}; double delta_max{delta_min};
for (size_t i : DartConstraints) { for (size_t i : DartConstraints) {
DartConstraint const &dc = m.dart_constraint(i); DartConstraint const &dc = m.dart_constraint(i);
double_t delta = std::abs(dc.S1 - dc.S0); double_t delta = std::abs(dc.S1 - dc.S0);
if (delta > delta_min) { if (delta > delta_min) {
m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart())); //m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
//m.add_dart_to_queue(dc.dart()); m.add_dart_to_queue(dc.dart());
delta_max = std::max(delta_max, delta); delta_max = std::max(delta_max, delta);
} }
} }
...@@ -48,8 +49,8 @@ void BoundaryConstraint::add_to_queue(Mesh &m, size_t dart) const { ...@@ -48,8 +49,8 @@ void BoundaryConstraint::add_to_queue(Mesh &m, size_t dart) const {
if (UniformDiscretization) { if (UniformDiscretization) {
queue_uniform(m, smallest_parametric_edge_length(m) / 2.0); queue_uniform(m, smallest_parametric_edge_length(m) / 2.0);
} else { } else {
m.add_dart_to_queue(dart); //m.add_dart_to_queue(dart);
//m.add_to_queue(std::make_unique<MidpointQueuer>(dart)); m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
} }
} }
......
...@@ -906,6 +906,18 @@ void Mesh::save_as(std::string path, std::string file_name) const { ...@@ -906,6 +906,18 @@ void Mesh::save_as(std::string path, std::string file_name) const {
fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n'; fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n';
} }
/*
for (auto t : Triangles) {
size_t e{t.edge()};
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(); fs.close();
} }
...@@ -1594,7 +1606,6 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) { ...@@ -1594,7 +1606,6 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) {
Edges[e].add_to_queue(*this); Edges[e].add_to_queue(*this);
} }
return AddToQueueResult::Midpoint; return AddToQueueResult::Midpoint;
} else if (result == LocateTriangleResult::Interior) { } else if (result == LocateTriangleResult::Interior) {
Queue.push_back(std::make_unique<CircumcenterQueuer>(vc, ei)); Queue.push_back(std::make_unique<CircumcenterQueuer>(vc, ei));
return AddToQueueResult::Circumcenter; return AddToQueueResult::Circumcenter;
......
...@@ -135,6 +135,19 @@ public: ...@@ -135,6 +135,19 @@ public:
return MappedBoundaries.back(); return MappedBoundaries.back();
} }
template<class T>
std::vector<std::shared_ptr<MappedBoundaryPair>> add_mapped_boundary_pair(std::vector<T> const &cbp_vec) {
/*
* Constructs a MappedBoundaryPair for each ContinuousBoundaryPair in cbp_vec
*/
for(T const &cbp : cbp_vec) {
MappedBoundaries.push_back(std::make_shared<MappedBoundaryPair>(*this, cbp));
}
return std::vector<std::shared_ptr<MappedBoundaryPair>>(MappedBoundaries.end() - cbp_vec.size(), MappedBoundaries.end());
}
DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; }; DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
DartTriangle const &triangle(size_t i) const { return Triangles[i]; }; DartTriangle const &triangle(size_t i) const { return Triangles[i]; };
......
...@@ -21,8 +21,11 @@ set(SOURCE_FILES ...@@ -21,8 +21,11 @@ set(SOURCE_FILES
./src/Physics.h ./src/Physics.cpp ./src/Physics.h ./src/Physics.cpp
./src/PhysicsConstants.h
./src/Triangle.h ./src/Triangle.cpp ./src/Triangle.h ./src/Triangle.cpp
)
src/Solution.cpp src/Solution.h)
add_library(physics SHARED ${SOURCE_FILES}) add_library(physics SHARED ${SOURCE_FILES})
......
...@@ -315,6 +315,45 @@ public: ...@@ -315,6 +315,45 @@ public:
fs.close(); fs.close();
} }
void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, 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 + ".oeve", 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 << ',' << Fx.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 != Fx.size(); ++i) {
fs << Fx(i) << ',' << Fy(i) << std::endl;
}
fs.close();
}
protected: protected:
std::vector<XY> Nodes; std::vector<XY> Nodes;
std::vector<Triangle<Order>> Triangles; std::vector<Triangle<Order>> Triangles;
......
...@@ -50,6 +50,78 @@ protected: ...@@ -50,6 +50,78 @@ protected:
std::shared_ptr<MagneticMaterialInterface> MagneticProperties; std::shared_ptr<MagneticMaterialInterface> MagneticProperties;
}; };
class NonlinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
NonlinearIsotropicMagneticMaterial() {
dB = 1.0;
BSaturation = 2.0;
double_t a0 = 0.0;
double_t b0 = 0.999 / mu_0;
double_t c0 = 0.0;
BHSpline.push_back(std::array<double_t, 3>{a0,b0,c0});
double_t a1 = -b0 / 2.0;
double_t b1 = -4.0 * a1;
double_t c1 = b0 - a1 - b1;
BHSpline.push_back(std::array<double_t, 3>{a1,b1,c1});
MSaturation = (a1 * BSaturation + b1) * BSaturation + c1;
};
void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) override {
for(size_t const &i : index) {
const double_t &Bx = Fx(i);
const double_t &By = Fy(i);
double_t B = sqrt(Bx * Bx + By * By);
double_t M, dMdB;
if (B > BSaturation) {
M = MSaturation;
dMdB = 0.0;
} else {
size_t j = (size_t)(B / dB);
const auto &s = BHSpline[j];
M = (s[0] * B + s[1]) * B + s[2];
dMdB = (2.0 * s[0] * B + s[1]);
}
if (B > 0) {
double_t Mx = Bx * (M / B);
double_t My = By * (M / B);
double_t dMxdBx = M / B + (Bx / B) * (Bx / B) * (dMdB - M / B);
double_t dMydBy = M / B + (By / B) * (By / B) * (dMdB - M / B);
double_t dMxdBy = (Bx / B) * (By / B) * (dMdB - M / B);
dFxdx(i) = 1.0 / mu_0 - dMxdBx;
dFydy(i) = 1.0 / mu_0 - dMydBy;
dFxdy(i) = -dMxdBy;
Fx(i) = Bx / mu_0 - Mx;
Fy(i) = By / mu_0 - My;
} else {
dFxdx(i) = 1.0 / mu_0 - dMdB;
dFydy(i) = 1.0 / mu_0 - dMdB;
dFxdy(i) = 0.0;
Fx(i) = 0.0;
Fy(i) = 0.0;
}
}
}
protected:
std::vector<std::array<double_t,3>> BHSpline;
double_t dB;
double_t MSaturation;
double_t BSaturation;
};
inline MaterialProperties Air() { inline MaterialProperties Air() {
return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0)); return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
} }
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "BoundaryCondition.h" #include "BoundaryCondition.h"
#include "FiniteElementMesh.h" #include "FiniteElementMesh.h"
#include "Forcing.h" #include "Forcing.h"
#include "Solution.h"
#include "Node.h" #include "Node.h"
#include "PhysicsConstants.h" #include "PhysicsConstants.h"
#include "Triangle.h" #include "Triangle.h"
...@@ -110,6 +111,32 @@ public: ...@@ -110,6 +111,32 @@ public:
Eigen::DiagonalMatrix<double, Eigen::Dynamic> init_element_matrix() const { return Eigen::DiagonalMatrix<double, Eigen::Dynamic>(Elements); } Eigen::DiagonalMatrix<double, Eigen::Dynamic> init_element_matrix() const { return Eigen::DiagonalMatrix<double, Eigen::Dynamic>(Elements); }
std::shared_ptr<Solution> initialize() {
auto s = std::make_shared<Solution>();
s->J = init_unknown_matrix(); // Jacobian
s->v = init_unknown_vector(); // solution
s->v.setZero();
s->r = init_unknown_vector(); // residual
s->f = init_unknown_vector(); // forcing
s->Fx = init_element_array();
s->Fy = init_element_array();
s->dFxdGx = init_element_array();
s->dFydGy = init_element_array();
s->dFxdGy = init_element_array();
return s;
}
void assemble() {
build_matrices();
apply_conditions();
}
Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); }; Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); };
Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); }; Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); };
...@@ -118,8 +145,34 @@ public: ...@@ -118,8 +145,34 @@ public:
Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(Elements); }; Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(Elements); };
void solve(std::shared_ptr<Solution> s) {
// TODO: Template on solver type
// TODO: Save J,v,r,f,Fx,Fy,dFxdx,dFydy,dFxdy to reduce number of memory allocations?
// TODO: This method is fairly general
calculate_forcing(s->f);
linearize(s->J, s->v, s->r, s->f, s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);
for (size_t iter = 0; iter != 20; ++iter) {
s->Solver.compute(s->J);
if (s->Solver.info() != Eigen::Success) {
std::cerr << "Factorization of Jacobian failed";
}
s->v -= 1.0 * s->Solver.solve(s->r);
if (s->Solver.info() != Eigen::Success) {
std::cerr << "Solution of linear system failed";
}
linearize(s->J, s->v, s->r, s->f, s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);
std::cout << s->r.norm() << std::endl;
}
}
void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f, void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f,
Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) { Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxGx, Eigen::ArrayXd &dFyGy, Eigen::ArrayXd &dFxGy) {
r = -f; r = -f;
J.setZero(); J.setZero();
...@@ -128,9 +181,9 @@ public: ...@@ -128,9 +181,9 @@ public:
Fx = Derivatives.dy(i).transpose() * v; Fx = Derivatives.dy(i).transpose() * v;
Fy = -Derivatives.dx(i).transpose() * v; Fy = -Derivatives.dx(i).transpose() * v;
// Calculate polarization // Calculate field intensity
for (auto &dr : Domain.regions()) { for (auto &dr : Domain.regions()) {
dr->material().h_from_b(dr->elements(), Fx, Fy, dFxdx, dFydy, dFxdy); dr->material().h_from_b(dr->elements(), Fx, Fy, dFxGx, dFyGy, dFxGy);
} }
// Calculate residual // Calculate residual
...@@ -138,10 +191,10 @@ public: ...@@ -138,10 +191,10 @@ public:
- Derivatives.dx(i) * (ElementWeights(i) * Fy).matrix(); // note: weak form introduces a negative sign into double-curl operator - Derivatives.dx(i) * (ElementWeights(i) * Fy).matrix(); // note: weak form introduces a negative sign into double-curl operator
// Calculate jacobian // Calculate jacobian
J += Derivatives.dx(i) * (ElementWeights(i) * dFydy).matrix().asDiagonal() * Derivatives.dx(i).transpose() J += Derivatives.dx(i) * (ElementWeights(i) * dFyGy).matrix().asDiagonal() * Derivatives.dx(i).transpose()
+ Derivatives.dy(i) * (ElementWeights(i) * dFxdx).matrix().asDiagonal() * Derivatives.dy(i).transpose() + Derivatives.dy(i) * (ElementWeights(i) * dFxGx).matrix().asDiagonal() * Derivatives.dy(i).transpose()
+ Derivatives.dx(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dy(i).transpose() - Derivatives.dx(i) * (ElementWeights(i) * dFxGy).matrix().asDiagonal() * Derivatives.dy(i).transpose()
+ Derivatives.dy(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dx(i).transpose(); - Derivatives.dy(i) * (ElementWeights(i) * dFxGy).matrix().asDiagonal() * Derivatives.dx(i).transpose();
} }
}; };
......
#include "Solution.h"
\ No newline at end of file
#ifndef OERSTED_SOLUTION_H
#define OERSTED_SOLUTION_H
#include "Eigen"
#include "Eigen/IterativeLinearSolvers"
#include "Eigen/SparseCholesky"
class Solution {
public:
Eigen::SparseMatrix<double_t> J;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double_t>> Solver;
Eigen::VectorXd v;
Eigen::VectorXd r;
Eigen::VectorXd f;
Eigen::ArrayXd Fx;
Eigen::ArrayXd Fy;
Eigen::ArrayXd dFxdGx;
Eigen::ArrayXd dFydGy;
Eigen::ArrayXd dFxdGy;
};
#endif //OERSTED_SOLUTION_H
...@@ -54,7 +54,6 @@ MirrorCopy::MirrorCopy(std::vector<std::shared_ptr<Curve const>> &input, std::sh ...@@ -54,7 +54,6 @@ MirrorCopy::MirrorCopy(std::vector<std::shared_ptr<Curve const>> &input, std::sh
//Verticies.push_back(new Vertex(px, py)); //Verticies.push_back(new Vertex(px, py));
Verticies.push_back(std::make_shared<Vertex>(px, py)); Verticies.push_back(std::make_shared<Vertex>(px, py));
++v; ++v;
} else { } else {
v = input_vlist.erase(v); v = input_vlist.erase(v);
......
...@@ -20,6 +20,8 @@ public: ...@@ -20,6 +20,8 @@ public:
void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) const override {}; void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) const override {};
std::vector<std::shared_ptr<Curve>> const &curves() const { return Curves; };
protected: protected:
std::vector<std::shared_ptr<Curve const>> Input; std::vector<std::shared_ptr<Curve const>> Input;
......
...@@ -2,6 +2,146 @@ ...@@ -2,6 +2,146 @@
std::string SAVE_DIR = "./test/output/Integration/"; std::string SAVE_DIR = "./test/output/Integration/";
TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
// Create Sketch
Sketch sk;
auto v0 = sk.new_element<Vertex>(0.0, 0.0);
auto v1 = sk.new_element<Vertex>(1.0, 0.0);
auto v2 = sk.new_element<Vertex>(1.0 * std::cos(M_PI * 2.0 / 3.0), 1.0 * std::sin(M_PI * 2.0 / 3.0));
auto v3 = sk.new_element<Vertex>(1.0 * std::cos(-M_PI * 2.0 / 3.0), 1.0 * std::sin(-M_PI * 2.0 / 3.0));
auto c0 = sk.new_element<CircularArc>(v1, v2, v0, 1.0);
auto c1 = sk.new_element<CircularArc>(v2, v3, v0, 1.0);
auto c2 = sk.new_element<CircularArc>(v3, v1, v0, 1.0);
auto f0 = sk.new_element<Fixation>(v0);
auto f1 = sk.new_element<Fixation>(v1);
auto f2 = sk.new_element<Fixation>(v2);
auto f3 = sk.new_element<Fixation>(v3);
double_t tol = sk.solve();
EXPECT_LE(tol, FLT_EPSILON);
bool result = sk.build();
EXPECT_TRUE(result);
EXPECT_EQ(sk.size_contours(), 1);
sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("full_circle_sketch"));
// Create Refineable Mesh
Mesh me{sk};
me.create();
me.MinimumElementSize = 0.01;
me.MaximumElementSize = 0.1;
me.MinimumElementQuality = 0.5;
me.refine();
me.save_as(SAVE_DIR, std::string("fulL_circle_mesh"));
// Convert to FiniteElementMesh
FiniteElementMesh<2, 1> fem{me};
for (std::shared_ptr<DiscreteBoundary<2>> const &b : fem.boundaries()) {
double_t a0{-2 * M_PI};
for (size_t i : b->nodes()) {
XY const &p = fem.node(i);
// Test radius of boundary
EXPECT_NEAR(1.0, std::hypot(p.x(), p.y()), FLT_EPSILON);
// Test ordering of boundary nodes
double_t a1{std::atan2(p.y(), p.x())};
if (a1 < 0 || (a1 == 0 && a0 > M_PI)) {
a1 += 2 * M_PI;
}
EXPECT_TRUE(a1 > a0);
a0 = a1;
}
}
// Create magnetostatic physics
Magnetostatic<2, 1, 1, FieldVariable::A> msph{fem};
msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(0)});
msph.add_magnetic_insulation({fem.boundary(0), fem.boundary(1), fem.boundary(2)});
msph.build_matrices();
msph.apply_conditions();
// Initialize matrices
auto J = msph.init_unknown_matrix();
auto v = msph.init_unknown_vector();
auto r = msph.init_unknown_vector();