function u = gmfem(brep, mesh, conduct, source, userdata);
%   u = gmfem(brep, mesh {, conduct { , source { , userdata}}});
% This function solves the boundary value problem
%    div (c grad u) = -f
% The arguments are: 
%   brep: the domain on which the bvp is specified.
%   mesh: a mesh of the domain, generated by gmmeshgen.
%     Quadratic elements not supported; all elements interpreted
%     as linear.
%   conduct: A function that computes conductivity of the domain.  
%      Overridden by the region's conductivity property.  Default is 1.
%   source: A function that returns f, the source term.  Overridden by
%     the regions's source property. Default is 0.
%  Return value u is a vector of nodal values of the solution.  It
%  is a zba.

if nargin < 2
  error('gmfem must have at least two arguments');
end
if nargin < 3
  conduct = '(const 1.0)';
end

if nargin < 4
  source = '(const 0.0)';
end

if nargin < 5
   userdata = '';
end


%% This routine solves the finite element problem using
%% the default ordering (min degree).

if ~isa(brep,'zba')
  error('First argument passed to gmfem not a zba');
end

if length(brep) < 1
  error('First argument passed to gmfem not a brep');
end
global GM_BREP_TYPE_CODE
typecode = double(brep{0});
if ~strcmpi(typecode, GM_BREP_TYPE_CODE)
  error('First argument passed to gmfem not a brep');
end

if ~isa(mesh,'zba')
  error('Second argument passed to gmfem not a zba');
end

if length(mesh) < 1
  error('Second argument passed to gmfem not a mesh');
end

global GM_SIMPCOMP_TYPE_CODE
typecode = double(mesh{0});
if ~strcmpi(typecode, GM_SIMPCOMP_TYPE_CODE)
  error('Second argument passed to gmfem not a mesh');
end

global GM_GEO_ID_PROP;
if ~strcmpi(gm_lookup_objprop(brep, GM_GEO_ID_PROP), ...
            gm_lookup_objprop(mesh, GM_GEO_ID_PROP))
  error('Brep-id code of mesh does not match with brep');
end

di = double(brep{2});
if di ~= 2 & di ~= 3
  error('Only dims 2 and 3 supported.')
end

if brep{1} ~= di
  error('Brep must be full dimensional.')
end
 
if mesh{2} ~= di
  error('Dimension mismatch between mesh and brep')
end

if mesh{1} ~= di
  error('Mesh must be full dimensional')
end

[scrap,numvtx] = size(mesh{4});
[scrap,numtoplev] = size(mesh{5+di});
[scrap,numtoplev2] = size(brep{5+di});
if numtoplev ~= numtoplev2
  error('Number of regions disagree between mesh and brep');
end
[scrap,numfacet] = size(mesh{4+di});
[scrap,numfacet2] = size(brep{4+di});
if numfacet ~= numfacet2
  error('Number of facets disagree between mesh and brep');
end

numtri = 0;
tri = zba([]);
for faceind = 0 : numtoplev - 1
  [scrap, numtri1] = size(mesh{5+di}{1,faceind});
  tri = [tri;mesh{5+di}{1,faceind}'];  
  numtri = numtri + numtri1;
end


[scrap,numtri] = size(mesh{5});
[scrap,numtoplev] = size(brep{5+di});


% Evaluate the conductivity and source.  Loop over top-level
% faces of the brep and find the conductivity and source for each.

con = zba(-ones(numtri,1));
sourceterm = zba(zeros(numvtx,1));
vtx_global_to_ordinal = ...
   sparse((double(mesh{4}(0,:)))'+1,ones(numvtx,1), (0:numvtx-1)');

vtx = mesh{4}(1:di,:)';

tri = zba(full(vtx_global_to_ordinal(double(tri)+1)));


pos = 0;

for fa = 0 : numtoplev - 1
  conduct1 = gm_lookup_prop(brep, di, fa, 'conductivity');
  if length(conduct1) == 0
    conduct1 = conduct;
  end
  facename = double(brep{5+di}{0,fa});
  disp(sprintf('Conductivity on face %s is %s', facename, conduct1))
  source1 = gm_lookup_prop(brep, di, fa, 'source');
  if length(source1) == 0
    source1 = source;
  end
  disp(sprintf('Source on face %s is %s', facename, source1))


%  Evaluate conductivity at the midpoints of the mesh elements
% (one-point integration).


% tri1 has the vertices of the triangles in
% the current toplev face, converted to ordinal numbering.

  tri1 = full(vtx_global_to_ordinal(double(mesh{5+di}{1,fa})+1))';
  [numtri1,scrap] = size(tri1);

  if di == 2
    midpoints = (vtx(tri1(:,1),:) + vtx(tri1(:,2),:) + ...
                 vtx(tri1(:,3),:)) / 3;
  else
    midpoints = (vtx(tri1(:,1),:) + vtx(tri1(:,2),:) + ...
                 vtx(tri1(:,3),:) + vtx(tri1(:,4),:)) / 4;
  end
  con1 = gm_conductivity(conduct1, userdata, midpoints);
  [m,n] = size(con1);
  if m ~= numtri1 | n ~= 1
    error('Conductivity call returned a vector of the wrong size.')
  end
  if any(con1 <= 0) 
    error('Conductivity routine returned nonpositive conductivity');
  end
  con(pos:pos+numtri1 - 1) = con1;

  % Evaluate the conductivity at all the mesh nodes of the
  % top level face.  First, let is_on_face be a mask for
  % vertices of the top-level face.

  is_on_face0 = sparse([], [], [], numvtx, 1);
  for i = 1 : di + 1
    is_on_face0(tri1(:,i)+1) = ones(numtri1,1);
  end
  is_on_face = find(is_on_face0) - 1;
  num_on_face = length(is_on_face);
  
  % Evaluate source term at those vertices.
  sourcec1 = gm_source(source1, userdata, vtx(is_on_face,:));
  
  [m,n] = size(sourcec1);
  if m ~= num_on_face | n ~= 1
    error('Source-term call returned a vector of the wrong size.');
  end
  sourcec = zba(sparse(is_on_face+1,ones(num_on_face,1),double(sourcec1),numvtx,1));

  % Compute simplex volumes as determinants.

  if di == 2
    trivol1 = (vtx(tri1(:,2),0) - vtx(tri1(:,1),0)) .* ...
              (vtx(tri1(:,3),1) - vtx(tri1(:,1),1));
    trivol2 = (vtx(tri1(:,3),0) - vtx(tri1(:,1),0)) .* ...
              (vtx(tri1(:,2),1) - vtx(tri1(:,1),1));
    trivol = trivol1 - trivol2;
  else
    trivol1 = (vtx(tri1(:,2),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,3),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,4),2) - vtx(tri1(:,1),2));
                
    trivol2 = (vtx(tri1(:,2),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,4),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,3),2) - vtx(tri1(:,1),2));

    trivol3 = (vtx(tri1(:,3),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,2),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,4),2) - vtx(tri1(:,1),2));

    trivol4 = (vtx(tri1(:,3),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,4),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,2),2) - vtx(tri1(:,1),2));

    trivol5 = (vtx(tri1(:,4),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,2),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,3),2) - vtx(tri1(:,1),2));

    trivol6 = (vtx(tri1(:,4),0) - vtx(tri1(:,1),0)) .*...
              (vtx(tri1(:,3),1) - vtx(tri1(:,1),1)) .*...
              (vtx(tri1(:,2),2) - vtx(tri1(:,1),2));

    trivol = trivol1 - trivol2 - trivol3 + trivol4 + trivol5 - trivol6;
  end

  % The volumes should all be the same sign.
  
  if (any(trivol * trivol(0)) <= 0) 
    error('Nonpositive simplex volume');
  end
  if trivol(0) < 0
    trivol = -trivol;
  end

  
% integrate the source term.  Use a low-order integration rule
% that is exact for constant and linear source terms.

  if di == 2
    fac1 = 1/24;
  else
    fac1 = 1/120;
  end
  [numtriface, scrap] = size(tri1);
  integ0 = zba(zeros(numtriface,1));
  trivol1 = trivol * fac1;
  for k = 1 : di + 1
   integ0 = integ0 + trivol1 .* (sourcec(tri1(:,k)));
  end
  
  for k = 1 : di + 1
    integ1 = integ0 + trivol1 .* (sourcec(tri1(:,k)));
    
    % add integ1(j) to sourceterm(tri1(j,k)).  Do this vectorized
    % over j with a sparse matrix trick.
    ss = zba(sparse(1:numtri1, tri1(:,k)+1, double(integ1), numtri1, numvtx));
    sourceterm = sourceterm + sum(ss)';
  end
  pos = pos + numtri1;
end

if any(con <= 0)
  error('Some conductivities were not evaluated');
end

% Set up dirichlet and Neuman BC's.

dirval = zba(zeros(numvtx,1));
isdir = zba(zeros(numvtx,1));
dircount = zba(zeros(numvtx,1));
neusum = zba(zeros(numvtx,1));

% Loop over faces of dimension di-1 and find out
% their boundary values.

for fa = 0 : numfacet - 1
  bc = gm_strim(gm_lookup_prop(brep, di - 1, fa, 'bc'));
  if length(bc) == 0
    btype = 'n';
    bfunc = '(const 0.0)';
  else

    %parse the boundary specification.

    if bc(1) ~= '(' | bc(end) ~= ')'
      error('Boundary condition must be an ordered pair (type func)');
    end
    bc = gm_strim(bc(2:length(bc)-1));
    [btype,bfunc] = strtok(bc);
    btype = lower(btype(1));
    if btype ~= 'n' & btype ~= 'd'
      error('Bc type must be either d or n');
    end
    bfunc = gm_strim(bfunc);
    if length(bfunc) == 0
      error('No boundary condition function specified')
    end
  end

  facename = double(brep{4+di}{0,fa});
  disp(sprintf('facet %s btype = %s bfunc = %s', ...
	facename, btype, bfunc))

  % Get the facets lying on this face, converted to
  % ordinal indices.

  facetlist = ...
    zba(full(vtx_global_to_ordinal(double(mesh{4+di}{1,fa})+1))');
  [numfacet1, scrap] = size(facetlist);
  
  % String together the vertices of these facets and
  % deduplicate for dirichlet bc.

  vtxselect = sparse([],[],[],numvtx,1);
  vtxselect(double(facetlist(:)) + 1) = ones(numfacet1 * di, 1);
  vtxlist = find(vtxselect) - 1;
  vtx1 = vtx(vtxlist,:);

  [numvtx1, scrap] = size(vtx1);
  if numvtx1 == 0
    disp(sprintf('Warning: no vertices on facet %d:%d  %s', ...
        di-1, fa, brep{4+di}{0,fa}))
  else
    bcval = gm_bdrycond(bfunc, userdata, vtx1);	
    [m,n] = size(bcval);
    if m ~= numvtx1 | n ~= 1
      error('Bfunc return value is the wrong size');
    end
    if btype == 'd'
      dirval(vtxlist) = dirval(vtxlist) + bcval;
      dircount(vtxlist) = dircount(vtxlist) + ones(numvtx1,1);
      isdir(vtxlist) = ones(numvtx1,1);
    else
      neusum  = neusum + gm_neumann(vtx, facetlist, bcval);
    end
  end
end

%% dirval holds the sums of the dirichlet boundary values;
%% divide by dircount to make them average dirichlet values.

dirval(find(isdir)) = dirval(find(isdir)) ./ dircount(find(isdir));


%% At least one dirichlet node per connected component to make
%% system nonsingular.

compo_label = gm_mcompo(mesh);
numcompo = max(double(compo_label)) + 1;
for componum = 0 : numcompo - 1
  if all(isdir(compo_label == componum) == 0)
    disp(sprintf('Note: all vertices of component %d have Neumann condition',...
          componum));
    disp(sprintf('Boundary integral of Neumann data = %e',...
       sum(neusum(compo_label == componum))));
    disp('(Should be close to zero)');
    ff = find(compo_label == componum);
    if length(ff) > 0
      disp(sprintf('Changing node %d to Dirichlet condition to make the problem well-posed', ff(1)));
      isdir(ff(0)) = 1;
      dirval(ff(0)) = 0;
    end
  end
end
  

% Assemble the fe matrix and solve it.
[K,f0] = gm_assemble(vtx, tri, isdir, dirval, con);
neusum1 = neusum(isdir == 0);
sourceterm1 = sourceterm(isdir == 0);
f = f0 + neusum1 + sourceterm1;
[m,n] = size(K);
disp(sprintf('Assembled K is %d-by-%d and contains %d nonzero entries',...
 m,n,nnz(K)));


startflops = flops;
u0 = K \ f;
endflops = flops;
disp(sprintf('%d flops to solve linear system', endflops-startflops))

% u0 has a value for all non-dirichlet boundary nodes.
% Insert dirichlet nodes.

u = zba(zeros(numvtx,1));
u(find(isdir == 0)) = u0;
u(find(isdir == 1)) = dirval(isdir == 1);


% ------------------------------------------------------------------
% Copyright (c) 1999 by Cornell University.  All rights reserved.
% See the accompanying file 'Copyright' for authorship information,
% the terms of the license governing this software, and disclaimers
% concerning this software.
% ------------------------------------------------------------------
% This file is part of the QMG software.  
% Version 2.0 of QMG, release date September 3, 1999
% ------------------------------------------------------------------