##Discretize:
Discretizing a mathematical structure.
21 lines (14 sloc) 320 Bytes
%RectangularDiscretization:
clear all counter = 1;
for(X = 0.05:0.05:1) for(Z = 0.05:0.05:2) x(counter) = -X^3; z(counter) = Z; counter = counter + 1; end end
Connectivity = delaunay(x,z); Vertices = [x',z']
figure trisurf(Connectivity, x,x*0,z)
##Mesh:
Final compilation of one of my Masters degree projects.
275 lines (272 sloc) 7530 Bytes
%Circle Geometry
Vertices = [ 0.354078223658165 -0.766695555351536 -0.298231892827986 0.372870950777673 0.297986811811660 -0.959324251681120 -0.212485229682307 -0.980997388468482 -0.224818474827896 0.140198925104464 0.841199824526692 -0.547087944789703 0.625479155779087 0.783224677636050 0.229859417206465 0.642635182263987 0.390468696828730 0.490410260635561 0.238166773991895 -0.614649349112947 -0.660720113005497 0.485063713088005 -0.403245531692523 0.213599624445504 -0.138199279919743 0.791505007313432 -0.534631205125817 0.353969834110798 0.614013734082178 0.521798768998622 -0.636076834531567 -0.032653119038256 -0.816620666467475 -0.142009839458299 -0.912511619520552 0.416026781542313 -0.469523231869694 -0.702160225580048 -0.075043317877789 -0.832392416359616 -0.361657265597155 -0.549631641488656 0.680256106516730 -0.057283776890683 -0.332340123936992 0.729557433580277 -0.704095891396357 0.715770548526409 0.699206069279343 -0.717745625706996 0.844831753599322 -0.091293768514259 -0.019169470908630 0.122728885287606 -1.002904300029686 -0.024109065305038 -0.415580703585440 0.520474979242920 -0.127148606389532 0.315397379840905 0.158866238978901 -0.420118791485332 -0.451920777949139 0.036343744744883 0.041596225057031 -1.003040501906917 0.456926240538339 0.894225648074915 -0.524666023496280 -0.390210597967998 -0.651442451929566 -0.764228369460929 -0.207615694271768 -0.221573949920743 -0.168169122200641 -0.425538318278553 -0.090792317219241 -0.058138735631788 0.096077925561974 -0.072397712115510 -0.827617928316233 0.567747742268008 -0.339339418726908 0.942923482104062 0.532861875838247 -0.850202071777010 -0.980685837877796 0.209788419940782 0.765937514331123 -0.294917048414618 0.493603043079254 -0.046377916720548 0.298667493857445 -0.051217615338435 0.078899444400723 0.817269475871917 -0.438899892179220 -0.171312358572774 0.002265290252134 1.002588538785814 0.981654705468083 -0.211564257172202 0.012558227632303 0.641679693402978 -0.029471727969854 0.463866058436891 -0.618917445740263 0.167655107615183 -0.199698142726620 0.565075216188424 0.301169823282179 0.305728453028272 0.489775190111495 -0.622515495113159 -0.279678083625676 -0.774456227063590 -0.170850603917301 0.988346757057512 0.922052325363514 -0.396339660125565 -0.803237974855761 -0.600572466408068 0.805613545935361 0.127109016373742 0.742086164564964 0.674764103404251 0.734451878556294 0.331057626600557 -0.743965065717372 -0.367215652902759 -0.625051189972950 -0.220353365680011 -0.283453167168401 -0.039833067270852 0.293843484579637 0.807247835831146 -0.606481310486604 -0.549905553509522 0.458779297138542 0.670055215887352 -0.763852770253396 0.295121054781917 0.207091657711513 -0.244260494835011 0.010494200435293 -0.493885472008037 0.609282014793535 0.145338152144908 0.405584586837439 0.139612949351751 0.672043384413539 -0.488679448089427 0.975609196668914 0.230450957033015 -0.348799782154616 -0.348669887927972 0.397735107242810 -0.237771780624337 -0.817028129455280 0.079430987698179 -0.907089424726309 -0.430870237128993 0.519665506047953 -0.411777185457409 0.088051017286548 0.289608438637462 -0.006487422836260 -0.269430144840844 0.927191469853954 0.382518350910049 0.216881354117256 0.980686245776811 0.173413903739520 0.470499197541963 0.599118508475341 -0.240935317618663 0.520030946326333 0.330166521065558 0.149859766447007 -0.829056778670696 -0.519332090055107 0.858504518485319 -0.529725242147416 0.654726105280486 0.846553908098404 0.535944245154634 -0.970104902703809 -0.255080238223342 0.192850234882940 0.117845481402531 -0.438676608346221 -0.903403777825454 1.003394518502078 0.027488276532593 0.344450338887639 -0.447022594690313 0.044690988771908 -0.671432554675671 -0.162259902428279 -0.627117142435382];
Connectivity = [ 11 14 29 17 28 94 100 38 21 54 14 71 69 19 21 78 21 38 37 38 84 19 36 96 65 69 35 65 17 94 32 49 67 78 49 35 66 49 16 49 78 37 35 21 78 17 80 28 17 66 16 44 80 71 2 55 29 91 92 23 11 71 14 5 12 32 29 92 11 92 29 23 62 22 26 23 13 42 13 23 55 30 2 5 95 27 40 37 78 38 79 47 72 27 5 39 64 74 62 100 21 58 31 84 73 39 67 37 90 33 3 57 10 1 73 38 100 98 10 57 36 19 69 58 19 96 21 35 69 58 96 4 20 58 4 21 19 58 20 4 33 90 20 33 99 73 100 1 90 3 99 100 20 20 100 58 20 90 99 10 99 90 10 98 31 32 16 49 66 35 49 40 39 84 67 49 37 84 39 37 40 27 39 10 73 99 84 38 73 14 54 12 80 17 16 65 81 61 61 36 69 61 69 65 81 65 94 66 65 35 66 17 65 28 80 44 44 71 18 80 16 54 80 54 71 16 32 54 5 32 67 12 54 32 18 71 11 41 11 24 24 92 91 24 11 92 41 18 11 12 2 14 23 42 91 55 23 29 13 59 42 55 52 13 48 8 68 13 50 59 13 52 48 53 83 87 13 48 50 52 87 8 14 2 29 53 52 55 5 2 12 30 55 2 39 5 67 27 83 30 27 30 5 53 55 30 56 83 95 53 30 83 9 87 56 52 53 87 47 95 40 83 27 95 89 56 75 87 83 56 89 75 74 56 95 75 70 68 8 48 52 8 48 86 50 68 34 86 48 68 86 8 9 70 70 9 15 8 87 9 7 70 15 34 68 70 70 7 34 47 75 95 15 63 7 15 93 63 15 64 93 93 64 85 15 9 89 15 89 64 9 56 89 46 74 75 64 89 74 64 62 77 26 22 45 64 77 85 45 51 26 62 74 22 62 26 97 88 82 76 26 51 97 62 97 77 88 46 79 22 74 46 84 72 40 46 75 47 31 72 84 47 40 72 88 79 82 46 47 79 22 88 45 22 46 88 79 98 82 79 72 98 10 31 73 98 72 31 1 3 43 57 1 43 10 90 1 76 57 25 82 98 57 57 43 25 57 76 82 25 6 76 76 45 88 76 6 45 45 6 60 45 60 51];
trisurf(Connectivity, Vertices(:,1), Vertices(:,2), Vertices(:,2)) view([0 0 1]) axis equal
##ABAQUS:
Some mathematics will help ya solve it.
32 lines (24 sloc) 558 Bytes
% Make a triangle close all clear all V1 = [0 0] V2 = [1 0] V3 = [0 2]
plot([V1(1) V2(1) V3(1) V1(1)], [V1(2) V2(2) V3(2) V1(2)],'r-o') hold on % basis function associated with node 1
Plane_Mat = [V1(1) V1(2) 1; V2(1) V2(2) 1; V3(1) V3(2) 1]; RHS = [0 1 0]';
abc = Plane_Mat\RHS;
a = abc(1); b = abc(2); c = abc(3);
for(x=0:0.1:1) for(y=0:0.1:(2-2x)) z = ax + b*y + c; plot3(x, y, z,'.')
end
end
% slope of basis function Centroid = (V1+V2+V3)/3;
quiver(Centroid(1), Centroid(2), a, b)
##Solve:
Final compilation of one of my Masters degree projects.
Now these files we should be able to compare to part A and B... ;) So we've got 4 pictures to check. :D
105 lines (77 sloc) 2970 Bytes
% Finite elements on our circle:
% discretizing the domain Vertices Connectivity
NNodes = length(Vertices); NElem = length(Connectivity);
% Initialize A and f A = zeros(NNodes, NNodes);%spalloc(NNodes, NNodes,6*NNodes); F = zeros(1, NNodes)
% Setup the A-matrix & RHS for(i=1:NElem) % Figure out what the triangle is!!!!! ConnectedNodes = Connectivity(i,:);
V1 = Vertices(ConnectedNodes(1),:);
V2 = Vertices(ConnectedNodes(2),:);
V3 = Vertices(ConnectedNodes(3),:);
plot([V1(1) V2(1) V3(1) V1(1)], [V1(2) V2(2) V3(2) V1(2)],'r-o')
%pause
hold on
% First calculate the area -- Oh No!!!!!!!!!!!!!!!!!!!!!
L1 = [V2-V1,0];
L2 = [V3-V1,0];
Area = norm(1/2*cross( L1, L2));
% Solve the matrix of basis function nodes
Plane_Mat = [V1(1) V1(2) 1; V2(1) V2(2) 1; V3(1) V3(2) 1];
RHS1 = [1 0 0]';
RHS2 = [0 1 0]';
RHS3 = [0 0 1]';
abc1 = Plane_Mat\RHS1;
abc2 = Plane_Mat\RHS2;
abc3 = Plane_Mat\RHS3;
a1 = abc1(1);
b1 = abc1(2);
c1 = abc1(3);
a2 = abc2(1);
b2 = abc2(2);
c2 = abc2(3);
a3 = abc3(1);
b3 = abc3(2);
c3 = abc3(3);
% Interactions between basis functions (gradient * gradient) OH
% NO!!!!!!!!!!!!!
ElementMatrix = Area*[a1*a1 + b1*b1 a2*a1 + b2*b1 a3*a1 + b3*b1;
a1*a2 + b1*b2 a2*a2 + b2*b2 a3*a2 + b3*b2;
a1*a3 + b1*b3 a2*a3 + b2*b3 a3*a3 + b3*b3];
CN = ConnectedNodes;
A(CN(1),CN(1)) = A(CN(1),CN(1)) + ElementMatrix(1,1);
A(CN(1),CN(2)) = A(CN(1),CN(2)) + ElementMatrix(1,2);
A(CN(1),CN(3)) = A(CN(1),CN(3)) + ElementMatrix(1,3);
A(CN(2),CN(1)) = A(CN(2),CN(1)) + ElementMatrix(2,1);
A(CN(2),CN(2)) = A(CN(2),CN(2)) + ElementMatrix(2,2);
A(CN(2),CN(3)) = A(CN(2),CN(3)) + ElementMatrix(2,3);
A(CN(3),CN(1)) = A(CN(3),CN(1)) + ElementMatrix(3,1);
A(CN(3),CN(2)) = A(CN(3),CN(2)) + ElementMatrix(3,2);
A(CN(3),CN(3)) = A(CN(3),CN(3)) + ElementMatrix(3,3);
F(CN(1)) = F(CN(1))+Area/3;
F(CN(2)) = F(CN(2))+Area/3;
F(CN(3)) = F(CN(3))+Area/3;
end
% Boundary conditions for(i = 1:NNodes) V1 = Vertices(i,1); V2 = Vertices(i,2);
if( (V1^2 + V2^2)>.99)
A(i,:) = 0;
A(i,i) = 1;
F(i) = 0;
end
% if( (V1^2 + V2^2)<0.06) & (V1>0.95|V1<-0.95|V2>0.95|V2<-0.95)) % A(i,:) = 0; % A(i,i) = 1; % F(i) = 0; % end
end
Solution = A-F';
figure trisurf(Connectivity, Vertices(:,1), Vertices(:,2), Solution)