function [element,A,b]=eleass2(element,mesh,quad,physics,varargin)

varargin

dimen=element(1).geometry.dimension;
degree=element(1).geometry.degree;

for i=1:length(varargin)
    if isa(varargin{i},'char') == 1
      switch varargin{i}
       case 'solution'
        SOLFLAG=true;
        usol=varargin{i+1};
        solptr=varargin{i+2};
       case 'adjoint'
        PHIFLAG=true;
        phisol=varargin{i+1};
        phiptr=varargin{i+2};
       case 'parameters'
        parameters=varargin{i+1};
       case 'mcpoints' %tell it to do monte carlo
        MCFLAG=true;
        mcpoints=varargin{i+1};
      end
    end
end

if exist('MCFLAG')==0 
    MCFLAG=false;
    intpts=quad;
    basisip=mesh.psiev;
    gbasisip=mesh.dpsiev;
end

MCFLAG
pts=size(mesh.p,1);
A=spalloc(pts,pts,pts*3*factorial(dimen+1)*degree);
b=zeros(size(mesh.p,1),1);
for geom=1:length(mesh.t)
  k=mesh.ptr(geom);
  if MCFLAG
    intpts=rand(dimen,mcpoints);
    switch element(k).geometry.shape 
      case 'tri'
       tempsum=sum(intpts);
       % if not in reference triangle REFLECT!
       for loop=1:mcpoints
         if tempsum(loop) > 1
           intpts(:,loop)=[1 - intpts(2,loop); 1 - intpts(1,loop)];
         end
       end
     case 'tet'
    end
    basisip=psiquadeval(mesh.cmat,intpts);
  end
  numbasis=length(element(k).geometry.nodenum);
  ei=element(k).geometry.nodenum;
  emass=zeros(numbasis,numbasis);
  erhs=zeros(numbasis,1);
  switch element(k).geometry.shape
   case {'lagrange','tet','tri'}
    vert=element(k).geometry.nodes(1:dimen+1,:);
   case {'quad','hex'}
    vert=element(k).geometry.nodes(1:2*dimen,:);
  end
  emass=zeros(numbasis);
  if isfield(physics,'a')==1 | isfield(physics,'b') == 1 | isfield(physics,'divb')==1
    for i=1:numbasis
      for j=1:dimen
        dpsi(j,:)=mesh.dpsiev(i,:,j);
      end
      gpsi(:,:,i)=element(k).geometry.jacinv*dpsi;
    end
  end
  if isfield(physics,'a')==1
    for i=1:numbasis
      for j=1:i
        emass(i,j)=integrator(vert,quad,physics.a,gpsi(:,:,i)', ...
                         gpsi(:,:,j)');
        emass(j,i)=emass(i,j);
      end
    end
  end
  if isfield(physics,'b')==1
    for i=1:numbasis
      for j=1:numbasis
        temp(j,i)=integrator(vert,quad,physics.b,gpsi(:,:,i)',mesh.psiev(j,:)');
      end
    end
    emass=emass+temp;
  end
  if isfield(physics,'divb')==1
      for i=1:numbasis
        for j=1:i
         temp(i,j)=integrator(vert,quad,physics.divb,mesh.psiev(i,:)', mesh.psiev(j,:)');
         temp(j,i)=temp(i,j);
       end
    end
    emass=emass+temp;
  end
  if isfield(physics,'c')==1
    for i=1:numbasis
      for j=1:i
        temp(i,j)=integrator(vert,quad,physics.c,mesh.psiev(i,:)', mesh.psiev(j,:)');
        temp(j,i)=temp(i,j);
      end
    end
    emass=emass+temp;
  end
  if isfield(physics,'f')==1
    if isa(physics.f,'function_handle')
      if exist('SOLFLAG')==1 
        if exist('PHIFLAG')==1
          for i=1:numbasis
            erhs(i)=integrator(vert,intpts,physics.f,'u+phi',usol,solptr,phisol, ...
                          phiptr,k,basisip(i,:)');
          end
        else
          for i=1:numbasis
            erhs(i)=integrator(vert,intpts,physics.f,'u',usol,solptr,k, ...
                            basisip(i,:)');
          end
        end
      elseif exist('PHIFLAG')==1
        for i=1:numbasis
          erhs(i)=integrator(vert,intpts,physics.f,'phi',phisol, ...
                        phiptr,k,basisip(i,:)');
        end
      else
        % user defined rhs
        if exist('parameters')==1
          erhs=feval(physics.f,vert,mesh.cmat,element(k).geometry.jacinv,parameters);
        else
          erhs=feval(physics.f,vert,mesh.cmat,element(k).geometry.jacinv);
        end
      end
    else
      for i=1:numbasis
        erhs(i)=integrator(vert,mesh.quad,physics.f,mesh.psiev(i,:)');
      end
    end
  end
  element(k).physics.mass=emass;
  element(k).physics.rhs=erhs;
  A(ei,ei)=A(ei,ei)+emass;
  b(ei)=b(ei)+erhs;
end
if exist('A')==0
  A=[];
end