From 894840dde599798aed3565770f3280e158601317 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 25 Nov 2019 23:11:58 +0100
Subject: [PATCH] simpler geometry + comments
---
SlidingSurface3D/rfpm.geo | 517 ++++++++++---------------------
SlidingSurface3D/rfpm.pro | 2 +-
SlidingSurface3D/rfpm_common.pro | 17 +-
3 files changed, 169 insertions(+), 367 deletions(-)
diff --git a/SlidingSurface3D/rfpm.geo b/SlidingSurface3D/rfpm.geo
index aa7b9d1..a168ace 100644
--- a/SlidingSurface3D/rfpm.geo
+++ b/SlidingSurface3D/rfpm.geo
@@ -1,413 +1,230 @@
Include "rfpm_common.pro";
+SetFactory("OpenCASCADE");
-// Tolerance for the geometric identification.
-tol = 1.1*SlidingSurfaceShift;
-
-// Auxiliary Gmsh functions for OpenCascade solid modelling
-
-Function CheckNumbering
-// Check whether BooleanFragments has preserved numbering as expected
-For num In {0:#Volumes()-1}
- If( Volumes[num] != vv(num) )
- Error("BooleanFragment changed volume tag %g -> %g",
- Volumes[num], vv(num));
- Abort;
- EndIf
-EndFor
-/* EXPLAIN
-For num In {0:#Surfaces()-1}
- If( Surfaces(num) != vv(num+#Volumes()) )
- Error("BooleanFragment changed surface tag %g -> %g",
- Surfaces(num), vv(num+#Volumes()));
- Abort;
- EndIf
-EndFor
-*/
-Return
+Mesh.CharacteristicLengthMin = lc;
+Mesh.CharacteristicLengthMax = lc;
+tol = 1e-6;
+
+AnglePole = (ModelAngleMax - ModelAngleMin) * deg;
+AngleMagnet = 0.8 * AnglePole;
+
+// radial dimensions
+lRotor = 44 * mm;
+lMagnet = 10 * mm;
+lGap = 20 * mm;
+lAirRing = 1 * mm;
+lStator = 30 * mm;
+lInf = 30 * mm;
+
+// heights
+hRotor = 80 * mm;
+hMagnet = 40 * mm;
+hTot = 130 * mm;
+
+// initial rotor bloc
+Cylinder(1) = {0,0,0, 0,0,hRotor, lRotor, AnglePole};
+
+// magnet
+Cylinder(2) = {0,0,0, 0,0,hMagnet, lRotor + lMagnet / 2, AngleMagnet};
+Cylinder(3) = {0,0,0, 0,0,hMagnet, lRotor - lMagnet / 2, AngleMagnet};
+BooleanDifference(4) = { Volume{2}; Delete; }{ Volume{3}; Delete; };
+Rotate{ {0,0,1}, {0,0,0}, (AnglePole - AngleMagnet) / 2 }{ Volume{4}; }
+Physical Volume("Magnet", 2) = 4;
+
+// rotor
+BooleanDifference(5) = { Volume{1}; Delete; }{ Volume{4}; };
+Physical Volume("Rotor", 1) = 5;
+
+// rotor air
+Cylinder(6) = {0,0,0, 0,0,hTot, lRotor + lGap / 2, AnglePole};
+BooleanDifference(7) = { Volume{6}; Delete; }{ Volume{4, 5}; };
+Physical Volume("AirRotor", 4) = 7;
+
+// make rotor parts conformal
+VolumesRotor() = {4,5,7};
+BooleanFragments{ Volume{VolumesRotor()}; Delete; }{}
+
+// air ring
+Cylinder(8) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 , AnglePole};
+Cylinder(9) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 + lAirRing, AnglePole};
+BooleanDifference(10) = { Volume{9}; Delete; }{ Volume{8}; Delete; };
+Physical Volume("AirRing", 6) = 10;
+
+// stator
+Cylinder(11) = {0,0,0, 0,0,hRotor, lRotor + lGap, AnglePole};
+Cylinder(12) = {0,0,0, 0,0,hRotor, lRotor + lGap + lStator, AnglePole};
+BooleanDifference(13) = { Volume{12}; Delete; }{ Volume{11}; Delete; };
+Physical Volume("Stator", 3) = 13;
+
+// stator air
+Cylinder(14) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 + lAirRing, AnglePole};
+Cylinder(15) = {0,0,0, 0,0,hTot, lRotor + lGap + lStator + lInf, AnglePole};
+BooleanDifference(16) = { Volume{15}; Delete; }{ Volume{14}; Delete; };
+BooleanDifference(17) = { Volume{16}; Delete; }{ Volume{13}; };
+Physical Volume("AirStator", 5) = 17;
+
+// make stator parts and air ring conformal - but don't glue them to the rotor:
+// we need a "double" surface to act as "SlidingMaster" and "SlidingSlave",
+// which will be matched programmatically with a Link constraint in getdp after
+// rotating the rotor alone
+VolumesStator() = {10,13,17};
+BooleanFragments{ Volume{VolumesStator()}; Delete; }{}
+
+// We're done with the volumes. Due to internal relabelling during
+// BooleanFragments that cannot be controlled by the user, the consistent way to
+// identify the surface tags is to proceed geometrically. The surface tags that
+// are needed are in general pieces of the outer boundary of the solid where
+// specific boundary conditions are to be imposed. Very often, the volume
+// occupied by a FE model is a parallelepipedic block in some coordinate system
+// (Cartesian, cylindrical, ...). The boundary conditions are thus associated
+// with the faces of that block and can therefore be identified by means of thin
+// bounding boxes in the coordinate space.
Function ChangeCoordinates
- If( CoordinateSystem == 1) // cylindrical coordinates
- coord() = {Hypot(point(0),point(1)),Atan2(point(1),point(0)),point(2)};
+ // return coordinates of a point depending on CoordinateSystem
+ //
+ // input: point() and CoordinateSystem
+ // output: coord()
+ If(CoordinateSystem == 1) // cylindrical coordinates
+ coord() = { Hypot(point(0), point(1)), Atan2(point(1), point(0)), point(2) };
Else // Cartesian coordinates
coord() = point();
EndIf
Return
-/*
- Due to internal relabelling that cannot be controlled by the user,
- the consistent way to identify the surface labels (tags)
- of a solid modelled with the OpenCascade factory is to proceed
- geometrically.
- The surface tags that are needed are in general
- pieces of the outer boundary of the solid
- where specific boundary conditions are to be imposed.
- Very often, the volume occupied by a FE model
- is a parallelepipedic block in some coordinate system
- (Cartesian, cylindrical, ...)
- The boundary conditions are thus associated
- with the faces of that block and can therefore be identified
- by means of thin bounding boxes in the coordinate space.
- The function 'Volume2Cube' returns
- the parallelepipedic coordinate block of a list of volumes
- 'Volumes'.
- The function 'Cube2Face' returns the list of faces of the model
- that have all their points located inside a thin box encosing
- the 'FaceId'th face of that coordinate block.
-*/
Function Volume2Cube
- //
+ // return the parallelepipedic coordinate block of a list of volumes
+ //
+ // input: Volumes()
+ // output: cube()
points() = PointsOf{ Volume{ Volumes() }; };
- point() = Point { points[0] };
+ point() = Point { points(0) };
Call ChangeCoordinates;
- cube() = {coord(0),coord(1),coord(2),coord(0),coord(1),coord(2)};
- For p In { 1:#points()-1 }
- point() = Point { points(p) };
+ cube() = {coord(0), coord(1), coord(2), coord(0), coord(1), coord(2)};
+ For p In {1 : #points() - 1}
+ point() = Point{ points(p) };
Call ChangeCoordinates;
- cube(0) = (coord(0) < cube(0) ? coord(0) : cube(0) );
- cube(1) = (coord(1) < cube(1) ? coord(1) : cube(1) );
- cube(2) = (coord(2) < cube(2) ? coord(2) : cube(2) );
- cube(3) = (coord(0) > cube(3) ? coord(0) : cube(3) );
- cube(4) = (coord(1) > cube(4) ? coord(1) : cube(4) );
- cube(5) = (coord(2) > cube(5) ? coord(2) : cube(5) );
+ cube(0) = (coord(0) < cube(0) ? coord(0) : cube(0));
+ cube(1) = (coord(1) < cube(1) ? coord(1) : cube(1));
+ cube(2) = (coord(2) < cube(2) ? coord(2) : cube(2));
+ cube(3) = (coord(0) > cube(3) ? coord(0) : cube(3));
+ cube(4) = (coord(1) > cube(4) ? coord(1) : cube(4));
+ cube(5) = (coord(2) > cube(5) ? coord(2) : cube(5));
EndFor
Return
Function Cube2Face
+ // return the list of surfaces of the model that have all their points located
+ // inside a thin box encosing the FaceId'th face of that coordinate block
+ //
+ // input: cube() and FaceId
+ // output: face()
bbox() = cube();
- bbox(Modulo(FaceId+3,6)) = bbox(FaceId);
- //Printf("bbox=", bbox());
- face()={};
+ bbox(Modulo(FaceId + 3, 6)) = bbox(FaceId);
+ face() = {};
surfaces() = Boundary{ Volume{ Volumes() }; };
For s In { 0:#surfaces()-1 }
NbIn = 0;
points() = PointsOf { Surface { surfaces(s) } ; };
- For p In { 0:#points()-1 }
+ For p In { 0 : #points() - 1 }
point() = Point { points(p) };
Call ChangeCoordinates;
- If( (CoordinateSystem == 1) && (coord(0)==0) )
+ If((CoordinateSystem == 1) && (coord(0)==0))
coord(1) = bbox(1);
EndIf
- If( (coord(0)>bbox(0)-tol) &&
- (coord(1)>bbox(1)-tol) &&
- (coord(2)>bbox(2)-tol) &&
- (coord(0)<bbox(3)+tol) &&
- (coord(1)<bbox(4)+tol) &&
- (coord(2)<bbox(5)+tol) )
+ If((coord(0)>bbox(0)-tol) &&
+ (coord(1)>bbox(1)-tol) &&
+ (coord(2)>bbox(2)-tol) &&
+ (coord(0)<bbox(3)+tol) &&
+ (coord(1)<bbox(4)+tol) &&
+ (coord(2)<bbox(5)+tol))
NbIn++;
EndIf
EndFor
- If( NbIn == #points() )
+ If(NbIn == #points())
face() += surfaces(s);
EndIf
EndFor
Return
-
-SetFactory("OpenCASCADE");
-
-Solver.AutoMesh = 2;
-Mesh.CharacteristicLengthMin = lc;
-Mesh.CharacteristicLengthMax = lc;
-// Hide mesh in post-processing views
-Mesh.VolumeEdges = 0;
-Mesh.SurfaceEdges = 0;
-
-L_r = 44*mm;
-L_m = 40*mm;
-L_s = 44*mm;
-pp = 3;
-anglePole = Pi/pp;
-e_in = 0;
-e_r = 3.5*mm;
-b_depth = 0.5*mm;
-e_m = 5.3*mm;
-e_sl = 0;
-e_ar = 1*mm;
-e_ss = 0.1*mm; // air ring width
-e_a = 1*mm; // width stator air gap
-e_s = 5*mm; // width stator core
-sRM = 0.75;
-A_infinity = 66.75*mm;
-R_infinity = 21.75*mm;
-
-// Radii
-R_in = e_in;
-R_r = R_in + e_r;
-R_m_in = R_r - b_depth;
-R_m = R_m_in + e_m;
-R_sl = R_m + e_sl;
-R_ar = R_sl + e_ar;
-R_ss = R_ar + e_ss;
-R_s = R_ar + e_a;
-R_o = R_s + e_s;
-
-// R O T O R
-
-p1r = newp; Point(p1r) = {R_in,0,0};
-p2r = newp; Point(p2r) = {R_in,0,L_r/2};
-p3r = newp; Point(p3r) = {R_r,0,L_r/2};
-p4r = newp; Point(p4r) = {R_r,0,0};
-
-l1r = newl; Line(l1r) = {p1r,p2r};
-l2r = newl; Line(l2r) = {p2r,p3r};
-l3r = newl; Line(l3r) = {p3r,p4r};
-l4r = newl; Line(l4r) = {p4r,p1r};
-
-ll_r = newll; Line Loop(ll_r) = {l1r, l2r, l3r, l4r};
-sRotor = news; Plane Surface(sRotor) = {ll_r};
-vRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
- Surface{sRotor};
-};
-
-
-// Magnets
-
-p1m = newp; Point(p1m) = {R_m_in,0,0};
-p2m = newp; Point(p2m) = {R_m_in,0,L_m/2};
-p3m = newp; Point(p3m) = {R_m,0,L_m/2};
-p4m = newp; Point(p4m) = {R_m,0,0};
-
-l1m = newl; Line(l1m) = {p1m,p2m};
-l2m = newl; Line(l2m) = {p2m,p3m};
-l3m = newl; Line(l3m) = {p3m,p4m};
-l4m = newl; Line(l4m) = {p4m,p1m};
-
-ll_m = newll; Line Loop(ll_m) = {l1m, l2m, l3m, l4m};
-sMagnet = news; Plane Surface(sMagnet) = {ll_m};
-
-Rotate{ {0,0,1}, {0,0,0}, anglePole*(1-sRM)/2 } {
- Surface{sMagnet};
-}
-vMagnetExt() = Extrude {{0,0,1}, {0,0,0}, anglePole*sRM} {
- Surface{sMagnet};
-};
-
-// Air rotor
-
-p1ar = newp; Point(p1ar) = {R_in,0,0};
-p2ar = newp; Point(p2ar) = {R_in,0,A_infinity/2};
-p3ar = newp; Point(p3ar) = {R_ar,0,A_infinity/2};
-p4ar = newp; Point(p4ar) = {R_ar,0,0};
-
-l1ar = newl; Line(l1ar) = {p1ar,p2ar};
-l2ar = newl; Line(l2ar) = {p2ar,p3ar};
-l3ar = newl; Line(l3ar) = {p3ar,p4ar};
-l4ar = newl; Line(l4ar) = {p4ar,p1ar};
-
-ll_ar = newll; Line Loop(ll_ar) = {l1ar, l2ar, l3ar, l4ar};
-sAirRotor = news; Plane Surface(sAirRotor) = {ll_ar};
-vAirRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
- Surface{sAirRotor};
-};
-
-// Create volumes with boolean operations
-
-vMagnet = vMagnetExt[1];
-vRotor = BooleanDifference{
- Volume{vRotorExt[1]}; Delete;
-}{
- Volume{vMagnet};
-};
-vAirRotor = BooleanDifference{
- Volume{vAirRotorExt[1]}; Delete;
-}{
- Volume{vRotor}; Volume{vMagnet};
-};
-
-// Glue up rotor pieces with 'BooleanFragment'.
-// This operation is necessary, as it connects the different
-// parts of the model together,
-// but it entails an internal renumbering of surfaces, lines, etc...
-// Hence the need for a geometrical identification afterwards (see below).
-
-Volumes() = { vMagnet, vRotor, vAirRotor };
-VolumesRotor() = BooleanFragments{
- Volume { Volumes() }; Delete; }{};
-
-vv() = VolumesRotor(); Call CheckNumbering;
-
-// S T A T O R
-
-// Air ring
-// This region is specific to the sliding surface mechanism
-// to represent rotor motion.
-// It is a thin volumic region of width 'e_ss' in contact
-// with the 'SlidingMaster' surface.
-
-p1ss = newp; Point(p1ss) = {R_ar,0,SlidingSurfaceShift};
-p2ss = newp; Point(p2ss) = {R_ar,0,SlidingSurfaceShift+A_infinity/2};
-p3ss = newp; Point(p3ss) = {R_ss,0,A_infinity/2};
-p4ss = newp; Point(p4ss) = {R_ss,0,0};
-
-l1ss = newl; Line(l1ss) = {p1ss,p2ss};
-l2ss = newl; Line(l2ss) = {p2ss,p3ss};
-l3ss = newl; Line(l3ss) = {p3ss,p4ss};
-l4ss = newl; Line(l4ss) = {p4ss,p1ss};
-
-ll_ss = newll; Line Loop(ll_ss) = {l1ss, l2ss, l3ss, l4ss};
-sAirRing = news; Plane Surface(sAirRing) = {ll_ss};
-vAirRingExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
- Surface{sAirRing};
-};
-
-vAirRing = vAirRingExt[1] ;
-
-// Air stator
-
-p1as = newp; Point(p1as) = {R_ss,0,0};
-p2as = newp; Point(p2as) = {R_ss,0,A_infinity/2};
-p3as = newp; Point(p3as) = {R_infinity,0,A_infinity/2};
-p4as = newp; Point(p4as) = {R_infinity,0,0};
-
-l1as = newl; Line(l1as) = {p1as,p2as};
-l2as = newl; Line(l2as) = {p2as,p3as};
-l3as = newl; Line(l3as) = {p3as,p4as};
-l4as = newl; Line(l4as) = {p4as,p1as};
-
-ll_as = newll; Line Loop(ll_as) = {l1as, l2as, l3as, l4as};
-sAirStator = news; Plane Surface(sAirStator) = {ll_as};
-vAirStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
- Surface{sAirStator};
-};
-
-// Stator core
-
-p1s = newp; Point(p1s) = {R_s,0,0};
-p2s = newp; Point(p2s) = {R_s,0,L_s/2};
-p3s = newp; Point(p3s) = {R_o,0,L_s/2};
-p4s = newp; Point(p4s) = {R_o,0,0};
-
-l1s = newl; Line(l1s) = {p1s,p2s};
-l2s = newl; Line(l2s) = {p2s,p3s};
-l3s = newl; Line(l3s) = {p3s,p4s};
-l4s = newl; Line(l4s) = {p4s,p1s};
-
-ll_s = newll; Line Loop(ll_s) = {l1s, l2s, l3s, l4s};
-sStator = news; Plane Surface(sStator) = {ll_s};
-vStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
- Surface{sStator};
-};
-
-// Create volumes by boolean operations
-
-vStator = vStatorExt[1];
-vAirStator = BooleanDifference{
- Volume{vAirStatorExt[1]}; Delete;
-}{
- Volume{vStator};
-};
-
-// Glue up stator pieces
-Volumes() = { vAirRing, vStator, vAirStator };
-VolumesStator() = BooleanFragments{
- Volume { Volumes() }; Delete; }{};
-
-Printf("VolumesStator=", VolumesStator());
-Printf("VolumesRotor=", VolumesRotor());
-
-// Recover Physical Surfaces geometrically
-
CoordinateSystem = 1; // cylindrical coordinates
+// get useful rotor surfaces
Volumes() = VolumesRotor(); Call Volume2Cube;
-Printf("cube rotor", cube());
-FaceId = 1; Call Cube2Face; // no statement here FIXME
+Printf("cube rotor: ", cube());
+FaceId = 1; Call Cube2Face; // FIXME: cannot have statement here
RotorPerMaster() = face();
-FaceId = 2; Call Cube2Face;
+FaceId = 2; Call Cube2Face;
RotorBottom() = face();
-FaceId = 3; Call Cube2Face;
+FaceId = 3; Call Cube2Face;
SlidingSlave() = face();
-FaceId = 4; Call Cube2Face;
+FaceId = 4; Call Cube2Face;
RotorPerSlave() = face();
-FaceId = 5; Call Cube2Face;
+FaceId = 5; Call Cube2Face;
RotorTop() = face();
-/* control if the geometrical identification succeeded
-Physical Surface("Test1", 101) = { RotorPerMaster() };
-Physical Surface("Test2", 102) = { RotorTop() };
-Physical Surface("Test3", 103) = { SlidingSlave() };
-Physical Surface("Test4", 104) = { RotorBottom() };
-Physical Surface("Test5", 105) = { RotorPerSlave() };
-*/
-
+// get useful stator surfaces
Volumes() = VolumesStator(); Call Volume2Cube;
-Printf("cube stator", cube());
-FaceId = 0; Call Cube2Face;
+Printf("cube stator: ", cube());
+FaceId = 0; Call Cube2Face;
SlidingMaster() = face();
-FaceId = 1; Call Cube2Face;
+FaceId = 1; Call Cube2Face;
StatorPerMaster() = face();
-FaceId = 2; Call Cube2Face;
+FaceId = 2; Call Cube2Face;
StatorBottom() = face();
-FaceId = 3; Call Cube2Face;
+FaceId = 3; Call Cube2Face;
StatorOuter() = face();
-FaceId = 4; Call Cube2Face;
+FaceId = 4; Call Cube2Face;
StatorPerSlave() = face();
-FaceId = 5; Call Cube2Face;
+FaceId = 5; Call Cube2Face;
StatorTop() = face();
-/* control if the geometrical identification succeeded
-Physical Surface("Test11", 111) = { SlidingMaster() };
-Physical Surface("Test12", 112) = { StatorPerMaster() };
-Physical Surface("Test13", 113) = { StatorBottom() };
-Physical Surface("Test14", 114) = { StatorOuter() };
-Physical Surface("Test15", 115) = { StatorPerSlave() };
-Physical Surface("Test16", 116) = { StatorTop() };
-*/
+// Make AirRing a transfinite region. The width of the air ring is arbitrary,
+// but it is in general thin, especially in machine models where it has to fit
+// into the air gap. To improve the quality of the mesh, it is thus appropriate
+// to have it transfinite in practice.
-
-// Make AirRing a transfinite region
-// The width of the air ring is arbitrary, but it is in general thin,
-// especially in machine models where it has to fit into the air gap.
-// To improve the quality of the mesh,
-// it is thus appropriate to have it transfinite in practice.
-
-surfaces() = Boundary { Volume{ vAirRing }; };
-//Printf("%g %g", vAirRing, #surfaces());
+surfaces() = Boundary { Volume{10}; };
lines() = Boundary { Surface{ SlidingSlave() }; };
For num In { 0:#surfaces()-1 }
- lines() += Boundary { Surface{ surfaces[num] }; };
- Printf("%g %g", surfaces[num], #lines());
+ lines() += Boundary { Surface{ surfaces(num) }; };
EndFor
For num In { 0:#lines()-1 }
points() = PointsOf { Line { lines(num) }; };
- If( #points() != 2 )
+ If(#points() != 2)
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
- // Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
- If( Fabs( Atan2( ptA(1), ptA(0) ] - Atan2[ ptB(1), ptB(0) ) ) > tol )
+ If(Fabs(Atan2(ptA(1), ptA(0)) - Atan2(ptB(1), ptB(0))) > tol)
// line along theta-axis
- Transfinite Line{ lines[num] } = NbrDiv+1;
- ElseIf( Fabs[ ptA(2) - ptB(2) ] > tol )
+ Transfinite Line{ lines(num) } = NbrDiv+1;
+ ElseIf(Fabs(ptA(2) - ptB(2)) > tol)
// line along z-axis
- Transfinite Line{ lines(num) } = Fabs( ptA(2) - ptB(2) )/lc-1;
- Else
+ Transfinite Line{ lines(num) } = Fabs(ptA(2) - ptB(2))/lc-1;
+ Else
// line along r-axis
Transfinite Line{ lines(num) } = 2;
EndIf
EndFor
-Transfinite Surface{ surfaces() } ;
-Transfinite Surface{ SlidingSlave() } Right;
-Transfinite Volume{ vAirRing };
+Transfinite Surface{ surfaces() } ;
+Transfinite Surface{ SlidingSlave() };
+Transfinite Volume{ 10 };
// Identify 'SlidingSubmaster' Curve
-AngleSubSurfaces = (ModelAngleMax-ModelAngleMin)*deg;
-
lines() = Boundary { Surface{ SlidingMaster() }; };
For l In { 0:#lines()-1 }
points() = PointsOf { Line { lines(l) }; };
- If( #points() != 2 )
+ If(#points() != 2)
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
- // Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
- If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces) < tol &&
- Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces) < tol )
+ If(Fabs(Atan2(ptA(1), ptA(0)) - AnglePole) < tol &&
+ Fabs(Atan2(ptB(1), ptB(0)) - AnglePole) < tol)
SlidingSubmaster = lines(l);
EndIf
EndFor
@@ -417,30 +234,20 @@ EndFor
lines() = Boundary { Surface{ SlidingSlave() }; };
For l In { 0:#lines()-1 }
points() = PointsOf { Line { lines(l) }; };
- If( #points() != 2 )
+ If(#points() != 2)
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
- If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces ) < tol &&
- Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces ) < tol )
+ If(Fabs(Atan2(ptA(1), ptA(0)) - AnglePole) < tol &&
+ Fabs(Atan2(ptB(1), ptB(0)) - AnglePole) < tol)
SlidingSubslave = lines(l);
EndIf
EndFor
-
-// P H Y S I C A L R E G I O N S
-
-Physical Volume("Rotor" , 1) = {vRotor};
-Physical Volume("Magnet" , 2) = {vMagnet};
-Physical Volume("Stator" , 3) = {vStator};
-Physical Volume("AirRotor" , 4) = {vAirRotor};
-Physical Volume("AirStator", 5) = {vAirStator};
-Physical Volume("AirRing" , 6) = {vAirRing};
-
Physical Surface("Outer" , 10) = {StatorOuter()};
-Physical Surface("Top" , 11) = {RotorTop(),StatorTop()};
-Physical Surface("Bottom" , 12) = {RotorBottom(),StatorBottom()};
+Physical Surface("Top" , 11) = {RotorTop(), StatorTop()};
+Physical Surface("Bottom" , 12) = {RotorBottom(), StatorBottom()};
Physical Surface("RotorPerMaster" , 13) = {RotorPerMaster()};
Physical Surface("RotorPerSlave" , 14) = {RotorPerSlave()};
Physical Surface("StatorPerMaster", 15) = {StatorPerMaster()};
@@ -448,24 +255,19 @@ Physical Surface("StatorPerSlave" , 16) = {StatorPerSlave()};
Physical Surface("SlidingMaster" , 17) = {SlidingMaster()};
Physical Surface("SlidingSlave" , 18) = {SlidingSlave()};
-Physical Line("Sliding_Submaster", 20) = {SlidingSubmaster};
-Physical Line("Sliding_Subslave", 21) = {SlidingSubslave};
-
+Physical Line("Sliding_Submaster" , 20) = {SlidingSubmaster};
+Physical Line("Sliding_Subslave" , 21) = {SlidingSubslave};
// Ensure identical meshes on the periodic faces
-// The 'Air ring' being meshed transfinite with prisms,
-// the periodic surface routine issues a harmless warning
-// Info : No periodic vertices in surface xxx
For num In {0:#RotorPerMaster()-1}
- Periodic Surface { RotorPerSlave(num) }
- = { RotorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
+ Periodic Surface { RotorPerSlave(num) } =
+ { RotorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, AnglePole };
EndFor
For num In {0:#StatorPerMaster()-1}
- Periodic Surface { StatorPerSlave(num) }
- = { StatorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
+ Periodic Surface { StatorPerSlave(num) } =
+ { StatorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, AnglePole };
EndFor
-
treeLines() = CombinedBoundary { Physical Surface { 10 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 11 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 12 }; } ;
@@ -477,4 +279,3 @@ treeLines() += CombinedBoundary { Physical Surface { 17 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 18 }; } ;
Physical Line("treeLines", 22) = { treeLines() };
-
diff --git a/SlidingSurface3D/rfpm.pro b/SlidingSurface3D/rfpm.pro
index c7fc558..cbc46fb 100644
--- a/SlidingSurface3D/rfpm.pro
+++ b/SlidingSurface3D/rfpm.pro
@@ -18,7 +18,7 @@ DefineConstant[
Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
RotorPosition = deg*
- DefineNumber[0, Name "Options/Rotor position [deg]",
+ DefineNumber[20, Name "Options/Rotor position [deg]",
Visible !Flag_QuasiStatic];
AngularStep = deg*
DefineNumber[1, Name "Options/1Angular step [deg]",
diff --git a/SlidingSurface3D/rfpm_common.pro b/SlidingSurface3D/rfpm_common.pro
index c4054fd..76b3f4f 100644
--- a/SlidingSurface3D/rfpm_common.pro
+++ b/SlidingSurface3D/rfpm_common.pro
@@ -1,18 +1,19 @@
mm = 1.e-3 ;
deg = Pi/180.;
-lc = mm * DefineNumber[2, Name "Options/Global mesh size (mm)"];
-If(lc == 0)
- Error("Glabal mesh size cannot be zero. Reset to default."); lc=2*mm;
+lc = mm * DefineNumber[10, Name "Options/Global mesh size (mm)"];
+If(lc <= 0)
+ Error("Glabal mesh size cannot be zero. Reset to default.");
+ lc = 10 * mm;
EndIf
// Number of elements on the sliding surface in rotation direction
-NbrDiv = 30*mm/lc;
+NbrDiv = 50*mm/lc;
-Gauge = DefineNumber[1, Name "Options/Gauging",
- Choices{ 0="MUMPS", 1="Tree"}];
+Gauge = DefineNumber[1, Name "Options/Gauging", Choices{ 0="MUMPS", 1="Tree"}];
-ModelAngleMin = 0. ;
+ModelAngleMin = 0. ;
ModelAngleMax = 60. ;
-SlidingSurfaceShift = 1.e-5;
+//SlidingSurfaceShift = 1.e-5;
+SlidingSurfaceShift = 0;
--
GitLab