Commit 8964ade8 authored by Pries, Jason's avatar Pries, Jason
Browse files

Work on Synchronous Reluctance Machine template meshing issues

parent 5fbe4a1a
......@@ -4,7 +4,7 @@ project(Oersted)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --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/)
......
import sys
import os
import numpy as np
import matplotlib.pyplot as plt
......@@ -86,19 +87,19 @@ def plot_vector(fname):
plt.savefig(fname+ext)
plt.clf()
root = os.getcwd()
root += '/build/test/output/'
#root += '/build/test/output/'
for path, _, files in os.walk(root):
for name in files:
name, fext = os.path.splitext(name)
fpath = os.path.join(path, name)
if fext == '.oesk':
plot_sketch(fpath)
elif fext == '.oeme':
plot_mesh(fpath)
elif fext == '.oesc':
plot_scalar(fpath)
elif fext == '.oeve':
plot_vector(fpath)
if len(sys.argv) == 1 or name.startswith(sys.argv[1]):
fpath = os.path.join(path, name)
if fext == '.oesk':
plot_sketch(fpath)
elif fext == '.oeme':
plot_mesh(fpath)
elif fext == '.oesc':
plot_scalar(fpath)
elif fext == '.oeve':
plot_vector(fpath)
......@@ -763,8 +763,11 @@ void Mesh::insert_internal_boundaries() {
// TODO: This method may fail for some situations with poor aspect ratios (small airgaps)
// TODO: Something to try: search for internal boundaries with two existing end points and insert those first by decreasing length
// TODO: Then repeat for internal boundaries with one existing end point
// TODO: Or just add a generic bounding box + insert
sort_constraints_by_length();
sort_constraints_by_length(); // TODO: Method is sensitive to sorting of constraints
// TODO: Method is also sensitive to lines connecting boundary verticies intersecting with internal features
// TODO: I.E., straight line approximations of the geometry hull must not intersect internal boundaries
// Find interior curves
std::vector<size_t> interior_index;
......@@ -818,7 +821,32 @@ void Mesh::insert_internal_boundaries() {
LocateTriangleResult result = locate_triangle(p0, ei);
if (result != LocateTriangleResult::Point) {
throw std::exception();
/*
std::cerr << "Could not locate triangle referencing point, performing brute force search" << std::endl;
for (auto p : Points) {
double_t d = hypot(p.X-p0.X,p.Y-p0.Y);
//std::cout << d << std::endl;
}
bool point_found = false;
for (size_t i = 0; i != Edges.size(); ++i) {
Point pe = point(Edges[i]);
double_t d = hypot(pe.X-p0.X,pe.Y-p0.Y);
if (d < DBL_MIN * 10.0) {
std::cout << d << std::endl;
ei = i;
point_found = true;
break;
}
}
if(!point_found) {
std::cerr << "Could not locate edge referencing point, aborting mesh" << std::endl;
return;
}
*/
std::cerr << "Could not locate edge referencing point, aborting mesh" << std::endl;
return;
}
if (find_attached(p1, ei)) {
......@@ -906,7 +934,7 @@ void Mesh::mark_triangles() {
}
}
void Mesh::refine_once(std::vector<size_t> index, std::vector<double> radii, std::vector<double> quality) {
void Mesh::refine_once(std::vector<size_t> index, std::vector<double_t> radii, std::vector<double_t> quality) {
for (size_t i = 0; i != Triangles.size(); ++i) {
size_t j = index[i];
if ((edge_from_triangle_index(j).Mark) && ((radii[j] > MaximumElementSize) || (radii[j] > MinimumElementSize && quality[j] < MinimumElementQuality))) {
......@@ -1055,8 +1083,8 @@ void Mesh::sort_constraints_by_length() {
* Sorts BoundaryConstraints by length for internal boundary insertion
*/
//auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); };
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); };
//auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); }; //ascending
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); }; //descending
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(), comp);
};
......@@ -1184,7 +1212,7 @@ void Mesh::add_encroached_edges_to_queue_and_insert() {
void Mesh::triangulate_boundary_polygon() {
/*
* Performs intial triangulation of the mesh boundary
* Performs initial triangulation of the mesh boundary
*/
Edges.reserve(3 * num_points());
......
......@@ -4,7 +4,10 @@ set(SOURCE_FILES
./include/ModelTemplates.hpp
./src/DistributedWindingStator.h ./src/DistributedWindingStator.cpp
src/SynchronousReluctanceRotor.cpp src/SynchronousReluctanceRotor.h src/ModelTemplate.h)
./src/SynchronousReluctanceRotor.cpp ./src/SynchronousReluctanceRotor.h
./src/ModelTemplate.h ./src/ModelTemplate.cpp
./src/PolarModelTemplate.cpp ./src/PolarModelTemplate.h
)
add_library(model_templates SHARED ${SOURCE_FILES})
......
......@@ -87,8 +87,10 @@ void DistributedWindingStator::add_to_sketch(Sketch &s, std::shared_ptr<Vertex>
std::shared_ptr<RotateCopy> rc = s.new_element<RotateCopy>(copy_curves,origin,at,ModeledTeeth-1,true);
if (Convexify) {
auto V0p = s.select_periodic_vertex(V[0],origin, (360.0 * ModeledTeeth) / TotalTeeth);
s.new_element<LineSegment>(origin, V[0]);
s.new_element<LineSegment>(origin, V0p);
convexify(s,origin,V[0], ModeledTeeth * (360.0 / TotalTeeth));
}
if (AddBoundingBox) {
add_bounding_box(s, origin, V[1], ro, ModeledTeeth * (360.0 / TotalTeeth));
}
}
\ No newline at end of file
#ifndef OERSTED_DISTRIBUTEDWINDINGSTATOR_H
#define OERSTED_DISTRIBUTEDWINDINGSTATOR_H
#include "ModelTemplate.h"
#include "PolarModelTemplate.h"
class DistributedWindingStator : public ModelTemplate {
class DistributedWindingStator : public PolarModelTemplate {
public:
void add_to_sketch(Sketch &s, std::shared_ptr<Vertex> origin) override;
......
#include "ModelTemplate.h"
\ No newline at end of file
#include "PolarModelTemplate.h"
void PolarModelTemplate::convexify(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp0, double_t angle) {
auto vp1 = s.select_periodic_vertex(vp0,origin, angle);
s.new_element<LineSegment>(origin, vp0);
s.new_element<LineSegment>(origin, vp1);
}
void PolarModelTemplate::add_bounding_box(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp00, double_t radius, double_t angle) {
double_t rb = radius * std::min(2.0, std::abs(sqrt(8.0 / (1.0 + std::cos(angle * M_PI / 180.0))) - 1.0));
auto vp10 = s.select_periodic_vertex(vp00, origin, angle);
auto vp01 = s.new_element<Vertex>(rb, 0.0);
auto vp11 = s.new_element<Vertex>(rb * std::cos(angle * M_PI / 180.0), rb * std::sin(angle * 180.0 / M_PI));
auto l0 = s.new_element<LineSegment>(vp00, vp01);
auto l1 = s.new_element<LineSegment>(vp10, vp11);
auto co = s.new_element<CircularArc>(vp01, vp11, origin, rb);
s.new_element<Radius>(co, rb);
s.new_element<Horizontal>(l0);
s.new_element<Angle>(l0, l1, angle);
s.new_element<Coincident<LineSegment>>(origin, l1);
}
\ No newline at end of file
#ifndef OERSTED_POLARMODELTEMPLATE_H
#define OERSTED_POLARMODELTEMPLATE_H
#include "ModelTemplate.h"
class PolarModelTemplate : public ModelTemplate {
public:
void convexify(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp, double_t angle);
void add_bounding_box(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp, double_t radius, double_t angle);
};
#endif //OERSTED_POLARMODELTEMPLATE_H
\ No newline at end of file
......@@ -171,26 +171,10 @@ void SynchronousReluctanceRotor::add_to_sketch(Sketch &s, std::shared_ptr<Vertex
tol = s.solve();
if (Convexify) {
auto V0p = s.select_periodic_vertex(V[0],origin, 360.0 / Poles);
s.new_element<LineSegment>(origin, V[0]);
s.new_element<LineSegment>(origin, V0p);
convexify(s,origin,V[0],360.0/Poles);
}
if (AddBoundingBox) {
double_t rb = ro * std::min(2.0, std::abs(sqrt(8.0 / (1.0 + std::cos(M_PI * 2.0 / Poles))) - 1.0));
auto v1p = s.select_periodic_vertex(V[1],origin,360.0/Poles);
auto v11 = s.new_element<Vertex>(rb,0.0);
auto v11p = s.new_element<Vertex>(rb * std::cos(M_PI * 2.0 / Poles), rb * std::sin(M_PI * 2.0 / Poles));
auto l0 = s.new_element<LineSegment>(V[1], v11);
auto l1 = s.new_element<LineSegment>(v1p, v11p);
auto co = s.new_element<CircularArc>(v11,v11p,origin,rb);
s.new_element<Radius>(co,rb);
s.new_element<Horizontal>(l0);
s.new_element<Angle>(l0,l1,360.0/Poles);
s.new_element<Coincident<LineSegment>>(origin,l1);
add_bounding_box(s, origin, V[1], ro, 360.0 / Poles);
}
}
\ No newline at end of file
#ifndef OERSTED_SYNCHRONOUSRELUCTANCEROTOR_H
#define OERSTED_SYNCHRONOUSRELUCTANCEROTOR_H
#include "ModelTemplate.h"
#include "PolarModelTemplate.h"
class SynchronousReluctanceRotor : public ModelTemplate {
class SynchronousReluctanceRotor : public PolarModelTemplate {
public:
void add_to_sketch(Sketch &s, std::shared_ptr<Vertex> origin) override;
......
......@@ -8,7 +8,7 @@
#include "Coincident.h"
#include "Angle.h"
// TODO: Move to a library for higher order functions
// TODO: Move to a library for higher order functions (Move to model templates?, using add_to_sketch)
// TODO: Need a better meshing strategy for small gaps
template <size_t Nl>
......@@ -16,13 +16,45 @@ class CylindricalAirgap {
public:
CylindricalAirgap() {};
CylindricalAirgap(Sketch &sketch, std::shared_ptr<const Vertex> vc, std::shared_ptr<const Vertex> vi0, std::shared_ptr<const Vertex> vo0, double_t angle) {
CylindricalAirgap(Sketch &sketch, std::shared_ptr<const Vertex> center, std::shared_ptr<const Vertex> vi0, std::shared_ptr<const Vertex> vo0, double_t angle) {
sketch.solve();
// Get matching periodic verticies
auto vi1 = sketch.select_periodic_vertex(vi0, vc, angle);
auto vo1 = sketch.select_periodic_vertex(vo0, vc, angle);
auto vi1 = sketch.select_periodic_vertex(vi0,center,angle);
auto vo1 = sketch.select_periodic_vertex(vo0,center,angle);
construct_airgap(sketch,center,vi0,vi1,vo0,vo1,angle);
}
CylindricalAirgap(Sketch &sketch, std::shared_ptr<const Vertex> center, double_t xi, double_t yi, double_t xo, double_t yo, double_t angle) {
auto vi0 = sketch.select_nearest_vertex(xi,yi);
auto vo0 = sketch.select_nearest_vertex(xo,yo);
auto vi1 = sketch.select_periodic_vertex(vi0,center,angle);
auto vo1 = sketch.select_periodic_vertex(vo0,center,angle);
construct_airgap(sketch,center,vi0,vi1,vo0,vo1,angle);
}
std::shared_ptr<Vertex> vertex0(size_t i) const { return Vertex0[i]; };
std::shared_ptr<Vertex> vertex1(size_t i) const { return Vertex1[i]; };
std::shared_ptr<LineSegment> line0(size_t i) const { return Line0[i]; };
std::shared_ptr<LineSegment> line1(size_t i) const { return Line1[i]; };
std::shared_ptr<CircularArc> interface(size_t i) const { return Interface[i]; };
protected:
std::array<std::shared_ptr<Vertex>,Nl> Vertex0;
std::array<std::shared_ptr<Vertex>,Nl> Vertex1;
std::array<std::shared_ptr<LineSegment>,Nl+1> Line0;
std::array<std::shared_ptr<LineSegment>,Nl+1> Line1;
std::array<std::shared_ptr<CircularArc>,Nl> Interface;
void construct_airgap(Sketch &sketch, std::shared_ptr<const Vertex> vc, std::shared_ptr<const Vertex> vi0, std::shared_ptr<const Vertex> vi1, std::shared_ptr<const Vertex> vo0, std::shared_ptr<const Vertex> vo1, double_t angle) {
// Determine radii
double_t ri = vc->hypot(vi0);
double_t ro = vc->hypot(vo0);
......@@ -75,25 +107,6 @@ public:
// TODO: Should constraints be added?
}
std::shared_ptr<Vertex> vertex0(size_t i) const { return Vertex0[i]; };
std::shared_ptr<Vertex> vertex1(size_t i) const { return Vertex1[i]; };
std::shared_ptr<LineSegment> line0(size_t i) const { return Line0[i]; };
std::shared_ptr<LineSegment> line1(size_t i) const { return Line1[i]; };
std::shared_ptr<CircularArc> interface(size_t i) const { return Interface[i]; };
protected:
std::array<std::shared_ptr<Vertex>,Nl> Vertex0;
std::array<std::shared_ptr<Vertex>,Nl> Vertex1;
std::array<std::shared_ptr<LineSegment>,Nl+1> Line0;
std::array<std::shared_ptr<LineSegment>,Nl+1> Line1;
std::array<std::shared_ptr<CircularArc>,Nl> Interface;
};
#endif //OERSTED_CYLINDRICALAIRGAP_H
......@@ -320,4 +320,19 @@ std::vector<std::shared_ptr<Contour const>> Sketch::select_contours(std::vector<
}
return output;
}
std::shared_ptr<Vertex> Sketch::select_nearest_vertex(double_t x, double_t y) const {
std::shared_ptr<Vertex> nearest;
double_t dmin = DBL_MAX;
for (auto v : Verticies) {
double_t d = hypot(x - v->x(), y - v->y());
if (d < dmin) {
nearest = v;
dmin = d;
}
}
return nearest;
}
\ No newline at end of file
......@@ -83,6 +83,9 @@ public:
std::vector<std::shared_ptr<Curve const>> const curves() { return std::vector<std::shared_ptr<Curve const>>{Curves.begin(), Curves.end()}; };
// Selection methods
std::shared_ptr<Vertex> select_nearest_vertex(double_t x, double_t y) const;
std::vector<PeriodicBoundaryPair> select_periodic_boundary_pairs(std::shared_ptr<Vertex const> v0, double_t angle) const;
std::vector<std::shared_ptr<Curve const>> select_radial_boundary(std::shared_ptr<Vertex const> v0, double_t radius) const;
......
......@@ -43,7 +43,7 @@ set(SOURCE_FILES
Model\ Templates/test_ModelTemplates.hpp
Model\ Templates/test_DistributedWindingStator.cpp
Model\ Templates/test_SynchronousReluctanceRotor.cpp "Model Templates/test_ModelTemplates.hpp")
Model\ Templates/test_SynchronousReluctanceRotor.cpp "Model Templates/test_ModelTemplates.hpp" "Model Templates/test_SynchronousReluctanceMachine.cpp")
add_executable(run_tests ${SOURCE_FILES})
......
#include "test_ModelTemplates.hpp"
class Synchronous_Reluctance_Machine: public ::testing::Test {
public:
virtual void SetUp() {
}
};
TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
SynchronousReluctanceRotor srr;
srr.Poles = 8;
srr.InnerRadius = 25.4e-3;
srr.OuterRadius = 89e-3-0.03*25.4e-3;
srr.AngularBridgeThickness = {0.5e-3,0.5e-3,0.5e-3};
srr.RadialBridgeThickness = {0.5e-3,0.5e-3,0.5e-3};
srr.RadialThickness = {5.0e-3,5.0e-3,5.0e-3,5.0e-3,5.0e-3,5.0e-3};
srr.AngularThickness = {3.2,3.2,3.2,3.2,3.2,3.2};
srr.Convexify = true;
srr.AddBoundingBox = false;
DistributedWindingStator dws;
dws.InnerRadius = 89e-3;
dws.ToothFaceThickness = 0.75e-3;
dws.SlotDepth = 17.4e-3;
dws.BackironThickness = 16.8e-3;
dws.SlotWidth = 6.05e-3;
dws.SlotOpening = 1.74e-3;
dws.TotalTeeth = 48;
dws.ModeledTeeth = 6;
dws.Convexify = false;
dws.AddBoundingBox = true;
Sketch sketch;
auto origin = sketch.new_element<Vertex>(0.0,0.0);
auto fo = sketch.new_element<Fixation>(origin);
srr.add_to_sketch(sketch,origin);
dws.add_to_sketch(sketch,origin);
CylindricalAirgap<1> cag{sketch,origin,srr.OuterRadius,0.0,dws.InnerRadius,0.0,360.0/8.0};
EXPECT_LE(sketch.solve(), 100e-3 * FLT_EPSILON);
EXPECT_TRUE(sketch.build());
sketch.save_as<SaveMethod::Rasterize>(SAVE_DIR,"Synchronous_Reluctance_Machine");
Mesh mesh{sketch};
mesh.create();
EXPECT_TRUE(mesh.edges_are_valid());
EXPECT_TRUE(mesh.edges_are_optimal());
mesh.MinimumElementQuality = 0.5;
EXPECT_TRUE(mesh.refine());
mesh.save_as(SAVE_DIR,"Synchronous_Reluctance_Machine");
}
\ 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