Skip to content

Commit a3d7032

Browse files
riccardotorchioriccardotorchio
committed
added analytical test case (100x100x1 m^3 square)
1 parent d469b9f commit a3d7032

16 files changed

+1084
-6
lines changed

‎FFT_VI.m

+10-3
Original file line numberDiff line numberDiff line change
@@ -9,19 +9,20 @@
99
%% BEGIN USER SETTINGS
1010
%%
1111
%% Directory
12-
name_dir='test1';
12+
name_dir='test2';
1313
%% Frequency
1414
freq = 50; %[Hz]
1515
%% Selections
1616
plot_vectorsJ_flag = 1; %quiver plot of real and imag of J
1717
plot_potential_flag = 1; %color plot of phi real and imag
1818
paraview_export_flag = 1; % export to paraviw
1919
refine.flag = 0; refine.x=1; refine.y=1; refine.z=1; % refine
20-
Integration_flag = 'NumNum'; %'NumAn'; 'NumNum' (Integration: NumericalNumerical or AnalyticalNumerical)
20+
Integration_flag = 'NumAn'; %'NumAn'; 'NumNum' (Integration: NumericalNumerical or AnalyticalNumerical)
2121
ext_field_flag = 1; % exernal field
2222
% below you can write the external electric field as a function of x,y,z
2323
% and omega. Active only if ext_field_flag=1
24-
Ex_ext = @(x,y,z,omega) -1j*0.5*omega*y; Ey_ext = @(x,y,z,omega) 1j*0.5*omega*x; Ez_ext = @(x,y,z,omega) 0*z; % external field
24+
scal=0.991/sqrt(2);
25+
Ex_ext = @(x,y,z,omega) -scal*1j*0.5*omega*y; Ey_ext = @(x,y,z,omega) scal*1j*0.5*omega*x; Ez_ext = @(x,y,z,omega) 0*z; % external field, this electric field gives a vertical magnetic field of magnitude "scal"
2526
%% Solver parameters
2627
tol = 1e-6;
2728
inner_it = 40;
@@ -276,4 +277,10 @@
276277
fun_for_ParaView_vec_HEXA(...
277278
jjR(idxV,:),jjI(idxV,:),P0,VP,dad,[modelname,'J']);
278279
warning on
280+
end
281+
%%
282+
Plosses=real(vsol(1:num_curr)'*([z_realx_loc*dx/(dy*dz);z_realy_loc*dy/(dx*dz);z_realz_loc*dz/(dy*dx)].*vsol(1:num_curr))) % losses
283+
if strcmp(name_dir,'test2') && freq==50
284+
Plosses_reference_value=226.97
285+
disp('Analytical results given in: https://doi.org/10.1103/PhysRevSTAB.14.062401 ')
279286
end

‎data/test2/data.mat

38.2 KB
Binary file not shown.

‎data/test2/mainVOXELISE.m

+205
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,205 @@
1+
clear
2+
clear global
3+
close all
4+
clc
5+
global meshXmin meshXmax meshYmin meshYmax meshZmin meshZmax
6+
%% BEGIN USER SETTINGS
7+
paraview_export_flag = 1;
8+
x_ray_flag = 1;
9+
model_name='plate';
10+
%
11+
stl_files(1).name = 'plate.stl';
12+
stl_files(1).tag = 'cond';
13+
stl_files(1).rho=1/1.33e6;
14+
% to scale a stl file from any unit to meters
15+
scal_geomery.x=1/1; scal_geomery.y=1/1; scal_geomery.z=1/1;
16+
% Box
17+
% number of voxels in the x y z directions
18+
Nx=100;
19+
Ny=100;
20+
Nz=1;
21+
% corners
22+
flag_auto=1; % if 1, user_data below are ignored
23+
% user_data
24+
meshXmin = 0; % (m)
25+
meshXmax = 100e-3; % (m)
26+
meshYmin = 0;% (m)
27+
meshYmax = 100e-3; % (m)
28+
meshZmin = 0;% (m)
29+
meshZmax = 1e-3; % (m)
30+
%% END USER SETTINGS
31+
how_many_stl=size(stl_files,2);
32+
%%
33+
dad=pwd;
34+
cd ..; cd ..; cd('fun'); addpath(genpath(pwd)); cd(dad)
35+
%% Plot the original STL mesh
36+
ccolor=distinguishable_colors(how_many_stl);
37+
figure
38+
hold on
39+
xmin=[];
40+
xmax=[];
41+
ymin=[];
42+
ymax=[];
43+
zmin=[];
44+
zmax=[];
45+
for ii = 1:how_many_stl
46+
[stlcoords] = READ_stl(stl_files(ii).name);
47+
xco = squeeze( stlcoords(:,1,:) )';
48+
yco = squeeze( stlcoords(:,2,:) )';
49+
zco = squeeze( stlcoords(:,3,:) )';
50+
[hpat] = patch(xco*scal_geomery.x,yco*scal_geomery.y,zco*scal_geomery.z,ccolor(ii,:));
51+
alpha(0.3) % for trasparency
52+
axis equal
53+
xlabel('x');
54+
ylabel('y');
55+
zlabel('z');
56+
view(3)
57+
title('stl')
58+
drawnow
59+
%
60+
xmin=min([xmin;xco(:)]);
61+
xmax=max([xmax;xco(:)]);
62+
ymin=min([ymin;yco(:)]);
63+
ymax=max([ymax;yco(:)]);
64+
zmin=min([zmin;zco(:)]);
65+
zmax=max([zmax;zco(:)]);
66+
end
67+
%%
68+
if flag_auto
69+
meshXmin=xmin;
70+
meshXmax=xmax;
71+
meshYmin=ymin;
72+
meshYmax=ymax;
73+
meshZmin=zmin;
74+
meshZmax=zmax;
75+
else
76+
meshXmin = meshXmin/scal_geomery.x;
77+
meshXmax = meshXmax/scal_geomery.x;
78+
meshYmin = meshYmin/scal_geomery.y;
79+
meshYmax = meshYmax/scal_geomery.y;
80+
meshZmin = meshZmin/scal_geomery.z;
81+
meshZmax = meshZmax/scal_geomery.z;
82+
end
83+
%%
84+
disp('===================================================================')
85+
disp('voxelize...')
86+
for ii = 1:how_many_stl
87+
[o(ii).OUTPUTgrid,...
88+
o(ii).gridCOx,...
89+
o(ii).gridCOy,...
90+
o(ii).gridCOz] = VOXELISE_mod2(Nx,Ny,Nz,stl_files(ii).name,'xyz');
91+
o(ii).gridCOx=o(ii).gridCOx*scal_geomery.x;
92+
o(ii).gridCOy=o(ii).gridCOy*scal_geomery.y;
93+
o(ii).gridCOz=o(ii).gridCOz*scal_geomery.z;
94+
end
95+
xyz= grid3dRT2(o(1).gridCOx,o(1).gridCOy,o(1).gridCOz);
96+
xd=xyz(:,:,:,1);
97+
yd=xyz(:,:,:,2);
98+
zd=xyz(:,:,:,3);
99+
%
100+
idx=zeros(Nx,Ny,Nz);
101+
for ii = 1:how_many_stl
102+
idx=idx+(o(ii).OUTPUTgrid);
103+
end
104+
idx=find(idx);
105+
% plot
106+
xidx=xd(idx);
107+
yidx=yd(idx);
108+
zidx=zd(idx);
109+
disp(' ')
110+
%%
111+
disp('===================================================================')
112+
disp('x_ray...')
113+
if x_ray_flag
114+
for ii = 1:how_many_stl
115+
figure
116+
subplot(1,3,1);
117+
title(['xray object',num2str(ii),' ZY'])
118+
hold on
119+
imagesc(squeeze(sum(o(ii).OUTPUTgrid,1)));
120+
colormap(gray(256));
121+
xlabel('Z-direction');
122+
ylabel('Y-direction');
123+
axis equal
124+
subplot(1,3,2);
125+
title(['xray object',num2str(ii),' ZX'])
126+
hold on
127+
imagesc(squeeze(sum(o(ii).OUTPUTgrid,2)));
128+
colormap(gray(256));
129+
xlabel('Z-direction');
130+
ylabel('X-direction');
131+
axis equal
132+
subplot(1,3,3);
133+
title(['xray object',num2str(ii),' YX'])
134+
hold on
135+
imagesc(squeeze(sum(o(ii).OUTPUTgrid,3)));
136+
colormap(gray(256));
137+
xlabel('Y-direction');
138+
ylabel('X-direction');
139+
axis equal
140+
end
141+
end
142+
drawnow
143+
%%
144+
if Nx>1
145+
dx=abs(o(1).gridCOx(2)-o(1).gridCOx(1));
146+
else
147+
dx=(max(xco(:))-min(xco(:)))*scal_geomery.x;
148+
end
149+
if Ny>1
150+
dy=abs(o(1).gridCOy(2)-o(1).gridCOy(1));
151+
else
152+
dy=(max(yco(:))-min(yco(:)))*scal_geomery.y;
153+
end
154+
if Nz>1
155+
dz=abs(o(1).gridCOz(2)-o(1).gridCOz(1));
156+
else
157+
dz=(max(zco(:))-min(zco(:)))*scal_geomery.z;
158+
end
159+
xidx=xd([idx]);
160+
yidx=yd([idx]);
161+
zidx=zd([idx]);
162+
disp(' ')
163+
%% paraview
164+
disp('===================================================================')
165+
disp('paraview...')
166+
if paraview_export_flag
167+
P0=[...
168+
[xidx-dx/2,yidx-dy/2,zidx+dz/2];...
169+
[xidx-dx/2,yidx-dy/2,zidx-dz/2];...
170+
[xidx+dx/2,yidx-dy/2,zidx-dz/2];...
171+
[xidx+dx/2,yidx-dy/2,zidx+dz/2];...
172+
[xidx-dx/2,yidx+dy/2,zidx+dz/2];...
173+
[xidx-dx/2,yidx+dy/2,zidx-dz/2];...
174+
[xidx+dx/2,yidx+dy/2,zidx-dz/2];...
175+
[xidx+dx/2,yidx+dy/2,zidx+dz/2]];
176+
nVox=length([idx]);
177+
VP=[1:nVox;...
178+
nVox+1:2*nVox;...
179+
2*nVox+1:3*nVox;...
180+
3*nVox+1:4*nVox;...
181+
4*nVox+1:5*nVox;...
182+
5*nVox+1:6*nVox;...
183+
6*nVox+1:7*nVox;...
184+
7*nVox+1:8*nVox];
185+
val=2*ones(length([idx]),1);
186+
val(1:length(idx))=1;
187+
dad=pwd;
188+
[Ricc] = fun_for_ParaView_sca_HEXA(...
189+
val,val,P0,VP,dad,model_name);
190+
end
191+
%%
192+
L=length(o(1).gridCOx);
193+
M=length(o(1).gridCOy);
194+
N=length(o(1).gridCOz);
195+
nVoxel=L*M*N;
196+
smeshx=dx;smeshy=dy;smeshz=dz;
197+
%%
198+
for ii = 1:how_many_stl
199+
Ind(ii).ind= find(o(ii).OUTPUTgrid);
200+
Ind(ii).tag=stl_files(ii).tag;
201+
Ind(ii).rho=stl_files(ii).rho;
202+
end
203+
%%
204+
Nmat = how_many_stl;
205+
save data.mat Ind L M N Nmat nVoxel smeshx smeshy smeshz xyz -v7.3

‎data/test2/plate.stl

+86
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
solid COMSOL rendering object blk1
2+
facet normal 0 0 -1
3+
outer loop
4+
vertex 0 0.100000001490116 0
5+
vertex 0.100000001490116 0.100000001490116 0
6+
vertex 0.100000001490116 0 0
7+
endloop
8+
endfacet
9+
facet normal 0 0 -1
10+
outer loop
11+
vertex 0.100000001490116 0 0
12+
vertex 0 0 0
13+
vertex 0 0.100000001490116 0
14+
endloop
15+
endfacet
16+
facet normal -1 0 0
17+
outer loop
18+
vertex 0 0 0.00100000004749745
19+
vertex 0 0.100000001490116 0.00100000004749745
20+
vertex 0 0.100000001490116 0
21+
endloop
22+
endfacet
23+
facet normal -1 0 0
24+
outer loop
25+
vertex 0 0.100000001490116 0
26+
vertex 0 0 0
27+
vertex 0 0 0.00100000004749745
28+
endloop
29+
endfacet
30+
facet normal 0 -1 0
31+
outer loop
32+
vertex 0.100000001490116 0 0
33+
vertex 0.100000001490116 0 0.00100000004749745
34+
vertex 0 0 0.00100000004749745
35+
endloop
36+
endfacet
37+
facet normal 0 -1 0
38+
outer loop
39+
vertex 0 0 0.00100000004749745
40+
vertex 0 0 0
41+
vertex 0.100000001490116 0 0
42+
endloop
43+
endfacet
44+
facet normal 0 0 1
45+
outer loop
46+
vertex 0.100000001490116 0 0.00100000004749745
47+
vertex 0.100000001490116 0.100000001490116 0.00100000004749745
48+
vertex 0 0.100000001490116 0.00100000004749745
49+
endloop
50+
endfacet
51+
facet normal 0 0 1
52+
outer loop
53+
vertex 0 0.100000001490116 0.00100000004749745
54+
vertex 0 0 0.00100000004749745
55+
vertex 0.100000001490116 0 0.00100000004749745
56+
endloop
57+
endfacet
58+
facet normal 1 0 0
59+
outer loop
60+
vertex 0.100000001490116 0.100000001490116 0
61+
vertex 0.100000001490116 0.100000001490116 0.00100000004749745
62+
vertex 0.100000001490116 0 0.00100000004749745
63+
endloop
64+
endfacet
65+
facet normal 1 0 0
66+
outer loop
67+
vertex 0.100000001490116 0 0.00100000004749745
68+
vertex 0.100000001490116 0 0
69+
vertex 0.100000001490116 0.100000001490116 0
70+
endloop
71+
endfacet
72+
facet normal 0 1 0
73+
outer loop
74+
vertex 0 0.100000001490116 0.00100000004749745
75+
vertex 0.100000001490116 0.100000001490116 0.00100000004749745
76+
vertex 0.100000001490116 0.100000001490116 0
77+
endloop
78+
endfacet
79+
facet normal 0 1 0
80+
outer loop
81+
vertex 0.100000001490116 0.100000001490116 0
82+
vertex 0 0.100000001490116 0
83+
vertex 0 0.100000001490116 0.00100000004749745
84+
endloop
85+
endfacet
86+
endsolid COMSOL rendering object blk1

‎data/test2/res_para/plate.hdf

163 KB
Binary file not shown.

‎data/test2/res_para/plate.xmf

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
<?xml version="1.0" ?>
2+
<Xdmf>
3+
<Domain>
4+
<Grid Name="Mesh_Grid">
5+
<Topology NumberOfElements="10000" TopologyType="Hexahedron">
6+
<DataItem DataType="Int" Dimensions="10000 8" Format="HDF">
7+
plate.hdf:/mesh/hex
8+
</DataItem>
9+
</Topology>
10+
<Geometry GeometryType="XYZ">
11+
<DataItem DataType="Float" Dimensions="80000 3" Format="HDF">
12+
plate.hdf:/mesh/coord
13+
</DataItem>
14+
</Geometry>
15+
<Attribute Center="Cell" Name="Scalar_re">
16+
<DataItem DataType="Float" Dimensions="10000" Format="HDF">
17+
plate.hdf:/mesh/sca_re
18+
</DataItem>
19+
</Attribute>
20+
<Attribute Center="Cell" Name="Scalar_im">
21+
<DataItem DataType="Float" Dimensions="10000" Format="HDF">
22+
plate.hdf:/mesh/sca_im
23+
</DataItem>
24+
</Attribute>
25+
</Grid>
26+
</Domain>
27+
</Xdmf>

‎data/test2/res_para/sphere_shell.hdf

2.04 MB
Binary file not shown.

‎data/test2/res_para/sphere_shell.xmf

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
<?xml version="1.0" ?>
2+
<Xdmf>
3+
<Domain>
4+
<Grid Name="Mesh_Grid">
5+
<Topology NumberOfElements="131846" TopologyType="Hexahedron">
6+
<DataItem DataType="Int" Dimensions="131846 8" Format="HDF">
7+
sphere_shell.hdf:/mesh/hex
8+
</DataItem>
9+
</Topology>
10+
<Geometry GeometryType="XYZ">
11+
<DataItem DataType="Float" Dimensions="1054768 3" Format="HDF">
12+
sphere_shell.hdf:/mesh/coord
13+
</DataItem>
14+
</Geometry>
15+
<Attribute Center="Cell" Name="Scalar_re">
16+
<DataItem DataType="Float" Dimensions="131846" Format="HDF">
17+
sphere_shell.hdf:/mesh/sca_re
18+
</DataItem>
19+
</Attribute>
20+
<Attribute Center="Cell" Name="Scalar_im">
21+
<DataItem DataType="Float" Dimensions="131846" Format="HDF">
22+
sphere_shell.hdf:/mesh/sca_im
23+
</DataItem>
24+
</Attribute>
25+
</Grid>
26+
</Domain>
27+
</Xdmf>

‎fortran/computeGREEN_f90_mexed.mexw64

0 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)