Skip to content

Commit 7615e6c

Browse files
mntalhamntalha
authored
Update Plot Related Files
1 parent 8f40ac8 commit 7615e6c

15 files changed

+605
-0
lines changed

‎Simulator/CheckMesh.m

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
function CheckMesh(mesh, vector)
2+
crd = mesh.crd;
3+
eqv = mesh.eqv;
4+
%
5+
ncrd = size(crd, 2);
6+
if (isempty(eqv))
7+
maxeqv = ncrd;
8+
else
9+
maxeqv = max(eqv(1, :));
10+
end
11+
%
12+
if (ncrd ~= maxeqv)
13+
error('Mesh coordinates inconsistent with equivalence array.')
14+
end
15+
%
16+
if (nargin == 2)
17+
if (isempty(eqv))
18+
nreduced = ncrd;
19+
else
20+
nreduced = min(eqv(1, :)) - 1;
21+
end
22+
lenvec = length(vector);
23+
if (lenvec ~= nreduced)
24+
error ('Vector does not match mesh.');
25+
end
26+
end

‎Simulator/CubPolytope.m

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
function cubp = CubPolytope()
2+
b1 = tan(pi/8);
3+
b2 = tan(pi/6);
4+
%
5+
n111 = UnitVector([1 1 1; 1 -1 1; -1 1 1; -1 -1 1]');
6+
matrix = [eye(3); n111'];
7+
matrix = [matrix; -matrix];
8+
%
9+
rhs = [b1 b1 b1 b2 b2 b2 b2];
10+
rhs = [rhs'; rhs'];
11+
%
12+
% Compute Vertices.
13+
%
14+
x = b1; z = b1^2;
15+
vertices = [...
16+
x x z; x z x; z x x;
17+
-x x z; -x z x; -z x x;
18+
-x -x z; -x -z x; -z -x x;
19+
x -x z; x -z x; z -x x;
20+
x x -z; x z -x; z x -x;
21+
-x x -z; -x z -x; -z x -x;
22+
-x -x -z; -x -z -x; -z -x -x;
23+
x -x -z; x -z -x; z -x -x;
24+
]';
25+
%
26+
% Form polygons for faces.
27+
%
28+
f100 = {...
29+
[ 1 2 11 10 22 23 14 13], ...
30+
[ 4 5 8 7 19 20 17 16], ...
31+
[ 1 3 6 4 16 18 15 13], ...
32+
[ 7 9 12 10 22 24 21 19], ...
33+
[ 2 3 6 5 8 9 12 11], ...
34+
[14 15 18 17 20 21 24 23] ...
35+
};
36+
%
37+
f111 = {[1 2 3], [4 5 6], [7 8 9], [10 11 12], ...
38+
[13 14 15], [16 17 18], [19 20 21], [22 23 24]};
39+
%
40+
faces = {f100{:}, f111{:}};
41+
%
42+
% Form output structure.
43+
%
44+
cubp = PolytopeStructure(matrix, rhs, vertices, faces);

‎Simulator/CycleIndices.m

+7
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
function cycle = CycleIndices(n)
2+
ident = 1:n;
3+
%
4+
cols = repmat(ident, n, 1);
5+
rows = cols';
6+
%
7+
cycle = 1 + mod(cols + rows - 2, n);

‎Simulator/MeshFaces.m

+39
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
function [faces, mult] = MeshFaces(con)
2+
d = size(con, 1); % dimension plus one
3+
n = size(con, 2); % number of elements
4+
dface = d - 1; % dimension of face
5+
nallf = d*n; % total number of faces, including interior ones
6+
%
7+
% efaces is the connectivity for faces within an element
8+
%
9+
efaces = CycleIndices(d);
10+
efaces = efaces(1:dface, :);
11+
%
12+
f = [];
13+
for i=1:d;
14+
f = [f con(efaces(:, i), :)];
15+
end
16+
%
17+
allfaces = sortrows(sort(f)')';
18+
[u, f_of, a_of] = unique(allfaces', 'rows');
19+
%
20+
% Evaluate output args.
21+
%
22+
faces = u';
23+
%
24+
% Multiplicity if requested.
25+
%
26+
% Note: can get rid of for-loop as in SymEquiv
27+
% (or use sparse matrix to accumulate)
28+
%
29+
if (nargout >= 2)
30+
nfaces = size(faces, 2);
31+
wunz = ones(1, nallf);
32+
%
33+
mult = sparse(wunz, a_of, wunz, 1, nfaces);
34+
mult = full(mult); % note sparse accumulates entries
35+
%
36+
% for i=1:nallf
37+
% mult(a_of(i)) = mult(a_of(i)) + 1;
38+
% end
39+
end

‎Simulator/MeshStructure.m

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
function mesh = MeshStructure(crd, con, eqv, varargin)
2+
3+
MyName = mfilename;
4+
%
5+
%-------------------- Defaults and Keyword Arguments
6+
%
7+
optcell = {...
8+
'ElementType', '', ...
9+
'Symmetries', [] ...
10+
};
11+
%
12+
options = OptArgs(optcell, varargin);
13+
%
14+
%-------------------- Main Code
15+
%
16+
if (nargin == 0)
17+
crd = [];
18+
con = [];
19+
eqv = [];
20+
elseif (nargin == 2)
21+
eqv = [];
22+
end
23+
%
24+
mesh = struct('crd', crd, 'con', con, 'eqv', eqv);
25+
%
26+
% Process optional arguments.
27+
%
28+
if nargin > 3
29+
%
30+
if options.ElementType
31+
mesh.etype = ElementTypeStruct(options.ElementType);
32+
end
33+
%
34+
if options.Symmetries
35+
mesh.symmetries = options.Symmetries;
36+
end
37+
%
38+
end

‎Simulator/OnOrOff.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function toggle = OnOrOff(string)
2+
3+
string = lower(string); % to lower case
4+
%
5+
if (strcmp(string, 'on'))
6+
toggle = 1;
7+
elseif (strcmp(string, 'off'))
8+
toggle = 0;
9+
else
10+
error('string is not either ''on'' or ''off''')
11+
end

‎Simulator/OptArgs.m

+16
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function opts = OptArgs(optkeys, optargs)
2+
3+
opts = struct(optkeys{:});
4+
optargs = reshape(optargs, [2 length(optargs)/2]);
5+
%
6+
for k=1:size(optargs, 2)
7+
%
8+
key = optargs{1, k};
9+
val = optargs{2, k};
10+
%
11+
if (isfield(opts, key))
12+
opts.(key) = val;
13+
else
14+
warning('OptArgs:NoKey', ['no such option: ', key])
15+
end
16+
end

‎Simulator/PlotFR.m

+112
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
function PlotFR(mesh, odf, varargin)
2+
3+
CheckMesh(mesh, odf) % check sizes
4+
%
5+
crd = mesh.crd; % Unpack input structures.
6+
con = mesh.con;
7+
eqv = mesh.eqv;
8+
%
9+
%odf = ToAllNodes(odf, eqv);
10+
%
11+
%--------------------Defaults and Options-------------------------------
12+
%
13+
optcell = {...
14+
'Symmetries', 'cubic', ...
15+
'ShowMesh', 'off', ...
16+
'Colormap', [], ...
17+
'BrightenColormap', 0, ...
18+
'NumberOfColors', 64 ...
19+
};
20+
%
21+
opts = OptArgs(optcell, varargin);
22+
%
23+
symtype = opts.Symmetries;
24+
%
25+
% Colormap.
26+
%
27+
cmap = opts.Colormap;
28+
if (isempty(cmap)) % default colormap
29+
cmap = jet(opts.NumberOfColors);
30+
end
31+
bright = opts.BrightenColormap;
32+
cmap = bright + (1 - bright)*cmap;
33+
%
34+
% Surface options.
35+
%
36+
plotsurfopts = {'ShowMesh', opts.ShowMesh};
37+
%
38+
%-------------------- Build figure.
39+
%
40+
f = figure;
41+
%
42+
colormap(cmap);
43+
%
44+
figscale = 1.0; % figure scale
45+
%
46+
% Reshape figure to reflect the 1x2 array of subplots.
47+
%
48+
vertscale = 0.6;
49+
%
50+
pos = get(f, 'Position');
51+
pos(4) = pos(4)*vertscale;
52+
%
53+
set(f, 'Position', pos);
54+
%-------------------- *** Subplot 1 (surface)
55+
%
56+
subplot(1, 2, 1)
57+
%
58+
% The first plot shows the surface, and the range of data
59+
% is complete since the surface plot uses the same nodal point
60+
% array, but only the surface elements.
61+
%
62+
[faces, multiplicity] = MeshFaces(con);
63+
scon = faces(:, (multiplicity == 1));
64+
surfmesh = MeshStructure(crd, scon, eqv);
65+
%
66+
PlotSurface(surfmesh, odf, plotsurfopts{:});
67+
%
68+
% Axes.
69+
%
70+
axis off
71+
%
72+
a11 = gca;
73+
clim = get(a11, 'CLIM'); % save color range for next plot
74+
%
75+
% Colorbar.
76+
%
77+
% This positions the colorbar in the middle of the figure.
78+
% The vertical size is 90% that of the axes, while the default
79+
% width remains the same.
80+
%
81+
% Position = [left bottom width height]
82+
%
83+
cb = colorbar;
84+
pos_a11 = get(a11, 'Position'); % save position after creation of colorbar
85+
poscb = get(cb, 'Position');
86+
%
87+
poscb(4) = 0.9*pos_a11(4);
88+
poscb(1) = 0.5 - 0.5*poscb(3);
89+
poscb(2) = 0.5 - 0.5*poscb(4);
90+
%
91+
set(cb, 'Position', poscb);
92+
%
93+
PlotFRPerimeter(symtype);
94+
%
95+
%-------------------- *** Subplot 2 (slices)
96+
%
97+
subplot(1, 2, 2)
98+
%
99+
% The second plot shows the slices and uses the same data
100+
% as the first plot, so that the colorbar is identical.
101+
%
102+
PlotFRSlices(mesh, ToAllNodes(odf, eqv), symtype)
103+
%
104+
% Adjust axes.
105+
%
106+
a12 = gca;
107+
%
108+
% Make same size as other axes.
109+
%
110+
pos_a12 = pos_a11;
111+
pos_a12(1) = 1 - pos_a12(1) - pos_a12(3);
112+
set(a12, 'CLIM', clim, 'Position', pos_a12);

‎Simulator/PlotFRPerimeter.m

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
function perim = PlotFRPerimeter(symtype)
2+
3+
if (strcmpi(symtype, 'cubic'))
4+
ptope = CubPolytope;
5+
elseif (strcmpi(symtype, 'hexagonal'))
6+
ptope = HexPolytope;
7+
elseif (strcmpi(symtype, 'orthorhombic'))
8+
ptope = OrtPolytope;
9+
else
10+
wid = 'PlotFRPerimeter:symmetry';
11+
msg = 'symmetry type not recognized';
12+
warning(wid, msg)
13+
return
14+
end
15+
%
16+
verts = ptope.vertices;
17+
faces = ptope.faces;
18+
nfaces = length(faces);
19+
%
20+
perim = [];
21+
for i=1:nfaces
22+
list = faces{i};
23+
pts = verts(:, [list list(1)]);
24+
prm = plot3(pts(1, :), pts(2, :), pts(3, :), 'ko-');
25+
perim = [perim; prm];
26+
hold on
27+
end

‎Simulator/PlotFRSlices.m

+57
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
function PlotFRSlices(mesh, odf, type, planes)
2+
3+
if (nargin < 3)
4+
type = 'cubic';
5+
end
6+
%
7+
if (nargin < 4)
8+
planes = SymPlanes(type);
9+
end
10+
%
11+
if isempty(planes)
12+
planes = SymPlanes('cubic');
13+
end
14+
%
15+
for p=planes
16+
[sm, e, ecrd] = SliceMesh(mesh, p.Point, p.Normal);
17+
%[e, ecrd] = tsearchn(mesh.crd', mesh.con', sm.crd');
18+
eodf = odf(mesh.con(:, e));
19+
slice_odf = dot(ecrd, eodf, 1);
20+
%
21+
PlotSurface(sm, slice_odf);
22+
%
23+
end
24+
%
25+
%--------------------Private Functions----------------------------------
26+
%
27+
function planes = SymPlanes(type)
28+
% CUBPLANES -
29+
%
30+
Z = [0; 0; 0];
31+
%
32+
if (strcmpi(type, 'cubic'))
33+
%
34+
E = eye(3);
35+
E1 = E(:, 1);
36+
E2 = E(:, 2);
37+
E3 = E(:, 3);
38+
planes = struct('Point', Z, 'Normal', {E1, E2, E3});
39+
%
40+
elseif (strcmpi(type, 'hexagonal'))
41+
%
42+
n_z = [0 0 1]';
43+
theta = (0:5)*2*pi/5;
44+
n_xy = [cos(theta); sin(theta); zeros(1, 6)];
45+
planes = struct('Point', Z, 'Normal', num2cell([n_z, n_xy], 1));
46+
%
47+
else
48+
%
49+
% Default planes to display.
50+
%
51+
E = eye(3);
52+
E1 = E(:, 1);
53+
E2 = E(:, 2);
54+
E3 = E(:, 3);
55+
planes = struct('Point', Z, 'Normal', {E1, E2, E3});
56+
%
57+
end

0 commit comments

Comments
 (0)