Need help with this MATLAB HW
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/589887/image.png)
function Example641()
%
%
% example for 2D steady-state diffusion equation on \Omega=(0,1)x(0,1)
%
% Neumann BC on all boundaries
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (~isdeployed)
addpath('MMPDElab/src_MMPDElab');
end
% set the basic parameters
jmax = 41;
a = 0;
b = 1;
npde = 1;
moving_mesh = false;
mmpde_tau = 1e-1;
mmpde_ncycles = 3;
mmpde_alpha = [];
nn = 10;
% set the initial meshes, find the indices of the corner points and fix them
kmax = jmax;
[X,tri] = MovMesh_rect2tri(linspace(a,b,jmax), linspace(a,b,kmax), 1);
TR = triangulation(tri,X);
tri_bf = freeBoundary(TR);
Nbf = length(tri_bf);
[Nv,d] = size(X);
N = size(tri, 1);
Xi_ref = X;
% find the indices of the corner points and fix them
corners = [a, a; b, a; b, b; a, b];
[~,nodes_fixed] = ismembertol(corners,Xi_ref,1e-10,'ByRows',true);
% nodes_fixed = unique(tri_bf); % for fixing all boundary nodes
% set initial conditions
% set the initial solution
U = zeros(Nv,npde);
figure(1)
clf
triplot(tri,X(:,1),X(:,2),'Color','r')
axis([a b a b]);
axis square;
% define PDE system and BCs
pdedef.bfMark = ones(Nbf,1); % for y = 0 (b1)
% Xbfm = (X(tri_bf(:,1),:)+X(tri_bf(:,2),:))*0.5;
% pdedef.bfMark(Xbfm(:,1)<1e-8) = 4; % for x = 0 (b4)
% pdedef.bfMark(Xbfm(:,1)>1-1e-8) = 2; % for x = 1 (b2)
% pdedef.bfMark(Xbfm(:,2)>1-1e-8) = 3; % for y = 1 (b3)
% define boundary types
pdedef.bftype = zeros(Nbf,npde);
% pdedef.bftype = ones(Nbf,npde);
% % for neumann bcs:
% pdedef.bftype(pdedef.bfMark==2|pdedef.bfMark==3,npde) = 0;
pdedef.volumeInt = @pdedef_volumeInt;
pdedef.boundaryInt = @pdedef_boundaryInt;
pdedef.dirichletRes = @pdedef_dirichletRes;
% perform integration (MP)
tcpu = cputime;
if (~moving_mesh)
nn = 1;
end
for n=1:nn
fprintf('--- n = %d\n', n);
% move the mesh
if (moving_mesh)
M = MovMesh_metric(U,X,tri,tri_bf,mmpde_alpha);
M = MovMesh_metric_smoothing(M,mmpde_ncycles,X,tri);
Xnew = MovMesh([0,1],Xi_ref,X,M,mmpde_tau,tri,tri_bf,nodes_fixed);
else
Xnew = X;
end
% solve physical PDEs
Unew = MovFEM_bvp(U,Xnew,tri,tri_bf,pdedef,'newtons');
%Unew = MovFEM_bvp(U,Xnew,tri,tri_bf,pdedef,'fsolve');
% update
X = Xnew;
U = Unew;
% figure(2)
% clf
% triplot(tri,X(:,1),X(:,2),'Color','r')
% axis([a b a b]);
% axis square;
% drawnow;
end
tcpu = cputime-tcpu;
fprintf('\n --- total elapsed cpu time = %e \n\n', tcpu);
% output
figure(3)
clf
trisurf(tri,X(:,1),X(:,2),U(:,1))
xlabel('x')
ylabel('y')
% err = MovFEM_Error_P1L2(@uexact,0,X,U,tri,tri_bf);
% Ue = uexact(0, X);
% fprintf('(Nv, N, max err, L2 err) = %d %d %e %e\n',Nv,N,norm(Ue-U,Inf),err);
% [Qgeo,Qeq,Qali] = MovMesh_MeshQualMeasure(X,tri,M);
% fprintf(' Mesh quality measures (Qgeo, Qeq, Qali) = %e %e %e\n', ...
% Qgeo, Qeq, Qali);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u = uexact(t, x)
u = sin(2*pi*x(:,1)).*sin(3*pi*x(:,2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = pdedef_volumeInt(du, u, ut, dv, v, x, t, ipde)
k1 = 1.0;
k2 = 0.1;
F = sin(2*pi*x(:,1)).*sin(3*pi*x(:,2));
F = k1*du(:,1).*dv(:,1) + k2*du(:,2).*dv(:,2) - F.*v(:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G = pdedef_boundaryInt(du, u, v, x, t, ipde, bfMark)
G = zeros(size(x,1),1);
% ID = find(bfMark==2);
% G(ID) = -2*pi*cos(2*pi*x(ID,1)).*sin(3*pi*x(ID,2)).*v(ID);
% ID = find(bfMark==3);
% G(ID) = -3*pi*sin(2*pi*x(ID,1)).*cos(3*pi*x(ID,2)).*v(ID);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Res = pdedef_dirichletRes(u, x, t, ipde, bfMark)
Res = u(:,1);
% Res = zeros(size(x,1),1);
% ID = find(bfMark==1|bfMark==4);
% Res(ID) = u(ID,1) - uexact(t,x(ID,:));
end
1 comentario
Respuestas (0)
Ver también
Categorías
Más información sobre Boundary Conditions en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!