[Up]
3.5.1999

Modifications caused by the boundaries , assemb

Follow the way PDEtool does the assembly. (Manual Ch. 4)

Suggestion for an outline of lassemb

We will explain the meaning of the boundary-b-matrix. It looks a bit mysterious in the beginning. Of course you don't have to follow the design of PDEtool , you can perfectly well design your own way, for instance just writing an m-file for each boundary segment for defining your boundary conditions.

On the other hand we try to give simplified instructions of using the b-matrix.

A third point of view is that the data structures in new versions of Matlab allow for easier and more natural design (than the toolbox b-matrix). Femlab does it the new way.

The b-matrix

In our example:
b =

     1     1     1     1
     1     1     0     1
     1     1     1     1
     1     1     1     1
     1     1    48     1
     3     3    49     1
    48    48    48    48
    48    48    48    48
    49    49    49    49
    49    49    48    49
    43    43     0     0
   120   121     0     0

>> char(b)

ans =

````
`` `
````
````
``0`
bb1`
0000
0000
1111
1101
++  
xy  

Meaning of the rows in b, consider only scalar case (N=1)

(see Manual p. 5-14). (Requires some "reading between lines")
  • Row 1: 1 in the scalar case (the dim. of the system in general)
  • Row 2: 1 for Dirichlet, 0 for Neumann
  • Row 3: Lengths of strings representring q (Neumann)
  • Row 4: Lengths of strins repr. g (Neumann)
  • Row 5: Lengths of strins repr. h (Dirichlet)
  • Row 6: Lengths of strins repr. r (Dirichlet)
  • Next rows: The stings representing the boundary cond. The above has to be interpreted so that in case of Dirichlet boundaries the string for r starts at row 10 and in case of Neumann, the string for q starts at row 5.
    
    
    

    Dirichlet-boundaries

    The form is

    h u = r

    On row 2 we see that columns 1,2,4 represent Dirichlet boundaries. We can agree to take h=1 always, thus we need to pick r which is obtained by

    for j=[1,2,4]
      char(b(10:10+b(6,1)-1,j)) 
    end;
    ans =
    
    1
    +
    x
    
    
    ans =
    
    1
    +
    y
    
    
    ans =
    
    1
    
    

    (Generalized) Neumann boundaries

    The form is

    n.(c grad(u)) + q u = g

    On row 2 of b we see that the 3rd colmumn represents a Neumann boundary.

    
    j=3;
    lq=b(3,j)  % length of string reprsenting q
    lg=b(4,j)  % length of string reprsenting g
    q=char(b(5:5+lq-1,j))  % The string for q starts at row 5 (a bit hard to
                           % read from the general explanation in the manual).
    g=char(b(5+lq:5+lq+lg-1,j))
    
    >> lq=b(3,j)
    
    lq =
    
         1
    
    >> lg=b(4,j)
    
    lg =
    
         1
    
    >> q=char(b(5:5+lq-1,j))
    
    q =
    
    0
    
    >> g=char(b(5+lq:5+lq+lg-1,j))
    
    g =
    
    1
    
    Note that c is a coefficient of the PDE and has been directly exported into the matlab workspace.

    A suggestion for the structure of lassemb

    function [G,R]=lassemb(b,p,e)  % 1 st. version, forget Q
    %function [Q,G,R]=lassemb(b,p,e)
    %QASSEMB Assembles boundary condition contributions in a PDE problem.
    %
    %	[Q,G,R]=LASSEMB(B,P,E) assembles the matrix Q and the vectors
    %	G and R.
    %	Q should be added to the system matrix and contains contributions
    %	from Mixed boundary conditions.
    %	G should be added to the right-hand side and contains contributions
    %	from Neumann and Mixed boundary conditions.
    %	R represents the Dirichlet type boundary conditions; first column
    %	contains indices of Dirichlet nodes and second column contains the
    %	values of the solution.
    %
    % 	The input parameter B is a Boundary Condition matrix exported from
    %	pdetool. 
    %
    np=size(p,2); % Number of points
    % Q=sparse(np,np);  % later
    G=sparse(np,1);
    R=sparse(np,1);
    Ri=sparse(np,1);
    
    ie=pdesde(e);  % Indices of outer boundary edges see help pdesde
    % Nothing to do?
    if length(ie)==0
      % Done
      return;
    end
    
    % Get rid of unwanted edges
    e=e(:,ie);
    ne=size(e,2);  % Number of edges
    
    % Diriclet conditions
    % Try this first:
    id=find(b(2,:)==1);  % Dirichlet cond. on these boundary segments
    for i=id
      str=char(b(10:10+b(6,i)-1,i))
    end;
    %
    % Embed this in something like:
    
    id=find(b(2,:)==1);  % Dirichlet cond. on these boundary segments
    for i=id
      
       for ...
      str=char(b(10:10+b(6,i)-1,i))  % See this, helps in using b-matrix
    
        x= ...
        y= ...
        r=eval(str);  % Remember the treatment of f above
     
       end
    end;
    
      
    % Generalized Neumann conditions:
    % Try this first:
    in=find(b(2,:)==0); % Neumann cond. on these boundary segments
    for j=in
      lq=b(3,j);
      lg=b(4,j);
      q=char(b(5:5+lq-1,j))
      g=char(b(5+lq:5+lq+lg-1,j))
    end
    
    
    % Then imbed it in
    for i=in
      for ...
    
    
      end;
    end;