GMSH Meshing updates for EIDORS

call_gmsh.m

function status= call_gmsh(geo_file,dim)
% call_gmsh: call Gmsh to create a vol_file from a geo_file
% status= call_netgen(geo_file, vol_file, msz_file, finelevel)
%  staus = 0 -> success , negative -> failure
%
% geo_file = geometry file (input)
% vol_file = FEM mesh file (output)
% msz_file = Meshsize file in netgen format
%
% Finelevel controls the fineness of the mesh
%   default is '' -> coarse
%   valid values are 'fine' or 'veryfine'
%
% $Id: call_gmsh.m 2628 2011-07-11 13:13:50Z aadler $
% (C) 2009 Bartosz Sawicki. Licensed under GPL V2

% default to 2-D model
if nargin<2
    dim = 2;
end

% Gmsh executable filename
gmsh_name = 'gmsh';

while( 1 )
   ldpath='';
    if exist('OCTAVE_VERSION') && strcmp(uname.sysname,'Linux')
     islinux =1;
   elseif  strfind(system_dependent('getos'),'Linux')
     islinux =1;
     s=version; ff= find(s=='.');
      if str2double(s(1:ff(2)-1))>=7
        %Version 7 under linux sets the LD_LIBRARY_PATH and that breaks netgen    
          ldpath ='LD_LIBRARY_PATH=;';
      end      
   else
     islinux =0;
   end    
   switch dim
       case 2
           status= system(sprintf( '%s %s %s -2 -v 2', ldpath, gmsh_name, geo_file));
       case 3
           status= system(sprintf( '%s %s %s -3 -v 2', ldpath, gmsh_name, geo_file));
   end
   

   if status==0; break; end

   if islinux
      disp(['It seems you are running Linux and Gmsh has not worked. ' ...
           'Check that it is installed and on the path. \n' ...
           'Perhaps LD_LIBRARY_PATH needs to be set?' ]);
       break;
   else
      fprintf([ ...
       'Gmsh call failed. Is Gmsh installed and on the search path?\n' ...
       'You are running under windows, I can attempt to create\n' ...
       'a batch file to access gmsh.\n' ...
       'Please enter the directory in which to find gmsh.\n' ...
       'If you dont have a copy, download it from' ...
       'http://www.geuz.org/gmsh/\n\n' ...
       'Note that you *MUST* use names without spaces. Thus\n' ...
       'instead of C:/Program Files/ write C:/Progra~1/\n\n' ]);
%      gmsh_path = input('gmsh_path? ','s');
%      if exist( sprintf('%s/gmsh.exe',gmsh_path) , 'file' ) 
%         disp('Found gmsh.');
%       fid= fopen('ng.bat','w');
%       fprintf(fid,'set TCL_LIBRARY=%s/lib/tcl8.3\n', netgen_path);
%       fprintf(fid,'set TIX_LIBRARY=%s/lib/tix8.1\n', netgen_path);
%       fprintf(fid,'%s/netgen.exe %%*\n', netgen_path);
%       fclose(fid);
%       elseif exist( sprintf('%s/ng431.exe',netgen_path) , 'file' ) 
%          disp('Found netgen version 4.3.1');
%          fid= fopen('ng.bat','w');
%          fprintf(fid,'set TCL_LIBRARY=%s/lib/tcl8.3\n', netgen_path);
%          fprintf(fid,'set TIX_LIBRARY=%s/lib/tcl8.2\n', netgen_path);
%          fprintf(fid,'%s/ng431.exe %%*\n', netgen_path);
%          fclose(fid);
%       else
%          warning(['cannot find a version of netgen that I know about\n' ...
%                   'Install netgen 4.4 or 4.3.1 or check the path\n']);
%       end
        break;
   end
end

create_2d2l_mesh_gmsh.m

function [mdl, mat_indices] = create_2d2l_mesh_gmsh(basename, in_rad, ext_rad, no_elec)
% Create a 2D Circular FEM using GMSH
% mdl= CREATE_GMSH_2D_CIRCLE(rad, n_elec)
%
% mdl - EIDORS forward model
% rad - model radius

% (C) 2006 Andy Adler. License: GPL version 2 or version 3
% $Id: create_2d2l_mesh_gmsh.m 2628 2011-07-11 13:13:50Z aadler $


geo_filename = sprintf('%s.geo', basename);
geo_fid= fopen(geo_filename,'w');

% Electrodes points on external boundary
theta= linspace(0,2*pi, no_elec+1); theta(end)=[];
point_no = 1;
for th=theta
    x=ext_rad*sin(th);
    y=ext_rad*cos(th);
    z=0;
    fprintf(geo_fid,'Point(%d) = {%f,%f,%f,%f};\n',point_no, x,y,z, ext_rad/100);
    fprintf(geo_fid,'Physical Point("%s%d") = {%d};\n','electrode-',point_no, point_no);
    point_no = point_no + 1;
end

% Points to define internal layer
center_no = point_no;
fprintf(geo_fid,'Point(%d) = {%f,%f,%f,%f};\n',center_no, 0, 0, 0, 1);
point_no = point_no + 1;
inpoint1_no = point_no;
fprintf(geo_fid,'Point(%d) = {%f,%f,%f,%f};\n',inpoint1_no, in_rad, 0,0, ext_rad/30);
point_no = point_no + 1;
inpoint2_no = point_no;
fprintf(geo_fid,'Point(%d) = {%f,%f,%f,%f};\n',inpoint2_no, -in_rad, 0,0, ext_rad/30);
point_no = point_no + 1;

% External lines / electrodes
line_no = 1;
for i = 1:no_elec
    start_point = i;
    end_point= i+1;
    if end_point > no_elec
        end_point = 1;
    end
    fprintf(geo_fid,'Line(%d) = {%d, %d};\n', line_no, start_point, end_point);
    line_no = line_no + 1;
end

extloop_no = line_no;
line_no = line_no + 1;
fprintf(geo_fid,'Line Loop(%d) = {', extloop_no);
for i = 1:(no_elec-1)
    fprintf(geo_fid,'%d,', i);
end
fprintf(geo_fid, '%d};\n', no_elec);

% Internal circle and loop
circle1_no = line_no;
line_no = line_no + 1;
fprintf(geo_fid,'Circle(%d) = {%d, %d, %d};\n', circle1_no, inpoint1_no, ...
    center_no, inpoint2_no);
circle2_no = line_no;
line_no = line_no + 1;
fprintf(geo_fid,'Circle(%d) = {%d, %d, %d};\n', circle2_no, inpoint2_no, ...
    center_no, inpoint1_no);

inloop_no = line_no;
line_no = line_no + 1;
fprintf(geo_fid,'Line Loop(%d) = {%d,%d};\n', inloop_no, circle1_no, ...
    circle2_no);

fprintf(geo_fid, 'Plane Surface(%d) = {%d, %d};\n', line_no, extloop_no, ...
    inloop_no);
fprintf(geo_fid, 'Physical Surface("extcircle") = {%d};\n', line_no);
line_no = line_no + 1;

fprintf(geo_fid, 'Plane Surface(%d) = {%d};\n', line_no, inloop_no);
fprintf(geo_fid, 'Physical Surface("incircle") = {%d};\n', line_no);
line_no = line_no + 1;

fclose(geo_fid);

% Call Gmsh 
disp('Calling Gmsh. Please wait ...');
call_gmsh( geo_filename);

msh_filename = sprintf('%s.msh', basename);

disp(['Now reading back data from file: ' msh_filename])
[mdl, mat_indices]= gmsh_mk_fwd_model( msh_filename, ...
    'Gmsh based 2 layer circural model' );

gmsh_mk_fwd_model.m

function [fwd_mdl, mat_indices]= ...
             gmsh_mk_fwd_model( vol_filename, name, eprefix, z_contact)
% GMSH_MK_FWD_MODEL: create a fwd_model object from a Gmsh file
% [fwd_mdl, mat_indices]= ...
%      gmsh_mk_fwd_model( vol_filename, centres, ...
%                       name, stim_pattern, z_contact)
%
%  vol_filename:      filename output from Gmsh
%  name:              name for object (if [] use vol_filename)
%  eprefix:           prefix used for names of electrodes
%                     (if [] or omitted use 'electrode-')
%  stim_pattern:      a stimulation pattern structure
%                     empty ([]) if stim_pattern is not available
%  z_contact:         vector or scalar electrode contact impedance
%
%  fwd_mdl:           eidors format fwd_model
%  mat_indices:       cell array of material indices from eidors 

% Gmsh mesher for EIRODS was based on Netgen interface.
% (C) 2009 Bartosz Sawicki. License: GPL version 2 or version 3
% Modified by James Snyder

if isempty(name); 
   name = ['fwd_mdl based on ', vol_filename];
end

if nargin<4
   eprefix = 'electrode-';
end

if isempty(eprefix); 
   eprefix = 'electrode-';
end

if nargin<5
   z_contact=0.01; % singular if z_contact=0
end



stim_pattern=[];
% Model Geometry
[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(vol_filename);
fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
                             stim_pattern, eprefix, z_contact, fc,phys_names);
mat_indices= mk_mat_indices( mat_ind);

% build fwd_model structure
function fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
    stim_pattern, eprefix, z_contact, fc, phys_names)
    mdl.nodes    = vtx;
    mdl.elems    = simp;
    mdl.boundary = srf;
    mdl.boundary_numbers=fc; 
    mdl.gnd_node = 1;
    mdl.np_fwd_solve.perm_sym = '{n}';
    mdl.name = name;
    
    % Model Stimulation
    if ~isempty(stim_pattern)
        mdl.stimulation= stim_pattern;
    end
    
    % Electrodes
    electrodes = find_elec(phys_names,eprefix,z_contact);

    mdl.electrode =     electrodes;
    mdl.solve=      @np_fwd_solve;
    mdl.jacobian=   @np_calc_jacobian;
    mdl.system_mat= @np_calc_system_mat;

    fwd_mdl= eidors_obj('fwd_model', mdl);

% Output cell array of indices into each material type
%   array order is sorted by length of material type
function mat_indices= mk_mat_indices( mat_ind);
    % find length of mat_indices 
    % test example: mat_ind=[10 12 14 14 12 12 14 12];
    sort_mi= sort(mat_ind(:));
    find_mi= find( diff([-1e8;sort_mi]) );
    len_mi = diff([find_mi;length(sort_mi)+1]);
    [jnk,idxs]= sort(-len_mi); %reverse sort
    l_idxs= length(idxs);
    mat_indices= cell(1, l_idxs);
    for i= 1:l_idxs;
        mat_idx_i= sort_mi(find_mi(idxs(i)));
        mat_indices{i}= find(mat_ind == mat_idx_i);
    end

% Assumes that electrodes are numbered starting at 1, with prefix provided
function electrodes = find_elec(phys_names,prefix,z_contact)
phys_elecs = find(arrayfun(@(x)strncmp(x.name,prefix,length(prefix)),phys_names));
for i = 1:length(phys_elecs)
    cur_elec = arrayfun(@(x)strcmp(sprintf('%s%d',prefix,i),x.name),phys_names(phys_elecs));
    electrodes(i).nodes = unique(phys_names(phys_elecs(cur_elec)).nodes(:));
    electrodes(i).z_contact = z_contact;
end

gmsh_read_mesh.m

function[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(filename)
%[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(filename)
% Function to read in a mesh model from Gmsh and saves it in
% five arrays; surface (srf), veritices (vtx), face no. (fc)
% volume (simp) and edges (edg)
%
% srf        = The surfaces indices into vtx
% simp       = The volume indices into vtx
% vtx        = The vertices matrix
% fc         = A one column matrix containing the face numbers
% edg        = Edge segment information
% filename   = Name of file containing NetGen .vol information
% mat_ind    = Material index
% phys_names = Structure of "Physical" entities in the mesh
%              .dim   = dimension
%              .name  = name (string)
%              .tag   = physical tag
%              .nodes = N-x-dim array of indices into vtx

% $Id: gmsh_read_mesh.m 2628 2011-07-11 13:13:50Z aadler $
% (C) 2009 Bartosz Sawicki. Licensed under GPL V2
% Modified by James Snyder <[email protected]>

fid = fopen(filename,'r');
phys_names = [];
while 1
    tline = fgetl(fid);
    if ~ischar(tline); fclose(fid); break; end

    if strcmp(tline,'$Elements')
       elements= parse_elements( fid );
    elseif strcmp(tline,'$Nodes')
       nodes= get_lines_with_nodes( fid );
    elseif strcmp(tline,'$PhysicalNames')
       phys_names= parse_names( fid );
    end
end

if ~isempty(phys_names)
    for i = 1:length(phys_names)
        tmpelements = find(arrayfun(@(x)x.phys_tag==phys_names(i).tag,elements));
        phys_names(i).nodes = cat(1,elements(tmpelements).simp);
    end
end

edg = [];
bc = [];
mat_ind = [];

% Select 2d vs 3d model by checking if Z is all the same
if length( unique( nodes(:,4) ) ) > 1 
    vtx = nodes(:,2:4);
    % Type 2: 3-node triangle
    tri = find(arrayfun(@(x)x.type==2,elements));
    % Type 4: 4-node tetrahedron
    tet = find(arrayfun(@(x)x.type==4,elements));
    simp = cat(1,elements(tet).simp);
    srf = cat(1,elements(tri).simp);
else
    vtx = nodes(:,2:3);
    tri = find(arrayfun(@(x)x.type==2,elements));
    simp = cat(1,elements(tri).simp);
    srf = [];
end

elemtags = cat(1,elements.phys_tag);
fc = elemtags(tri,1);
end


function mat = get_lines_with_nodes( fid )
% Line Format:
% node-number x-coord y-coord z-coord
tline = fgetl(fid);
n_rows = sscanf(tline,'%d');
mat= fscanf(fid,'%f',[4,n_rows])';
end

function names = parse_names( fid )
% Line Format:
% physical-dimension physical-number "physical-name"
tline = fgetl(fid);
n_rows = sscanf(tline,'%d');
names = struct('tag',{},'dim',{},'name',{});
for i = 1:n_rows
    tline = fgetl(fid);
    if exist('OCTAVE_VERSION')
        parts = strsplit(tline,' ');
    else
        parts = regexp(tline,' ','split');
    end
    nsz = size(names,2)+1;
    names(nsz).dim = str2double( parts(1) );
    names(nsz).tag = str2double( parts(2) );
    tname = parts(3);
    names(nsz).name = strrep(tname{1},'"','');
end
end

function elements = parse_elements( fid )
% Line Format:
% elm-number elm-type number-of-tags < tag > ... node-number-list
tline = fgetl(fid);
n_rows = sscanf(tline,'%d');
elements = struct('simp',{},'phys_tag',{},'geom_tag',{});
for i = 1:n_rows
    tline = fgetl(fid);
    n = sscanf(tline, '%d')';
    nsz = size(elements,2)+1;
    elements(nsz).simp = n(n(3) + 4:end);
    % 
    elements(nsz).type = n(2);
    if n(3) > 0 % get tags if they exist
        % By default, first tag is number of parent physical entity
        % second is parent elementary geometrical entity
        % third is number of parent mesh partitions followed by
        % partition ids
        tags = n(4:3+n(3));
        if length(tags) >= 1
            elements(nsz).phys_tag = tags(1);
            if length(tags) >= 2
                elements(nsz).geom_tag = tags(2);
            end
        end
    end
end
end