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');endif nargin < 3  conduct = '(const 1.0)';endif nargin < 4  source = '(const 0.0)';endif 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');endif length(brep) < 1  error('First argument passed to gmfem not a brep');endglobal GM_BREP_TYPE_CODEtypecode = double(brep{0});if ~strcmpi(typecode, GM_BREP_TYPE_CODE)  error('First argument passed to gmfem not a brep');endif ~isa(mesh,'zba')  error('Second argument passed to gmfem not a zba');endif length(mesh) < 1  error('Second argument passed to gmfem not a mesh');endglobal GM_SIMPCOMP_TYPE_CODEtypecode = double(mesh{0});if ~strcmpi(typecode, GM_SIMPCOMP_TYPE_CODE)  error('Second argument passed to gmfem not a mesh');endglobal 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');enddi = double(brep{2});if di ~= 2 & di ~= 3  error('Only dims 2 and 3 supported.')endif brep{1} ~= di  error('Brep must be full dimensional.')end if mesh{2} ~= di  error('Dimension mismatch between mesh and brep')endif 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');endnumtri = 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;   sourceterm1 = gm_source(source1, userdata, midpoints);  [m,n] = size(sourceterm1);  if m ~= numtri1 | n ~= 1    error('Source-term call returned a vector of the wrong size.');  end  sourceterm(pos:pos+numtri1 - 1) = sourceterm1;  pos = pos + numtri1;endif 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  endend%% 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  endend  % 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% ------------------------------------------------------------------