diff --git a/.gitignore b/.gitignore
index 41379149..5068d811 100644
--- a/.gitignore
+++ b/.gitignore
@@ -198,3 +198,7 @@ matlab/MagTense_matlab_project.prj
*.vsidx
*.lock
+
+python/examples/micromagnetism/.idea/inspectionProfiles/
+
+python/examples/micromagnetism/.idea/
diff --git a/buildMagTenseMEX.m b/buildMagTenseMEX.m
index 2c772a6d..aeb73e8e 100644
--- a/buildMagTenseMEX.m
+++ b/buildMagTenseMEX.m
@@ -56,7 +56,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
if (ispc)
build_str_NO_CUDA_MagTenseMicroMag = [build_str '_no_CUDA'];
if (USE_CUDA)
- CUDA_str = '''-Lc:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v11.6/lib/x64/'' -lcublas -lcudart -lcuda -lcusparse';
+ CUDA_str = '''-Lc:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v12.1/lib/x64/'' -lcublas -lcudart -lcuda -lcusparse';
build_str_MagTenseMicroMag = build_str;
else
CUDA_str = '';
@@ -115,7 +115,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
%% MagTenseLandauLifshitzSolver_mex
Source_str = 'source/MagTenseMEX/MagTenseMEX/MagTenseLandauLifshitzSolver_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagTenseMicroMag_str ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' CUDA_str ' ' CVODE_str ' ' MKL_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb']);
@@ -127,7 +127,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
MagTenseMicroMag_str = ['-Lsource/MagTenseMicroMag/' build_str_NO_CUDA_MagTenseMicroMag '/ ' MagTenseMicroMag_lib_str ' -Isource/MagTenseMicroMag/' build_str_MagTenseMicroMag '/'];
Source_str = 'source/MagTenseMEX/MagTenseMEX/MagTenseLandauLifshitzSolver_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagTenseMicroMag_str ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' CVODE_str ' ' MKL_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['MagTenseLandauLifshitzSolver_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/MagTenseLandauLifshitzSolverNoCUDA_mex.mex' MEX_str '64.pdb']);
@@ -139,7 +139,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
%% IterateMagnetization_mex
Source_str = 'source/MagTenseMEX/MagTenseMEX/IterateMagnetization_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['IterateMagnetization_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/IterateMagnetization_mex.mex' MEX_str '64.pdb']);
@@ -150,7 +150,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
%% getHFromTiles_mex
Source_str = 'source/MagTenseMEX/MagTenseMEX/getHFromTiles_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' TileDemagTensor_str ' ' NumericalIntegration_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['getHFromTiles_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getHFromTiles_mex.mex' MEX_str '64.pdb']);
@@ -161,7 +161,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
%% getNFromTile_mex
Source_str = 'source/MagTenseMEX/MagTenseMEX/getNFromTile_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' DemagField_str ' ' NumericalIntegration_str ' ' TileDemagTensor_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['getNFromTile_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getNFromTile_mex.mex' MEX_str '64.pdb']);
@@ -172,7 +172,7 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
%% getMagForce_mex
Source_str = 'source/MagTenseMEX/MagTenseMEX/getMagForce_mex.f90';
mex_str = ['mex ' compiler_str ' ' Debug_flag ' ' MagneticForceIntegrator_str ' ' DemagField_str ' ' NumericalIntegration_str ' ' TileDemagTensor_str ' ' Source_str ' ' Options_str];
-eval(mex_str)
+eval_MEX(mex_str)
if ~USE_RELEASE
pause(pause_time)
movefile(['getMagForce_mex.mex' MEX_str '64.pdb'],['matlab/MEX_files/getMagForce_mex.mex' MEX_str '64.pdb']);
@@ -181,3 +181,26 @@ function buildMagTenseMEX(USE_RELEASE,USE_CUDA,USE_CVODE)
movefile(['getMagForce_mex.mex' MEX_str '64'],['matlab/MEX_files/getMagForce_mex.mex' MEX_str '64']);
end
+
+function eval_MEX(mex_str)
+
+try
+ eval(mex_str);
+catch ME
+ if (strcmp(ME.message(91:117),'mt : general error c101008d'))
+ fail_mex = true;
+ while fail_mex
+ try
+ disp('Microsoft manifest tool error - retrying')
+ eval(mex_str);
+ fail_mex = false;
+ catch
+ continue
+ end
+ end
+ else
+ disp(ME.message)
+ rethrow(ME)
+ end
+end
+end
diff --git a/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Standard_problem_4_unstructured_cart.m b/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Standard_problem_4_unstructured_cart.m
index 0eedefd1..34048514 100644
--- a/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Standard_problem_4_unstructured_cart.m
+++ b/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Standard_problem_4_unstructured_cart.m
@@ -81,7 +81,6 @@
addpath('../../../MEX_files');
addpath('../../../util');
-addpath('../../../micromagnetism_matlab_only_implementation');
%% --------------------------------------------------------------------------------------------------------------------------------------
%% ------------------------------------------------------------------- MAGTENSE ---------------------------------------------------------
diff --git a/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Std_prob_4_unstructured_mesh_grains_6_res_80_20_ref_2.mat b/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Std_prob_4_unstructured_mesh_grains_6_res_80_20_ref_2.mat
index 9f492128..40b097b2 100644
Binary files a/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Std_prob_4_unstructured_mesh_grains_6_res_80_20_ref_2.mat and b/matlab/examples/Micromagnetism/mumag_micromag_Std_problem_4/Std_prob_4_unstructured_mesh_grains_6_res_80_20_ref_2.mat differ
diff --git a/matlab/util/cartesianUnstructuredMeshPlot.m b/matlab/util/cartesianUnstructuredMeshPlot.m
index 2a8c00a3..f3352f54 100644
--- a/matlab/util/cartesianUnstructuredMeshPlot.m
+++ b/matlab/util/cartesianUnstructuredMeshPlot.m
@@ -24,6 +24,7 @@ function cartesianUnstructuredMeshPlot(pos, dims, GridInfo, iIn, fnamesave,hF)
%%
cols = [hsv(numel(iIn)-1);[1,1,1].*0.4];
iOne = sum(abs(GridInfo.TheSigns),1) == 1 ;
+% plot3(GridInfo.Xf(iOne),GridInfo.Yf(iOne),GridInfo.Zf(iOne),'.')
iOneInd = find(iOne) ;
TheLW = 1.5 ;
for ik=1:numel(iOneInd)
@@ -48,21 +49,34 @@ function cartesianUnstructuredMeshPlot(pos, dims, GridInfo, iIn, fnamesave,hF)
TheLS = '-' ;
end
end
+
if abs(GridInfo.fNormX(k))
- hP(ik) = patch(pos(n,1)+[1,1,1,1].*GridInfo.fNormX(k).*dims(n,1)./2,...
+ if isfield(GridInfo,'DimsF')
+ hP(ik) = patch(GridInfo.Xf(k)+[0,0,0,0],GridInfo.Yf(k)+([0,0,1,1]-1/2).*GridInfo.DimsF(k,2),GridInfo.Zf(k)+([0,1,1,0]-1/2).*GridInfo.DimsF(k,3),ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ else
+ hP(ik) = patch(pos(n,1)+[1,1,1,1].*GridInfo.fNormX(k).*dims(n,1)./2,...
pos(n,2)+([0,0,1,1]-1/2).*dims(n,2),...
pos(n,3)+([0,1,1,0]-1/2).*dims(n,3),ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ end
end
if abs(GridInfo.fNormY(k))
- hP(ik) = patch(pos(n,1)+([0,0,1,1]-1/2).*dims(n,1),...
+ if isfield(GridInfo,'DimsF')
+ hP(ik) = patch(GridInfo.Xf(k)+([0,1,1,0]-1/2).*GridInfo.DimsF(k,1),GridInfo.Yf(k)+[0,0,0,0],GridInfo.Zf(k)+([0,0,1,1]-1/2).*GridInfo.DimsF(k,3),ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ else
+ hP(ik) = patch(pos(n,1)+([0,0,1,1]-1/2).*dims(n,1),...
pos(n,2)+[1,1,1,1].*GridInfo.fNormY(k).*dims(n,2)./2,...
pos(n,3)+([0,1,1,0]-1/2).*dims(n,3),ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ end
end
if abs(GridInfo.fNormZ(k))
- hP(ik) = patch(pos(n,1)+([0,1,1,0]-1/2).*dims(n,1),...
+ if isfield(GridInfo,'DimsF')
+ hP(ik) = patch(GridInfo.Xf(k)+([0,0,1,1]-1/2).*GridInfo.DimsF(k,1),GridInfo.Yf(k)+([0,1,1,0]-1/2).*GridInfo.DimsF(k,2),GridInfo.Zf(k)+[0,0,0,0],ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ else
+ hP(ik) = patch(pos(n,1)+([0,1,1,0]-1/2).*dims(n,1),...
pos(n,2)+([0,0,1,1]-1/2).*dims(n,2),...
pos(n,3)+[1,1,1,1].*GridInfo.fNormZ(k).*dims(n,3)./2,ik,ColProperty,thisCol,'linestyle',TheLS,'linewidth',TheLW) ;
+ end
end
'' ;
% drawnow ;
diff --git a/matlab/util/cartesianUnstructuredMeshPlotSolution.m b/matlab/util/cartesianUnstructuredMeshPlotSolution.m
new file mode 100644
index 00000000..7445b8d3
--- /dev/null
+++ b/matlab/util/cartesianUnstructuredMeshPlotSolution.m
@@ -0,0 +1,228 @@
+function [hs,hq] = cartesianUnstructuredMeshPlotSolution(GeomTRI,mesh,IntPlane,Mx_arr,My_arr,Mz_arr,hF)
+%
+% CartesianUnstructuredMeshPlotSolution(GeomTRI,mesh,IntPlane,Mx_arr,My_arr,Mz_arr,hF)
+%
+% Intersects the geometry with a plane, and shows a cutaway view of the
+% solution.
+%
+% GeomTri is the triangulation of the geometry (can be obtained using
+% convhulln. For this function to work, the geometry is assumed to be
+% convex anyway.
+% mesh is the Cartesian Unstructured Mesh corresponding to the solution
+% IntPlane is a structure with two fields (each of which is a 1 by 3
+% array): IntPlane.Normal is the normal to the plane
+% IntPlane.Point is a point passing by the plane
+% Mx_arr, My_arr, and Mz_arr are the arrays of the solution
+%
+% Requires the PDE toolbox (used for creating a tetrahedral mesh -
+% and the correponding triangular surface mesh on which
+% the solution is plotted)
+
+%%
+X = [min(GeomTRI.Points(:,1)),max(GeomTRI.Points(:,1))] ;
+Y = [min(GeomTRI.Points(:,2)),max(GeomTRI.Points(:,2))] ;
+Z = [min(GeomTRI.Points(:,3)),max(GeomTRI.Points(:,3))] ;
+IntPlaneN = norm(IntPlane.Normal) ;
+IntPlane.Normal = IntPlane.Normal./IntPlaneN ;
+
+%%
+ik = [1,2,3,1] ;
+
+% keep only the points on the negative side of the plane
+iIn = sum((GeomTRI.Points-repmat(IntPlane.Point,size(GeomTRI.Points,1),1)).*repmat(IntPlane.Normal,size(GeomTRI.Points,1),1),2)<0 ;
+% Cycle over each segment
+jN = 0 ; clear xNew
+for j=1:size(GeomTRI.ConnectivityList,1)
+ for k=1:3
+ % find the intersection point between the segment and the
+ % intersecting plane
+ xA = GeomTRI.Points(GeomTRI.ConnectivityList(j,ik(k)),:) ;
+ xB = GeomTRI.Points(GeomTRI.ConnectivityList(j,ik(k+1)),:) ;
+ t = -sum((xA-IntPlane.Point).*IntPlane.Normal)/sum((xB-xA).*IntPlane.Normal) ;
+ xInt = xA + t.*(xB-xA) ;
+ if (t>=0) & (t<=1) ; plot3(xInt(1),xInt(2),xInt(3),'.') ;
+ jN = jN + 1 ;
+ xNew(jN,:) = xInt ;
+
+
+ end
+ end
+end
+
+BisectGeomTRI.Points = GeomTRI.Points ;
+BisectGeomTRI.Points(find(~iIn),:) = [] ;
+BisectGeomTRI.Points = [BisectGeomTRI.Points;xNew] ;
+
+BisectGeomTRI.ConnectivityList = convhulln(BisectGeomTRI.Points) ;
+
+
+BiSectVorStruct.pos = sum(BisectGeomTRI.Points,1) ;
+BiSectVorStruct.vorvx{1} = BisectGeomTRI.Points ;
+BiSectVorStruct.ThatK{1} = BisectGeomTRI.ConnectivityList ;
+
+VoronoiStruct.pos = sum(GeomTRI.Points,1) ;
+VoronoiStruct.vorvx{1} = GeomTRI.Points ;
+VoronoiStruct.ThatK{1} = GeomTRI.ConnectivityList ;
+
+
+% PlotVoronoiGeometry(BiSectVorStruct,X,Y,Z,'BisectedHexagonalPlatelet01',hF,[.5,.5,.5])
+BisectGeomTRI = triangulation(BisectGeomTRI.ConnectivityList,BisectGeomTRI.Points) ;
+stlwrite(BisectGeomTRI,'BisectGeomTRI.stl')
+model = createpde(1);
+importGeometry(model,'BisectGeomTRI.stl');
+delete('BisectGeomTRI.stl') ;
+generateMesh(model,'GeometricOrder','linear','Hmax',diff(X)/50) ;
+hF0 = figure ;
+h = pdeplot3D(model) ;
+SurfMesh = triangulation(h.Faces,h.Vertices) ;
+delete(hF0)
+
+%% Mesh for quiver plots
+generateMesh(model,'GeometricOrder','linear','Hmin',diff(X)/10) ;
+hF0 = figure ;
+h = pdeplot3D(model) ;
+QuiverMesh = triangulation(h.Faces,h.Vertices) ;
+inQuiv = sum((QuiverMesh.Points-repmat(IntPlane.Point,size(QuiverMesh.Points,1),1)).*repmat(IntPlane.Normal,size(QuiverMesh.Points,1),1),2)> (-.01*diff(X)) ;
+delete(hF0)
+%%
+
+Fmx = scatteredInterpolant(mesh.pos_out(:,1),mesh.pos_out(:,2),mesh.pos_out(:,3),Mx_arr) ;
+Fmy = scatteredInterpolant(mesh.pos_out(:,1),mesh.pos_out(:,2),mesh.pos_out(:,3),My_arr) ;
+Fmz = scatteredInterpolant(mesh.pos_out(:,1),mesh.pos_out(:,2),mesh.pos_out(:,3),Mz_arr) ;
+
+Mx = Fmx(SurfMesh.Points(:,1),SurfMesh.Points(:,2),SurfMesh.Points(:,3)) ;
+My = Fmy(SurfMesh.Points(:,1),SurfMesh.Points(:,2),SurfMesh.Points(:,3)) ;
+Mz = Fmz(SurfMesh.Points(:,1),SurfMesh.Points(:,2),SurfMesh.Points(:,3)) ;
+
+
+MxQ = Fmx(QuiverMesh.Points(:,1),QuiverMesh.Points(:,2),QuiverMesh.Points(:,3)) ;
+MyQ = Fmy(QuiverMesh.Points(:,1),QuiverMesh.Points(:,2),QuiverMesh.Points(:,3)) ;
+MzQ = Fmz(QuiverMesh.Points(:,1),QuiverMesh.Points(:,2),QuiverMesh.Points(:,3)) ;
+MxQ(~inQuiv) = 0 ; MyQ(~inQuiv) = 0 ; MzQ(~inQuiv) = 0 ;
+QuiverMesh2.Points = QuiverMesh.Points + (.01*diff(X)).*repmat(IntPlane.Normal,size(QuiverMesh.Points,1),1) ;
+ThisColTriangle = ColorFromHorPsiTheta01(Mx,My,Mz) ;
+%%
+
+
+if ~exist('hF','var')
+ hF = figure('position',[0 0 600 600],'Color',[1 1 1]);
+ ppsz = 1.*[20,19] ;
+ ppps = 1.*[0,0,20,19] ;
+ set(gcf,'PaperUnits','centimeters','PaperSize',ppsz,'PaperPosition',ppps) ;
+else
+ ppsz=get(gcf,'PaperSize');
+ ppps=get(gcf,'PaperPosition');
+ if isequal(get(hF,'type'),'figure')
+ set(hF,'PaperUnits','centimeters')
+
+ figure(hF)
+% clf(hF)
+ else
+ axes(hF) ;
+ end
+end
+
+hold on
+
+axis('equal')
+
+%%
+
+
+% PlotVoronoiGeometry(VoronoiStruct,X,Y,Z,'HexagonalPlatelet01',hF,[.5,.5,.5])
+PlotOnlyRealEdges(VoronoiStruct.pos,VoronoiStruct.vorvx,VoronoiStruct.ThatK,1)
+PlotOnlyRealEdges(BiSectVorStruct.pos,BiSectVorStruct.vorvx,BiSectVorStruct.ThatK,2)
+
+
+hs = trisurf(SurfMesh.ConnectivityList,SurfMesh.Points(:,1),SurfMesh.Points(:,2),SurfMesh.Points(:,3),'FaceVertexCData',ThisColTriangle,'linestyle','none','facelighting','flat','SpecularStrength',0,'facealpha',1) ;
+hq = quiver3(QuiverMesh2.Points(:,1),QuiverMesh2.Points(:,2),QuiverMesh2.Points(:,3),MxQ,MyQ,MzQ,'color','k','LineWidth',1.5) ;
+
+view(90,30) ;
+light('position',[0,0,2].*[X(end),Y(end),Z(end)])
+light('position',[0,-2,0].*[X(end),Y(end),Z(end)],'color',.3.*[1,1,1])
+
+set(gca,'xlim',[X(1),X(end)],...
+ 'ylim',[Y(1),Y(end)],...
+ 'zlim',[Z(1),Z(end)],'visible','off')
+
+% xlabel('x') ; ylabel('y') ; zlabel('z') ;
+
+end
+
+function PlotOnlyRealEdges(pos,vorvx,ThatK,LW)
+
+TotFaces = 0 ;
+if ~exist('col','var')
+col = [hsv(size(pos,1));[1,1,1].*0.4];
+end
+for i = 1:size(pos,1)
+ K = ThatK{i} ;
+ % K = convhulln(vorvx{i});
+ % [K,vorvx{i}] = CleanUp(K,vorvx{i}) ;
+ TheSharedEdges{i} = zeros(size(K,1),size(K,1)) ;
+% trisurf(K,vorvx{i}(:,1),vorvx{i}(:,2),vorvx{i}(:,3),'FaceColor',col(i,:),'FaceAlpha',.5,'EdgeAlpha',0,'facelighting','flat','SpecularStrength',0)
+ hold on;
+ % stlwrite(TR,['testSTLwrite',num2str(i),'.stl']) ;
+ for k=1:size(K,1) % cycle over triangles
+ nn{i}{k} = cross((vorvx{i}(K(k,2),:)-vorvx{i}(K(k,1),:)),(vorvx{i}(K(k,3),:)-vorvx{i}(K(k,1),:))) ;
+ nn{i}{k} = nn{i}{k}./norm( nn{i}{k}) ;
+ % line(sigma.*(vorvx{i}(K(k,:),1)-pos(i,1))+pos(i,1),sigma.*(vorvx{i}(K(k,:),2)-pos(i,2)) + pos(i,2),sigma.*(vorvx{i}(K(k,:),3)-pos(i,3)) + pos(i,3))
+ end
+ ThatDot = zeros(size(K,1),size(K,1)) ;
+ for k1=1:size(K,1)
+ for k2=(k1+1):size(K,1)
+ ThatDot(k1,k2) = dot( nn{i}{k1}, nn{i}{k2} ) ;
+ end
+ end
+ ThatDot = ThatDot+ThatDot.' ; % triangles belonging to the same planes
+ ijk = [1,2,3,1] ;
+ CommonEdges = zeros(size(K,1),3) ;
+ for k1=1:size(K,1)
+ for k2=1:size(K,1)
+ [C,ia,ib] = intersect( K(k1,:),K(k2,:)) ;
+ if numel(C)==2 & abs(abs(ThatDot(k1,k2))-1)<1e-9
+ ia = sort(ia) ;
+ ib = sort(ib) ;
+ TheSide1 = isequal(ia,[1;2]).*1 +isequal(ia,[2;3]).*2 + isequal(ia,[1;3]).*3 ;
+ TheSide2 = isequal(ib,[1;2]).*1 +isequal(ib,[2;3]).*2 + isequal(ib,[1;3]).*3 ;
+ CommonEdges(k1,TheSide1) = 1 ;
+ CommonEdges(k2,TheSide2) = 1 ;
+ TheSharedEdges{i}(k1,k2) = 1 ; TheSharedEdges{i}(k2,k1) = 1 ;
+ % TheSharedEdges{i}(k2,TheSide2) = k1 ;
+ '' ;
+ end
+ end
+ end
+
+ '' ;
+ % FindPolygons01(TheSharedEdges{i})
+ CCC = zeros(max(K(:)),max(K(:))) ;
+ for k=1:size(K,1)
+ for j= 1:3
+ if ~CommonEdges(k,j)
+ CCC(K(k,ijk(j)),K(k,ijk(j+1))) = 1 ;
+ CCC(K(k,ijk(j+1)),K(k,ijk(j))) = 1 ;
+ % line(vorvx{i}(K(k,ijk(j:j+1)),1),vorvx{i}(K(k,ijk(j:j+1)),2),vorvx{i}(K(k,ijk(j:j+1)),3),'color','k','linewidth',2)
+ end
+ end
+ end
+ bins = conncomp(graph(TheSharedEdges{i})) ;
+ for kb = 1:max(bins)
+ BB = bins==kb ; kB = K(BB,:) ;
+ [ukB,ikB,ckB] = unique(kB(:)) ;
+
+ GG = graph(CCC(ukB,ukB)) ;
+ nOrd = dfsearch(GG,1) ;
+ iiord = kB(ikB(nOrd)) ;
+ iiord = [iiord(:);iiord(1)];
+ AllFaces{i,kb} = iiord ;
+ TotFaces = TotFaces + 1 ;
+ AllSingleFaces{TotFaces} = iiord ;
+ AllSingleFacesZones(TotFaces) = i ;
+ line(vorvx{i}(iiord,1),vorvx{i}(iiord,2),vorvx{i}(iiord,3),'color','k','linewidth',LW)
+ end
+ '' ;
+
+
+end
+end
\ No newline at end of file
diff --git a/readme_CUDA.txt b/readme_CUDA.txt
index 8039d95f..384f2a8c 100644
--- a/readme_CUDA.txt
+++ b/readme_CUDA.txt
@@ -2,46 +2,29 @@ How to compile changes in the CUDA code in Intel Fortran in MagTense on Windows:
Intel's Fortran Compiler and the PG Fortran compiler (which is used for doing CUDA directly in Fortran) are not compatible. One simply cannot link object files from the two compilers as there is no standard for this.
-The strategy then is to make the CUDA GPU kernel in C++, compile this with the nvcc compiler (that uses MSVC at the core), output to an object file and then compile a C++ wrapper with Intel C++ compiler (icl) that includes and uses the nvcc compiled file. The output here should be an obj file that can be called from Fortran (using the Intel Fortran compiler) via the standard way (iso_c_binding) for calling C++ functions from Fortran.
+The strategy then is to make the CUDA GPU kernel in C++, compile this with the nvcc compiler (that uses MSVC at the core), output to an object file and then compile a C++ wrapper with Intel C++ compiler (icx) that includes and uses the nvcc compiled file. The output here should be an obj file that can be called from Fortran (using the Intel Fortran compiler) via the standard way (iso_c_binding) for calling C++ functions from Fortran.
The following details the required steps if changes has been made to the CUDA code.
Step 0.
Install VS 2019 with Intel's C++ and Fortran compilers (and MKL as this is used also by MagTense) and with MSVC compiler
Install CUDA SDK: https://developer.nvidia.com/cuda-downloads
+ Check that the PATH variable is pointing to the latest version of CUDA. This can be set using the registry at Computer\HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Session Manager\Environment\Path
Step 1.
- From the MS Visual Studio x64 native prompt navigate to the MagTense dir and compile with nvcc:
+ From a standard command prompt, navigate to the MagTense dir "MagTense/source/MagTenseFortranCUDA/cuda" and compile with nvcc, using a specific CUDA version and the latest version of Visual Studio:
- cd MagTense/source/MagTenseFortranCUDA/cuda
- nvcc -c MagTenseCudaBlas.cu
+ "C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.1\bin\nvcc.exe" -c MagTenseCudaBlas.cu -ccbin "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Tools\MSVC\14.36.32532\bin\Hostx64\x64"
Step 2.
- From the Intel 64 prompt (Start Menu -> Intel Parallel Studio -> latest x64 cmd) compile the C++ wrapper with icl including the cuda stuff*:
+ From the Intel 64 prompt (Start Menu -> Intel oneAPI 20XX -> Intel oneAPI command prombt for Intel 64 for Visual Studio XX compile the C++ wrapper with icx including the cuda stuff*:
- icl -c MagTenseCudaBlasICLWrapper.cxx MagTenseCudaBlas.obj /link "c:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v11.6\lib\x64\cublas.lib"
+ icx -c MagTenseCudaBlasICLWrapper.cxx
Step 3.
Compile MagTense as usual (there is a project in the MagTense solution that deals with CUDA)
-
- * CUDA v10.2 is supported (for now) -- simply rename the file MagTenseCudaBlasV10_2.cu -> MagTenseCudaBlas.cu, and point to v10.2 cublas lib when compiling the wrapper.
-
-Old comment - no longer applies:
-Previously there was a problem getting cusparse BLAS to compile (in Windows)
-
-There one needed to need to hack the file cusparse.h (in the folder c:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\include\)
-Find the pre-processor tags:
-
-#if !defined(_WIN32) (line 7346)
-
-#endif // !defined(_WIN32) (line 7705)
-
-Out-comment them (in lack of a better way to fix this - nvcc does compile with x64 as default, so I am clueless as to why this win32 thing is an issue)
-->
-
-//#if !defined(_WIN32) (line 7346)
-
-//#endif // !defined(_WIN32) (line 7705)
+Step 4. (optional)
+ Check the dependencies of the compiled MEX-file using the Visual Studio dumpbin utility:
-At which point one could proceed with the regular compilation
\ No newline at end of file
+ "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Tools\MSVC\14.36.32532\bin\Hostx64\x64\dumpbin.exe" /dependents MagTenseLandauLifshitzSolver_mex.mexw64
\ No newline at end of file
diff --git a/source/DemagField/DemagField/IterateMagnetSolution.f90 b/source/DemagField/DemagField/IterateMagnetSolution.f90
index 7a68e5e3..7ddae80e 100644
--- a/source/DemagField/DemagField/IterateMagnetSolution.f90
+++ b/source/DemagField/DemagField/IterateMagnetSolution.f90
@@ -204,6 +204,12 @@ subroutine iterateMagnetization( tiles, n, stateFunction, n_stf, T, err_max, max
pts = pts + tiles(i)%offset !< Finally add the offset to translate the circ piece
+ case(tileTypeTetrahedron)
+ !! No rotation is needed as the center of the tetrahedron is with respect to the vertices
+ pts(1) = (tiles(i)%vert(1,1) + tiles(i)%vert(1,2) + tiles(i)%vert(1,3) + tiles(i)%vert(1,4))/4.0
+ pts(2) = (tiles(i)%vert(2,1) + tiles(i)%vert(2,2) + tiles(i)%vert(2,3) + tiles(i)%vert(2,4))/4.0
+ pts(3) = (tiles(i)%vert(3,1) + tiles(i)%vert(3,2) + tiles(i)%vert(3,3) + tiles(i)%vert(3,4))/4.0
+
case default
diff --git a/source/MagTenseFortranCuda/cuda/MagTenseCudaBlas.obj b/source/MagTenseFortranCuda/cuda/MagTenseCudaBlas.obj
index 34a07281..fd6cfc90 100644
Binary files a/source/MagTenseFortranCuda/cuda/MagTenseCudaBlas.obj and b/source/MagTenseFortranCuda/cuda/MagTenseCudaBlas.obj differ
diff --git a/source/MagTenseFortranCuda/cuda/MagTenseCudaBlasICLWrapper.obj b/source/MagTenseFortranCuda/cuda/MagTenseCudaBlasICLWrapper.obj
index d4f0dbbb..809fdddc 100644
Binary files a/source/MagTenseFortranCuda/cuda/MagTenseCudaBlasICLWrapper.obj and b/source/MagTenseFortranCuda/cuda/MagTenseCudaBlasICLWrapper.obj differ
diff --git a/source/MagTenseMEX/MagTenseMEX/VTune Amplifier Results/MagTenseMEX/MagTenseMEX.amplxeproj b/source/MagTenseMEX/MagTenseMEX/VTune Amplifier Results/MagTenseMEX/MagTenseMEX.amplxeproj
deleted file mode 100644
index 4912cb35..00000000
--- a/source/MagTenseMEX/MagTenseMEX/VTune Amplifier Results/MagTenseMEX/MagTenseMEX.amplxeproj
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-
- 1570471019
- DTU-18419184102.win.dtu.dk
- windows
- Intel® VTune™ Amplifier 2019 Update 6
- 602217
- 16
- 8
- 1
- 3700000000
- 6
- 85
- 4
- Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz
- avx512
-
-
diff --git a/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90 b/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90
index 1c2ecc2b..e657769a 100644
--- a/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90
+++ b/source/MagTenseMicroMag/LandauLifshitzEquationSolver.f90
@@ -699,7 +699,7 @@ subroutine updateDemagfield( problem, solution)
beta = 0.0
!Hmx = Kxx * Mx
call gemv( problem%Kxx, solution%Mx_s, solution%HmX, alpha, beta )
-
+
beta = 1.0
!Hmx = Hmx + Kxy * My
call gemv( problem%Kxy, solution%My_s, solution%HmX, alpha, beta )
@@ -914,8 +914,9 @@ subroutine ComputeDemagfieldTensor( problem )
CALL SYSTEM_CLOCK(c1)
- call mkl_set_num_threads(problem%nThreadsMatlab)
- call omp_set_num_threads(problem%nThreadsMatlab)
+ !call mkl_set_num_threads(problem%nThreadsMatlab)
+ !call omp_set_num_threads(problem%nThreadsMatlab)
+ call omp_set_num_threads(1)
if ( problem%grid%gridType .eq. gridTypeUniform ) then
@@ -984,10 +985,10 @@ subroutine ComputeDemagfieldTensor( problem )
call displayGUIMessage( 'Averaging the N_tensor not supported for this tile type' )
endif
- !$OMP PARALLEL shared(problem)
- !$omp do private(ind, tile, H, Nout)
-
+ !$OMP PARALLEL shared(problem) private(ind, indx_ele, tile, H, Nout, pts_arr, i)
+
!for each element find the tensor for all evaluation points (i.e. all elements)
+ !$omp do
do i=1,nx
!Setup template tile
tile(1)%tileType = 5 !(for tetrahedron)
@@ -1001,7 +1002,13 @@ subroutine ComputeDemagfieldTensor( problem )
allocate(Nout(1,ntot,3,3))
allocate(H(ntot,3))
- call getFieldFromTiles( tile, H, problem%grid%pts, 1, ntot, Nout, .false. )
+ allocate(pts_arr(ntot,3))
+ pts_arr(:,1) = problem%grid%pts(:,1)
+ pts_arr(:,2) = problem%grid%pts(:,2)
+ pts_arr(:,3) = problem%grid%pts(:,3)
+
+ !call getFieldFromTiles( tile, H, problem%grid%pts, 1, ntot, Nout, .false. )
+ call getFieldFromTiles( tile, H, pts_arr, 1, ntot, Nout, .false. )
!Copy Nout into the proper structure used by the micro mag model
ind = i
@@ -1022,13 +1029,15 @@ subroutine ComputeDemagfieldTensor( problem )
problem%Kzz(:,ind) = sngl(Nout(1,:,3,3))
!Clean up
+ deallocate(pts_arr)
deallocate(Nout)
deallocate(H)
enddo
-
!$omp end do
+
!$OMP END PARALLEL
-
+
+
elseif ( problem%grid%gridType .eq. gridTypeUnstructuredPrisms ) then
!$OMP PARALLEL shared(problem)
@@ -1129,8 +1138,7 @@ subroutine ComputeDemagfieldTensor( problem )
!call displayGUIMessage( 'Kxx(1,1):' )
!write (prog_str,'(f10.3)') problem%Kxx(1,1)
- !call displayGUIMessage( prog_str )
-
+ !call displayGUIMessage( prog_str )
!Write the demag tensors to disk if asked to do so
if ( problem%demagTensorReturnState .gt. DemagTensorReturnMemory ) then
diff --git a/source/MagTenseMicroMag/MagTenseMicroMag.vfproj b/source/MagTenseMicroMag/MagTenseMicroMag.vfproj
index 90f2a59b..101a5c8c 100644
--- a/source/MagTenseMicroMag/MagTenseMicroMag.vfproj
+++ b/source/MagTenseMicroMag/MagTenseMicroMag.vfproj
@@ -6,7 +6,7 @@
-
+
@@ -15,7 +15,7 @@
-
+
@@ -24,7 +24,7 @@
-
+
@@ -34,7 +34,7 @@
-
+
@@ -45,7 +45,7 @@
-
+
@@ -57,7 +57,7 @@
-
+
@@ -68,7 +68,7 @@
-
+
@@ -79,7 +79,7 @@
-
+
@@ -88,8 +88,8 @@
-
-
+
+
diff --git a/source/MagTenseMicroMag/MicroMagParameters.f90 b/source/MagTenseMicroMag/MicroMagParameters.f90
index 9a5d5409..10465f35 100644
--- a/source/MagTenseMicroMag/MicroMagParameters.f90
+++ b/source/MagTenseMicroMag/MicroMagParameters.f90
@@ -105,7 +105,7 @@ module MicroMagParameters
integer :: demag_approximation !> Flag for how to approximate the demagnetization tensor as specified in the parameters below
integer :: demagTensorReturnState !> Flag describing how or if the demag tensor should be returned
integer :: demagTensorLoadState !> Flag describing how or if to load the demag tensor (from disk e.g.)
- integer*4 :: nThreadsMatlab !> Number of threads to use in the OpenMP demag tensor allocation
+ integer :: nThreadsMatlab !> Number of threads to use in the OpenMP demag tensor allocation
integer,dimension(3) :: N_ave !> Number of points to average the demag tensor in in the recieving tile, N_ave(1) = N_x etc
character*256 :: demagTensorFileOut, demagTensorFileIn !> Filename (including path) for output (input) of demag tensor if it is to be returned as a file (demagTensorReturnState >2 and the value is equal to the length of the file including path)
diff --git a/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90 b/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90
index 06d26086..676b5bed 100644
--- a/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90
+++ b/source/TileDemagTensor/TileDemagTensor/TileTriangle.f90
@@ -1,5 +1,6 @@
module TileTriangle
use SPECIALFUNCTIONS
+ use, intrinsic :: ieee_arithmetic
implicit none
real,parameter :: pi=3.14159265359
@@ -169,8 +170,8 @@ function Nxz( r, l, h )
! Returns the Nxz tensor component in the local coordinate system
Nxz = -1./(4.*pi) * ( F_Nxz(r,h,l,h) - F_Nxz(r,0.,l,h) - ( G_Nxz(r,h) - G_Nxz(r,0.) ) )
-
- end function Nxz
+
+ end function Nxz
function F_Nxz( r, yp, l, h )
real :: F_Nxz
@@ -189,6 +190,10 @@ function G_Nxz( r, yp )
real,dimension(3),intent(in) :: r
real,intent(in) :: yp
G_Nxz = atanh( ( r(2)-yp )/ (sqrt(r(1)**2+(r(2)-yp)**2+r(3)**2)) )
+
+ if (.not.ieee_is_finite(G_Nxz)) then ! variable is not finite
+ G_Nxz = 0;
+ endif
end function G_Nxz
!---------------------------------------------------------------------------
@@ -223,6 +228,10 @@ function L_Nyz( r, xp )
real,dimension(3),intent(in) :: r
real,intent(in) :: xp
L_Nyz = atanh( (r(1) - xp) / (sqrt((r(1)-xp)**2+r(2)**2+r(3)**2)) )
+
+ if (.not.ieee_is_finite(L_Nyz)) then ! variable is not finite
+ L_Nyz = 0;
+ endif
end function L_Nyz
!---------------------------------------------------------------------------