MATLAB: Geometry M-File Troubleshooting

[ephedyn]

I've highlighted the areas which might have problems in red, especially boundary segment 1. I know

x(ii)=interp1([dl(1,1),dl(2,1)],[0 2],s(ii));

should be wrong, but I have no idea what to replace it with such that it works. There shouldn't be a problem with my parameter values I've highlighted it just in case. Does anyone know what's wrong with my code?

Thanks in advance!

function [x,y]=expg2(bs,s)
%EXPG2 Creates a geometry file for an enclosed region of 3 line segments and 1 exp curve.

% Number of boundary segments
nbs=4;

if nargin==0 % Number of boundary segments
    x=nbs;
return
end

dl=[
[COLOR="Red"]  0 0 0 0 % start parameter value
  2 1 1 1 % end parameter value
  0 0 0 0 % left hand region
  1 1 1 1 % right hand region[/COLOR]
];
    
    bs1=bs(:)';
    
if find(bs1<1 | bs1>nbs),
  error('PDE:expg2:InvalidBs', 'Non existent boundary segment number.')
end
    
if nargin==1,
  x=dl(:,bs1);
  return
end

x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
if m==1 && n==1,
  bs=bs*ones(size(s)); % expand bs
elseif m~=size(s,1) || n~=size(s,2),
  error('PDE:expg2:SizeBs', 'bs must be scalar or of same size as s.');
end

if ~isempty(s),

[COLOR="Red"]% boundary segment 1
ii=find(bs==1);
if length(ii)
x(ii)=interp1([dl(1,1),dl(2,1)],[0 2],s(ii));
y(ii)=exp(-s(ii));
end[/COLOR]
    
% boundary segment 2
ii=find(bs==2);
if length(ii)
x(ii)=interp1([dl(1,2),dl(2,2)],[2 2],s(ii));
y(ii)=interp1([dl(1,2),dl(2,2)],[exp(-2) 0],s(ii));
end

% boundary segment 3
ii=find(bs==3);
if length(ii)
x(ii)=interp1([dl(1,3),dl(2,3)],[2 0],s(ii));
y(ii)=interp1([dl(1,3),dl(2,3)],[0 0],s(ii));
end

% boundary segment 4
ii=find(bs==4);
if length(ii)
x(ii)=interp1([dl(1,4),dl(2,4)],[0 0],s(ii));
y(ii)=interp1([dl(1,4),dl(2,4)],[0 2],s(ii));;
end

end

 

[Jon_]

Well,for starters.I copied your code into my version of MATLAB 2007b and ran it. There were two errors neither of which really addresses your topic.
First error was an extra semicolon on line 66.

Second error(s) were in the % boundary segment N, your "if" statements "length(ii)". If instead, you used ~isempty(ii) in it's place.

Also, when running these scripts in editor,you may want to use the "show M-Lint Report" after each run to see what additional errors might exist. But as far as coding is considered, I'm an amateur at best. Sorry, I can't be of much more help.

 

[ephedyn]

Hey,

Thanks for helping. I don't think you're amateurish - surely, you're better than me at this.

I must have accidentally typed an extra colon while I was copying it over. And yep, I replaced the parts with ~isempty(ii) as suggested so it looks like this now (see below).

function [x,y]=expg2(bs,s)
%EXPG2 Creates a geometry file for an enclosed region.

% Number of boundary segments
nbs=4;

if nargin==0 % Number of boundary segments
    x=nbs;
return
end

[COLOR="Red"]dl=[
  0 0 0 0 % start parameter value
  2 1 1 1 % end parameter value
  0 0 0 0 % left hand region
  1 1 1 1 % right hand region
];[/COLOR]
    
    bs1=bs(:)';
    
if find(bs1<1 | bs1>nbs),
  error('PDE:expg2:InvalidBs', 'Non existent boundary segment number.')
end
    
if nargin==1,
  x=dl(:,bs1);
  return
end

x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
if m==1 && n==1,
  bs=bs*ones(size(s)); % expand bs
elseif m~=size(s,1) || n~=size(s,2),
  error('PDE:expg2:SizeBs', 'bs must be scalar or of same size as s.');
end

if ~isempty(s),

[COLOR="Red"]% boundary segment 1
ii=find(bs==1);
if ~isempty(ii)
x(ii)=interp1([dl(1,1),dl(2,1)],[0 2],s(ii));
y(ii)=exp(-s(ii));
end[/COLOR]
    
% boundary segment 2
ii=find(bs==2);
if ~isempty(ii)
x(ii)=interp1([dl(1,2),dl(2,2)],[2 2],s(ii));
y(ii)=interp1([dl(1,2),dl(2,2)],[exp(-2) 0],s(ii));
end

% boundary segment 3
ii=find(bs==3);
if ~isempty(ii)
x(ii)=interp1([dl(1,3),dl(2,3)],[2 0],s(ii));
y(ii)=interp1([dl(1,3),dl(2,3)],[0 0],s(ii));
end

% boundary segment 4
ii=find(bs==4);
if ~isempty(ii)
x(ii)=interp1([dl(1,4),dl(2,4)],[0 0],s(ii));
y(ii)=interp1([dl(1,4),dl(2,4)],[0 2],s(ii));
end

end

But the geometry still isn't working. I'm sorry that I didn't elaborate earlier, but I'm trying to initiate a mesh over it using the following commands:

pdegplot('expg2');
[p,e,t]=initmesh('expg2');
pdemesh(p,e,t), axis equal

which shows the following error:

? Error using ==> pdevoron
Geometry error.

Error in ==> pderespe at 32
  [p,t,c,h]=pdevoron(p,t,c,h,x,y,tol,Hmax,Hgrad);

Error in ==> initmesh at 181
[p,e,t,c]=pderespe(g,p,e,t,c,Hmax,tol,tol2,Hmax,Hgrad);

Do you (or anyone else) know what's wrong?

 

[Jon_]

Well, I looked over it again and the only conclusion thus far is,somewhere, you are suppressing output.

Anytime you get an output as "ans=" you have placed a semicolon where there ought not to be one. Please be advised, I have only one semester exposure to MatLab 2007 and to this date haven't used it since.

I will ask around, to see if any of my collegues know the answer to your question.
Jon_

 

[Jon_]

Ok, when in your command window, type in help pdemesh you will get info regarding your function. If you denote the PDEMESH with another variable say , H, then this may help:

>> help pdevoron
PDEVORON Add points to triangular mesh.

>> help pdemesh
PDEMESH Plot a PDE triangular mesh.

PDEMESH(P,E,T) plots the mesh specified by P, E, and T.

PDEMESH(P,E,T,U) plots the solution column vector U using
a mesh plot. If U is a column vector, node data is
assumed. If U is a row vector, triangle data is assumed.
This command plots the solution substantially faster than the
PDESURF command.

H=PDEMESH(P,E,T) additionally returns handles to the plotted
axes objects.

You probably already know this. This was another attempt (by me) without asking for assistance.

Jon_

 

[Jon_]

The point is,(not to sound condesending), you need to return a value,scalar,etc.
You might want to check your nodes, do you have the proper number? I had written a program some time ago,in mathcad and a similar problem occurred when I was doing array calculations and didn't have the proper dimensions. For every grid you will have
N+1 nodes.

Jon_

 

[Jon_]

Here is more:

>> help pdevoron
PDEVORON Add points to triangular mesh.

>> help initmesh
INITMESH Build an initial PDE triangular mesh.

[P,E,T]=INITMESH(G) returns a triangular mesh using the
geometry specification function G. It uses a Delaunay triangulation
algorithm. The mesh size is determined from the shape of the geometry.

G describes the geometry of the PDE problem. See PDEGEOM for details.

The outputs P, E, and T are the mesh data.

In the point matrix P, the first and second rows contain
x- and y-coordinates of the points in the mesh.

In the edge matrix E, the first and second rows contain indices
of the starting and ending point, the third and fourth rows contain
the starting and ending parameter values, the fifth row contains
the boundary segment number, and the sixth and seventh row
contain the left- and right-hand side subdomain numbers.

In the triangle matrix T, the first three rows contain indices to
the corner points, given in counter clockwise order, and the fourth
row contains the subdomain number.

[P,E,T]=INITMESH(G,'Property',Value,...) can take property-value pairs to
be used by the meshing algorithm. Valid property/value pairs include:

Property Value/{Default} Description
-----------------------------------------------------------
Hmax numeric {estimate} Maximum edge size
Hgrad numeric {1.3} Triangle growth rate
Box on|{off} Preserve bounding box
Init on|{off} Boundary triangulation
Jiggle off|{mean}|min Call JIGGLEMESH
JiggleIter numeric {10} Maximum iterations

The Hmax parameter controls the size of the triangles on the
mesh. INITMESH creates a mesh where no triangle side exceeds
Hmax.

Both the Box and Init parameters are related to the way the
mesh algorithm works. By turning on Box you can get a good idea
of how the mesh generation algorithm works within the bounding box.
By turning on Init you can see the initial triangulation of
the boundaries.

I hope this helps you.
Jon_

 

[Jon_]

I have a question: Is the "triangulation of boundaries" related in any way to "TRI-DIAGONAL" matrices?

 

[ephedyn]

All right, just got up after a horrible day and was amazed that I'm almost late for school again. I must say I really appreciate your help. I'm sorry that I haven't had time to read through everything, I'll get back to your posts when lessons' over.

But to round things up (or rather, down), the reason why I'm suspecting %boundary segment 1 rather than any of the commands etc. which you have pointed out is because everything works if I simply change the lines

% boundary segment 1
ii=find(bs==1);
if ~isempty(ii)
[COLOR="Blue"]x(ii)=interp1([dl(1,1),dl(2,1)],[0 2],s(ii));
y(ii)=exp(-s(ii));[/COLOR]
end

to

% boundary segment 1
ii=find(bs==1);
if ~isempty(ii)
[COLOR="Blue"]x(ii)=interp1([dl(1,1),dl(2,1)],[0 2],s(ii));
y(ii)=interp1([dl(1,1),dl(2,1)],[2 exp(-2)],s(ii));[/COLOR]
end

for a region bounded by straight lines i.e. a polygon (but the problem is I need the curve instead of a straight line for boundary segment 1). That's why it's unlikely that anything else is wrong, but at the same time, I can't think of what to replace those 2 lines with (I've been trying).

But I'll still go through your posts in detail when I'm back. Thanks for extending your help, I'll be really grateful if you let me know if you found anything else.

Be right back.

EDIT:

Oh wait, I think I got it working! It was actually segment 4.

% boundary segment 4
ii=find(bs==4);
if ~isempty(ii)
x(ii)=interp1([dl(1,4),dl(2,4)],[0 0],s(ii));
y(ii)=interp1([dl(1,4),dl(2,4)],[0 [COLOR="Red"]1[/COLOR]],s(ii));
end

Just realized it because I had to change a few things after going through with you, so thanks a lot! Now I hope nothing goes wrong when I solve the PDE.

OK now I'm really running late.

 

[Jon_]

By the way, I have only solved PDE's by hand and never used MatLab to do numerical calculations. It's really cool what you can do, if you know a language. I hope all goes well for you.


Q: Did you right that entire script yourself or did you copy it from a text?
( I ask, 'cause if you wrote yourself, I want to know how I can write programs this complicated as well, for myself and other projects for future use or reference)

I wish the university I currently attend would offer a course in BVP/PDEs as a numerical analysis course in the undergraduate program, but they don't!

 

[Jon_]

If you do get it working, show the entire script/program in this forum, I'd like to copy/paste on my MatLab and see what it does, is that ok??

Jon_

 

[Jon_]

If the script still doesn't work, look at this option:

>> help PDEGEOM
PDEGEOM Geometry M-file.

NE=PDEGEOM is the number of boundary segments

D=PDEGEOM(BS) is a matrix with one column for each boundary segment
specified in BS.
Row 1 contains the start parameter value.
Row 2 contains the end parameter value.
Row 3 contains the label of the left-hand region.
Row 4 contains the label of the right-hand region.

[X,Y]=PDEGEOM(BS,S) produces coordinates of boundary points.
bs specifies the boundary segments and s the corresponding
parameter values. BS can be a scalar.

See also initmesh, refinemesh

I think you already have this, if so, disregard.

Jon_

 

[Jon_]

yes, you already have it it is your "dl= [...];"
ok.

 

[ephedyn]

Yup, I don't mind sharing everything I've got once I'm done with the research. But basically that's the basic geometry m-file.

You also have to write an m-file for the desired boundary conditions. After that all you've got to do is to use. By varying the 4 parameters of assempde, you are able to specify the PDE.

[p,e,t]=initmesh('geometry.file'); %initialise a mesh
[p,e,t]=refinemesh('geometry.file',p,e,t); %refine it
u=assempde('boundaryconds.file',...); %u is your solution of the PDE
pdeplot(...); %this gives you the plot of the solution in 2D

or

pdesurf(...); %this gives you the plot of the solution in 3D

It's nothing really special if you snoop around here and there and patch different tutorials together, but the real difficult part was to figure out what the values of each of the parameters meant, e.g. dl=[...] and b1=[...] (the matrix of the boundary conditions) in my opinion. I did figure out the code from scratch, but the more I worked on it the more I figured out that existing codes are a notch up on perfection, and eventually my code looked so similar to existing ones, I decided to lift the existing ones wholesale ('circleg.m' and 'circleb1.m' were particularly useful for illustrating the idea), so I daren't say it's any creative work, it's just a matter of time spent studying the existing codes.

And yes, I've never had a single course on programming, Matlab included. And it's great to know you're keen about this. I share my mentor's opinion on this matter: when you don't know something, you're more like "unfamiliar with the code" rather than "amateurish". And he was really encouraging on the idea of working on it without any training: he heads a department but when he brought up the use of Matlab, he said a crash course wasn't necessary because "...to be honest, [he] didn't know the code himself," and can't teach us much (I'm sure he's joking to an extent), and even he himself has to refer to the guide.

Actually, I appreciate the analytical solution more (challenge-wise), but both have their practical advantages, so I wish I could solve everything by hand; I'm looking at flow uniformity in this case, which was why I adopted a computational approach.

 

[Jon_]

That is so true. Now, when you access the graphical capablilties of MatLab, you will strike ahead with animations that will make your presentations nicely packaged in a professional format. Man, I wish I had your insight and PAITENCE, to probe into deeper constructs and coding "To boldly go, where no man has gone before". I have had one programming course in ANSI C and to tell you the truth, I haven't found that it hasn't helped me one iota.

Perhaps, it would have been better to take on MatLab without a prior programming course. I feel that C programming has confounded my understanding of "piece-wise" programming, like MatLab coding.

Jon_

 

[Jon_]

... "I have found that it hasn't helped me one iota".

In place of the previous thread

 

[Jon_]

You know, there is another software package you might consider, if you haven't already been redirected to it by now. It's called PDEflex. It's a program that follows closely to the structure of "C" programming, that is if you want to learn the elements of a programming language.

You can download the program for FREE, from: www.pdesolutions.com
I have used this program in the past, for one project on heat conduction. It's basically an interactive CAD program that incorporates the FEM approach.

I have used CAD programs in the past and this is similar to a basic CAD. You can create your own boundaries by following the coded samples from previous work of others.
It's really cool too.It computes and animates before your eyes. There is the student version demo and the professional version for a few thousand dollars.

I used it to model the heat conduction from a physical example of a 2N3904 power(modeled as a circular heat sink) transistor, mounted on a(n) aluminum(?)spelling(?) cooling fin. I played around with the code like you said in your earlier remarks. I was able to interchange region1 to region 2 by doing so, instead of having a "circular hole in the geometry" and create the cooling fins in a similar manner, after all said and done, I compiled and ran the code.

You can model just about anything and not only can you model two dimensional you can model 3 D as well. Actually, using the student version, you can use the "EXTRUDE" to add a third dimension. It's so cool!

Depending upon your skill level, you can model the "flow between two interfaces" be it through heat conduction, ground water, chemical, etc.

The student version provides examples to show you how it works. You can even change some of the parameters of the demos and create your own problem, as I did.

Just go too: www.pdesolutions.com
I have directed a lot of people to this site, similar to the ones you are currently faced with.

Best of success, and from the sound of it, you are doing just that.
Jon_

 

[Jon_]

"I was able to interchange region1 with region 2 by so doing, instead of having an empty circular heat source, ie: a hole in my geometrical figure, I had a a circular heat source bounded by the cooling fin. You could clearly see the circular boundary bounded by the cooling fin."

Sorry, I can't communicate that well via keyboard

 

[Jon_]

The cooling fin was rectangular. The kind you would see on a cooling fin on a computer or some other appliance. By the way, I did manage to get it to work. I displayed the heat flow with a vector field and the scalar field with contours.

Some coded examples are real easy to use, like: contour( ) vector( , ). just take a look at what you can do. If anything, this can be used as a learning tool supplemental to the language you currently use. I think there is a MatLab/PDEflex partnership where you can write code in either and execute in either environment, not sure though. Check it out.

Jon_

 

[Jon_]

I also believe you can create your own scripts as well!

Jon_

 

[Jon_]

From your earlier remark:"I'm looking at flow uniformity in this case, which was why I adopted a computational approach."

would bode well with the pdesolutions.com

I do not want to overburden you, just trying to help.

Jon_