Need help with this MATLAB HW

Question

Dhinesh Nagarajan 2021-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/807367-need-help-with-this-matlab-hw

评论： Walter Roberson 2021-4-20

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