Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % Poisson FEM based on MIT lectures by Qiqi Wang
- % Export the mesh data [p,e,t] in the wokspace
- % Compute the gradient of the basis functions
- grad = zeros(size(t,2),3,2);
- % loop over number of triangles
- for k = 1: size(t,2)
- % Store the coordinates of the triangles
- p1 = p(:,t(1,k));
- p2 = p(:,t(2,k));
- p3 = p(:,t(3,k));
- % Compute the gradient
- % the coefficient matrix: differences of coordinates
- C1 = [(p1-p2)';(p1-p3)'];
- % The RHS. Note that the basis function has the value 1 at p1 (say)
- r1 = [1;1];
- grad(k,1,:) = C1 r1;
- % Do it for the other nodes in a similar fashion
- C2 = [(p2-p1)';(p2-p3)'];
- grad(k,2,:) = C2 r1;
- C3 = [(p3-p1)';(p3-p2)'];
- grad(k,3,:) = C3 r1;
- end
- % Define the coefficient matrix
- C = zeros(size(p,2));
- % Note: the size of A must be the number of nodes. Ideally...
- % number of interior nodes.
- % Compute the area of each finite elements (triangles)
- % loop over number of triangles
- A = zeros (size(t,2),1);
- for k = 1:size(t,2)
- % Store the coordinates of the triangles
- p1 = p(:,t(1,k));
- p2 = p(:,t(2,k));
- p3 = p(:,t(3,k));
- % Area(k) = det (coefficient matrix of basis) /2
- A(k) = abs(det([(p1-p2)';(p1-p3)']))/2;
- end
- % Compute the stiffness matrix
- for k = 1:size(t,2)
- cr = t(1:3,k)';
- for i = 1:3
- for j = 1:3
- C(cr(i) , cr(j)) = C(cr(i) , cr(j)) +...
- A(k) * dot(grad(k,i,:),grad(k,j,:));
- end
- end
- end
- % Compute the global load vector (RHS)
- b = zeros(size(p,2),1);
- % Define the RHS function
- f=1;
- for i=1:size(t,2)
- xcell=num2cell(t(1:3,i));
- [i1,i2,i3] = deal(xcell{:});
- % Compute the centroid of the triangle
- centroid = (p(:,i1)+p(:,i2)+p(:,i3))/3;
- % The integral is appoximated by f(centroid)*area/3
- b(i1) = b(i1)+A(i)/3;
- b(i2) = b(i2)+A(i)/3;
- b(i3) = b(i3)+A(i)/3;
- end
- M=size(p,2);
- % Define the solution matrix
- sol = diag(zeros(M));
- % Get the number of interior nodes
- dof=setdiff(1:M,e);
- % Get the solution (Finally !)
- sol(dof) = C(dof,dof)b(dof);
- % Plot the solution
- figure
- a=trisurf(t',p(1,:)',p(2,:)',sol);
Add Comment
Please, Sign In to add comment