code

Signed-off-by: 刘明宏 <liuminghong@mail.sdu.edu.cn>

example2/ano5-20-5-aD/Input.m

@@ -0,0 +1,48 @@
+%%  SET UP PARAMETER
+
+mesh.name ='Mesh_Kalscheuer2008_1ano_v2';% Contains grid information for the model, output by GeoMesh
+[mesh.node2coord,mesh.elem2node,mesh.entity] = load_comsol(mesh.name);
+% load_mesh;
+%% Set the conductivity
+P.airID     = [2  4];      %  Domain ID of air
+P.sig_air   = 1e-8;        % Conductivity of air S/m
+P.sigp      = 1./[10000];  % Background (residual domain) conductivity
+P.sig_anoID = [5];         % Domain number of the abnormal body
+P.sigax    = 1./[10000];   % Conductivity of the abnormal body 
+P.sigay    = 1./[10000] ;
+P.sigaz    = 1./[10000];
+
+P.sig = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.sigax,P.sigay,P.sigaz,P.sig_air,P.sigp) ;
+
+% Set the Euler rotation Angle of the abnormal body 
+P.aS   = [aS];     
+P.aD   = [aD];
+P.aL   =  [0];
+P.AniAngle = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.aS,P.aD,P.aL,0,0) ;
+P.sig  = get_sig_ani(P.sig(:,1),P.sig(:,2),P.sig(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));
+
+%% Set magnetic permitivity
+P.mu0 = 4*pi*1e-7 ;  % V·s/ (A·m)
+P.mu_anoID    = [];        %: anomalous domain ID for comsol mesh,
+P.mura        = [];
+P.mur = Set_material_3Ani(mesh.entity,P.mu_anoID,P.airID,P.mura,P.mura,P.mura,1,1) ;
+P.mur = P.mur(:,1);
+P.mu  = P.mur.*P.mu0 ;
+
+%% Set dielectric permittivity
+P.epsr_air   = 1;    
+P.epsrp      = 1;
+P.eps_anoID  = [5];        % Domain ID of the abnormal body    
+P.epsax      = [epsrx];   % dielectric permittivity of the abnormal body
+P.epsay      = [epsry];   
+P.epsaz      = [epsrz];
+P.epsr = Set_material_3Ani(mesh.entity,P.eps_anoID,P.airID,P.epsax,P.epsay,P.epsaz,P.epsr_air,P.epsrp) ;
+P.epsr = get_sig_ani(P.epsr(:,1),P.epsr(:,2),P.epsr(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));
+P.eps0 = 8.854187817620389*10^-12 ;
+P.eps  = P.epsr *P.eps0 ;
+
+%% Set frequency
+freq  = logspace(4,6,11);
+omega = 2*pi*freq;
+nfreq = size(freq,2);
+

example2/ano5-20-5-aD/Mapdata_xyplane.m

@@ -0,0 +1,56 @@
+function Mapdata_xyplane(data,rec,comp,F,name)
+
+% load rec.mat
+pointsize =100;
+freq  = unique(data.Frequency);
+nfreq =size(freq,1);
+
+for ifreq =9%6:nfreq
+    d            = data(data.Frequency==freq(ifreq),:);
+    d_comp = F(d(:,comp).Variables);
+
+    d_comp = reshape(smoothdata(reshape(d_comp,41,[]), "lowess",5),[],1); %"gaussian"  "movmedian" "movmean" "lowess"
+
+    % figure
+    scatter(rec(d.RX,1),rec(d.RX,2),pointsize, d_comp,'s',"filled");
+    % contourf(reshape(rec(:,2),41,[]),reshape(rec(:,3),41,[]),reshape( d_comp,41,[]),50,'LineStyle','none');
+
+
+    hold on;plot([-50 -50 50 50 -50],[50 -50 -50 50 50],'--w',LineWidth=2)
+    % hold on ; plot([-50 -50 50 50 -50],[20 65 65 20 20],'--w',LineWidth=2)
+    % Eul = plot([0 0],[-80 80],'g',[-80 80],[0 0],'r',LineWidth=2);
+    % rotate(Eul,[0 0 1],30)
+
+    if ismember(comp,{'phaXX' 'phaXY' 'phaYX' 'phaYY' } )
+        clim([20,45]);
+    elseif ismember(comp,{ 'rhoXY' 'rhoYX'  } )
+        clim([2000,11000]);
+    elseif ismember(comp,{'rhoXX' ;'rhoYY' } )
+        clim([0,200]);
+    elseif ismember(comp,{'Tzx' ;'Tzy' } )
+        clim([-0.12,0.12]);
+    end
+
+    % colorbar;
+    colormap(jet );
+    set(gca,'FontSize',14,'FontWeight','bold','LineWidth',1.5)
+
+    set(gca,"YDir","reverse")
+    box on;
+    axis equal
+    % title([num2str(round(freq(ifreq)/1000)) ' kHz ' comp])
+    % ylabel('Y(m)'); xlabel('X(m)')
+    set(gca,'XTickLabel',[],'YTickLabel',[])
+
+    ylim([-210,210])
+    xlim([-210,210])
+
+    set(gca,'LooseInset',get(gca,"TightInset"))
+    % % set(gca,'LooseInset',[0.02 0.02 0.05 0.02])
+    % set(gcf,'Position',[680 578 560 283])
+    % set(gcf,'Position', [799 578 441 420])
+
+    % saveas(gcf, [num2str(round(freq(ifreq)/1000)) ' kHz ' func2str(F) ' ' comp name],'fig')
+    % saveas(gcf, [num2str(round(freq(ifreq)/1000)) ' kHz ' func2str(F) ' ' comp name],'tif')
+
+end