Prev: ODE Event- Run on first iteration
Next: compiling C program in Matlab using windows
From: Rui Silva on 18 Jun 2008 22:24
Someone pls reply...any help or hint would be great...

"Rui Silva" <el_tuka(a)iol.pt> wrote in message
<g3bhnq$1sa$1(a)fred.mathworks.com>...
> Hi ppl!
>
> I have been trying to mesh a circular cross section, divided
> in 2 shells ou layers, and i believe it's easier to use
> polar coordinates, to specify the position of my 16
electrodes.
>
> I am using this code:
>
> num_elec=16;
> num_shells=2;
> %element numbers
> nel=num_elec*(num_shells+1);
> %nodes per element
> nnel=3;
> %DOF per node
> ndof=1;
> %Total number of nodes
> nnode=num_elec*num_shells+1;
> %Total DOF
> sdof=nnode*ndof;
> %
> %Nodal coordinates
> angle_min=-pi;
> angle_inc=2*pi/(num_elec);
> angle_max=pi;
> %
> theta=(angle_min:angle_inc:angle_max)';
> %
> radius_min=0;
> radius_max=1;
> radius_inc=-(radius_max/num_shells);
> %
> r=(radius_max:radius_inc:radius_min)';
> %
> lr=length(r);
> lt=length(theta);
> k=1;
> for i=1:lr
> for j=1:lt
> nodesxy(k,:)=[rr(i)*cos(tt(j)),rr(i)*sin(tt(j))];
> k=k+1;
> end
> end
> %
> ....
> and when i try to plot it so i get a circular cross section
> using delaunay, i just get a quadrangular structure, with
> the nodes all out of order...
>
> And also, i am meshing from the outside to the inside, and i
> can't figure out how to tell matlab just to create the
> origin when the radius=0...
>
> If someone can help, i would appreciate :P
>
> Cheers!
>
>

From: John D'Errico on 18 Jun 2008 22:54
"Rui Silva" <el_tuka(a)iol.pt> wrote in message
<g3cg02$59h$1(a)fred.mathworks.com>...
> Someone pls reply...any help or hint would be great...
>
> "Rui Silva" <el_tuka(a)iol.pt> wrote in message
> <g3bhnq$1sa$1(a)fred.mathworks.com>...
> > Hi ppl!
> >
> > I have been trying to mesh a circular cross section, divided
> > in 2 shells ou layers, and i believe it's easier to use
> > polar coordinates, to specify the position of my 16
> electrodes.
> >
> > I am using this code:
> >
> > num_elec=16;
> > num_shells=2;
> > %element numbers
> > nel=num_elec*(num_shells+1);
> > %nodes per element
> > nnel=3;
> > %DOF per node
> > ndof=1;
> > %Total number of nodes
> > nnode=num_elec*num_shells+1;
> > %Total DOF
> > sdof=nnode*ndof;
> > %
> > %Nodal coordinates
> > angle_min=-pi;
> > angle_inc=2*pi/(num_elec);
> > angle_max=pi;
> > %
> > theta=(angle_min:angle_inc:angle_max)';
> > %
> > radius_min=0;
> > radius_max=1;
> > radius_inc=-(radius_max/num_shells);
> > %
> > r=(radius_max:radius_inc:radius_min)';
> > %
> > lr=length(r);
> > lt=length(theta);
> > k=1;
> > for i=1:lr
> > for j=1:lt
> > nodesxy(k,:)=[rr(i)*cos(tt(j)),rr(i)*sin(tt(j))];
> > k=k+1;
> > end
> > end
> > %
> > ....
> > and when i try to plot it so i get a circular cross section
> > using delaunay, i just get a quadrangular structure, with
> > the nodes all out of order...
> >
> > And also, i am meshing from the outside to the inside, and i
> > can't figure out how to tell matlab just to create the
> > origin when the radius=0...
> >
> > If someone can help, i would appreciate :P

Well, after I removed the replicated points
using my consolidator, what is wrong with
what delaunay did to mesh it?

John
From: Rui Silva on 19 Jun 2008 05:10

> Well, after I removed the replicated points
> using my consolidator, what is wrong with
> what delaunay did to mesh it?
>
> John


Well ,i'll try to explain better, i want a mesh to be a
circular one, and with the code i inserted here i just get a
quadrangular structure.

And with some manipulation, i get sort of an circular
structure but the points are all over the place, they don't
follow the connectivity that i coded...if necessary i''l
upload the m-file.

And, i'm sorry, i am kinda new to matlab, but what do you
mean with replicated points and consolidator??

thks for the reply! Cheers!
From: us on 19 Jun 2008 05:27
"Rui Silva":
<SNIP wants to meshing his/her mush...

> I have been trying to mesh a circular cross section...

this is an excellent source for a wide variety of 2d/3d
meshing problems:

www-math.mit.edu/~persson/mesh/persson04mesh.pdf

and the authors, <per-olof persson> and <gilbert strang>,
generously offer their ML package here

http://www-math.mit.edu/~persson/mesh

us
From: Rui Silva on 19 Jun 2008 06:00
lol yes, i already know that one! :) that's a very good
project, one of the most simple an intuitive i have seen so
far...however i really wanted to build my own code, and just
use someone else's project if that fails....sorry...

Anyway, this is the code, i compressed it the post isn't huge:

function V2
num_elec=16;
num_shells=2;
nel=num_elec*(num_shells+1);
nnel=3;
ndof=1;
nnode=num_elec*num_shells+1;
sdof=nnode*ndof;
angle_min=-pi;
angle_inc=2*pi/(num_elec);
angle_max=pi;
theta=(angle_min:angle_inc:angle_max)';
radius_min=0;
radius_max=1;
radius_inc=-(radius_max/num_shells);
r=(radius_max:radius_inc:radius_min)';
lr=length(r);
ltheta=length(theta)-1;
j=1;
for radius=1:lr
for angle=1:ltheta
node(j,:)=[radius*cos(angle) radius*sin(angle)];
j=j+1;
end
end
gcoord(1,1)=node(1,1);gcoord(1,2)=node(1,2);
gcoord(2,1)=node(2,1);gcoord(2,2)=node(2,2);
gcoord(3,1)=node(3,1);gcoord(3,2)=node(3,2);
gcoord(4,1)=node(4,1);gcoord(4,2)=node(4,2);
gcoord(5,1)=node(5,1);gcoord(5,2)=node(5,2);
gcoord(6,1)=node(6,1);gcoord(6,2)=node(6,2);
gcoord(7,1)=node(7,1);gcoord(7,2)=node(7,2);
gcoord(8,1)=node(8,1);gcoord(8,2)=node(8,2);
gcoord(9,1)=node(9,1);gcoord(9,2)=node(9,2);
gcoord(10,1)=node(10,1);gcoord(10,2)=node(10,2);
gcoord(11,1)=node(11,1);gcoord(11,2)=node(11,2);
gcoord(12,1)=node(12,1);gcoord(12,2)=node(12,2);
gcoord(13,1)=node(13,1);gcoord(13,2)=node(13,2);
gcoord(14,1)=node(14,1);gcoord(14,2)=node(14,2);
gcoord(15,1)=node(15,1);gcoord(15,2)=node(15,2);
gcoord(16,1)=node(16,1);gcoord(16,2)=node(16,2);
gcoord(17,1)=node(17,1);gcoord(17,2)=node(17,2);
gcoord(18,1)=node(18,1);gcoord(18,2)=node(18,2);
gcoord(19,1)=node(19,1);gcoord(19,2)=node(19,2);
gcoord(20,1)=node(20,1);gcoord(20,2)=node(20,2);
gcoord(21,1)=node(21,1);gcoord(21,2)=node(21,2);
gcoord(22,1)=node(22,1);gcoord(22,2)=node(22,2);
gcoord(23,1)=node(23,1);gcoord(23,2)=node(23,2);
gcoord(24,1)=node(24,1);gcoord(24,2)=node(24,2);
gcoord(25,1)=node(25,1);gcoord(25,2)=node(25,2);
gcoord(26,1)=node(26,1);gcoord(26,2)=node(26,2);
gcoord(27,1)=node(27,1);gcoord(27,2)=node(27,2);
gcoord(28,1)=node(28,1);gcoord(28,2)=node(28,2);
gcoord(29,1)=node(29,1);gcoord(29,2)=node(29,2);
gcoord(30,1)=node(30,1);gcoord(30,2)=node(30,2);
gcoord(31,1)=node(31,1);gcoord(31,2)=node(31,2);
gcoord(32,1)=node(32,1);gcoord(32,2)=node(32,2);
gcoord(33,1)=node(33,1);gcoord(33,2)=node(33,2);
%
nodes(1,1)=1;nodes(1,2)=2;nodes(1,3)=18;
nodes(2,1)=2;nodes(2,2)=3;nodes(2,3)=19;
nodes(3,1)=3;nodes(3,2)=4;nodes(3,3)=20;
nodes(4,1)=4;nodes(4,2)=5;nodes(4,3)=21;
nodes(5,1)=5;nodes(5,2)=6;nodes(5,3)=22;
nodes(6,1)=6;nodes(6,2)=7;nodes(6,3)=23;
nodes(7,1)=7;nodes(7,2)=8;nodes(7,3)=24;
nodes(8,1)=8;nodes(8,2)=9;nodes(8,3)=25;
nodes(9,1)=9;nodes(9,2)=10;nodes(9,3)=26;
nodes(10,1)=10;nodes(10,2)=11;nodes(10,3)=27;
nodes(11,1)=11;nodes(11,2)=12;nodes(11,3)=28;
nodes(12,1)=12;nodes(12,2)=13;nodes(12,3)=29;
nodes(13,1)=13;nodes(13,2)=14;nodes(13,3)=30;
nodes(14,1)=14;nodes(14,2)=15;nodes(14,3)=31;
nodes(15,1)=15;nodes(15,2)=16;nodes(15,3)=33;
nodes(16,1)=16;nodes(16,2)=1;nodes(16,3)=17;
nodes(17,1)=1;nodes(17,2)=18;nodes(17,3)=17;
nodes(18,1)=2;nodes(18,2)=19;nodes(18,3)=18;
nodes(19,1)=3;nodes(19,2)=20;nodes(19,3)=19;
nodes(20,1)=4;nodes(20,2)=21;nodes(20,3)=20;
nodes(21,1)=5;nodes(21,2)=22;nodes(21,3)=21;
nodes(22,1)=6;nodes(22,2)=23;nodes(22,3)=22;
nodes(23,1)=7;nodes(23,2)=24;nodes(23,3)=23;
nodes(24,1)=8;nodes(24,2)=25;nodes(24,3)=24;
nodes(25,1)=9;nodes(25,2)=26;nodes(25,3)=25;
nodes(26,1)=10;nodes(26,2)=27;nodes(26,3)=26;
nodes(27,1)=11;nodes(27,2)=28;nodes(27,3)=27;
nodes(28,1)=12;nodes(28,2)=29;nodes(28,3)=28;
nodes(29,1)=13;nodes(29,2)=30;nodes(29,3)=29;
nodes(30,1)=14;nodes(30,2)=31;nodes(30,3)=30;
nodes(31,1)=15;nodes(31,2)=32;nodes(31,3)=31;
nodes(32,1)=16;nodes(32,2)=17;nodes(32,3)=32;
nodes(33,1)=17;nodes(33,2)=18;nodes(33,3)=33;
nodes(34,1)=18;nodes(34,2)=19;nodes(34,3)=33;
nodes(35,1)=19;nodes(35,2)=20;nodes(35,3)=33;
nodes(36,1)=20;nodes(36,2)=21;nodes(36,3)=33;
nodes(37,1)=21;nodes(37,2)=22;nodes(37,3)=33;
nodes(38,1)=22;nodes(38,2)=23;nodes(38,3)=33;
nodes(39,1)=23;nodes(39,2)=24;nodes(39,3)=33;
nodes(40,1)=24;nodes(40,2)=25;nodes(40,3)=33;
nodes(41,1)=25;nodes(41,2)=26;nodes(41,3)=33;
nodes(42,1)=26;nodes(42,2)=27;nodes(42,3)=33;
nodes(43,1)=27;nodes(43,2)=28;nodes(43,3)=33;
nodes(44,1)=28;nodes(44,2)=29;nodes(44,3)=33;
nodes(45,1)=29;nodes(45,2)=30;nodes(45,3)=33;
nodes(46,1)=30;nodes(46,2)=31;nodes(46,3)=33;
nodes(47,1)=31;nodes(47,2)=32;nodes(47,3)=33;
nodes(48,1)=32;nodes(48,2)=17;nodes(48,3)=33;
%
N=size(gcoord,1);
t=delaunayn(gcoord);
trimesh(t,gcoord(:,1),gcoord(:,2),zeros(N,1))
view(2),axis equal,axis off,drawnow
%
bcdof(1)=1;
bcval(1)=0;
bcdof(2)=2;
bcval(2)=0;
bcdof(3)=3;
bcval(3)=0.031;
bcdof(4)=4;
bcval(4)=0.008;
bcdof(5)=5;
bcval(5)=0.003;
bcdof(6)=6;
bcval(6)=0.003;
bcdof(7)=7;
bcval(7)=0.001;
bcdof(8)=8;
bcval(8)=0.001;
bcdof(9)=9;
bcval(9)=0.001;
bcdof(10)=10;
bcval(10)=0.001;
bcdof(11)=11;
bcval(11)=0.001;
bcdof(12)=12;
bcval(12)=0.001;
bcdof(13)=13;
bcval(13)=0.002;
bcdof(14)=14;
bcval(14)=0.004;
bcdof(15)=15;
bcval(15)=0.176;
bcdof(16)=16;
bcval(16)=0.192;
%
ff=zeros(sdof,1);
kk=zeros(sdof,sdof);
index=zeros(nnel*ndof,1);
%
for iel=1:nel
nd(1)=nodes(iel,1);
nd(2)=nodes(iel,2);
nd(3)=nodes(iel,3);
x1=gcoord(nd(1),1);y1=gcoord(nd(1),2);
x2=gcoord(nd(2),1);y2=gcoord(nd(2),2);
x3=gcoord(nd(3),1);y3=gcoord(nd(3),2);
[index]=feeldof(nd,nnel,ndof);
[k]=felp2dt3(x1,y1,x2,y2,x3,y3);
[kk]=feasmbl(kk,k,index);
end
%
[kk,ff]=feaplyc2d(kk,ff,bcdof,bcval);
femsol=kk\ff;
num=1:1:sdof;
results=[num' femsol];
%%
function [k]=felp2dt3(x1,y1,x2,y2,x3,y3)
%
A=0.5*(x2*y3+x1*y2+x3*y1-x2*y1-x1*y3-x3*y2);
k(1,1)=((x3-x2)*(x3-x2)+(y2-y3)*(y2-y3))/(4*A);
k(1,2)=((x3-x2)*(x1-x3)+(y2-y3)*(y3-y1))/(4*A);
k(1,3)=((x3-x2)*(x2-x1)+(y2-y3)*(y1-y2))/(4*A);
k(2,1)=k(1,2);
k(2,2)=((x1-x3)*(x1-x3)+(y3-y1)*(y3-y1))/(4*A);
k(2,3)=((x1-x3)*(x2-x1)+(y3-y1)*(y1-y2))/(4*A);
k(3,1)=k(1,3);
k(3,2)=k(2,3);
k(3,3)=((x2-x1)*(x2-x1)+(y1-y2)*(y1-y2))/(4*A);
%%
function [index]=feeldof(nd,nnel,ndof)
edof=nnel*ndof;
k=0;
for i=1:edof
start=(nd(i)-1)*ndof;
for j=1:ndof
k=k+1;
index(k)=start+j;
end
end
%%
function [kk]=feasmbl(kk,k,index)
edof=length(index);
%
for i=1:edof
ii=index(i);
for j=1:edof
jj=index(j);
kk(ii,jj)=kk(ii,jj)+k(i,j);
end
end
%%
function [kk,ff]=feaplyc2d(kk,ff,bcdof,bcval)
n=length(bcdof);
sdof=size(kk);
%
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j)=0;
end
kk(c,c)=1;
ff(c)=bcval(i);
end

 |  Next  |  Last
Pages: 1 2
Prev: ODE Event- Run on first iteration
Next: compiling C program in Matlab using windows