[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;