dm_mk_fwd_model

 DM_MK_FWD_MODEL: create a fwd_model object using distmesh
 fwd_mdl= dm_mk_fwd_model( fd, fh, h0, bbox,...
                          elec_nodes, stim_pattern, z_contact, name);

  fd:        Shape function - +ve for space inside area, -ve outside
  fh:        Mesh refinement function - use [] for uniform
  nnodes:    Estimate of number of nodes in model
  bbox:      Bounding box [xmin,ymin, {zmin}; xmax,ymax,{zmax}]

  elec_nodes:        cell of matrix N x [x,y,{z}] for each electrode
  refine_nodes:      vector of fixed nodes to add to model to refine it
  z_contact:         vector or scalar electrode contact impedance
  name:              name for eidors object

  fwd_mdl:           eidors format fwd_model

0001 function [fwd_mdl]= dm_mk_fwd_model( fd, fh, nnodes, bbox, elec_nodes, ...
0002                      refine_nodes, z_contact, name)
0003 % DM_MK_FWD_MODEL: create a fwd_model object using distmesh
0004 % fwd_mdl= dm_mk_fwd_model( fd, fh, h0, bbox,...
0005 %                          elec_nodes, stim_pattern, z_contact, name);
0006 %
0007 %  fd:        Shape function - +ve for space inside area, -ve outside
0008 %  fh:        Mesh refinement function - use [] for uniform
0009 %  nnodes:    Estimate of number of nodes in model
0010 %  bbox:      Bounding box [xmin,ymin, {zmin}; xmax,ymax,{zmax}]
0011 %
0012 %  elec_nodes:        cell of matrix N x [x,y,{z}] for each electrode
0013 %  refine_nodes:      vector of fixed nodes to add to model to refine it
0014 %  z_contact:         vector or scalar electrode contact impedance
0015 %  name:              name for eidors object
0016 %
0017 %  fwd_mdl:           eidors format fwd_model
0018 
0019 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0020 % $Id: dm_mk_fwd_model.m 5505 2017-06-05 14:25:08Z aadler $
0021 
0022 if nargin <7
0023    error('dm_mk_fwd_model requires 7 or 8 parameters');
0024 else
0025    name = 'MDL from dm_mk_fwd_model';
0026 end
0027 
0028 if isempty(fh); fh= @huniform; end
0029 
0030 h0= estimate_h0(bbox, nnodes);
0031 fwd_mdl= create_refined_model(name, fd, fh, h0, bbox, elec_nodes, ...
0032                         refine_nodes, z_contact); 
0033 
0034 fwd_mdl = mdl_normalize( fwd_mdl, 0);
0035 
0036 % estimate initial edge length to get nnodes
0037 function  h0= estimate_h0(bbox, nnodes);
0038    dims= size(bbox,2);
0039    area_est= prod(abs(diff(bbox,[],1)));
0040    h0 = (area_est/nnodes)^(1/dims);
0041 
0042 function fmdl= create_refined_model(name, fd, fh, h0, bbox, elec_nodes, ...
0043                         refine_nodes, z_contact); 
0044 
0045    % fixed_node= [elec_nodes{:}]; - wish we could do it like this - matlab bug!!
0046    fixed_node= [];
0047    for i= 1:prod(size(elec_nodes)) 
0048       fixed_node= [fixed_node; elec_nodes{i}];
0049    end
0050    n_elec_nodes= size(fixed_node,1);
0051    fixed_node= [fixed_node; refine_nodes];
0052    [vtx,simp,srf] = call_distmesh(fd,fh, h0,bbox,fixed_node, n_elec_nodes);
0053 
0054    fmdl = construct_fwd_model(srf,vtx,simp, name, ...
0055                           elec_nodes, z_contact);
0056 
0057 
0058 % build fwd_model structure
0059 function mdl= construct_fwd_model(srf,vtx,simp, name, ...
0060                        elec_nodes, z_contact)
0061    mdl= eidors_obj('fwd_model', name);
0062 
0063    mdl.nodes    = vtx;
0064    mdl.elems    = simp;
0065    mdl.boundary = srf;
0066    mdl.gnd_node=           1;
0067    mdl.name = name;
0068 
0069    % Electrodes and z_contact
0070 
0071    % set the z_contact
0072    n_elec= prod(size(elec_nodes));
0073    z_contact= z_contact.*ones(n_elec,1);
0074    curr_e_node=0;
0075    for i= 1:n_elec
0076       this_elec= size(elec_nodes{i},1);
0077 
0078       electrodes(i).nodes    = curr_e_node + (1:this_elec);
0079       electrodes(i).z_contact= z_contact(i);
0080 
0081       curr_e_node= curr_e_node + this_elec;
0082    end
0083 
0084 
0085    mdl.electrode =     electrodes;
0086    mdl.solve=          'eidors_default';
0087    mdl.jacobian=       'eidors_default';
0088    mdl.system_mat=     'eidors_default';
0089 
0090 function [vtx,simp,srf] = call_distmesh(fd,fh,h0,bbox, ...
0091                                  fixed_node, n_elec_nodes);
0092    [vtx,simp] = distmeshnd(fd,fh,h0,bbox,fixed_node);
0093    srf= find_boundary(simp);
0094    % Test if distmesh puts extra unneeded nodes on boundary
0095    rmnode= [];
0096    % Get srf_nodes which aren't in electrodes
0097    srf_nodes= unique( srf(:));
0098    srf_nodes= srf_nodes( srf_nodes > n_elec_nodes);
0099    for nd = srf_nodes(:)'
0100        ff= find( any(srf == nd, 2) );
0101        this_bdy= srf(ff,:);
0102        this_bdy= unique(this_bdy(:));
0103        this_bdy( find(this_bdy==nd) ) = [];
0104        if all( this_bdy < n_elec_nodes);
0105            rmnode= [rmnode, nd];
0106 %          disp([nd, this_bdy(:)']);
0107        end
0108    end
0109 
0110    if ~isempty(rmnode)
0111       vtx(rmnode,:) = [];
0112       simp = delaunayn( vtx);
0113       srf= find_boundary(simp);
0114    end
0115 
0116    % Test and remove if distmesh creates degenerate elements
0117    vol = get_elem_volume(struct('nodes',vtx,'elems',simp,'type','fwd_model'));
0118    mvol = mean(vol);
0119    degen = find(vol/mvol < 1e-10); 
0120    if length(degen)>0
0121       warning(['distmesh: degenerate (zero volume) elements created ', ...
0122               'and removed. Your mesh may contain holes.']);
0123       simp(degen,:) = [];
0124    end
0125 
0126 function h= huniform(p);
0127    h= ones(size(p,1),1);