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 6926 2024-05-31 15:34:13Z bgrychtol $
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, 'default'); % also orders fields
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    
0068    % Electrodes and z_contact
0069 
0070    % set the z_contact
0071    n_elec= prod(size(elec_nodes));
0072    z_contact= z_contact.*ones(n_elec,1);
0073    curr_e_node=0;
0074    for i= 1:n_elec
0075       this_elec= size(elec_nodes{i},1);
0076 
0077       electrodes(i).nodes    = curr_e_node + (1:this_elec);
0078       electrodes(i).z_contact= z_contact(i);
0079 
0080       curr_e_node= curr_e_node + this_elec;
0081    end
0082 
0083 
0084    mdl.electrode =     electrodes;
0085    mdl.solve=          'eidors_default';
0086    mdl.jacobian=       'eidors_default';
0087    mdl.system_mat=     'eidors_default';
0088 
0089 function [vtx,simp,srf] = call_distmesh(fd,fh,h0,bbox, ...
0090                                  fixed_node, n_elec_nodes);
0091    [vtx,simp] = distmeshnd(fd,fh,h0,bbox,fixed_node);
0092    srf= find_boundary(simp);
0093    % Test if distmesh puts extra unneeded nodes on boundary
0094    rmnode= [];
0095    % Get srf_nodes which aren't in electrodes
0096    srf_nodes= unique( srf(:));
0097    srf_nodes= srf_nodes( srf_nodes > n_elec_nodes);
0098    for nd = srf_nodes(:)'
0099        ff= find( any(srf == nd, 2) );
0100        this_bdy= srf(ff,:);
0101        this_bdy= unique(this_bdy(:));
0102        this_bdy( find(this_bdy==nd) ) = [];
0103        if all( this_bdy < n_elec_nodes);
0104            rmnode= [rmnode, nd];
0105 %          disp([nd, this_bdy(:)']);
0106        end
0107    end
0108 
0109    if ~isempty(rmnode)
0110       vtx(rmnode,:) = [];
0111       simp = delaunayn( vtx);
0112       srf= find_boundary(simp);
0113    end
0114 
0115    % Test and remove if distmesh creates degenerate elements
0116    vol = get_elem_volume(struct('nodes',vtx,'elems',simp,'type','fwd_model'));
0117    mvol = mean(vol);
0118    degen = find(vol/mvol < 1e-10); 
0119    if length(degen)>0
0120       warning(['distmesh: degenerate (zero volume) elements created ', ...
0121               'and removed. Your mesh may contain holes.']);
0122       simp(degen,:) = [];
0123    end
0124 
0125 function h= huniform(p);
0126    h= ones(size(p,1),1);