function [ks,km,fm] = assema2DImpl(thePde, p, t, c, a, f,  u, time, sdl)
%ASSEMA2DIMPL Calculate global FE matrices for 2D models
% This undocumented function may be changed or removed in a future release.

%       Copyright 2014-2016 The MathWorks, Inc.

gotu=false;
gottime=false;
gotsdl=false;
narginchk(6,9);

% if nargin==6
  % No action
if nargin==7
  if size(u,1)>1
    gotu=true;
  else
    time=u;
    gottime=true;
  end
elseif nargin==8
  if size(u,1)>1
    gotu=true;
    gottime=true;
  else
    sdl=time;
    time=u;
    gottime=true;
    gotsdl=true;
  end
elseif nargin==9
  gotu=true;
  gottime=true;
  gotsdl=true;
end

if ~gottime
  time=[];
end

if ~gotu
  u = [];
end

if ~gotsdl
  sdl = [];
end

N = thePde.PDESystemSize;

% Choose triangles to assemble
if ~isempty(sdl)
  it=pdesdt(t,sdl);
  t=t(:,it);
end

np=size(p,2); % Number of points

% Corner point indices
it1=t(1,:);
it2=t(2,:);
it3=t(3,:);

% Triangle geometries:
[ar,g1x,g1y,g2x,g2y,g3x,g3y]=pdetrg(p,t);

% Find midpoints of triangles
x=(p(1,it1)+p(1,it2)+p(1,it3))/3;
y=(p(2,it1)+p(2,it2)+p(2,it3))/3;

sd=t(4,:);

gotu = ~isempty(u);

if ~gotu
  uu=[];
  u=[];
  ux=[];
  uy=[];
else
  uu=u;
  Nu=length(uu)/np;
  if(rem(length(uu),np))
    error(message('pde:assema:invalidLengthU'));
  end
  if(N ~= Nu)
    error(message('pde:assema:ULength', length(uu), N*np));
  end
  uu=reshape(uu,np,N);
  u=(uu(it1,:).'+uu(it2,:).'+uu(it3,:).')/3;
  ux=uu(it1,:).'.*(ones(N,1)*g1x)+uu(it2,:).'.*(ones(N,1)*g2x)+ ...
     uu(it3,:).'.*(ones(N,1)*g3x);
  uy=uu(it1,:).'.*(ones(N,1)*g1y)+uu(it2,:).'.*(ones(N,1)*g2y)+ ...
     uu(it3,:).'.*(ones(N,1)*g3y);
  uu=reshape(uu,np*N,1);
end

f=pdetfxpd(p,t,uu,time,f);
if(size(f,1) ~= N)
  error(message('pde:assema:FLength', size(f,1), N));
end


% Stiffness matrix
c=pdetfxpd(p,t,uu,time,c);

if(any(c(:)) || any(isnan(c(:))))
  ks=sparse(N*np,N*np);

  nrc=size(c,1);

  if nrc>=1 && nrc<=4 % Block scalar c
    ks1=pdeasmc(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,g1x,g1y,g2x,g2y,g3x,g3y,c);
    [ii,jj,kss]=find(ks1);
    for k=1:N
      ks=ks+sparse(ii+(k-1)*np,jj+(k-1)*np,kss,N*np,N*np);
    end
  elseif nrc==N || nrc==2*N || nrc==3*N || nrc==4*N % Block diagonal c
    nb=nrc/N;
    m1=1;
    m2=nb;
    for k=1:N
      ks1=pdeasmc(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,g1x,g1y,g2x,g2y,g3x,g3y,c(m1:m2,:));
      [ii,jj,kss]=find(ks1);
      ks=ks+sparse(ii+(k-1)*np,jj+(k-1)*np,kss,N*np,N*np);
      m1=m1+nb;
      m2=m2+nb;
    end
  elseif nrc==2*N*(2*N+1)/2 % Symmetric c
    m1=1;
    m2=4;
    for l=1:N
      for k=1:l-1
        ks1=pdeasmc(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,g1x,g1y,g2x,g2y,g3x,g3y,c(m1:m2,:));
        [ii,jj,kss]=find(ks1);
        ks=ks+sparse(ii+(k-1)*np,jj+(l-1)*np,kss,N*np,N*np);
        m1=m1+4;
        m2=m2+4;
      end
      m1=m1+3;
      m2=m2+3;
    end
    ks=ks+ks.';
    m1=1;
    m2=3;
    for k=1:N
      ks1=pdeasmc(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,g1x,g1y,g2x,g2y,g3x,g3y,c(m1:m2,:));
      [ii,jj,kss]=find(ks1);
      ks=ks+sparse(ii+(k-1)*np,jj+(k-1)*np,kss,N*np,N*np);
      m1=m1+3+4*k;
      m2=m2+3+4*k;
    end
  elseif nrc==4*N*N % General (unsymmetric) c
    m1=1;
    m2=4;
    for l=1:N
      for k=1:N
        ks1=pdeasmc(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,g1x,g1y,g2x,g2y,g3x,g3y,c(m1:m2,:));
        [ii,jj,kss]=find(ks1);
        ks=ks+sparse(ii+(k-1)*np,jj+(l-1)*np,kss,N*np,N*np);
        m1=m1+4;
        m2=m2+4;
      end
    end
  else
    error(message('pde:assema:NumRowsC'));
  end % nrc
  clear c ks1;
else
  clear c;
  ks=sparse(N*np,N*np);
end

% Mass matrix
a=pdetfxpd(p,t,uu,time,a);

if(any(a(:)) || any(isnan(a(:))))
  km=sparse(N*np,N*np);

  nra=size(a,1);

  if nra==1 % Scalar a
    km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a);
    [ii,jj,kmm]=find(km1);
    for k=1:N
      km=km+sparse(ii+(k-1)*np,jj+(k-1)*np,kmm,N*np,N*np);
    end
  elseif nra==N % Diagonal a
    for k=1:N
      km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a(k,:));
      [ii,jj,kmm]=find(km1);
      km=km+sparse(ii+(k-1)*np,jj+(k-1)*np,kmm,N*np,N*np);
    end
  elseif nra==N*(N+1)/2 % Symmetric a
    m=1;
    for l=1:N
      for k=1:l-1
        km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a(m,:));
        [ii,jj,kmm]=find(km1);
        km=km+sparse(ii+(k-1)*np,jj+(l-1)*np,kmm,N*np,N*np);
        m=m+1;
      end
      m=m+1;
    end
    km=km+km.';
    m=1;
    for k=1:N
      km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a(m,:));
      [ii,jj,kmm]=find(km1);
      km=km+sparse(ii+(k-1)*np,jj+(k-1)*np,kmm,N*np,N*np);
      m=m+1+k;
    end
  elseif nra==N*N % General (unsymmetric) a
    m=1;
    for l=1:N
      for k=1:N
        km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a(m,:));
        [ii,jj,kmm]=find(km1);
        km=km+sparse(ii+(k-1)*np,jj+(l-1)*np,kmm,N*np,N*np);
        m=m+1;
      end
    end
  else
    error(message('pde:assema:NumRowsA'));
  end % nra
  clear a km1;
else
  clear a;
  km=sparse(N*np,N*np);
end

% RHS

if(any(f(:)) || any(isnan(f(:))))
  fm=zeros(N*np,1);

  nrf=size(f,1);

  if nrf==1 % Scalar f
    fm1=pdeasmf(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,f);
    for k=1:N
      fm((k-1)*np+1:k*np,:)=fm1;
    end
  elseif nrf==N % Vector f
    for k=1:N
      fm1=pdeasmf(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,f(k,:));
      fm((k-1)*np+1:k*np,:)=fm1;
    end
  else
    error(message('pde:assema:NumRowsF'));
  end % size(f,1)
  clear f fm1;
else
  clear f;
  fm=zeros(N*np,1);
end


end

