I'm struggling to get my time dependent boundary condition function working. I went back to the beginning and copied the example straight out of the documentation for "Boundary Conditions for Scalar PDE".
There appears to be two issues with this example. First, the critical step of generating the mesh is noticeably absent. The boundary condition function, as shown, only accounts for the eight edges of the basic geometry. However, once the initmesh function is performed, the edge matrix becomes much larger (all the edges of the individual triangular elements) and the edge boundary labels in the boundary condition function (1-8) are no longer valid.
Second, and probably most important, is that the "time" variable is an empty matrix causing errors in the calculations.
Has anyone had any success with this method?
Thanks – Michael
Here is the script I'm using:
% Time Variant Boundary Conditions
clearclc% Rectangle is code 3, 4 sides,
% followed by x-coordinates and then y-coordinates
R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]';% Circle is code 1, center (.5,0), radius .2
C1 = [1,.5,0,.2]';% Pad C1 with zeros to enable concatenation with R1
C1 = [C1;zeros(length(R1)-length(C1),1)];geom = [R1,C1];% Names for the two geometric objects
ns = (char('R1','C1'))';% Set formula
sf = 'R1-C1';% Create geometry
gd = decsg(geom,sf,ns);% View geometry
figure(1)pdegplot(gd,'edgeLabels','on')xlim([-1.1 1.1])axis equal[p, e, t] = initmesh(gd,'Hmax',Inf);% [p, e, t] = refinemesh(gd,p,e,t);
% p = jigglemesh(p,e,t);
figure(2)pdemesh(p,e,t)xlim([-1.1 1.1])axis equalb = @BC_fun;u0 = 0;c = 1;a = 0;f = 10;d = 1;tlist = 0:.01:1;u = parabolic(u0,tlist,b,p,e,t,c,a,f,d);
And here is the function file:
function [ qmatrix, gmatrix, hmatrix, rmatrix ] = BC_fun( p, e, u, time )% PDE Boundary Condition Function
ne = size(e,2); % number of edges
qmatrix = zeros(1,ne);gmatrix = qmatrix;hmatrix = zeros(1,2*ne);rmatrix = hmatrix;for k = 1:ne x1 = p(1,e(1,k)); % x at first point in segment
x2 = p(1,e(2,k)); % x at second point in segment
xm = (x1 + x2)/2; % x at segment midpoint
y1 = p(2,e(1,k)); % y at first point in segment
y2 = p(2,e(2,k)); % y at second point in segment
ym = (y1 + y2)/2; % y at segment midpoint
switch e(5,k) case {1,2,3,4} % rectangle boundaries
hmatrix(k) = 1; hmatrix(k+ne) = 1; rmatrix(k) = time*(x1 - y1); rmatrix(k+ne) = time*(x2 - y2); otherwise % same as case {5,6,7,8}, circle boundaries
qmatrix(k) = 1; gmatrix(k) = xm^2 + ym^2; endend
Here is the error returned:
In an assignment A(I) = B, the number of elements in B and I must be the same.Error in BC_fun (line 21) rmatrix(k) = time*(x1 - y1);Error in pdeefxpd (line 10) [q,g,h,r]=feval(bl,p,e,u,time);Error in pdeexpd (line 40)[q,g,h,r]=pdeefxpd(p,e,u,time,bl);Error in pdeODEInfo/setN (line 178) [q,g,h,r]=pdeexpd(self.p,self.e(:,1),self.b);Error in pdeODEInfo/checkFuncDepen (line 56) self=self.setN(u0);Error in pdeParabolicInfo (line 14) obj=obj.checkFuncDepen(u0, tlist);Error in parabolic (line 41) pdeprb=pdeParabolicInfo(u0,tlist,b,p,e,t,c,a,f,d);Error in TimeVaryBC (line 46)u = parabolic(u0,tlist,b,p,e,t,c,a,f,d);
Best Answer