Skip to content
Snippets Groups Projects
Commit 9d8fbc49 authored by Carolina Bergeling's avatar Carolina Bergeling
Browse files

orginial files

parents
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CompareWithCAREOptimalCutOut2.m 0 → 100644
+ 137
− 0
View file @ 9d8fbc49
% Comparison optimal gamma
%1. CARE in matlab
%3. Closed-form method
clear all;
close all;
refine = 0;
prtlevel = 0;
deltas = [0.3 0.25 0.2 0.15 0.1 0.05];
deltas = [0.18 0.13 0.11];
for ii = 1:3
% actuator half width
ep = 0.5;
% actuator location
al = [2 2]; % center
% disturbance location
dl = [1 1.3];
epd = 0.1;
% mesh size
delta = deltas(ii); % change this to change the order of the problem
% Debug level 0 - no messages
debuglevel = 0;
% Building the approximation.
[triangles, coordinates] = mymeshdata(delta, prtlevel);
%[E, A, B,~,~] = heat_2d_original_var2(refine, coordinates, triangles, al, ep, al, delta, 1e-12, prtlevel);
%[~, ~, H,~,~] = heat_2d_original_var2(refine, coordinates, triangles, dl, epd, dl, delta, 1e-12, prtlevel);
r = al;
s = [1 1];
nu = 0.25;
d =1;
[E, A, B,C,H] = heat_2d_original_var2C(refine, coordinates, triangles, r,s,dl,ep,nu,epd, delta, 1e-12, prtlevel);
% Generalized form matrices
% Edot{z} = Az+Bu+Hv
% zeta = Cz+Du
% y = z
[n,~] = size(A);
[~ , m] = size(B);
[~,q] = size(H);
[k,~] = size(C);
L = chol(E);
A = inv(L')*A*inv(L);
B = inv(L')*B;
H = inv(L')*H;
H = eye(n,n);
[~,q] = size(H);
C = C*inv(L);
C2 = [eye(n,n); zeros(m,n)];
D = [zeros(n,m); eye(m,m)];
G = ss(A,[H zeros(n,k) B],[C2; C], [[zeros(n+m,q+k); zeros(k,q) eye(k,k)] [D;zeros(k,m)]]);
% Save system size to "order" vector
order(ii) = n;
%---------------------------------------------------------------------%
% Algorithm 3: Closed-form solution
%---------------------------------------------------------------------%
tic;
a = 1.1;
ainv = inv(a);
gc = @(a)sqrt(norm(H'*inv(A*a)^2*H,2))
gsf = @(a)sqrt(norm(H'*inv(A^2+B*B'-(1-a)^2*A^2)*H,2))
gse = @(a)sqrt(norm(H'*inv(A^2+C'*C-(1-a)^2*A^2)*H,2))
g = max(gc(a),gsf(a));
g = max(g,gse(a));
gap = min(g-gsf(1),g-gse(1));
optimality(ii,:) = [g gap];
g = 1.1*g;
g2inv = 1/g^2;
A_K = A-ainv*g2inv*inv(A)+ainv*B*B'*inv(A)+ainv*inv(eye(n)-1/(g*a)^2*inv(A)^2)*inv(A)*C'*C;
B_K = -ainv*inv(eye(n)-1/(g*a)^2*inv(A)^2)*inv(A)*C';
C_K = ainv*B'*inv(A);
D_K = zeros(m,k);
K = ss(A_K,B_K,C_K,D_K);
timeCF(ii) = toc;
%cl = lft(G,K);
%clall(ii) = cl;
%hinfnorm(pck(cl.a,cl.b,cl.c,cl.d))
%%---------------------------------------------------------------------%
% Algorithm 1: Schur Method (through MATLABs CARE)
%---------------------------------------------------------------------%
% Tolerance in the optimal gamma.
gamtol = 1.1*gap;
% Starting values for largest and smallest value of gamma.
% How to choose these?
gmax = 10; % H2 solution?
gmin = 0; % 0?
tic;
[K,closed,gfin] = hinfsyn(pck(G.a,G.b,G.c,G.d),k,m,gmin,gmax,gamtol)
timeSchur(ii) = toc;
end;
%%
figure
semilogy(order,timeSchur,'ro')
hold on
semilogy(order,timeCF,'ko')
xlabel('problem order')
ylabel('comp. time')
legend('Schur','Closed-form','Location', 'SouthEast')
heat_2d_original_var2C.m 0 → 100644
+ 371
− 0
View file @ 9d8fbc49
function [M,N,B,C,D,status,triangles,coordinates] = heat_2d_original_var2C(refine, coordinates, triangles, r,s,dl,ep,nu,epd, delta,mytol,prtlevel)
%% [M,N] = myheat(r,s,ep,nu,delta)
% Written by Dhanaraja Kasinathan and Kirsten A Morris, 2014. Altered by
% Carolina Bergeling 2019.
% The following (MATLAB) files are mandatory:
% distmesh(directory): Mesh generator could be downloaded
% from http://www-math.mit.edu/~persson/mesh/
% Path of this directory needs
% to be added using addpath command
% distmesh2d.m: Generates the mesh. (present inside
% distmesh(directory))
% fixmesh.m: Fix for removing the redundant coordinates. (present inside
% distmesh(directory)).
% stima3.m: Computes the element stiffness matrix for each
% triangle.
%mycharacteristic.m: Performs double integration (quad twice) over the triangle that
% partially overlaps the square (actuator/sensor).
%etaj.m: returns the function that needs to be integrated
% over the square patch. This is used only when the
% actuator/sensor is completely sitting inside one
% triangle.
%overlap.m: Finds the triangles that intersect with the given
% square
%myverifyb.m: Calculates the error between the [b]_j and the
% b(x,y) on each triangle.
%show.m: Plots the given value on the triangular mesh.(optional)
%
%FEM2D_HEAT sets up the state-space matrices for the PDE described
%below
% d/dt z = div(grad(z)) + b(x,y) f in Omega
% z = 0 on the Dirichlet boundary
% z(x,y,0) = z_0 in Omega
% y = int_omega c(x,y)z d(x,y) (Integrating u over omega)
%on the geometry as described below, where ,with XI denoting the indicator function
%b(x,y) = 1/(2*ep) * XI_B(r,ep)(x,y) ;
%c(x,y) = 1/(2*nu) * XI_B(s,nu)(x,y) ;
%where B(r,ep) is the square with each side of width 2*ep centered at r
%
% The input arguments are the following
% r: 1 X 2 row vector,the 2D coordinates of the position of actuator r =
% [rx,ry]. Must be within the Geometry
% s: 1 X 2 row vector,the 2D coordinates of the position of sensor s =
% [sx,sy]. Must be within the Geometry
% ep: half width of the actuator( a square patch with each side of width '2*ep'
% centered at 'r'.)
% nu: half width of the sensor ( a square patch with each side of width '2*nu'
% centered at 's'.)
% delta: positive scalar, mesh size.
% prtlevel (default 0, no display messages), 1-debug messages.
%
% The output arguments are the corresponding system matrices
% such that the system above could be represented in the state-space form
% as follows
% M d/dt z = K z + F u
% y = C z
% M : Mass matrix.
% N : Stiffness matrix.
% F : Load matrix.
% C: Output operator (= C(sx,sy)) = {1/(2*nu), |x-sx| < nu & |y-sy| < nu; 0, otherwise}
% Cfull: Output operator with penaly on all states. (C = I, that is y=z(x,y,t))
% Using this representation, the standard system matrices can be derived
% as follows,
% A = M\ N: 2D laplacian operator
% B = M\ F: Input operator (= B(rx,ry)) = {1/(2*ep), |x-rx| < ep & |y-ry| < ep; 0, otherwise}
%
% Example (usage):
% [M,N,F,C,Cfull] = myheat([2,2],[2,2],0.2,0.2,0.1)
%
%
% Geometry :
% A Square of side 4, with main diagonal from (0,0) to (4,4).
% A circle of radius of 0.4 centered at (3,1) is cut out.
% A mesh is generated with the width of each edge being "delta" given by
% the user.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Book keeping : To remove warnings from delayunayn and quad.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s1 = warning('off','MATLAB:delaunayn:DuplicateDataPoints');
s2 = warning('off','MATLAB:quad:MinStepSize');
s3 = warning('off','MATLAB:quad:ImproperFcnValue');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setting the accuracy of the integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mytol=1e-12;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 'distmesh' is used to generate the required geometry
% % Geometry : A square with a hole
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % Square coordinates and circle's centre and radii needs to be changed here
% fd=inline('ddiff(drectangle(p,0,4,0,4),dcircle(p,3,1,0.4))','p');
% box=[0,0;4,4];
% fix = [0,0;0,4;4,0;4,4];
% if (prtlevel > 0) disp('Start the mesh generator...'); end;
% tic;
% [coordinates,triangles]=distmesh2d(fd,@huniform,delta,box,fix);
% tmesh = toc;
% [coordinates,triangles] = fixmesh(coordinates,triangles); % Need this to remove the redundant coordinate.
% if (prtlevel > 0) str = sprintf('myheat_refined ::Finished Meshing in %1.2e refining ...',tmesh); disp(str); end;
%
% % Retrieving the coordinates of the actuator and sensor
% rx=r(1,1);ry=r(1,2);sx=s(1,1);sy=s(1,2);
% % Retrieving the endpoints of the square patch for actuator and sensor.
% actleft = rx - ep; actright = rx+ ep; actbottom = ry - ep; acttop = ry + ep;
% senleft = sx - nu; senright = sx+ nu; senbottom = sy - nu; sentop = sy + nu;
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Plotting the actuator and sensor on top of the generated mesh.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% if (prtlevel > 0)
% hold on, text('Position',r,'String','A','Fontsize',12,'FontWeight','bold'), text('Position',s,'String','S','Fontsize',12,'FontWeight','bold'),
% plot([rx-ep rx-ep rx+ep rx+ep rx-ep],[ry-ep ry+ep ry+ep ry-ep ry-ep],'b','linewidth',2),
% plot([sx-nu sx-nu sx+nu sx+nu sx-nu],[sy-nu sy+nu sy+nu sy-nu sy-nu],'m','linewidth',2), hold off
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Removing the roundoff-error (zero-correction) generated by the mesh on the grid points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% for i=1:size(coordinates,1)
% if (abs(coordinates(i,1)) < 1e-8)
% coordinates(i,1) = 0;
% end
% if (abs(coordinates(i,2)) < 1e-8)
% coordinates(i,2) = 0;
% end
% if (abs(coordinates(i,1)-4) < 1e-8)
% coordinates(i,1) = 4;
% end
% if (abs(coordinates(i,2)-4) < 1e-8)
% coordinates(i,2) = 4;
% end
% end
tic;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refining the mesh on all triangles that lie in and around the actuator.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (refine == 1)
% Meshing thrice with delta/4 had a better value of error as delta goes
% to zero.
for i=1:2
tri = overlaptriwin(coordinates,triangles,r,ep+delta/4); % Finds the triangles that overlaps the square actuator.
[coordinates,triangles] = myrefine(coordinates,triangles,tri); % Refines those triangles.
end;
%for i=1:2
% tri = overlaptriwin(coordinates,triangles,d,ep+delta/4); % Finds the triangles that overlaps the square actuator.
% [coordinates,triangles] = myrefine(coordinates,triangles,tri); % Refines those triangles.
%end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Finding the boundary of the geometry described above. This includes the
%edges along the cut-out circle as well. Dirichlet condition needs to be
%applied on this boundary.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dirichlet = boundedges(coordinates,triangles);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Error-Check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (any(size(r) ~= [1,2])) || (any(r) > 4) || (any(r) < 0)
error('Invalid input arguments. see help fem2d_heat');
end
if ((r(1,1)-3)^2+(r(1,2)-1)^2 < 0.4+ep)
%warning('myheatrefineddata :: Actuator location is not inside the geometry.');
status = 0; M = []; N = []; B = []; D = []; C = [];
return;
else
status = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finding the triangles that partially overlap('actpartial') or completely
% inside ('acttriinsquare') the square. If the mesh size is coarse, then the
% last parameter('actsquareintri') provides the triangle in which the actuator/sensor is present.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Actuator
%[actpartial,acttriinsquare,actsquareintri] = overlap(coordinates,triangles,r,ep);
% Sensor
%[senpartial,sentriinsquare,sensquareintri] = overlap(coordinates,triangles,s,nu);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialising the global stiffness matrix K, mass matrix M, force matrix F,
% Output matrix C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = zeros(size(coordinates,1));
M = zeros(size(coordinates,1));
B = zeros(size(coordinates,1),1);
D = zeros(size(coordinates,1),1);
C = zeros(1,size(coordinates,1));
if (prtlevel > 0) disp('myheatrefineddata :: Finished refining. Building matrices ...'); end;
tic;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assembling the global stiffness matrix K with variable diffusion (stima3var routine does this) and
% mass matrix M (this is a direct formula)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 1:size(triangles,1)
N(triangles(j,:),triangles(j,:)) = N(triangles(j,:),triangles(j,:)) ...
+ stima3var_different(coordinates(triangles(j,:),:)); % to account for dx^2dy^2
M(triangles(j,:),triangles(j,:)) = M(triangles(j,:),triangles(j,:)) ...
+ det([1,1,1;coordinates(triangles(j,:),:)'])*[2,1,1;1,2,1;1,1,2]/24;
valb = mycharacteristic(coordinates(triangles(j,:),:),r,ep,mytol);
B(triangles(j,1),1) = B(triangles(j,1),1) + valb(1);
B(triangles(j,2),1) = B(triangles(j,2),1) + valb(2);
B(triangles(j,3),1) = B(triangles(j,3),1) + valb(3);
vald = mycharacteristic(coordinates(triangles(j,:),:),dl,epd,mytol);
D(triangles(j,1),1) = D(triangles(j,1),1) + vald(1);
D(triangles(j,2),1) = D(triangles(j,2),1) + vald(2);
D(triangles(j,3),1) = D(triangles(j,3),1) + vald(3);
val = mycharacteristic(coordinates(triangles(j,:),:),s,nu,mytol);
C(1,triangles(j,1)) = C(1,triangles(j,1)) + val(1);
C(1,triangles(j,2)) = C(1,triangles(j,2)) + val(2);
C(1,triangles(j,3)) = C(1,triangles(j,3)) + val(3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assembling the force matrix on triangles that lie partially inside the
% actuator and partially outside the actuator. This is done by the function
% 'mycharacteristic'. For every triangle that lie partially on the boundary
% I am calling this function with corresponding vertices, center and the half-width of the actuator.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for t=1:size(actpartial,1)
% val = mycharacteristic(coordinates(triangles(actpartial(t,1),:),:),r,ep,mytol);
% F(triangles(actpartial(t,1),1),1) = F(triangles(actpartial(t,1),1),1) + val(1);
% F(triangles(actpartial(t,1),2),1) = F(triangles(actpartial(t,1),2),1) + val(2);
% F(triangles(actpartial(t,1),3),1) = F(triangles(actpartial(t,1),3),1) + val(3);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assembling the force matrix on triangles that lie completely inside the
% actuator: The value for this is a direct formula = 1/(2*ep) * Je * (1/6).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for t=1:size(acttriinsquare,1)
% Je = abs(det([1 1 1;coordinates(triangles(acttriinsquare(t,1),:),:)']))/(12*ep);
% F(triangles(acttriinsquare(t,1),1),1) = F(triangles(acttriinsquare(t,1),1),1) + Je;
% F(triangles(acttriinsquare(t,1),2),1) = F(triangles(acttriinsquare(t,1),2),1) + Je;
% F(triangles(acttriinsquare(t,1),3),1) = F(triangles(acttriinsquare(t,1),3),1) + Je;
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If the actuator is very small when compared to the mesh size then it may
% lie completely inside one triangle. Assmebling the force matrix for that
% triangle. Limits of integration are the end points of the square.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if (size(actsquareintri,1) > 0)
% F(triangles(actsquareintri,1),1) = F(triangles(actsquareintri,1),1) + dblquad('etaj',actbottom,acttop,actleft,actright,1e-12,[],coordinates(triangles(actsquareintri,:),:),ep,1);
% F(triangles(actsquareintri,2),1) = F(triangles(actsquareintri,2),1) + dblquad('etaj',actbottom,acttop,actleft,actright,1e-12,[],coordinates(triangles(actsquareintri,:),:),ep,2);
% F(triangles(actsquareintri,3),1) = F(triangles(actsquareintri,3),1) + dblquad('etaj',actbottom,acttop,actleft,actright,1e-12,[],coordinates(triangles(actsquareintri,:),:),ep,3);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assembling the output matrix on triangles that lie partially inside the
% sensor and partially outside the sensor. This is done by the function
% 'mycharacteristic'. For every triangle that lie partially on the boundary
% I am calling this function with corresponding vertices, center
% and half-width of the sensor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for t=1:size(senpartial,1)
% val = mycharacteristic(coordinates(triangles(senpartial(t,1),:),:),s,nu,mytol);
% C(1,triangles(senpartial(t,1),1)) = C(1,triangles(senpartial(t,1),1)) + val(1);
% C(1,triangles(senpartial(t,1),2)) = C(1,triangles(senpartial(t,1),2)) + val(2);
% C(1,triangles(senpartial(t,1),3)) = C(1,triangles(senpartial(t,1),3)) + val(3);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assembling the force matrix on triangles that lie completely inside the
% sensor. The value for this is a direct formula = 1/(2*nu) * Je * (1/6).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for t=1:size(sentriinsquare,1)
% Je = abs(det([1 1 1;coordinates(triangles(sentriinsquare(t,1),:),:)']))/(12*nu);
% C(1,triangles(sentriinsquare(t,1),1)) = C(1,triangles(sentriinsquare(t,1),1)) + Je;
% C(1,triangles(sentriinsquare(t,1),2)) = C(1,triangles(sentriinsquare(t,1),2)) + Je;
% C(1,triangles(sentriinsquare(t,1),3)) = C(1,triangles(sentriinsquare(t,1),3)) + Je;
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If the sensor is very small when compared to the mesh size then it may
% lie completely inside only one triangle. Assmebling the force matrix for that
% triangle. Limits of integration are the end points of the square.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if (size(sensquareintri,1) > 0)
% C(1,triangles(sensquareintri,1)) = C(1,triangles(sensquareintri,1)) + dblquad('etaj',senbottom,sentop,senleft,senright,1e-12,[],coordinates(triangles(sensquareintri,:),:),nu,1);
% C(1,triangles(sensquareintri,2)) = C(1,triangles(sensquareintri,2)) + dblquad('etaj',senbottom,sentop,senleft,senright,1e-12,[],coordinates(triangles(sensquareintri,:),:),nu,2);
% C(1,triangles(sensquareintri,3)) = C(1,triangles(sensquareintri,3)) + dblquad('etaj',senbottom,sentop,senleft,senright,1e-12,[],coordinates(triangles(sensquareintri,:),:),nu,3);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Error-Check: If all entries of F( or C) are zeroes, then the actuator( or
%sensor) is not within the geometry.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (any(B) == 0)
warning('myheatdata :: Actuator location is not inside the geometry.');
status = 0;
else
status = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Incorporating the Dirichlet boundary condition in the global stiffness
%matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N(unique(dirichlet),:) = zeros(size(unique(dirichlet),1),size(coordinates,1));
N(:,unique(dirichlet)) = zeros(size(coordinates,1),size(unique(dirichlet),1));
N(unique(dirichlet),unique(dirichlet)) = eye(size(unique(dirichlet),1));
N=-N; %Negative sign appears while multiplying the laplacian operator with test functions and integrating by parts.
tbuild = toc;
if (prtlevel > 0) str = sprintf('myheatdata :: ..Done in %1.2e. System order-%d.',tbuild,size(N,2));disp(str); end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computing A = M\K , B = M\F
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A = M\N;
% Btilde = M\B;
% Dtilde = M\D;
% d = sort(eig(A));
% figure, plot(real(d),imag(d),'*');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plotting theoretical value of B on the triangles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure,
% show(triangles,[],coordinates,Btilde);
% figure,
% show(triangles,[],coordinates,Dtilde);
%disp('Computing error...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computing the error in B :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% acterror = 0;
% for i=1:size(triangles,1)
% acterror = acterror + myverifyb(coordinates(triangles(i,:),:),B(triangles(i,:),:),r,ep,mytol);
% end
% acterror
end
mymeshdata.m 0 → 100644
+ 148
− 0
View file @ 9d8fbc49
function [triangles,coordinates] = mymeshdata(delta,prtlevel)
%% [M,N] = mymeshdata(delta,prtlevel)
% Written by Dhanaraja Kasinathan and Kirsten A Morris, 2014. Altered by
% Carolina Bergeling 2019.
% Generates the mesh for the base mesh size.
% Populates the coordinates of the points and the nodes in each element.
% The following (MATLAB) files are mandatory:
% distmesh(directory): Mesh generator could be downloaded
% from http://www-math.mit.edu/~persson/mesh/
% Path of this directory needs
% to be added using addpath command
% distmesh2d.m: Generates the mesh. (present inside
% distmesh(directory))
% fixmesh.m: Fix for removing the redundant coordinates. (present inside
% distmesh(directory)).
% stima3.m: Computes the element stiffness matrix for each
% triangle.
%mycharacteristic.m: Performs double integration (quad twice) over the triangle that
% partially overlaps the square (actuator/sensor).
%etaj.m: returns the function that needs to be integrated
% over the square patch. This is used only when the
% actuator/sensor is completely sitting inside one
% triangle.
%overlap.m: Finds the triangles that intersect with the given
% square
%myverifyb.m: Calculates the error between the [b]_j and the
% b(x,y) on each triangle.
%show.m: Plots the given value on the triangular mesh.(optional)
%
%FEM2D_HEAT sets up the state-space matrices for the PDE described
%below
% d/dt z = div(grad(z)) + b(x,y) f in Omega
% z = 0 on the Dirichlet boundary
% z(x,y,0) = z_0 in Omega
% y = int_omega c(x,y)z d(x,y) (Integrating u over omega)
%on the geometry as described below, where ,with XI denoting the indicator function
%b(x,y) = 1/(2*ep) * XI_B(r,ep)(x,y) ;
%c(x,y) = 1/(2*nu) * XI_B(s,nu)(x,y) ;
%where B(r,ep) is the square with each side of width 2*ep centered at r
%
% The input arguments are the following
% r: 1 X 2 row vector,the 2D coordinates of the position of actuator r =
% [rx,ry]. Must be within the Geometry
% s: 1 X 2 row vector,the 2D coordinates of the position of sensor s =
% [sx,sy]. Must be within the Geometry
% ep: half width of the actuator( a square patch with each side of width '2*ep'
% centered at 'r'.)
% nu: half width of the sensor ( a square patch with each side of width '2*nu'
% centered at 's'.)
% delta: positive scalar, mesh size.
% prtlevel (default 0, no display messages), 1-debug messages.
%
% The output arguments are the corresponding system matrices
% such that the system above could be represented in the state-space form
% as follows
% M d/dt z = K z + F u
% y = C z
% M : Mass matrix.
% N : Stiffness matrix.
% F : Load matrix.
% C: Output operator (= C(sx,sy)) = {1/(2*nu), |x-sx| < nu & |y-sy| < nu; 0, otherwise}
% Cfull: Output operator with penaly on all states. (C = I, that is y=z(x,y,t))
% Using this representation, the standard system matrices can be derived
% as follows,
% A = M\ N: 2D laplacian operator
% B = M\ F: Input operator (= B(rx,ry)) = {1/(2*ep), |x-rx| < ep & |y-ry| < ep; 0, otherwise}
%
% Example (usage):
% [M,N,F,C,Cfull] = myheat([2,2],[2,2],0.2,0.2,0.1)
%
%
% Geometry :
% A Square of side 4, with main diagonal from (0,0) to (4,4).
% A circle of radius of 0.4 centered at (3,1) is cut out.
% A mesh is generated with the width of each edge being "delta" given by
% the user.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Book keeping : To remove warnings from delayunayn and quad.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s1 = warning('off','MATLAB:delaunayn:DuplicateDataPoints');
s2 = warning('off','MATLAB:quad:MinStepSize');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setting the accuracy of the integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mytol=1e-12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'distmesh' is used to generate the required geometry
% Geometry : A square with a hole
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Square coordinates and circle's centre and radii needs to be changed here
fd=inline('ddiff(drectangle(p,0,4,0,4),dcircle(p,3,1,0.4))','p');
box=[0,0;4,4];
fix = [0,0;0,4;4,0;4,4];
if (prtlevel > 0) disp('mymeshdata :: Start the mesh generator...'); end;
tic;
[coordinates,triangles]=distmesh2d(fd,@huniform,delta,box,fix);
tmesh = toc;
[coordinates,triangles] = fixmesh(coordinates,triangles); % Need this to remove the redundant coordinate.
if (prtlevel > 0) str = sprintf('mymeshdata ::Finished Meshing in %1.2e',tmesh); disp(str); end;
%
% % Retrieving the coordinates of the actuator and sensor
% rx=r(1,1);ry=r(1,2);sx=s(1,1);sy=s(1,2);
% % Retrieving the endpoints of the square patch for actuator and sensor.
% actleft = rx - ep; actright = rx+ ep; actbottom = ry - ep; acttop = ry + ep;
% senleft = sx - nu; senright = sx+ nu; senbottom = sy - nu; sentop = sy + nu;
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Plotting the actuator and sensor on top of the generated mesh.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% if (prtlevel > 0)
% hold on, text('Position',r,'String','A','Fontsize',12,'FontWeight','bold'), text('Position',s,'String','S','Fontsize',12,'FontWeight','bold'),
% plot([rx-ep rx-ep rx+ep rx+ep rx-ep],[ry-ep ry+ep ry+ep ry-ep ry-ep],'b','linewidth',2),
% plot([sx-nu sx-nu sx+nu sx+nu sx-nu],[sy-nu sy+nu sy+nu sy-nu sy-nu],'m','linewidth',2), hold off
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Removing the roundoff-error (zero-correction) generated by the mesh on the grid points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:size(coordinates,1)
if (abs(coordinates(i,1)) < 1e-8)
coordinates(i,1) = 0;
end
if (abs(coordinates(i,2)) < 1e-8)
coordinates(i,2) = 0;
end
if (abs(coordinates(i,1)-4) < 1e-8)
coordinates(i,1) = 4;
end
if (abs(coordinates(i,2)-4) < 1e-8)
coordinates(i,2) = 4;
end
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment