function MyPLHeatEqSolver() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Solves the pde % % du/dt + div q = f, q = - d grad u, % % with boundary condition % % n q = k (u - g) % % and initial condition % % u = u0 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % defining geometry and meshing geom = geometry; [p,e,t] = initmesh(geom, 'HMax', .1); % [p,e,t] = refinemesh(geom, p, e, t); pdemesh(p, e, t); shg pause(1) % defining data functions x1 = p(1,:); x2 = p(2,:); nn = size(p,2); d = ones(1,nn); m = ones(1,nn); f = zeros(1,nn); k = 100000 * ones(1,nn); g0 = x1.*x2.*(1-x2); % preliminary, stationary part of g u0 = zeros(nn,1); % initial values % initialising matrises D = sparse(zeros(nn,nn)); M = sparse(zeros(nn,nn)); K = sparse(zeros(nn,nn)); F = zeros(nn,1); G = zeros(nn,1); U = u0; % storing % assembling element contributions for el = 1 : size(t,2) elNodeInd = t(1: 3, el); elCoords = p(:,elNodeInd); elDiff = d(elNodeInd); dD = elDiffMatrix(elCoords, elDiff); D(elNodeInd, elNodeInd) = D(elNodeInd, elNodeInd) + dD; elMass = m(elNodeInd); dM = elMassMatrix(elCoords, elMass); M(elNodeInd,elNodeInd) = M(elNodeInd,elNodeInd) + dM; elLoad = f(elNodeInd); dF = elLoadVector(elCoords, elLoad); F(elNodeInd) = F(elNodeInd) + dF; end % assembling boundary element conributions for bel = 1 : size(e, 2) belNodeInd = e(1: 2, bel); belCoords = p(:, belNodeInd); belPerm = k(belNodeInd); dK = elBdryMatrix(belCoords, belPerm); K(belNodeInd, belNodeInd) = K(belNodeInd, belNodeInd) + dK; % belLoad = g(belNodeInd); % dG = elBdryVector(belCoords, belPerm, belLoad); % G(belNodeInd) = G(belNodeInd) + dG; end % initializations time = 0; u = u0; FinalTime = 3; dt = .03; A = M + dt/2 * D + dt/2 * K; L = M - dt/2 * D - dt/2 * K; % time stepping h = waitbar(0,'Just a minute..'); while time < FinalTime G = zeros(nn,1); for bel = 1 : size(e, 2) belNodeInd = e(1: 2, bel); belCoords = p(:, belNodeInd); belPerm = k(belNodeInd); belLoad = g0(belNodeInd) * sin(2*pi*(time + dt/2)); dG = elBdryVector(belCoords, belPerm, belLoad); G(belNodeInd) = G(belNodeInd) + dG; end b = L * U + dt * F + dt * G; U = A \ b; u = [u U]; time = time + dt; waitbar(time/FinalTime) end delete(h) % plotting set(gcf,'DoubleBuffer','on') for i = 1 : size(u, 2) pdeplot(p,e,t,'xydata',U,'zdata',u(:,i),'mesh','on'); set(gca,'ZLim',[-.5 .5]) xlabel('x_1') ylabel('x_2') zlabel('y') title(['solution y = u(x_1, x_2, t), t = ' num2str((i-1)*dt,'%1.1f')]) drawnow end % subroutines function dD = elDiffMatrix(coords, diffusion) dx = area(coords); x = sum(coords, 2) / 3; Dphi = Dphi(coords); d = sum(diffusion)/length(diffusion); dD = Dphi' * (d * Dphi) * dx; function dM = elMassMatrix(coords, mass) dx = area(coords); x = sum(coords, 2) / 3; phi = [1/3 1/3 1/3]; m = sum(mass)/length(mass); dM = phi' * (m * phi) * dx; function dF = elLoadVector(coords, domainLoad) dx = area(coords); x = sum(coords, 2) / 3; phi = [1/3 1/3 1/3]; f = sum(domainLoad)/length(domainLoad); dF = phi' * f * dx; function dK = elBdryMatrix(coords, bdryPerm) v = coords(:, 2) - coords(:, 1); ds = sqrt(v' * v); p = 0.57735026918963; points = [.5-p/2; .5+p/2]; weights = [.5 .5]; x = sum(coords, 2) / 2; phi = [1/2 1/2]; phi1 = [1-points(1) points(1)]; phi2 = [1-points(2) points(2)]; k = sum(bdryPerm)/length(bdryPerm); dK = phi1' * (k * phi1) * ds / 2 ... + phi2' * (k * phi2) * ds / 2; function dG = elBdryVector(coords, bdryPerm, bdryLoad) v = coords(:, 2) - coords(:, 1); ds = sqrt(v' * v); x = sum(coords, 2) / 2; phi = [1/2 1/2]; k = sum(bdryPerm)/length(bdryPerm); g = sum(bdryLoad)/length(bdryLoad); dG = phi' * k * g * ds; function dx = area(coords) v = coords(:, 2) - coords(:, 1); w = coords(:, 3) - coords(:, 1); dx = (v(1) * w(2) - v(2) * w(1)) / 2; function Dphi = Dphi(coords) v = coords(:, 2) - coords(:, 1); w = coords(:, 3) - coords(:, 1); n2v = [-v(2) v(1)]'; % normal to v n2v = n2v / sqrt(n2v' * n2v); % normalizing pw = w' * n2v; % height of triangle Dphi(:, 3) = n2v / pw; n2w = [w(2) -w(1)]'; n2w = n2w / sqrt(n2w' * n2w); pv = v' * n2w; Dphi(:, 2) = n2w / pv; Dphi(:, 1) = - Dphi(:, 2) - Dphi(:, 3); % data functions function geom = geometry() e1 = [2 0 1 0 0 1 0]'; e2 = [2 1 1 0 1 1 0]'; e3 = [2 1 0 1 1 1 0]'; e4 = [2 0 0 1 0 1 0]'; geom = [e1 e2 e3 e4];