Add files via upload

Functions/FPE.prj

@@ -0,0 +1,2 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<MATLABProject xmlns="http://www.mathworks.com/MATLABProjectFile" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0"/>

Functions/Force_diff_xy.m

@@ -0,0 +1,13 @@
+function [f1,f2] = Force_diff_xy(x1,x2)
+
+%分化图景的动力对坐标矩阵的函数
+
+a =1;
+b=1;
+S=0.5;
+n=4;
+k=1;
+Sn = S^n;
+f1 = a*x1.^n./(Sn+x1.^n) + b*S^n./(Sn+x2.^n)-k*x1;
+f2 = a*x2.^n./(Sn+x2.^n) + b*S^n./(Sn+x1.^n)-k*x2;
+end

Functions/Force_ring_xy.m

@@ -0,0 +1,11 @@
+function [f1,f2] = Force_ring_xy(X,Y)
+%环流动力-坐标juzhen函数  X Y分别为网格的横纵坐标矩阵
+a = 0.1;
+b = 0.1;
+c = 100;
+eps = 0.1;
+tau = 5;
+f1 = (eps^2+X.^2)./((1+X.^2).*(1+Y))-a*X;
+f2 = (b-Y./(1+c*X.^2))/tau;
+end
+

Functions/deri_x1.m

@@ -0,0 +1,48 @@
+function [M] = deri_x1(L,Nx,D,bc,f)
+%Force = ring
+x = linspace(0,L,Nx);
+dx = x(2)-x(1);
+
+%DERI_X x方向一阶导数
+N = Nx^2;
+M = zeros(N,N);
+
+%内点
+for i= 2:Nx-1
+    for j = 2:Nx-1
+        n = Nx*(i-1);
+        M(n+j,n+j+1)=1/2;
+        M(n+j,n+j-1)=-1/2;
+    end
+end
+M = M/dx;
+
+
+%left and right
+for i = 1:Nx
+    [f1,f2]=Force_ring_xy(0,(i-1)*dx);
+    M((i-1)*Nx+1,(i-1)*Nx+1) = f1/D;
+    [f1,f2]=Force_ring_xy(L,(i-1)*dx);
+    M(i*Nx, i*Nx) = f1/D;
+end
+
+if bc == 1
+   n=Nx*(Nx-1);
+   for j=2:Nx-1
+       [f1,f2]=Force_ring_xy((j-1)*dx,0);
+       M(j,j)=f1/D;
+       [f1,f2]=Force_ring_xy((j-1)*dx,L);
+       M(n+j,n+j)=f1/D;
+   end
+else
+    for j=2:Nx-1
+       M(j,j-1)=-1/(2*dx);
+       M(j,j+1)=1/(2*dx);
+       M(n+j,n+j-1)=-1/(2*dx);
+       M(n+j,n+j+1)=1/(2*dx);
+   end
+
+
+
+end
+end

Functions/grid-laplace.asv

@@ -0,0 +1,32 @@
+function [A] = grid-laplace(L,Nx)
+%GRID-LAPLACE 均匀正方网格拉普拉斯算子
+%Grid
+x=linspace(0,L,Nx);
+dx = x(2)- x(1);
+
+%initialize Matrices
+ N = Nx^2;
+ M = zeros(N,N)
+
+
+ %内点
+ for i =2:Nx-1
+     for j = 2:Nx-1
+         n=i+(j-1)*Nx;
+         M(n,n) = -4;
+         M(n,n-1) = 1;
+         M(n,n+1) = 1;
+         M(n,n-Nx) = 1;
+         M(n,n+Nx) = 1;
+     end
+ end
+
+ %up
+ for j = 2:Nx-1
+
+
+
+
+
+end
+

Functions/grid_laplace.m

@@ -0,0 +1,126 @@
+function [M] = grid_laplace(L,Nx,BC,D)
+%GRID-LAPLACE 均匀正方网格拉普拉斯算子
+x = linspace(0,L,Nx);
+dx = x(2)-x(1);
+
+
+%initialize Matrices
+ N = Nx^2;
+ M = zeros(N,N);
+
+
+ %内点
+ for i =2:Nx-1
+     for j = 2:Nx-1
+         n=i+(j-1)*Nx;
+         M(n,n) = -4;
+         M(n,n-1) = 1;
+         M(n,n+1) = 1;
+         M(n,n-Nx) = 1;
+         M(n,n+Nx) = 1;
+     end
+ end
+
+
+
+
+ %down and up 内点
+ n = (Nx-1)*Nx;
+ for j = 2:Nx-1
+     M(j,j) = -2;
+     M(j,j-1)=1;
+     M(j,j+1)=1;
+     M(n+j,n+j)=-2;
+     M(n+j,n+j-1)=1;
+     M(n+j,n+j+1)=1;
+
+ end
+
+ %left and Right 内点
+ for i = 2:Nx-1
+     M(Nx*(i-1)+1,Nx*(i-1)+1)=-2;
+     M(Nx*(i-1)+1,Nx*(i-2)+1)=1;
+     M(Nx*(i-1)+1,Nx*i+1)=1;
+     M(Nx*i,Nx*i)=-2;
+     M(Nx*i,Nx*(i-1))=1;
+     M(Nx*i,Nx*(i+1))=1;
+ end
+
+
+
+ %down and up 内点
+ for j = 2:Nx-1
+    [f1,f2] = Force_ring_xy(x(j),0);
+     M(j,j) = -2*f2*dx/D-2;
+     M(j,j+Nx) = 2;
+    [f1,f2] = Force_ring_xy(x(j),L);
+     M(n+j,n+j)=2*f2*dx/D+2;
+     M(n+j,n+j-Nx)=-2;
+ end
+
+  %left and Right 内点
+ for i = 2:Nx-1
+     [f1,f2] = Force_ring_xy(0,x(i));
+     M(Nx*(i-1)+1,Nx*(i-1)+1)=-2*f1*dx/D-2;
+     M(Nx*(i-1)+1,Nx*(i-1)+2)=2;
+     [f1,f2] = Force_ring_xy(L,x(i));
+     M(Nx*i,Nx*i)=2*f1*dx/D+2;
+     M(Nx*i,Nx*i-1)=-2;  
+ end
+
+ %corner
+ %左上
+ [f1,f2] = Force_ring_xy(0,L);
+ n=Nx*(Nx-1);
+ M(n+1,n+1) = 2*(f2-f1)*D*dx -4;
+ M(n+1,n-Nx+1) = 2;
+ M(n+1,n+2)=2;
+
+ %右上
+ [f1,f2] = Force_ring_xy(L,L);
+ n = Nx^2;
+ M(n,n) = 2*dx*(f1+f2)/D-4;
+ M(n,n-1)=2;
+ M(n,n-Nx)=2;
+
+ %左下
+ [f1,f2] = Force_ring_xy(0,0);
+ M(1,1)=-4-2*dx*(f1+f2)/D;
+ M(1,2)=2;
+ M(1,1+Nx)=2;
+
+ %右下
+ [f1,f2] = Force_ring_xy(L,0);
+ M(Nx,Nx) = 2*(f1-f2)*dx/D-4;
+ M(Nx,Nx-1)=2;
+ M(Nx,2*Nx)=2;
+
+
+
+
+
+ if BC == 0 %Dirichlet边界 J=0
+ %up and down 内点
+ n = (Nx-1)*Nx;
+ for j = 2:Nx-1
+
+
+ end
+
+ %left and Right 内点
+ for i = 2:Nx-1
+
+ end
+
+
+
+
+
+
+
+
+
+
+end
+
+

Functions/main.m

@@ -0,0 +1,70 @@
+Lx = 4;
+Ly = 4;
+Nx = 100;
+Ny = 100;
+m = createMesh2D(Nx,Ny,Lx,Ly);
+
+
+%边界
+D_val = 0.005; %diffusion coeeficient
+BC = createBC(m); % all Neumann boundary condition structure
+a = BC.left.a;
+BC.left.a = -D_val*a;
+BC.right.a = -D_val*a;
+BC.top.a = -D_val*a;
+BC.bottom.a = -D_val*a;
+
+[xb,yb] = Force_diff_xy(0,m.cellcenters.y); %diff/ring
+BC.left.b = transpose(xb);
+[xb,yb] = Force_diff_xy(Lx,m.cellcenters.y);
+BC.right.b = transpose(xb);
+
+[xb,yb] = Force_diff_xy(m.cellcenters.x,0);
+BC.bottom.b = transpose(yb);
+[xb,yb] = Force_diff_xy(m.cellcenters.x,Ly);
+BC.top.b = transpose(yb);
+
+c = 0.01*ones(1,Nx);
+BC.left.c = -c;
+BC.right.c = c;
+BC.top.c = c;
+BC.bottom.c = -c;
+
+%diffussion term
+D = createCellVariable(m, D_val); % assign dif. coef. to all the cells
+Dave = harmonicMean(D); % convert a cell variable to face variable
+
+%convection term
+u_face = createFaceVariable(m, 1); % assign velocity value to cell faces
+[x,y] = ndgrid(m.facecenters.x, m.cellcenters.y);
+[xv,yv] = Force_diff_xy(x,y);
+u_face.xvalue = xv;
+
+[x,y] = ndgrid(m.cellcenters.x, m.facecenters.y)
+[xv,yv] = Force_diff_xy(x,y);
+u_face.yvalue = yv;
+
+%solve
+Mconv =  convectionTerm(u_face); % convection term, central, second order
+Mdiff = diffusionTerm(Dave); % diffusion term
+[Mbc, RHS] = boundaryCondition(BC); % boundary condition discretization
+M = Mconv-Mdiff+Mbc; % matrix of coefficient for central scheme
+c = solvePDE(m, M, RHS); % solve for the central scheme
+c.value = c.value-min(min(c.value));
+c.value = log(c.value);
+%c.value = 
+visualizeCells2D(c);
+
+
+
+
+
+
+
+
+
+
+
+
+
+

Functions/resources/project/EEtUlUb-dLAdf0KpMVivaUlztwA/fMVdNXBE2QxYwW97LIvhwXs6zQod.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info Ref="" Type="Relative" />

Functions/resources/project/EEtUlUb-dLAdf0KpMVivaUlztwA/fMVdNXBE2QxYwW97LIvhwXs6zQop.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="1ce60dc3-b3ec-43cb-aa8f-087d59134185" type="Reference" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/2kj09UetkV_lru3gvSPXnY6-nM4d.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="测试" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/2kj09UetkV_lru3gvSPXnY6-nM4p.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="test" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/KKyDJtbdIBOlaeHmIZd5VX6vqx8d.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="其他" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/KKyDJtbdIBOlaeHmIZd5VX6vqx8p.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="other" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/QWNDYJD5mGW1bWYvPx9DtKnxzw4d.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="便利" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/QWNDYJD5mGW1bWYvPx9DtKnxzw4p.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="convenience" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/R1RggVhA72agIvELiuhWPRS8F0Id.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="无" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/R1RggVhA72agIvELiuhWPRS8F0Ip.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="none" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/aEHSZBIY-yve10yGis12Zr5DLZod.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="派生" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/aEHSZBIY-yve10yGis12Zr5DLZop.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="derived" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/j4xwF_j8iFTVayUMfxLgMnTbencd.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="设计" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/j4xwF_j8iFTVayUMfxLgMnTbencp.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="design" type="Label" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/r8LR4nLmg9ai3oHrW1r_-KocQzkd.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="READ_ONLY" Name="工件" />

Functions/resources/project/NjSPEMsIuLUyIpr2u1Js5bVPsOs/r8LR4nLmg9ai3oHrW1r_-KocQzkp.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="artifact" type="Label" />

Functions/resources/project/Project.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info MetadataType="fixedPathV2" />

Functions/resources/project/fjRQtWiSIy7hIlj-Kmk87M7s21k/NjSPEMsIuLUyIpr2u1Js5bVPsOsd.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info ReadOnly="1" SingleValued="1" DataType="None" Name="分类" />

Functions/resources/project/fjRQtWiSIy7hIlj-Kmk87M7s21k/NjSPEMsIuLUyIpr2u1Js5bVPsOsp.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="FileClassCategory" type="Category" />

Functions/resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/AU57kxAjJtOFVtmwUqEp7goA-j8d.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

Functions/resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/AU57kxAjJtOFVtmwUqEp7goA-j8p.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="Force_ring_xy.m" type="File" />