Initial commit

MatlabCode/EllipseTwoHolesOnCat.m

+%This script is used to create figure 1 in the paper "Reconstructing 
+%anisotropic conductivities on two-dimensional Riemannian manifolds from
+%power densities"
+
+%The reconstruction problem was solved in Python and we use Matlab to
+%illustrate the reconstructions on a catenoid
+load('EllipseTwoHoles.mat')
+
+map = redblue(255);
+
+N = 128;
+h = 2/N;
+x = -1+h:h:1-h;
+xop = linspace(1-h,-1+h,N-1);
+xtest = 0.5:h:2;
+thetatest = -pi:h:pi;
+[X1,X2] = meshgrid(xtest,thetatest);
+
+xcord = xy2(:,1);
+ycord = xy2(:,2);
+
+Ncord = length(xcord);
+
+x = @(u,v) cosh(u).*cos(v);
+y = @(u,v) cosh(u).*sin(v);
+z = @(u,v) u;
+
+%Limits
+umin = -1.4;
+umax = 0.1;
+vmin = -pi;
+vmax = pi;
+
+
+MDc = 12;
+
+C_c = 'cyan';
+
+set(groot,'defaultAxesTickLabelInterpreter','latex');  
+set(groot,'defaulttextinterpreter','latex');
+set(groot,'defaultLegendInterpreter','latex');
+
+%%
+%The true function \xi in the fundamental domain
+
+figure,
+scatter(xcord,ycord,4,xit,'filled')
+colormap(map)
+set(gca,'fontsize',20)
+axis equal
+caxis([min(xit) max(xit)])
+xlim([-0.63 0.63])
+ylim([-0.95 0.95])
+xticks([-0.5 0 0.5])
+yticks([-0.5 0 0.5])
+
+%%
+%The true function \zeta in the fundamental domain
+
+figure,
+scatter(xcord,ycord,4,zetat,'filled')
+colormap(map)
+set(gca,'fontsize',20)
+axis equal
+caxis([min(zetat) max(zetat)])
+xlim([-0.63 0.63])
+ylim([-0.95 0.95])
+xticks([-0.5 0 0.5])
+yticks([-0.5 0 0.5])
+
+%%
+%The true function (\det \gamma)^{1/2} in the fundamental domain
+
+figure,
+scatter(xcord,ycord,4,detgamsqrtt,'filled')
+colormap(map)
+set(gca,'fontsize',20)
+axis equal
+caxis([min(detgamsqrtt) max(detgamsqrtt)])
+xlim([-0.63 0.63])
+ylim([-0.95 0.95])
+xticks([-0.5 0 0.5])
+yticks([-0.5 0 0.5])
+
+%% 
+%The reconstructed \xi on a catenoid
+
+fig = figure
+h=fsurf(x,y,z,[umin umax vmin vmax],'meshdensity',MDc,'facecolor',[0.7, 0.8, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umin-3e-2 umin vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umax umax+6e-2 vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+set(fig,'units', 'normalized', 'Position', [0.1 0.2 0.15*2 0.3*2])
+
+x2 = @(u,v) 1.03*cosh(log(sqrt(u.^2+v.^2))).*cos(angle(u+sqrt(-1)*v));
+y2 = @(u,v) 1.03*cosh(log(sqrt(u.^2+v.^2))).*sin(angle(u+sqrt(-1)*v));
+z2 = @(u,v) log(sqrt(u.^2+v.^2));
+
+H0g = xi;
+
+colormap(map)
+minc = 2;
+maxc = 3.5;
+caxis([minc maxc])
+num = 255;
+cmap = map;
+cvals = linspace(minc,maxc,num);
+indexH0 = zeros(size(H0g));
+for i = 1:length(H0g)
+    for j = 1:num
+        diffvec = abs(H0g(i)-cvals);
+        [m,mI] = min(diffvec);
+        indexH0(i) = mI;
+    end
+end
+
+scatter3(x2(xcord,ycord),y2(xcord,ycord),z2(xcord,ycord),9,cmap(indexH0,:),'filled')
+cbar = colorbar;
+cbar.Location = 'southoutside';
+cbar.FontSize = 50;
+cbar.TickLabelInterpreter = 'latex';
+haxes = gca;
+colourRange = caxis(haxes);
+cbar.Ticks = linspace(colourRange(1),colourRange(end),4);
+grid off
+axis off
+view(265,50)
+[X,Y] = meshgrid(linspace(-1.2,1.2,100),linspace(-1.2,1.2,100));
+
+%%
+%The reconstructed \zeta on a catenoid
+
+fig=figure
+h=fsurf(x,y,z,[umin umax vmin vmax],'meshdensity',MDc,'facecolor',[0.7, 0.8, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umin-3e-2 umin vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umax umax+6e-2 vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+set(fig,'units', 'normalized', 'Position', [0.1 0.2 0.15*2 0.3*2])
+
+x2 = @(u,v) 1.03*cosh(log(sqrt(u.^2+v.^2))).*cos(angle(u+sqrt(-1)*v));
+y2 = @(u,v) 1.03*cosh(log(sqrt(u.^2+v.^2))).*sin(angle(u+sqrt(-1)*v));
+z2 = @(u,v) log(sqrt(u.^2+v.^2));
+
+H0g = zeta;
+colormap(map)
+minc = min(zetat);
+maxc = max(zetat);
+caxis([minc maxc])
+num = 255;
+cmap = map;
+cvals = linspace(minc,maxc,num);
+indexH0 = zeros(size(H0g));
+for i = 1:length(H0g)
+    for j = 1:num
+        diffvec = abs(H0g(i)-cvals);
+        [m,mI] = min(diffvec);
+        indexH0(i) = mI;
+    end
+end
+
+scatter3(x2(xcord,ycord),y2(xcord,ycord),z2(xcord,ycord),9,cmap(indexH0,:),'filled')
+cbar = colorbar;
+cbar.Location = 'southoutside';
+cbar.FontSize = 50;
+cbar.TickLabelInterpreter = 'latex';
+haxes = gca;
+colourRange = caxis(haxes);
+cbar.Ticks = linspace(colourRange(1),colourRange(end),3);
+grid off
+axis off
+view(265,50)
+[X,Y] = meshgrid(linspace(-1.2,1.2,100),linspace(-1.2,1.2,100));
+
+%%
+%The reconstructed (\det \gamma)^{1/2} on a catenoid
+
+fig=figure
+h=fsurf(x,y,z,[umin umax vmin vmax],'meshdensity',MDc,'facecolor',[0.7, 0.8, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umin-3e-2 umin vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+h=fsurf(x,y,z,[umax umax+6e-2 vmin vmax],'meshdensity',MDc,'facecolor',[0.5, 0.5, 1])
+xlabel('X')
+ylabel('Y')
+zlabel('Z')
+grid off
+set(h,'edgecolor',[0.5, 0.5, 1])
+hold on
+set(fig,'units', 'normalized', 'Position', [0.1 0.2 0.15*2 0.3*2])
+
+hold on
+x2 = @(u,v) 1.05*cosh(log(sqrt(u.^2+v.^2))).*cos(angle(u+sqrt(-1)*v));
+y2 = @(u,v) 1.05*cosh(log(sqrt(u.^2+v.^2))).*sin(angle(u+sqrt(-1)*v));
+z2 = @(u,v) log(sqrt(u.^2+v.^2));
+
+H0g = detgamsqrt;
+
+colormap(map)
+minc = min(detgamsqrtt);
+maxc = 2;
+caxis([minc maxc])
+num = 255;
+cmap = map;
+cvals = linspace(minc,maxc,num);
+indexH0 = zeros(size(H0g));
+for i = 1:length(H0g)
+    for j = 1:num
+        diffvec = abs(H0g(i)-cvals);
+        [m,mI] = min(diffvec);
+        indexH0(i) = mI;
+    end
+end
+
+scatter3(x2(xcord,ycord),y2(xcord,ycord),z2(xcord,ycord),9,cmap(indexH0,:),'filled')
+hold on
+cbar = colorbar;
+cbar.Location = 'southoutside';
+cbar.FontSize = 50;
+cbar.TickLabelInterpreter = 'latex';
+haxes = gca;
+colourRange = caxis(haxes);
+cbar.Ticks = linspace(colourRange(1),colourRange(end),3);
+grid off
+axis off
+view(265,50)
+[X,Y] = meshgrid(linspace(-1.2,1.2,100),linspace(-1.2,1.2,100));
\ No newline at end of file

MatlabCode/redblue.m

+function c = redblue(m)
+%REDBLUE    Shades of red and blue color map
+%   REDBLUE(M), is an M-by-3 matrix that defines a colormap.
+%   The colors begin with bright blue, range through shades of
+%   blue to white, and then through shades of red to bright red.
+%   REDBLUE, by itself, is the same length as the current figure's
+%   colormap. If no figure exists, MATLAB creates one.
+%
+%   For example, to reset the colormap of the current figure:
+%
+%             colormap(redblue)
+%
+%   See also HSV, GRAY, HOT, BONE, COPPER, PINK, FLAG, 
+%   COLORMAP, RGBPLOT.
+
+%   Adam Auton, 9th October 2009
+
+if nargin < 1, m = size(get(gcf,'colormap'),1); end
+
+if (mod(m,2) == 0)
+    % From [0 0 1] to [1 1 1], then [1 1 1] to [1 0 0];
+    m1 = m*0.5;
+    r = (0:m1-1)'/max(m1-1,1);
+    g = r;
+    r = [r; ones(m1,1)];
+    g = [g; flipud(g)];
+    b = flipud(r);
+else
+    % From [0 0 1] to [1 1 1] to [1 0 0];
+    m1 = floor(m*0.5);
+    r = (0:m1-1)'/max(m1,1);
+    g = r;
+    r = [r; ones(m1+1,1)];
+    g = [g; 1; flipud(g)];
+    b = flipud(r);
+end
+
+c = [r g b]; 
\ No newline at end of file

PythonCode/modules.sh

+#!/bin/bash
+
+# Load modules
+#module load FEniCS/2017.2.0-with-petsc-and-slepc
+module load FEniCS/2018.1.0-with-petsc-and-slepc-and-scotch-and-newmpi
+module load scipy/0.19.1-python-3.6.2
+module load matplotlib/2.0.2-python-3.6.2
+module load ffmpeg/3.4
+module load numba/0.35.0-python-3.6.2