From 13227a1b996e51e3f381e18b48f07d9989391643 Mon Sep 17 00:00:00 2001
From: Francesco Carpanese <francesco.carpanese@epfl.ch>
Date: Sun, 8 Dec 2019 15:51:17 +0100
Subject: [PATCH] Moving to a unique version of LIUQE IDS to include all cases
---
IDS/gen_filament.m | 162 ++++++++++++++-------------------------------
1 file changed, 51 insertions(+), 111 deletions(-)
diff --git a/IDS/gen_filament.m b/IDS/gen_filament.m
index bee54bfb..8a2e09aa 100644
--- a/IDS/gen_filament.m
+++ b/IDS/gen_filament.m
@@ -1,181 +1,121 @@
-function [discretizationr, discretizationz, Twc] = eval_trap(vvdata, num, type, varargin )
-% Evaluate the trapezoid from vv data structure, and computes the
-% discretization
-% type = 2 -> rectangle description
-% type = 3 -> oblique description
-% type = 4 -> contour of the coils are specified varargin = [r, z ]
-% contours
-% vv data contains R, Z, dR, dZ for type = 2.
-% vv data contains R, Z, dR, dZ, AC, AC2 for type = 3 where AC, AC2 are the
-% angles of the parallelogram. Refer to GA document for specifications.
+function [discretizationr, discretizationz, Twc] = gen_filament(vvdata, num, type_coil_description,varargin)
+% Generate filamentary(windings) discretization from different geometrical coil descriptions
-if mod(num,1)~=0, warning('The number of turns must be an integer. Ceil operation on it'); end
-num = ceil(num);
+% type_coil_description = 2 -> rectangle description
+% type_coil_description = 3 -> oblique description
+% type_coil_description = 4 -> contour of the coils are specified varargin = [r, z ] contours
+% vvdata = R, Z, dR, dZ for type_coil_description = 2.
+% vvdata = R, Z, dR, dZ, AC, AC2 for type_coil_description = 3 where AC, AC2 are the angles of the parallelogram.
R = vvdata(1, :);
Z = vvdata(2, :);
dR = vvdata(3, :);
dZ = vvdata(4, :);
-
-switch type
+switch type_coil_description
case 1 % outline option
- disp('Not yet implemented. It will contain outline option IMAS')
- return
+ error('Filamentary discretization for ''outline'' coils description in ids type_coil_description not yet implmented')
case 2 % rectangle
AC = 0*ones(size(R));
AC2 = 90*ones(size(R));
-
- case 3 % oblique
+ case 3 % oblique (parallelogram)
AC = vvdata(5, :);
AC2 = vvdata(6, :);
-
case 4 % contour of the polygon
-
+ error('Filamentary discretizaion for ''polygon'' coils description in ids typ not yet implemented')
otherwise
- disp('case not defined')
- return
+ error('Coils description not supported for filamenteray discretization ')
end
% Initialize the extrema of the trapezoid
pointr = zeros(numel(AC),4);
pointz = zeros(numel(AC),4);
-
discretizationr = [];
discretizationz = [];
-
Twc = []; % Grouping matrix from windings to active coil
-for ii = 1:numel(AC) % loop all the coils
+for ii = 1:numel(AC) % loop all the coils given vvstructure
+ % Find 4 extrema of the parallelogram
+ % D----C
+ % / /
+ % A----B
if AC2(ii) ~= 0 && AC(ii) ==0
DH = dZ(ii)/tand(AC2(ii));
pointr(ii,1) = R(ii) -DH/2 -dR(ii)/2;
pointr(ii,2) = R(ii) +DH/2 -dR(ii)/2;
pointr(ii,3) = R(ii) +DH/2 +dR(ii)/2;
pointr(ii,4) = R(ii) -DH/2 +dR(ii)/2;
-
pointz(ii,1) = Z(ii) -dZ(ii)/2;
pointz(ii,2) = Z(ii) +dZ(ii)/2;
pointz(ii,3) = Z(ii) +dZ(ii)/2;
pointz(ii,4) = Z(ii) -dZ(ii)/2;
-
elseif AC(ii) ~= 0 && AC2(ii) ==0
DH = dR(ii)*tand(AC(ii));
-
pointr(ii,1) = R(ii) - dR(ii)/2;
pointr(ii,2) = R(ii) - dR(ii)/2;
pointr(ii,3) = R(ii) + dR(ii)/2;
pointr(ii,4) = R(ii) + dR(ii)/2;
-
pointz(ii,1) = Z(ii) - DH/2- dZ(ii)/2;
pointz(ii,2) = Z(ii) - DH/2 +dZ(ii)/2;
pointz(ii,3) = Z(ii) + DH/2 +dZ(ii)/2;
pointz(ii,4) = Z(ii) + DH/2 -dZ(ii)/2;
-
elseif AC2(ii) == 0 && AC(ii) ==0
pointr(ii,1) = R(ii) -dR(ii)/2;
pointr(ii,2) = R(ii) -dR(ii)/2;
pointr(ii,3) = R(ii) +dR(ii)/2;
pointr(ii,4) = R(ii) +dR(ii)/2;
-
pointz(ii,1) = Z(ii) -dZ(ii)/2;
pointz(ii,2) = Z(ii) +dZ(ii)/2;
pointz(ii,3) = Z(ii) +dZ(ii)/2;
pointz(ii,4) = Z(ii) -dZ(ii)/2;
end
-
-
+ % Define base vector of the 2D deformed space for the parallelogram
AD = [pointr(ii,4)-pointr(ii,1), pointz(ii,4)-pointz(ii,1)];
AB = [pointr(ii,2)-pointr(ii,1), pointz(ii,2)-pointz(ii,1)];
ex = [1,0];
ey = [0,1];
-
+ % Rotation matrix containing the scalar product of the deformed space
+ % and orthogonal space
S = ([ex*AD', ex*AB'; ey*AD', ey*AB']);
- %fixed_area = 0.01; %Decide the area represented by each single coils
- %area_of_element = abs(AB(1)*AD(2)-AB(2)*AD(1)); %Evaluate the closer area to the element
+ % Try to use a number of filaments that respect the shape of the parallelogram
+ ratio = sqrt((AD*AD'))/sqrt((AB*AB'));
+ num1 = ceil(sqrt(num(ii)*ratio));
+ num2 = ceil(num1/ratio);
- %num = round(area_of_element/fixed_area); %Find the number of integer
- %filaments closer to a fixed desired area for each element.
-
-
- if numel(num) == 1 %selfsimilar equally spaced
-
- pointx = zeros(num*num, 1);
- pointy = zeros(num*num, 1);
-
- % Mode with selfsimilar shape
- for kk=1:num
- for jj=1:num
- pointx((kk-1)*num+jj) = (2*kk-1)/(2*num); %uniform grid in the deformed space defined by the parallelogram
- pointy((kk-1)*num+jj) = (2*jj-1)/(2*num);
-
- tmp = S*[pointx((kk-1)*(num)+jj), pointy((kk-1)*(num)+jj)]'; % evaluate component in the r,z system of reference
- discretizationr(end+1) = pointr(ii,1)+ tmp(1) ; % translate to the original point
- discretizationz(end+1) = pointz(ii,1) + tmp(2);
- end
+ pointx = zeros(num1*num2, 1);
+ pointy = zeros(num1*num2, 1);
+ for kk=1:num1
+ for jj=1:num2
+ pointx((kk-1)*num1+jj) = (2*kk-1)/(2*num1); %uniform grid in the deformed space defined by the parallelogram
+ pointy((kk-1)*num2+jj) = (2*jj-1)/(2*num2);
+ tmp = S*[pointx((kk-1)*(num1)+jj), pointy((kk-1)*(num2)+jj)]'; % evaluate component in the r,z system of reference
+ discretizationr(end+1) = pointr(ii,1) + tmp(1) ; % translate to the original central location
+ discretizationz(end+1) = pointz(ii,1) + tmp(2); % translate to the original central location
end
-
- % Creating the grouping matrix
- tmp = zeros(num*num, numel(AC));
- tmp(:,ii) = 1;
-
- Twc = [Twc ; tmp];
-
- else %Depending on the turns and same similar aspect ration
- ratio = sqrt(dot(AB,AB)/dot(AD,AD));
-
- nx = ceil(sqrt(num(ii)/ratio));
- ny = ceil(nx*ratio);
-
- if ny>=2
- pointy = linspace(1/(ny-1)/2, 1- 1/(ny-1)/2, ny);
- else
- pointy = 0.5;
- end
-
- if nx>=2
- pointx = linspace(1/(nx-1)/2, 1- 1/(nx-1)/2, nx);
- else
- pointx = 0.5;
- end
-
-
- [ pointxx, pointyy ] = meshgrid(pointx, pointy);
-
- for kk=1:ny
- for jj=1:nx
-
- tmp = S*[pointxx(kk,jj); pointyy(kk,jj)]; % evaluate component in the r,z system of reference
- discretizationr(end+1) = pointr(ii,1)+ tmp(1) ; % translate to the original point
- discretizationz(end+1) = pointz(ii,1) + tmp(2);
- end
- end
-
- % Creating the grouping matrix
- tmp = zeros(nx*ny, numel(num));
- tmp(:,ii) = 1;
- Twc = [Twc ; tmp];
-
end
-
-
-
+ % Only the number of turns is specified but not the position of
+ % each individual winding. The function eval_trap generates a winding
+ % distrubution that nicely resamble the geometry of the coil section,
+ % but will not necessarily generate a number of windings = number of
+ % turns. The following normalization garanties that the total current
+ % flowing in the cross section of coil Iw*nw = Ia*nturns respect the input data.
+ Twc(end+1:end+num1*num2,end+1) = 1*sign(num(ii))/(num1*num2)*num(ii);
end
-
-%% Plot the coils discretization from the original data
-%{
-figure
-hold on
-axis equal
-for ii = 1:numel(AC) % loop all the coils
- patch(pointr(ii,:),pointz(ii,:),'red');
-end
-scatter(discretizationr, discretizationz, '.');
-%}
+% Reshape filaments(windings) into unique vector
discretizationr = reshape(discretizationr, numel(discretizationr),1);
discretizationz = reshape(discretizationz, numel(discretizationz),1);
+% Plot the coils discretization
+if ~isempty(varargin) && varargin{1} == 1
+ figure
+ hold on
+ axis equal
+ for ii = 1:numel(AC) % loop all the coils
+ patch(pointr(ii,:),pointz(ii,:),'red');
+ end
+ scatter(discretizationr, discretizationz, '.');
+end
end
--
GitLab