Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB - geometry M-file

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

    Code (Text):

    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!

    Code (Text):

    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
     
     
    Last edited: Jul 28, 2008
  2. jcsd
  3. Jul 29, 2008 #2
    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" statments "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.
     
  4. Jul 30, 2008 #3
    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).

    Code (Text):

    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:

    Code (Text):

    pdegplot('expg2');
    [p,e,t]=initmesh('expg2');
    pdemesh(p,e,t), axis equal
     
    which shows the following error:

    Code (Text):

    ??? 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?
     
  5. Jul 30, 2008 #4
    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_
     
    Last edited: Jul 30, 2008
  6. Jul 30, 2008 #5
    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_
     
    Last edited: Jul 30, 2008
  7. Jul 30, 2008 #6
    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 occured when I was doing array calculations and didn't have the proper dimensions. For every grid you will have
    N+1 nodes.

    Jon_
     
    Last edited: Jul 30, 2008
  8. Jul 30, 2008 #7
    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_
     
    Last edited: Jul 30, 2008
  9. Jul 30, 2008 #8
    I have a question: Is the "triangulation of boundaries" related in any way to "TRI-DIAGONAL" matrices?
     
  10. Jul 31, 2008 #9
    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

    Code (Text):

    % 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

    Code (Text):

    % 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). :frown:

    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.

    Code (Text):

    % 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.
     
    Last edited: Jul 31, 2008
  11. Aug 1, 2008 #10
    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! :cry:
     
    Last edited: Aug 1, 2008
  12. Aug 1, 2008 #11
    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_
     
  13. Aug 1, 2008 #12
    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_
     
    Last edited: Aug 1, 2008
  14. Aug 1, 2008 #13
    yes, you already have it it is your "dl= [...];"
    ok.
     
  15. Aug 3, 2008 #14
    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.
     
  16. Aug 3, 2008 #15
    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_
     
    Last edited: Aug 3, 2008
  17. Aug 3, 2008 #16
    ... "I have found that it hasn't helped me one iota".

    In place of the previous thread
     
  18. Aug 3, 2008 #17
    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_
     
    Last edited: Aug 3, 2008
  19. Aug 3, 2008 #18
    "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
     
    Last edited: Aug 3, 2008
  20. Aug 3, 2008 #19
    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_
     
  21. Aug 3, 2008 #20
    I also believe you can create your own scripts as well!!

    Jon_
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: MATLAB - geometry M-file
  1. MATLAB and M-Files (Replies: 4)

Loading...