Commit cc6dd272 authored by JasonPries's avatar JasonPries
Browse files

Commit due to fubar kernel

parent 0e76bba5
......@@ -75,11 +75,11 @@ def plot_vector(fname):
yc = (nodes[tris[:,0],1]+nodes[tris[:,1],1]+nodes[tris[:,2],1])/3.0
plt.figure()
plt.set_cmap('viridis')
plt.set_cmap('inferno')
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.quiver(xc,yc,vals[:,0],vals[:,1],mag,pivot='mid',edgecolors='k',linewidths=0.1)
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
for ext in imgtypes:
......
This diff is collapsed.
......@@ -53,28 +53,34 @@ protected:
class NonlinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
NonlinearIsotropicMagneticMaterial() {
dB = 1.0;
DeltaB = 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});
MBSpline.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});
MBSpline.push_back(std::array<double_t, 3>{a1,b1,c1});
MSaturation = (a1 * BSaturation + b1) * BSaturation + c1;
};
NonlinearIsotropicMagneticMaterial(std::vector<std::array<double_t,3>> mbspline, double_t db, double_t bsat) : MBSpline(mbspline), DeltaB(db), BSaturation(bsat) {
const auto &s = MBSpline.back();
MSaturation = (s[0] * BSaturation + s[1]) * BSaturation + s[2];
}
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;
......@@ -83,8 +89,8 @@ public:
M = MSaturation;
dMdB = 0.0;
} else {
size_t j = (size_t)(B / dB);
const auto &s = BHSpline[j];
size_t j = (size_t)(B / DeltaB);
const auto &s = MBSpline[j];
M = (s[0] * B + s[1]) * B + s[2];
dMdB = (2.0 * s[0] * B + s[1]);
......@@ -116,8 +122,8 @@ public:
}
protected:
std::vector<std::array<double_t,3>> BHSpline;
double_t dB;
std::vector<std::array<double_t,3>> MBSpline;
double_t DeltaB;
double_t MSaturation;
double_t BSaturation;
};
......@@ -126,4 +132,6 @@ inline MaterialProperties Air() {
return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
}
MaterialProperties JFE_35JN210();
#endif //OERSTED_MATERIALPROPERTIES_H
......@@ -152,22 +152,29 @@ public:
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) {
size_t iter = 0;
double_t norm_dv = 1.0;
while (iter < 30 && norm_dv > sqrt(DBL_EPSILON) * std::sqrt(s->v.size())) {
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);
auto dv = s->Solver.solve(s->r);
if (s->Solver.info() != Eigen::Success) {
std::cerr << "Solution of linear system failed";
}
s->v -= dv;
norm_dv = dv.norm() / s->v.norm();
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;
iter++;
}
}
......
......@@ -153,7 +153,7 @@ public:
v0 = sk.new_element<Vertex>(0.0, 0.0);
v1 = sk.new_element<Vertex>(1.0, 0.0);
v2 = sk.new_element<Vertex>(2.0, 0.0);
v3 = sk.new_element<Vertex>(sqrt(3.0)/2.0, 0.5);
v3 = sk.new_element<Vertex>(sqrt(3.0) / 2.0, 0.5);
v4 = sk.new_element<Vertex>(sqrt(3.0), 1.0);
f0 = sk.new_element<Fixation>(v0);
......@@ -261,9 +261,9 @@ public:
me.save_as(SAVE_DIR, std::string("sixth_circle_mirror_copy_mesh"));
// Create FiniteElementMesh
fem = FiniteElementMesh<2,1>(me);
fem = FiniteElementMesh<2, 1>(me);
for (auto & pb : periodic_boundary) {
for (auto &pb : periodic_boundary) {
EXPECT_EQ(fem.boundary(pb.curve0()).size(), 1);
EXPECT_EQ(fem.boundary(pb.curve1()).size(), 1);
......@@ -288,7 +288,7 @@ public:
Sketch sk;
Mesh me;
FiniteElementMesh<2,1> fem;
FiniteElementMesh<2, 1> fem;
// Sketch Variables
std::shared_ptr<Vertex> v0;
......@@ -299,7 +299,7 @@ public:
std::shared_ptr<Fixation> f0;
std::shared_ptr<LineSegment> l01 ;
std::shared_ptr<LineSegment> l01;
std::shared_ptr<LineSegment> l12;
std::shared_ptr<LineSegment> l04;
......@@ -375,7 +375,7 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Periodic)
double_t c_f = -2.0 * a_f * std::log(2.0) - a_f;
double_t a_h = 2 * a_f;
double_t c_h = -2.0 * a_f * log(2.0);
for(size_t i = 0; i!= fem.size_nodes(); ++i) {
for (size_t i = 0; i != fem.size_nodes(); ++i) {
XY const &n = fem.node(i);
double_t r = std::hypot(n.x(), n.y());
if (r <= 1.0) {
......@@ -412,7 +412,7 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Periodic)
EXPECT_NEAR(Bmag(i), val, std::abs(2 * a_f * 0.05));
} else {
double_t val = std::abs(a_h / r);
EXPECT_NEAR(Bmag(i), val, std::abs(a_h* 0.05));
EXPECT_NEAR(Bmag(i), val, std::abs(a_h * 0.05));
}
double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;
......@@ -476,9 +476,9 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Antieriod
double_t c_h = a_f / 321.0;
double_t d_h = -64.0 / 321.0 * a_f;
double_t r_max = - (2.0 * b_f) / (3.0 * a_f);
double_t r_max = -(2.0 * b_f) / (3.0 * a_f);
double_t tol_a = (a_f * r_max + b_f) * r_max * r_max * 0.02;
for(size_t i = 0; i!= fem.size_nodes(); ++i) {
for (size_t i = 0; i != fem.size_nodes(); ++i) {
XY const &n = fem.node(i);
double_t r = std::hypot(n.x(), n.y());
double_t a = std::atan2(n.y(), n.x());
......@@ -614,9 +614,9 @@ public:
me.save_as(SAVE_DIR, std::string("sixth_circle_mesh"));
// Create FiniteElementMesh
fem = FiniteElementMesh<2,1>(me);
fem = FiniteElementMesh<2, 1>(me);
for (auto & pb : periodic_boundary) {
for (auto &pb : periodic_boundary) {
EXPECT_EQ(fem.boundary(pb.curve0()).size(), 1);
EXPECT_EQ(fem.boundary(pb.curve1()).size(), 1);
......@@ -639,7 +639,7 @@ public:
Sketch sk;
Mesh me;
FiniteElementMesh<2,1> fem;
FiniteElementMesh<2, 1> fem;
// Sketch Variables
std::shared_ptr<Vertex> v0;
......@@ -650,7 +650,7 @@ public:
std::shared_ptr<Fixation> f0;
std::shared_ptr<LineSegment> l01 ;
std::shared_ptr<LineSegment> l01;
std::shared_ptr<LineSegment> l12;
std::shared_ptr<LineSegment> l03;
......@@ -954,7 +954,8 @@ public:
// Stator
rsi = rro + airgap;
rso = 132e-3;
rsb = (rso + rsi) / 2.0;
rsb = (1.0 * rso + 1.0 * rsi) / 2.0;
tooth_width = rsi * (2.0 * M_PI / (4.0 * 3.0 * 2.0)) * (2.0 / 3.0);
auto vs0 = sk.new_element<Vertex>(rsi, 0.0);
auto vs1 = sk.new_element<Vertex>(rso, 0.0);
......@@ -970,8 +971,10 @@ public:
auto ls54 = sk.new_element<LineSegment>(vs5, vs4);
auto h01 = sk.new_element<Horizontal>(ls01);
auto h54 = sk.new_element<Horizontal>(ls54);
std::cout << "//TODO: Fix distance between h54 and h01. h54 will not need the horizontal constraint" << std::endl;
auto d01_54 = sk.new_element<Distance<LineSegment>>(ls01, ls54, tooth_width);
//auto h54 = sk.new_element<Horizontal>(ls54);
//std::cout << "//TODO: Fix distance between l54 and l01. l54 will not need the horizontal constraint" << std::endl;
auto a63 = sk.new_element<Angle>(ls01, ls63, 360.0 / 24.0);
auto a32 = sk.new_element<Angle>(ls01, ls32, 360.0 / 24.0);
......@@ -1043,18 +1046,18 @@ public:
// Create initial mesh
me.create();
me.MaximumElementSize = rro;
me.MaximumElementSize = 5.0e-3;
me.MinimumElementSize = airgap;
me.MinimumElementQuality = 0.5;
// Refine and save
me.refine();
me.save_as(SAVE_DIR, std::string("saleint_pole_synchrel_mesh"));
me.save_as(SAVE_DIR, std::string("salient_pole_synchrel_mesh"));
// Create FiniteElementMesh
fem = FiniteElementMesh<2,1>(me);
fem = FiniteElementMesh<2, 1>(me);
for (auto & pb : periodic_boundary) {
for (auto &pb : periodic_boundary) {
EXPECT_EQ(fem.boundary(pb.curve0()).size(), 1);
EXPECT_EQ(fem.boundary(pb.curve1()).size(), 1);
......@@ -1118,12 +1121,14 @@ public:
double_t airgap;
double_t tooth_width;
std::shared_ptr<Vertex> origin;
std::shared_ptr<CircularArc> continuous_airgap;
Sketch sk;
Mesh me;
FiniteElementMesh<2,1> fem;
FiniteElementMesh<2, 1> fem;
std::shared_ptr<Contour const> phase_a_contour;
std::shared_ptr<Contour const> phase_b_contour;
......@@ -1140,18 +1145,17 @@ TEST_F(Salient_Pole_Synchrel, Test) {
msph.time(0.0);
// Set material properties
//MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
//fem.region(rotor_iron)->material(Linear1000);
//fem.region(stator_iron)->material(Linear1000);
//MaterialProperties electrical_steel = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
//MaterialProperties electrical_steel = MaterialProperties(std::make_shared<NonlinearIsotropicMagneticMaterial>());
MaterialProperties electrical_steel = JFE_35JN210();
MaterialProperties Nonlinear796 = MaterialProperties(std::make_shared<NonlinearIsotropicMagneticMaterial>());
fem.region(rotor_iron)->material(Nonlinear796);
fem.region(stator_iron)->material(Nonlinear796);
fem.region(rotor_iron)->material(electrical_steel);
fem.region(stator_iron)->material(electrical_steel);
// Conditions
msph.add_current_density([](double t) { return -30e6 * sin(2*M_PI*t - M_PI*2.0/3.0); }, {phase_a_contour});
msph.add_current_density([](double t) { return 30e6 * sin(2*M_PI*t); }, {phase_b_contour});
msph.add_current_density([](double t) { return -30e6 * sin(2*M_PI*t + M_PI*2.0/3.0); }, {phase_c_contour});
msph.add_current_density([](double t) { return -2e6 * sin(2 * M_PI * t - M_PI * 2.0 / 3.0); }, {phase_a_contour});
msph.add_current_density([](double t) { return 2e6 * sin(2 * M_PI * t); }, {phase_b_contour});
msph.add_current_density([](double t) { return -2e6 * sin(2 * M_PI * t + M_PI * 2.0 / 3.0); }, {phase_c_contour});
msph.add_magnetic_insulation(outer_boundary);
......@@ -1163,31 +1167,27 @@ TEST_F(Salient_Pole_Synchrel, Test) {
msph.assemble();
// Solve
std::cout << "//TODO: msph.time() += N * dt, *position += N (see below)" << std::endl;
auto solution = msph.initialize();
for(size_t iter = 0; iter != 1; ++iter){ //position->size(); ++iter) {
for (size_t iter = 0; iter != position->size() / 32; ++iter) {
msph.solve(solution);
// Verify equation is solved
//for (size_t i = 0; i != solution->r.size(); ++i) {
// EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
//}
//Verify equation is solved
EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
// Save image
Eigen::VectorXd vv = msph.recover_boundary(solution->v);
fem.write_scalar(vv, SAVE_DIR, std::string("salient_pole_syncrel_A_") + std::to_string(iter));
Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * solution->v;
Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * solution->v;
fem.write_vector(Bx, By, SAVE_DIR, std::string("salient_pole_synchrel_B_") + std::to_string(iter));
//std::cout << Bx << std::endl;
// Increment Position
for (size_t i = 0; i != 2; ++i) {
// TODO: msph.time() += N * dt, *position += N
for (size_t i = 0; i != 32; ++i) {
++*position;
msph.time(msph.time() + dt);
std::cout << "//TODO: msph.time() += N * dt, *position += N" << std::endl;
}
msph.assemble();
}
......
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