ng_mk_geometric_models

 NG_MK_GEOMETRIC_MODELS: create geometric mesh models using Netgen. Body
 and electrodes are defined as combinations of solid primitives. The 3-D
 surface intersection of the electrode and body volumes define the
 electrode surfaces.

[fmdl, mat_idx] = ng_mk_geometric_models(body_geometry, electrode_geometry)

 INPUT:
  body_geometry      - Structure whose fields describe body geometries as
                       combinations of unions, intersections and 
                       complements of solid primitives. Cell array
                       should be used when describing several body 
                       geometries.

  electrode_geometry - Structure whose fields describe electrode
                       geometries as combinations of unions, intersections
                       and complements of solid primitives. Cell array
                       should be used when describing several electrode 
                       geometries.

 The following field names are available for geometry descriptions. Each 
 field is followed by the available subfields whose default values if left 
 unspecified are indicated between parentheses. 

 complement_flag:    If true, the desired geometry description is the
                     complement of the geometry description. (false)

 body_of_revolution: A body of revolution is described with the following
                     subfields: axis_point_a ([0; 0; 0]), axis_point_b 
                     ([0; 0; 1]), points (1 1; 1 2; 2 2; 2 1]),
                     segments ([1 2; 2 3; 3 4; 4 1]), complement_flag
                     (false).    

 cone:               A cone is described with the following subfields:
                     bottom_center ([0; 0; 0]), bottom_radius (1),
                     top_center ([0; 0; 1]), top_radius (0.5),
                     complement_flag (false).                    

 cylinder:           A cylinder is described with the following subfields:
                     bottom_center ([0; 0; 0]), top_center ([0; 0; 1]),
                     radius (1), complement_flag (false).

 ellipsoid:          An ellipsoid is described with the following
                     subfields: center ([0; 0; 0]), axis_a ([1; 0; 0]),
                     axis_b ([0; 1; 0]), axis_c ([0; 0; 1]),
                     complement_flag (false).

 elliptic_cylinder:  An elliptic cylinder is described with the following
                     subfields: bottom_center ([0; 0; 0]),
                     top_center ([0; 0; 1]), axis_a ([1; 0; 0]),
                     axis_b ([0; 1; 0]), complement_flag (false). 

 enter_body_flag:    This flag can be used only for electrode geometry
                     descriptions to indicate that the associated
                     electrode solid enters the body solids. It can only
                     be defined at the top level of each geometry
                     description. If this flag is true, it means that the
                     volume of the electrode intersecting with any body is
                     part of the electrode, otherwise it is part of a
                     body. (false)

 half_space          A half-space is described by the following subfields:
                     point ([0; 0; 0]), outward_normal_vector ([0; 0; 1]),
                     complement_flag (false).

 intersection:       This fields indicates to perform the intersection of
                     all subfields. Subfields: complement_flag (false).

 keep_material_flag: This flag can be used only for electrode geometry 
                     descriptions to indicate that the associated
                     electrode material should be kept in the final mesh.
                     It can only be defined at the top level of each
                     geometry description. If true, it means that the
                     volume of the electrode is meshed. Volume elements
                     that are part of the mesh are indicated in mat_idx
                     output argument. (false)

 max_edge_length:    This parameter is used to adjust the maximum size of
                     the element composing the mesh. It can only be used
                     at the top level of each geometry description. (inf)

 name:               This parameter is used to name the geometry
                     description.

 ortho_brick:        An ortho-brick is described by the following
                     subfields: opposite_corner_a ([0; 0; 0]),
                     opposite_corner_b ([1; 1; 1]), complement_flag
                     (false).

 parallelepiped:     A parallelepiped is described by the following
                     subfields: vertex ([0; 0; 0]), vector_a ([1; 0; 0]),
                     vector_b ([0; 1; 0]), vector_c ([0; 0; 1]),
                     complement_flag (false).

 point:              This parameter describes a point. It can only be used
                     at the top level of an electrode geometry
                     description. It must be the only field in the
                     structure. ([])

 sphere:             A sphere is described with the following subfields:
                     center ([0; 0; 0]), radius (1), complement_flag
                     (false).

 union:              This will perform the union of all its subfields.
                     There is an implicit union operator at the top level
                     of the geometry structure. Subfields: complement_flag
                     (false)

 OUTPUT:
  fmdl               - EIDORS forward model object.

  mat_idx            - Vector indicating for each mesh element the indices
                       of materials corresponding as separately defined
                       by input argument body_geometry.

 USAGE EXAMPLES:
 % 3D cylinder with radius 1. One plane of 16 electrodes with radius 0.1
   body_geometry.cylinder = struct;
   n_elect = 16;
   theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
   for i = 1:n_elect
     electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
     electrode_geometry{i}.sphere.radius = 0.1;
   end
   fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);

0001 function [fmdl, mat_idx] = ng_mk_geometric_models(body_geometry, electrode_geometry)
0002 % NG_MK_GEOMETRIC_MODELS: create geometric mesh models using Netgen. Body
0003 % and electrodes are defined as combinations of solid primitives. The 3-D
0004 % surface intersection of the electrode and body volumes define the
0005 % electrode surfaces.
0006 %
0007 %[fmdl, mat_idx] = ng_mk_geometric_models(body_geometry, electrode_geometry)
0008 %
0009 % INPUT:
0010 %  body_geometry      - Structure whose fields describe body geometries as
0011 %                       combinations of unions, intersections and
0012 %                       complements of solid primitives. Cell array
0013 %                       should be used when describing several body
0014 %                       geometries.
0015 %
0016 %  electrode_geometry - Structure whose fields describe electrode
0017 %                       geometries as combinations of unions, intersections
0018 %                       and complements of solid primitives. Cell array
0019 %                       should be used when describing several electrode
0020 %                       geometries.
0021 %
0022 % The following field names are available for geometry descriptions. Each
0023 % field is followed by the available subfields whose default values if left
0024 % unspecified are indicated between parentheses.
0025 %
0026 % complement_flag:    If true, the desired geometry description is the
0027 %                     complement of the geometry description. (false)
0028 %
0029 % body_of_revolution: A body of revolution is described with the following
0030 %                     subfields: axis_point_a ([0; 0; 0]), axis_point_b
0031 %                     ([0; 0; 1]), points (1 1; 1 2; 2 2; 2 1]),
0032 %                     segments ([1 2; 2 3; 3 4; 4 1]), complement_flag
0033 %                     (false).
0034 %
0035 % cone:               A cone is described with the following subfields:
0036 %                     bottom_center ([0; 0; 0]), bottom_radius (1),
0037 %                     top_center ([0; 0; 1]), top_radius (0.5),
0038 %                     complement_flag (false).
0039 %
0040 % cylinder:           A cylinder is described with the following subfields:
0041 %                     bottom_center ([0; 0; 0]), top_center ([0; 0; 1]),
0042 %                     radius (1), complement_flag (false).
0043 %
0044 % ellipsoid:          An ellipsoid is described with the following
0045 %                     subfields: center ([0; 0; 0]), axis_a ([1; 0; 0]),
0046 %                     axis_b ([0; 1; 0]), axis_c ([0; 0; 1]),
0047 %                     complement_flag (false).
0048 %
0049 % elliptic_cylinder:  An elliptic cylinder is described with the following
0050 %                     subfields: bottom_center ([0; 0; 0]),
0051 %                     top_center ([0; 0; 1]), axis_a ([1; 0; 0]),
0052 %                     axis_b ([0; 1; 0]), complement_flag (false).
0053 %
0054 % enter_body_flag:    This flag can be used only for electrode geometry
0055 %                     descriptions to indicate that the associated
0056 %                     electrode solid enters the body solids. It can only
0057 %                     be defined at the top level of each geometry
0058 %                     description. If this flag is true, it means that the
0059 %                     volume of the electrode intersecting with any body is
0060 %                     part of the electrode, otherwise it is part of a
0061 %                     body. (false)
0062 %
0063 % half_space          A half-space is described by the following subfields:
0064 %                     point ([0; 0; 0]), outward_normal_vector ([0; 0; 1]),
0065 %                     complement_flag (false).
0066 %
0067 % intersection:       This fields indicates to perform the intersection of
0068 %                     all subfields. Subfields: complement_flag (false).
0069 %
0070 % keep_material_flag: This flag can be used only for electrode geometry
0071 %                     descriptions to indicate that the associated
0072 %                     electrode material should be kept in the final mesh.
0073 %                     It can only be defined at the top level of each
0074 %                     geometry description. If true, it means that the
0075 %                     volume of the electrode is meshed. Volume elements
0076 %                     that are part of the mesh are indicated in mat_idx
0077 %                     output argument. (false)
0078 %
0079 % max_edge_length:    This parameter is used to adjust the maximum size of
0080 %                     the element composing the mesh. It can only be used
0081 %                     at the top level of each geometry description. (inf)
0082 %
0083 % name:               This parameter is used to name the geometry
0084 %                     description.
0085 %
0086 % ortho_brick:        An ortho-brick is described by the following
0087 %                     subfields: opposite_corner_a ([0; 0; 0]),
0088 %                     opposite_corner_b ([1; 1; 1]), complement_flag
0089 %                     (false).
0090 %
0091 % parallelepiped:     A parallelepiped is described by the following
0092 %                     subfields: vertex ([0; 0; 0]), vector_a ([1; 0; 0]),
0093 %                     vector_b ([0; 1; 0]), vector_c ([0; 0; 1]),
0094 %                     complement_flag (false).
0095 %
0096 % point:              This parameter describes a point. It can only be used
0097 %                     at the top level of an electrode geometry
0098 %                     description. It must be the only field in the
0099 %                     structure. ([])
0100 %
0101 % sphere:             A sphere is described with the following subfields:
0102 %                     center ([0; 0; 0]), radius (1), complement_flag
0103 %                     (false).
0104 %
0105 % union:              This will perform the union of all its subfields.
0106 %                     There is an implicit union operator at the top level
0107 %                     of the geometry structure. Subfields: complement_flag
0108 %                     (false)
0109 %
0110 % OUTPUT:
0111 %  fmdl               - EIDORS forward model object.
0112 %
0113 %  mat_idx            - Vector indicating for each mesh element the indices
0114 %                       of materials corresponding as separately defined
0115 %                       by input argument body_geometry.
0116 %
0117 % USAGE EXAMPLES:
0118 % % 3D cylinder with radius 1. One plane of 16 electrodes with radius 0.1
0119 %   body_geometry.cylinder = struct;
0120 %   n_elect = 16;
0121 %   theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
0122 %   for i = 1:n_elect
0123 %     electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
0124 %     electrode_geometry{i}.sphere.radius = 0.1;
0125 %   end
0126 %   fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
0127 
0128 % (C) Herve Gagnon, 2015. Licenced under GPL v2 or v3
0129 % $Id: ng_mk_geometric_models.m 5836 2018-10-09 07:23:36Z aadler $
0130 
0131 % Check if function is called in UNIT_TEST mode.
0132 if (ischar(body_geometry) && strcmp(body_geometry, 'UNIT_TEST'))
0133     do_unit_test; 
0134     return; 
0135 end
0136 
0137 % Validate function parameters.
0138 if (isempty(body_geometry) || ~isstruct(body_geometry) && ~iscell(body_geometry))
0139    error('Parameter body_geometry must be a structure or a cell.');
0140 end
0141 
0142 % Check if parameter electrode_geometry is specified.
0143 if (nargin < 2 || isempty(electrode_geometry))
0144     electrode_geometry = {};
0145 end
0146 
0147 if (~isstruct(electrode_geometry) && ~iscell(electrode_geometry))
0148    error('Parameter electrode_geometry must be a structure or a cell.');
0149 end
0150 
0151 % If parameters are not of cell type, convert them to cell.
0152 if (~iscell(body_geometry))
0153     body_geometry = {body_geometry};
0154 end
0155 
0156 if (~iscell(electrode_geometry))
0157     electrode_geometry = {electrode_geometry};
0158 end
0159 
0160 % Check if result is already in cache. Otherwise, compute and store in cache.
0161 copt.cache_obj = {body_geometry, electrode_geometry};
0162 copt.fstr = 'ng_mk_geometric_models';
0163 args = {body_geometry, electrode_geometry};
0164 fmdl_mat_idx = eidors_cache(@mk_geometric_models, args, copt);
0165 
0166 % Reformat output arguments.
0167 fmdl    = fmdl_mat_idx{1};
0168 mat_idx = fmdl_mat_idx{2};
0169 
0170 function [fmdl_mat_idx] = mk_geometric_models(body_geometry, electrode_geometry)     
0171     % Find number of subdomains
0172     n_body      = numel(body_geometry);
0173     n_electrode = numel(electrode_geometry);
0174 
0175     % Allocate cell memory.
0176     body_solid_code       = cell(size(body_geometry));
0177     body_extra_code       = cell(size(body_geometry));
0178     body_extra_param      = cell(size(body_geometry));
0179     electrode_solid_code  = cell(size(electrode_geometry));
0180     electrode_extra_code  = cell(size(electrode_geometry));
0181     electrode_extra_param = cell(size(electrode_geometry));
0182     
0183     % Parse geometry for each body subdomain.
0184     for i = 1:n_body
0185         if (isfield(body_geometry{i}, 'point'))
0186             error('Field name "point" is not allowed for body geometry.');
0187         end
0188         if (isfield(body_geometry{i}, 'enter_body_flag'))
0189             error('Field name "enter_body_flag" is not allowed for body geometry.');
0190         end
0191         if (isfield(body_geometry{i}, 'keep_material_flag'))
0192             error('Field name "keep_material_flag" is not allowed for body geometry.');
0193         end
0194         [body_solid_code{i} body_extra_code{i} body_extra_param{i}] = parse_geometry(body_geometry{i});  
0195     end
0196     
0197     % Parse geometry for each body subdomain.
0198     for i = 1:n_electrode
0199         [electrode_extra_code{i} electrode_extra_param{i}] = parse_geometry_point(electrode_geometry{i});
0200         if (isempty(electrode_extra_param{i}.point))
0201             [electrode_solid_code{i} electrode_extra_code{i} electrode_extra_param{i}] = parse_geometry(electrode_geometry{i}); 
0202         end
0203     end
0204 
0205     % Define temporary unique filenames.
0206     fn_prefix = tempname;
0207     geo_fn    = [fn_prefix, '.geo'];
0208     vol_fn    = [fn_prefix, '.vol'];
0209    
0210      % Add body names if unspecified.
0211     for i = 1:numel(body_solid_code)
0212         if (isempty(body_extra_param{i}.name))
0213             body_extra_param{i}.name = sprintf('body%04d', i);
0214         end
0215     end
0216     
0217     % Add electrode names if unspecified.
0218     for i = 1:numel(electrode_solid_code)
0219         if (isempty(electrode_extra_param{i}.name))
0220             electrode_extra_param{i}.name = sprintf('electrode%04d', i);
0221         end
0222     end
0223 
0224     % Write geo file for netgen.
0225     write_geo_file(geo_fn, body_solid_code, electrode_solid_code, body_extra_code, electrode_extra_code, body_extra_param, electrode_extra_param);
0226    
0227     % Call netgen.
0228     call_netgen(geo_fn, vol_fn);
0229  
0230     % Read vol file generated by netgen.
0231     fmdl_mat_idx{1} = read_vol_file(vol_fn, electrode_extra_param);
0232     
0233     % Delete temporary files.
0234     if ~eidors_debug('query','ng_mk_geometric_models:keep_temp_files')
0235        delete(geo_fn);
0236        delete(vol_fn);
0237     end
0238 
0239     % Complete fmdl object.
0240     fmdl_mat_idx{1} = complete_fmdl(fmdl_mat_idx{1}, electrode_extra_param);
0241 
0242     % Assign mat_idx value from fmdl.
0243     fmdl_mat_idx{2} = fmdl_mat_idx{1}.mat_idx;
0244 
0245 function radius = assign_radius(struct, n_structs, struct_name, field_name, default_value)
0246 
0247     radius = default_value;
0248     
0249     for i = 1:n_structs
0250         value = struct(i).(field_name);
0251 
0252         if (~isempty(value))
0253             if (isscalar(value) && isnumeric(value) && isreal(value) && value > 0)
0254                 radius(i) = value;
0255             else 
0256                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0257             end
0258         end
0259     end
0260     
0261 function point = assign_point(struct, n_structs, struct_name, field_name, default_value)
0262         
0263     point = default_value;
0264     
0265     for i = 1:n_structs
0266         value =  struct(i).(field_name);
0267 
0268         if (~isempty(value))
0269             if (numel(value) == 3 && isnumeric(value) && isreal(value))
0270                 point(:, i) = value;
0271             else
0272                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0273             end
0274         end
0275     end
0276  
0277 function point_list = assign_list_of_3D_points(struct, n_structs, struct_name, field_name, default_value)
0278         
0279     point_list = cell(n_structs, 1);
0280     
0281     for i = 1:n_structs
0282         
0283         point_list{i} = default_value;
0284         
0285         value = struct(i).(field_name);
0286 
0287         if (~isempty(value))
0288             if (size(value, 2) == 3 && isnumeric(value) && isreal(value))
0289                 point_list{i} = value;
0290             else
0291                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0292             end
0293         end
0294     end
0295     
0296 function segment_list = assign_list_of_2D_points(struct, n_structs, struct_name, field_name, default_value)
0297         
0298     segment_list = cell(n_structs, 1);
0299     
0300     for i = 1:n_structs
0301         
0302         segment_list{i} = default_value;
0303         
0304         value = struct(i).(field_name);
0305 
0306         if (~isempty(value))
0307             if (size(value, 2) == 2 && isnumeric(value) && isreal(value))
0308                 segment_list{i} = value;
0309             else
0310                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0311             end
0312         end
0313     end 
0314     
0315  function segment_list = assign_list_of_2D_or_3D_points(struct, n_structs, struct_name, field_name, default_value)
0316         
0317     segment_list = cell(n_structs, 1);
0318     
0319     for i = 1:n_structs
0320         
0321         segment_list{i} = default_value;
0322         
0323         value = struct(i).(field_name);
0324 
0325         if (~isempty(value))
0326             if ((size(value, 2) == 2 || size(value, 2) == 3) && isnumeric(value) && isreal(value))
0327                 segment_list{i} = value;
0328             else
0329                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0330             end
0331         end
0332     end    
0333     
0334 function flag = assign_flag(struct, n_structs, struct_name, field_name, default_value)
0335 
0336     flag = default_value;
0337 
0338     for i = 1:n_structs
0339         value = struct(i).(field_name);
0340  
0341         if (~isempty(value))
0342             if (isscalar(value) && (islogical(value) || (isnumeric(value) && ...
0343                                      isreal(value) && (value == 0 || value == 1))))
0344                 flag(i) = value;
0345             else
0346 value
0347                 error('%s(%d).%s value is not valid.', struct_name, i, field_name);
0348             end
0349         end
0350     end
0351 
0352 function [extra_code extra_param] = parse_geometry_point(geometry)
0353 
0354     % Initialize extra param values to default values.
0355     extra_code                     = '';
0356     extra_param.point              = [];
0357     extra_param.max_edge_length    = inf;
0358     extra_param.enter_body_flag    = false;
0359     extra_param.keep_material_flag = false;
0360     extra_param.name               = '';
0361 
0362     % Check if geometry is a point, return otherwise.
0363     if (isfield(geometry, 'point'))
0364         % Check if a single point is defined.
0365         if (numel(geometry) ~= 1)
0366             error('Field name "point" must define only a single point.');
0367         end
0368         
0369         % Get structure field names.
0370         field_names = fieldnames(geometry);
0371         n_fields = numel(field_names);
0372         
0373         % Check if it is the only field names.
0374         if (n_fields ~= 1)
0375             if (isfield(geometry, 'name'))
0376                 extra_param.name = geometry.name;
0377             else
0378                 error('Field name "point" must be used as a single field.');
0379             end
0380         end
0381         
0382         % Check point is 3D.
0383         if (numel(geometry.point) ~= 3)
0384             error('geometry.point value is not valid.');
0385         end
0386         
0387         % Set returning values.
0388         extra_param.point = geometry.point(:);
0389         
0390         extra_code = sprintf('point(%g, %g, %g);\n', extra_param.point);
0391     end
0392     
0393 function [geo_code extra_code extra_param] = parse_geometry(geometry, field_operator_string, element_operator_string)
0394 
0395     % Extra code is initialized empty;
0396     extra_code = '';
0397     
0398     % Initialize extra param values to default values.
0399     extra_param.max_edge_length     = inf;
0400     extra_param.enter_body_flag     = false;
0401     extra_param.keep_material_flag  = false;
0402     extra_param.name                = '';
0403     
0404     % If called from top_level, operators default to or.
0405     if (nargin == 1)
0406         top_level_flag = 1;
0407         field_operator_string = ' or ';
0408         element_operator_string = ' or ';
0409     else
0410         top_level_flag = 0;
0411     end
0412 
0413     eidors_msg('@@@ called with (field: "%s", element: "%s").', field_operator_string, element_operator_string, 4);
0414     
0415     % Validate that geometry is a non-empty structure.
0416     if (~isstruct(geometry) || isempty(geometry))
0417         error('Parameter geometry must be a valid structure.');        
0418     else
0419         % Get number of geometries.
0420         n_geometries = numel(geometry);
0421         
0422         % Get structure field names.
0423         field_names = fieldnames(geometry);
0424         n_fields = numel(field_names);
0425 
0426         % Recursively parse all geometry fields.
0427         geo_code = '(';
0428         for i = 1:n_geometries
0429             % complement_flag field has to be processed first.
0430             if (isfield(geometry(i), 'complement_flag') && ~isempty(geometry(i).complement_flag) && geometry(i).complement_flag)
0431                  geo_code = [geo_code '(not('];
0432             else
0433                  geo_code = [geo_code '('];
0434             end
0435             % Process name field.
0436             if (isfield(geometry(i), 'name'))
0437                 if (~top_level_flag)
0438                     error('Field "name" can only be specified at the top level of the geometry description');
0439                 end
0440                 extra_param.name = geometry(i).name;
0441 
0442                 if (isempty(extra_param.name) || ~ischar(extra_param.name))
0443                     error('name value is not valid.');
0444                 end
0445             end
0446            % Process max_edge_length field.
0447             if (isfield(geometry(i), 'max_edge_length'))
0448                 if (~top_level_flag)
0449                     error('Field "max_edge_length" can only be specified at the top level of the geometry description');
0450                 end
0451                 extra_param.max_edge_length = geometry(i).max_edge_length;
0452 
0453                 if (isempty(extra_param.max_edge_length) || ~isscalar(extra_param.max_edge_length) || ~isnumeric(extra_param.max_edge_length) || ~isreal(extra_param.max_edge_length) || extra_param.max_edge_length <= 0)
0454                     error('max_edge_length value is not valid.');
0455                 end
0456             end
0457             % Process enter_body_flag field.
0458             if (isfield(geometry(i), 'enter_body_flag'))
0459                 if (~top_level_flag)
0460                     error('Field "enter_body_flag" can only be specified at the top level of the geometry description');
0461                 end
0462                 extra_param.enter_body_flag = geometry(i).enter_body_flag;
0463                 
0464                 if (isempty(extra_param.enter_body_flag) || ~isscalar(extra_param.enter_body_flag) || (~islogical(extra_param.enter_body_flag) && ...
0465                         (~isnumeric(extra_param.enter_body_flag) || ~isreal(extra_param.enter_body_flag) || (extra_param.enter_body_flag ~= 0 && extra_param.enter_body_flag ~= 1))))
0466                     error('Field "enter_body_flag value" is not valid.');
0467                 end
0468             end
0469             % Process keep_material_flag field.
0470             if (isfield(geometry(i), 'keep_material_flag'))
0471                 if (~top_level_flag)
0472                     error('Field "keep_material_flag" can only be specified at the top level of the geometry description');
0473                 end
0474                 extra_param.keep_material_flag = geometry(i).keep_material_flag;
0475                 
0476                 if (isempty(extra_param.keep_material_flag) || ~isscalar(extra_param.keep_material_flag) || (~islogical(extra_param.keep_material_flag) && ...
0477                         (~isnumeric(extra_param.keep_material_flag) || ~isreal(extra_param.keep_material_flag) || (extra_param.keep_material_flag ~= 0 && extra_param.keep_material_flag ~= 1))))
0478                     error('Field "keep_material_flag" value is not valid.');
0479                 end
0480             end  
0481             first_internal_term = 1;
0482             for j = 1:n_fields
0483                 if (~isempty(geometry(i).(field_names{j})) && ~strcmp(field_names{j}, 'complement_flag') && ~strcmp(field_names{j}, 'name') && ~strcmp(field_names{j}, 'max_edge_length') && ~strcmp(field_names{j}, 'keep_material_flag') && ~strcmp(field_names{j}, 'enter_body_flag'))
0484                     if (first_internal_term)
0485                         first_internal_term = 0;
0486                     else
0487                         geo_code = [geo_code field_operator_string];
0488                     end
0489                     switch (field_names{j})
0490                         case 'body_of_extrusion'
0491                             [geo_code_temp extra_code_temp] = ...
0492                                         parse_geometry_body_of_extrusion(geometry(i).(field_names{j}), field_operator_string);
0493                             geo_code   = [geo_code geo_code_temp];        
0494                             extra_code = [extra_code extra_code_temp];
0495                         case 'body_of_revolution'
0496                             [geo_code_temp extra_code_temp] = ...
0497                                         parse_geometry_body_of_revolution(geometry(i).(field_names{j}), field_operator_string);
0498                             geo_code   = [geo_code geo_code_temp];        
0499                             extra_code = [extra_code extra_code_temp];
0500                         case 'cone'
0501                             geo_code = [geo_code ...
0502                                         parse_geometry_cone(geometry(i).(field_names{j}), field_operator_string)];
0503                         case 'cylinder'
0504                             [geo_code_temp, extra_code_temp] = ...
0505                                         parse_geometry_cylinder(geometry(i).(field_names{j}), field_operator_string);
0506                             geo_code   = [geo_code geo_code_temp];        
0507                             extra_code = [extra_code extra_code_temp];
0508                         case 'ellipsoid'
0509                             geo_code = [geo_code ...
0510                                         parse_geometry_ellipsoid(geometry(i).(field_names{j}), field_operator_string)];
0511                         case 'elliptic_cylinder'
0512                             geo_code = [geo_code ...
0513                                         parse_geometry_elliptic_cylinder(geometry(i).(field_names{j}), field_operator_string)];
0514                         case 'half_space'
0515                             geo_code = [geo_code ...
0516                                         parse_geometry_half_space(geometry(i).(field_names{j}), field_operator_string)];    
0517                         case 'ortho_brick'
0518                             geo_code = [geo_code ...
0519                                         parse_geometry_ortho_brick(geometry(i).(field_names{j}), field_operator_string)];
0520                         case 'parallelepiped'
0521                             geo_code = [geo_code ...
0522                                         parse_geometry_parallelepiped(geometry(i).(field_names{j}), field_operator_string)];       
0523                         case 'sphere'
0524                             geo_code = [geo_code ...
0525                                         parse_geometry_sphere(geometry(i).(field_names{j}), field_operator_string)];
0526                         case 'intersection'
0527                             [geo_code_temp extra_code_temp] = ...
0528                                         parse_geometry(geometry(i).(field_names{j}), ' and ', field_operator_string);
0529                             geo_code   = [geo_code geo_code_temp];        
0530                             extra_code = [extra_code extra_code_temp];
0531                         case 'union'
0532                             [geo_code_temp extra_code_temp] = ...
0533                                         parse_geometry(geometry(i).(field_names{j}), ' or ', field_operator_string);
0534                             geo_code   = [geo_code geo_code_temp];        
0535                             extra_code = [extra_code extra_code_temp];
0536                         otherwise
0537                             error(['Field name "%s" is not valid for a geometry.\nAvailable field names for a geometry are: '...
0538                                    'complement_flag, intersection, union, body_of_extrusion, body_of_revolution, cone, cylinder, ellipsoid, elliptic_cylinder, half_space, ortho_brick, parallelepiped, point, sphere, keep_material_flag, enter_body_flag, name, and max_edge_length.'], field_names{j});
0539                     end
0540                 end
0541             end
0542             if (isfield(geometry(i), 'complement_flag') && ~isempty(geometry(i).complement_flag) && geometry(i).complement_flag)
0543                 geo_code = [geo_code '))'];
0544             else
0545                 geo_code = [geo_code ')'];  
0546             end
0547            
0548             if (i < n_geometries)
0549                 geo_code = [geo_code element_operator_string];         
0550             end           
0551         end
0552         geo_code = [geo_code ')'];  
0553     end
0554     
0555     eidors_msg('@@@ returned with (field: "%s", element: "%s").', field_operator_string, element_operator_string, 4);
0556  
0557 function [geo_code extra_code] = parse_geometry_body_of_extrusion(body_of_extrusion, operator_string)
0558 
0559     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0560 
0561     if (~isstruct(body_of_extrusion) || isempty(body_of_extrusion))
0562         error('Parameter body_of_extrusion must be a valid structure.');        
0563     else
0564         % Get number of body_of_extrusion.
0565         n_body_of_extrusions = numel(body_of_extrusion);
0566         
0567         % Get structure field names.
0568         field_names = fieldnames(body_of_extrusion);
0569         n_fields = numel(field_names);
0570         
0571         % Assign default values.
0572         vector_d            = [0; 1; 0]*ones(1, n_body_of_extrusions);
0573         profile_points{1}   = [1 1; 1 2; 2 2; 2 1];
0574         profile_segments{1} = [1 2; 2 3; 3 4; 4 1];
0575         path_points{1}      = [0 0 0; 0 0 1; 0 0 2; 0 0 3];
0576         path_segments{1}    = [1 2; 2 3; 3 4];
0577         complement_flag     = false(1, n_body_of_extrusions);
0578         
0579         % Parse all structure fields.
0580         for i = 1:n_fields
0581             switch (field_names{i})
0582                 case 'vector_d'
0583                     vector_d = assign_point(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, vector_d);
0584                 case 'profile_points'
0585                     profile_points = assign_list_of_2D_points(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, profile_points{1});
0586                 case 'profile_segments'
0587                     profile_segments = assign_list_of_2D_or_3D_points(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, profile_segments{1});
0588                 case 'path_points'
0589                     path_points = assign_list_of_3D_points(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, path_points{1});
0590                 case 'path_segments'
0591                     path_segments = assign_list_of_2D_or_3D_points(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, path_segments{1});
0592                 case 'complement_flag'
0593                     complement_flag = assign_flag(body_of_extrusion, n_body_of_extrusions, 'body_of_extrusion', field_names{i}, complement_flag);
0594                 otherwise
0595                     error(['Field name "%s" is not allowed for a body_of_extrusion!\nAllowed field names for a body_of_extrusion are: ' ...
0596                            'path_points, path_segments, profile_points, profile_segments, vector_d, and complement_flag.'], field_names{i});
0597             end
0598         end
0599         
0600         % Start geo code with an opening parenthesis.
0601         geo_code = '(';
0602         extra_code = '';
0603         
0604         % Add geo code for each body_of_extrusion.
0605         for i = 1:n_body_of_extrusions
0606             
0607             for j = 1:size(path_segments{i}, 1)
0608                 if (dot(vector_d(:,i), path_points{i}(path_segments{i}(j, 1), :) - path_points{i}(path_segments{i}(j, end), :)) ~= 0)
0609                     error('vector_d and path must be perpendicular for a body of extrusion.');
0610                 end
0611             end
0612             
0613             n_points = size(profile_points{i}, 1);
0614             n_segments = size(profile_segments{i}, 1);
0615             
0616             if (size(profile_segments{i}, 2) == 2)
0617                 extra_code = [extra_code sprintf('curve2d Extrusion2DProfileCurve%d = (%g ; ', i, n_points) ...
0618                                          sprintf('%g, %g ; ', profile_points{i}') ...
0619                                          sprintf(' %g ', n_segments) ...
0620                                          sprintf('; 2, %g, %g ', profile_segments{i}') ...
0621                                          sprintf(');\n\n')];
0622             else
0623                 extra_code = [extra_code sprintf('curve2d Extrusion2DProfileCurve%d = (%g ; ', i, n_points) ...
0624                                          sprintf('%g, %g ; ', profile_points{i}') ...
0625                                          sprintf(' %g ', n_segments) ...
0626                                          sprintf('; 3, %g, %g, %g ', profile_segments{i}') ...
0627                                          sprintf(');\n\n')];
0628             end
0629   
0630             n_points = size(path_points{i}, 1);
0631             n_segments = size(path_segments{i}, 1);
0632             
0633             if (size(path_segments{i}, 2) == 2)
0634                 extra_code = [extra_code sprintf('curve3d Extrusion3DPathCurve%d = (%g ; ', i, n_points) ...
0635                                          sprintf('%g, %g, %g ; ', path_points{i}') ...
0636                                          sprintf(' %g ', n_segments) ...
0637                                          sprintf('; 2, %g, %g ', path_segments{i}') ...
0638                                          sprintf(');\n\n')];
0639             else
0640                  extra_code = [extra_code sprintf('curve3d Extrusion3DPathCurve%d = (%g ; ', i, n_points) ...
0641                                          sprintf('%g, %g, %g ; ', path_points{i}') ...
0642                                          sprintf(' %g ', n_segments) ...
0643                                          sprintf('; 3, %g, %g, %g ', path_segments{i}') ...
0644                                          sprintf(');\n\n')];               
0645             end
0646                                        
0647             if (complement_flag(i))
0648                 geo_code = [geo_code 'not '];
0649             end
0650             
0651             % Check if path is closed.
0652             if (path_segments{i}(end) == path_segments{i}(1))
0653                 geo_code = [geo_code sprintf('extrusion(Extrusion3DPathCurve%d ; Extrusion2DProfileCurve%d ; %g, %g, %g)', i, i, vector_d(:,i))];
0654             else
0655                 %error('Unclosed path are not yet supported for body of extrusion!');
0656                 %for j = 1: [1 1; 1 2; 2 2; 2 1];
0657                 %polyhedron1 = 'polyhedron(-1, 1, 0; -1, 2, 0; -2, 2, 0; -2, 1, 0; -1.5, 1.5, 0; -1.5, 1.5, 1 ;; 2, 1, 5; 3, 2, 5; 4, 3, 5; 1, 4, 5; 1, 2, 6; 2, 3, 6; 3, 4, 6; 4, 1, 6)';
0658                 %polyhedron2 = 'polyhedron(-1, 1, 3; -1, 2, 3; -2, 2, 3; -2, 1, 3; -1.5, 1.5, 3; -1.5, 1.5, 2 ;; 1, 2, 5; 2, 3, 5; 3, 4, 5; 4, 1, 5; 2, 1, 6; 3, 2, 6; 4, 3, 6; 1, 4, 6)';
0659                 %polyhedron1 = 'plane(0, 0, 3; 0, 0, 1)';
0660                 %polyhedron2 = 'plane(0, 0, 0; 0, 0, -1)';
0661                 first_point  = path_points{i}(path_segments{i}(1, 1), :);
0662                 first_vector = first_point - path_points{i}(path_segments{i}(1, end), :);
0663                 last_point   = path_points{i}(path_segments{i}(end, end), :);
0664                 last_vector  = last_point - path_points{i}(path_segments{i}(end, 1), :);
0665                 geo_code = [geo_code sprintf('(extrusion(Extrusion3DPathCurve%d ; Extrusion2DProfileCurve%d ; %g, %g, %g) and plane(%g, %g, %g; %g, %g, %g) and plane(%g, %g, %g; %g, %g, %g))', ...
0666                             i, i, vector_d(:,i), first_point, first_vector, last_point, last_vector)];
0667                 %geo_code = [geo_code sprintf('(%s or %s)', polyhedron1, polyhedron2)];
0668             end
0669 
0670             if (i < n_body_of_extrusions)
0671                 geo_code = [geo_code operator_string];
0672             else
0673                 geo_code = [geo_code ')'];             
0674             end
0675         end
0676     end
0677     
0678     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
0679 
0680 function [geo_code extra_code] = parse_geometry_body_of_revolution(body_of_revolution, operator_string)
0681 
0682     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0683 
0684     if (~isstruct(body_of_revolution) || isempty(body_of_revolution))
0685         error('Parameter body_of_revolution must be a valid structure.');        
0686     else
0687         % Get number of body_of_revolution.
0688         n_body_of_revolutions = numel(body_of_revolution);
0689         
0690         % Get structure field names.
0691         field_names = fieldnames(body_of_revolution);
0692         n_fields = numel(field_names);
0693         
0694         % Assign default values.
0695         axis_point_a   = [0;0;0]*ones(1, n_body_of_revolutions);
0696         axis_point_b   = [0;0;1]*ones(1, n_body_of_revolutions);
0697         points{1}      = [1 1; 1 2; 2 2; 2 1];
0698         segments{1}    = [1 2; 2 3; 3 4; 4 1];
0699         complement_flag = false(1, n_body_of_revolutions);
0700         
0701         % Parse all structure fields.
0702         for i = 1:n_fields
0703             switch (field_names{i})
0704                 case 'axis_point_a'
0705                     axis_point_a = assign_point(body_of_revolution, n_body_of_revolutions, 'body_of_revolution', field_names{i}, axis_point_a);
0706                 case 'axis_point_b'
0707                     axis_point_b = assign_point(body_of_revolution, n_body_of_revolutions, 'body_of_revolution', field_names{i}, axis_point_b);
0708                 case 'points'
0709                     points = assign_list_of_2D_points(body_of_revolution, n_body_of_revolutions, 'body_of_revolution', field_names{i}, points{1});
0710                 case 'segments'
0711                     segments = assign_list_of_2D_or_3D_points(body_of_revolution, n_body_of_revolutions, 'body_of_revolution', field_names{i}, segments{1});
0712                 case 'complement_flag'
0713                     complement_flag = assign_flag(body_of_revolution, n_body_of_revolutions, 'body_of_revolution', field_names{i}, complement_flag);
0714                 otherwise
0715                     error(['Field name ''%s'' is not valid for a body_of_revolution.\Available field names for a body_of_revolution are: ' ...
0716                            'axis_point_a, axis_point_b, points, segments, and complement_flag.'], field_names{i});
0717             end
0718         end
0719         
0720         % Start geo code with an opening parenthesis.
0721         geo_code = '(';
0722         extra_code = '';
0723         
0724         % Add geo code for each body_of_revolution.
0725         for i = 1:n_body_of_revolutions
0726             
0727             n_points = size(points{i}, 1);
0728             n_segments = size(segments{i}, 1);
0729             
0730             if (size(segments{i}, 2) == 2)
0731                 extra_code = [extra_code sprintf('curve2d Revolution2DCurve%d = (%g ; ', i, n_points) ...
0732                                                sprintf('%g, %g ; ', points{i}') ...
0733                                                sprintf(' %g ', n_segments) ...
0734                                                sprintf('; 2, %g, %g ', segments{i}') ...
0735                                                sprintf(');\n\n')];
0736             else
0737                 extra_code = [extra_code sprintf('curve2d Revolution2DCurve%d = (%g ; ', i, n_points) ...
0738                                                sprintf('%g, %g ; ', points{i}') ...
0739                                                sprintf(' %g ', n_segments) ...
0740                                                sprintf('; 3, %g, %g, %g ', segments{i}') ...
0741                                                sprintf(');\n\n')];
0742             end
0743             
0744             if (complement_flag(i))
0745                 geo_code = [geo_code 'not '];
0746             end
0747             
0748             geo_code = [geo_code sprintf('revolution(%g, %g, %g ; %g, %g, %g ; Revolution2DCurve%d)', ...
0749                         axis_point_a(:, i), axis_point_b(:, i), i)];
0750 
0751             if (i < n_body_of_revolutions)
0752                 geo_code = [geo_code operator_string];
0753             else
0754                 geo_code = [geo_code ')'];             
0755             end
0756         end
0757     end
0758     
0759     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
0760     
0761 function geo_code = parse_geometry_cone(cone, operator_string)
0762 
0763     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0764 
0765     if (~isstruct(cone) || isempty(cone))
0766         error('Parameter cone must be a valid structure.');        
0767     else
0768         % Get number of cones.
0769         n_cones = numel(cone);
0770         
0771         % Get structure field names.
0772         field_names = fieldnames(cone);
0773         n_fields = numel(field_names);
0774         
0775         % Assign default values.
0776         top_radius      = 0.5*ones(1, n_cones);
0777         bottom_radius   = ones(1, n_cones);
0778         top_center      = [0;0;1]*ones(1, n_cones);
0779         bottom_center   = [0;0;0]*ones(1, n_cones);
0780         complement_flag = false(1, n_cones);
0781         
0782         % Parse all structure fields.
0783         for i = 1:n_fields
0784             switch (field_names{i})
0785                 case 'top_radius'
0786                     top_radius = assign_radius(cone, n_cones, 'cone', field_names{i}, top_radius);
0787                 case 'bottom_radius'
0788                     bottom_radius = assign_radius(cone, n_cones, 'cone', field_names{i}, bottom_radius);
0789                 case {'top_center', 'top_centre'}
0790                     top_center = assign_point(cone, n_cones, 'cone', field_names{i}, top_center);
0791                 case {'bottom_center', 'bottom_centre'}
0792                     bottom_center = assign_point(cone, n_cones, 'cone', field_names{i}, bottom_center);
0793                 case 'complement_flag'
0794                     complement_flag = assign_flag(cone, n_cones, 'cone', field_names{i}, complement_flag);
0795                 otherwise
0796                     error(['Field name ''%s'' is not valid for a cone.\Available field names for a cone are: ' ...
0797                            'bottom_center, bottom_radius, top_center, top_radius, and complement_flag.'], field_names{i});
0798             end
0799         end
0800         
0801         % Start geo code with an opening parenthesis.
0802         geo_code = '(';
0803 
0804         % Add geo code for each cone.
0805         for i = 1:n_cones
0806             if (complement_flag(i))
0807                 geo_code = [geo_code 'not '];
0808             end
0809             
0810             n_vector = top_center(:,i) - bottom_center(:,i); 
0811             
0812             geo_code = [geo_code sprintf('(cone(%g, %g, %g ; %g ; %g, %g, %g ; %g) and plane(%g, %g, %g ; %g, %g, %g) and plane(%g, %g, %g ; %g, %g, %g))', ...
0813                         bottom_center(:,i), bottom_radius(i), top_center(:,i), top_radius(i), bottom_center(:,i), -n_vector, top_center(:,i), n_vector)];
0814 
0815             if (i < n_cones)
0816                 geo_code = [geo_code operator_string];
0817             else
0818                 geo_code = [geo_code ')'];             
0819             end
0820         end
0821     end
0822     
0823     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
0824     
0825 function [geo_code,extra_code] = parse_geometry_cylinder(cylinder, operator_string)
0826 
0827     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0828 
0829     if (~isstruct(cylinder) || isempty(cylinder))
0830         error('Parameter cylinder must be a valid structure.');        
0831     else
0832         % Get number of cylinders.
0833         n_cylinders = numel(cylinder);
0834         
0835         % Get structure field names.
0836         field_names = fieldnames(cylinder);
0837         n_fields = numel(field_names);
0838         
0839         % Assign default values.
0840         radius          = ones(1, n_cylinders);
0841         top_center      = [0;0;1]*ones(1, n_cylinders);
0842         bottom_center   = [0;0;0]*ones(1, n_cylinders);
0843         complement_flag = false(1, n_cylinders);
0844         max_edge_length = inf(1,n_cylinders);;
0845         
0846         % Parse all structure fields.
0847         for i = 1:n_fields
0848             switch (field_names{i})
0849                 case 'radius'
0850                     radius = assign_radius(cylinder, n_cylinders, 'cylinder', field_names{i}, radius);     
0851                 case {'top_center', 'top_centre'}
0852                     top_center = assign_point(cylinder, n_cylinders, 'cylinder', field_names{i}, top_center);
0853                 case {'bottom_center', 'bottom_centre'}
0854                     bottom_center = assign_point(cylinder, n_cylinders, 'cylinder', field_names{i}, bottom_center);
0855                 case 'complement_flag'
0856                     complement_flag = assign_flag(cylinder, n_cylinders, 'cylinder', field_names{i}, complement_flag);
0857                 case 'max_edge_length'
0858                     for ii=1:length(cylinder);
0859                        mel = cylinder(ii).max_edge_length;
0860                        if isempty(mel)
0861                           max_edge_length(ii) = inf;
0862                        else
0863                           max_edge_length(ii) = mel;
0864                        end
0865                     end
0866                 otherwise
0867                     error(['Field name ''%s'' is not valid for a cylinder!\nAvailable field names for a cylinder are: '...
0868                            'bottom_center, top_center, radius, max_edge_length and complement_flag.'], field_names{i});
0869             end
0870         end
0871         
0872         % Start geo code with an opening parenthesis.
0873         extra_code = '';
0874         geo_code = '(';
0875 
0876         % Add geo code for each cylinder.
0877         for i = 1:n_cylinders
0878             if (complement_flag(i))
0879                 geo_code = [geo_code 'not '];
0880             end
0881             
0882             n_vector = top_center(:,i) - bottom_center(:,i); 
0883             
0884             if isinf(max_edge_length(i))
0885                geo_code = [geo_code sprintf('(cylinder(%g, %g, %g ; %g, %g, %g ; %g) and plane(%g, %g, %g ; %g, %g, %g) and plane(%g, %g, %g ; %g, %g, %g))', ...
0886                            bottom_center(:,i), top_center(:,i), radius(i), bottom_center(:,i), -n_vector, top_center(:,i), n_vector)];
0887             else
0888                tmpvar   = sprintf('V%015d',round(1e15*rand));
0889                extra_code = [extra_code, sprintf( ...
0890                   'solid %s = cylinder(%g, %g, %g ; %g, %g, %g ; %g) -maxh=%g;\n', ...
0891                   tmpvar, bottom_center(:,i), top_center(:,i), radius(i), max_edge_length(i))];
0892                geo_code = [geo_code sprintf('((%s) and plane(%g, %g, %g ; %g, %g, %g) and plane(%g, %g, %g ; %g, %g, %g))', ...
0893                            tmpvar, bottom_center(:,i), -n_vector, top_center(:,i), n_vector)];
0894             end
0895                     
0896             if (i < n_cylinders)
0897                 geo_code = [geo_code operator_string];
0898             else
0899                 geo_code = [geo_code ')'];             
0900             end
0901         end
0902     end
0903     
0904     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
0905     
0906 function geo_code = parse_geometry_ellipsoid(ellipsoid, operator_string)
0907 
0908     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0909 
0910     if (~isstruct(ellipsoid) || isempty(ellipsoid))
0911         error('Parameter ellipsoid must be a valid structure.');        
0912     else
0913         % Get number of ellipsoids.
0914         n_ellipsoids = numel(ellipsoid);
0915         
0916         % Get structure field names.
0917         field_names = fieldnames(ellipsoid);
0918         n_fields = numel(field_names);
0919   
0920         % Assign default values.
0921         axis_a          = [1;0;0]*ones(1, n_ellipsoids);
0922         axis_b          = [0;1;0]*ones(1, n_ellipsoids);
0923         axis_c          = [0;0;1]*ones(1, n_ellipsoids);
0924         center          = [0;0;0]*ones(1, n_ellipsoids);
0925         complement_flag = false(1, n_ellipsoids);
0926         
0927         % Parse all structure fields.
0928         for i = 1:n_fields
0929             switch (field_names{i})   
0930                 case {'center', 'centre'}
0931                     center = assign_point(ellipsoid, n_ellipsoids, 'ellipsoid', field_names{i}, center);
0932                 case 'axis_a'
0933                     axis_a = assign_point(ellipsoid, n_ellipsoids, 'ellipsoid', field_names{i}, axis_a);
0934                 case 'axis_b'
0935                     axis_b = assign_point(ellipsoid, n_ellipsoids, 'ellipsoid', field_names{i}, axis_b);
0936                 case 'axis_c'
0937                     axis_c = assign_point(ellipsoid, n_ellipsoids, 'ellipsoid', field_names{i}, axis_c);
0938                 case 'complement_flag'
0939                     complement_flag = assign_flag(ellipsoid, n_ellipsoids, 'ellipsoid', field_names{i}, complement_flag);
0940                 otherwise
0941                      error(['Field name ''%s'' is not valid for an ellipsoid!\nAvailable field names for an ellipsoid are: '...
0942                            'center, axis_a, axis_b, axis_c, and complement_flag.'], field_names{i});
0943             end
0944         end
0945         
0946         % Start geo code with an opening parenthesis.
0947         geo_code = '(';
0948 
0949         % Add geo code for each ellipsoid.
0950         for i = 1:n_ellipsoids
0951             if (dot(axis_a(:,i), axis_b(:,i)) ~= 0)
0952                 error('axis_a and axis_b have to be perpendicular for an ellipsoid.');
0953             elseif (dot(axis_a(:,i), axis_c(:,i)) ~= 0)
0954                 error('axis_a and axis_c have to be perpendicular for an ellipsoid.');
0955             elseif (dot(axis_b(:,i), axis_c(:,i)) ~= 0)
0956                 error('axis_b and axis_c have to be perpendicular for an ellipsoid.');
0957             end
0958             
0959             if (complement_flag(i))
0960                 geo_code = [geo_code 'not '];
0961             end
0962                 
0963             geo_code = [geo_code sprintf('ellipsoid(%g, %g, %g ; %g, %g, %g ; %g, %g, %g ; %g, %g, %g)', ...
0964                         center(:,i), axis_a(:, i), axis_b(:,i) , axis_c(:,i))];
0965                     
0966             if (i < n_ellipsoids)
0967                 geo_code = [geo_code operator_string];
0968             else
0969                 geo_code = [geo_code ')'];             
0970             end
0971         end
0972     end
0973     
0974     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
0975       
0976 function geo_code = parse_geometry_elliptic_cylinder(elliptic_cylinder, operator_string)
0977 
0978     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
0979 
0980     if (~isstruct(elliptic_cylinder) || isempty(elliptic_cylinder))
0981         error('Parameter elliptic_cylinder must be a valid structure.');        
0982     else
0983         % Get number of elliptic_cylinders.
0984         n_elliptic_cylinders = numel(elliptic_cylinder);
0985         
0986         % Get structure field names.
0987         field_names = fieldnames(elliptic_cylinder);
0988         n_fields = numel(field_names);
0989         
0990         % Assign default values.
0991         top_center      = [0;0;1]*ones(1, n_elliptic_cylinders);
0992         bottom_center   = [0;0;0]*ones(1, n_elliptic_cylinders);
0993         axis_a          = [1;0;0]*ones(1, n_elliptic_cylinders);
0994         axis_b          = [0;1;0]*ones(1, n_elliptic_cylinders);
0995         complement_flag = false(1, n_elliptic_cylinders);
0996         
0997         % Parse all structure fields.
0998         for i = 1:n_fields
0999             switch (field_names{i})   
1000                 case {'top_center', 'top_centre'}
1001                     top_center = assign_point(elliptic_cylinder, n_elliptic_cylinders, 'elliptic_cylinder', field_names{i}, top_center);
1002                 case {'bottom_center', 'bottom_centre'}
1003                     bottom_center = assign_point(elliptic_cylinder, n_elliptic_cylinders, 'elliptic_cylinder', field_names{i}, bottom_center);
1004                 case 'axis_a'
1005                     axis_a = assign_point(elliptic_cylinder, n_elliptic_cylinders, 'elliptic_cylinder', field_names{i}, axis_a);
1006                 case 'axis_b'
1007                     axis_b = assign_point(elliptic_cylinder, n_elliptic_cylinders, 'elliptic_cylinder', field_names{i}, axis_b);
1008                 case 'complement_flag'
1009                     complement_flag = assign_flag(elliptic_cylinder, n_elliptic_cylinders, 'elliptic_cylinder', field_names{i}, complement_flag);
1010                 otherwise
1011                     error(['Field name ''%s'' is not valid for an elliptic cylinder!\nAvailable field names for an elliptic cylinder are: '...
1012                            'bottom_center, top_center, axis_a, axis_b, and complement_flag.'], field_names{i});
1013             end
1014         end
1015         
1016         % Start geo code with an opening parenthesis.
1017         geo_code = '(';
1018 
1019         % Add geo code for each cylinder.
1020         for i = 1:n_elliptic_cylinders
1021             if (complement_flag(i))
1022                 geo_code = [geo_code 'not '];
1023             end
1024             
1025             central_axis = top_center(:,i) - bottom_center(:,i);
1026             
1027             if (dot(axis_a(:,i), axis_b(:,i)) ~= 0)
1028                 error('axis_a and axis_b have to be perpendicular for an elliptic cylinder.');
1029             elseif (dot(axis_a(:,i), central_axis(:,i)) ~= 0)
1030                 error('axis_a and the central axis have to be perpendicular for an elliptic cylinder.');
1031             elseif (dot(axis_b(:,i), central_axis(:,i)) ~= 0)
1032                 error('axis_b and the central axis have to be perpendicular for an elliptic cylinder.');
1033             end
1034             
1035             geo_code = [geo_code sprintf('(ellipticcylinder(%g, %g, %g ; %g, %g, %g ; %g, %g, %g) and plane(%g, %g, %g ; %g, %g, %g) and plane(%g, %g, %g ; %g, %g, %g))', ...
1036                         bottom_center(:,i), axis_a(:,i), axis_b(:,i), bottom_center(:,i), -central_axis, top_center(:,i), central_axis)];
1037                     
1038             if (i < n_elliptic_cylinders)
1039                 geo_code = [geo_code operator_string];
1040             else
1041                 geo_code = [geo_code ')'];             
1042             end
1043         end
1044     end
1045     
1046     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
1047     
1048 function geo_code = parse_geometry_half_space(half_space, operator_string)
1049 
1050     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
1051 
1052     if (~isstruct(half_space) || isempty(half_space))
1053         error('Parameter half_space must be a valid structure.');        
1054     else
1055         % Get number of half_spaces.
1056         n_half_spaces = numel(half_space);
1057         
1058         % Get structure field names.
1059         field_names = fieldnames(half_space);
1060         n_fields = numel(field_names);
1061         
1062         % Assign default values.
1063         point           = [0;0;0]*ones(1, n_half_spaces);
1064         outward_normal_vector  = [0;0;1]*ones(1, n_half_spaces);
1065         complement_flag = false(1, n_half_spaces);
1066         
1067         % Parse all structure fields.
1068         for i = 1:n_fields
1069             switch (field_names{i})  
1070                 case 'point'
1071                     point = assign_point(half_space, n_half_spaces, 'half_space', field_names{i}, point);
1072                 case 'outward_normal_vector'
1073                     outward_normal_vector = assign_point(half_space, n_half_spaces, 'half_space', field_names{i}, outward_normal_vector);
1074                 case 'complement_flag'
1075                     complement_flag = assign_flag(half_space, n_half_spaces, 'half_space', field_names{i}, complement_flag);
1076                 otherwise
1077                     error(['Field name ''%s'' is not valid for a half_space!\Available field names for a half_space are: ' ...
1078                            'point, outward_normal_vector, and complement_flag.'], field_names{i});
1079             end
1080         end
1081         
1082         % Start geo code with an opening parenthesis.
1083         geo_code = '(';
1084 
1085         % Add geo code for each half_space.
1086         for i = 1:n_half_spaces
1087             if (complement_flag(i))
1088                 geo_code = [geo_code 'not '];
1089             end
1090             
1091             geo_code = [geo_code sprintf('plane(%g, %g, %g ; %g, %g, %g)', ...
1092                         point(:,i), outward_normal_vector(:,i))];
1093                     
1094             if (i < n_half_spaces)
1095                 geo_code = [geo_code operator_string];
1096             else
1097                 geo_code = [geo_code ')'];             
1098             end
1099         end
1100     end
1101     
1102     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
1103 
1104 function geo_code = parse_geometry_ortho_brick(ortho_brick, operator_string)
1105 
1106     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
1107 
1108     if (~isstruct(ortho_brick) || isempty(ortho_brick))
1109         error('Parameter ortho_brick must be a valid structure.');        
1110     else
1111         % Get number of ortho_bricks.
1112         n_ortho_bricks = numel(ortho_brick);
1113         
1114         % Get structure field names.
1115         field_names = fieldnames(ortho_brick);
1116         n_fields = numel(field_names);
1117         
1118         % Assign default values.
1119         opposite_corner_a = [0;0;0]*ones(1, n_ortho_bricks);
1120         opposite_corner_b = [1;1;1]*ones(1, n_ortho_bricks);
1121         complement_flag   = zeros(1, n_ortho_bricks);
1122         
1123         % Parse all structure fields.
1124         for i = 1:n_fields
1125             switch (field_names{i})  
1126                 case {'opposite_corner_a'}
1127                     opposite_corner_a = assign_point(ortho_brick, n_ortho_bricks, 'ortho_brick', field_names{i}, opposite_corner_a);
1128                 case {'opposite_corner_b'}
1129                     opposite_corner_b = assign_point(ortho_brick, n_ortho_bricks, 'ortho_brick', field_names{i}, opposite_corner_b);
1130                 case 'complement_flag'
1131                     complement_flag = assign_flag(ortho_brick, n_ortho_bricks, 'ortho_brick', field_names{i}, complement_flag);
1132                 otherwise
1133                     error(['Field name "%s" is not allowed for an ortho_brick!\nAllowed field names for an ortho_brick are: ' ...
1134                            'opposite_corner_a, opposite_corner_b, and complement_flag.'], field_names{i});
1135             end
1136         end
1137         
1138         % Start geo code with an opening parenthesis.
1139         geo_code = '(';
1140 
1141         % Add geo code for each ortho_brick.
1142         for i = 1:n_ortho_bricks
1143             if (complement_flag(i))
1144                 geo_code = [geo_code 'not '];
1145             end
1146             
1147             geo_code = [geo_code sprintf('orthobrick(%g, %g, %g ; %g, %g, %g)', ...
1148                         min([opposite_corner_a(:, i) opposite_corner_b(:, i)], [], 2), ...
1149                         max([opposite_corner_a(:, i) opposite_corner_b(:, i)], [], 2))];
1150                     
1151             if (i < n_ortho_bricks)
1152                 geo_code = [geo_code operator_string];
1153             else
1154                 geo_code = [geo_code ')'];             
1155             end
1156         end
1157     end
1158     
1159     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
1160     
1161 function geo_code = parse_geometry_parallelepiped(parallelepiped, operator_string)
1162 
1163     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
1164 
1165     if (~isstruct(parallelepiped) || isempty(parallelepiped))
1166         error('Parameter parallelepiped must be a valid structure.');        
1167     else
1168         % Get number of parallelepiped.
1169         n_parallelepipeds = numel(parallelepiped);
1170         
1171         % Get structure field names.
1172         field_names = fieldnames(parallelepiped);
1173         n_fields = numel(field_names);
1174         
1175         % Assign default values.
1176         vertex          = [0;0;0]*ones(1, n_parallelepipeds);
1177         vector_a        = [1;0;0]*ones(1, n_parallelepipeds);
1178         vector_b        = [0;1;0]*ones(1, n_parallelepipeds);
1179         vector_c        = [0;0;1]*ones(1, n_parallelepipeds);
1180         complement_flag = zeros(1, n_parallelepipeds);
1181         
1182         % Parse all structure fields.
1183         for i = 1:n_fields
1184             switch (field_names{i})  
1185                 case 'vertex'
1186                     vertex = assign_point(parallelepiped, n_parallelepipeds, 'parallelepiped', field_names{i}, vertex);
1187                 case 'vector_a'
1188                     vector_a = assign_point(parallelepiped, n_parallelepipeds, 'parallelepiped', field_names{i}, vector_a);
1189                 case 'vector_b'
1190                     vector_b = assign_point(parallelepiped, n_parallelepipeds, 'parallelepiped', field_names{i}, vector_b);
1191                 case 'vector_c'
1192                     vector_c = assign_point(parallelepiped, n_parallelepipeds, 'parallelepiped', field_names{i}, vector_c);
1193                 case 'complement_flag'
1194                     complement_flag = assign_flag(parallelepiped, n_parallelepipeds, 'parallelepiped', field_names{i}, complement_flag);
1195                 otherwise
1196                     error(['Field name "%s" is not allowed for a parallelepiped!\nAllowed field names for a parallelepiped are: ' ...
1197                            'vertex, vector_a, vector_b, vector_c, and complement_flag.'], field_names{i});
1198             end
1199         end
1200         
1201         % Start geo code with an opening parenthesis.
1202         geo_code = '(';
1203 
1204         % Add geo code for each parallelepiped.
1205         for i = 1:n_parallelepipeds
1206             if (complement_flag(i))
1207                 geo_code = [geo_code 'not '];
1208             end
1209              
1210             % Make sure all vectors are not coplanar.
1211             if (abs(dot(vector_a(:,i), cross(vector_b(:,i), vector_c(:,i)))) < eps)
1212                 error('parallelepiped(%d) description includes coplanar vectors.', i);
1213             end
1214             
1215             % Compute opposite vertex.
1216             opposite_vertex = vertex(:,i) + vector_a(:,i) + vector_b(:,i) + vector_c(:,i);
1217             
1218             % Compute normal vectors.
1219             n_vector_ab = cross(vector_a(:,i), vector_b(:,i));
1220             n_vector_ac = cross(vector_a(:,i), vector_c(:,i));
1221             n_vector_bc = cross(vector_b(:,i), vector_c(:,i));
1222             
1223             % Check normal vectors directions.
1224             if (dot(n_vector_ab, vector_c(:,i)) < 0)
1225                 n_vector_ab = -n_vector_ab;
1226             end
1227             if (dot(n_vector_ac, vector_b(:,i)) < 0)
1228                 n_vector_ac = -n_vector_ac;
1229             end
1230             if (dot(n_vector_bc, vector_a(:,i)) < 0)
1231                 n_vector_bc = -n_vector_bc;
1232             end
1233             
1234             geo_code = [geo_code sprintf(['(plane(%g, %g, %g ; %g, %g, %g) and' ...
1235                                           ' plane(%g, %g, %g ; %g, %g, %g) and' ...
1236                                           ' plane(%g, %g, %g ; %g, %g, %g) and' ...
1237                                           ' plane(%g, %g, %g ; %g, %g, %g) and' ...
1238                                           ' plane(%g, %g, %g ; %g, %g, %g) and' ...
1239                                           ' plane(%g, %g, %g ; %g, %g, %g))'], ...
1240                         vertex(:,i), -n_vector_ab, ...
1241                         vertex(:,i), -n_vector_ac, ...
1242                         vertex(:,i), -n_vector_bc, ...
1243                         opposite_vertex, n_vector_ab, ...
1244                         opposite_vertex, n_vector_ac, ...
1245                         opposite_vertex, n_vector_bc)];
1246                     
1247             if (i < n_parallelepipeds)
1248                 geo_code = [geo_code operator_string];
1249             else
1250                 geo_code = [geo_code ')'];             
1251             end
1252         end
1253     end
1254     
1255     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
1256     
1257 function geo_code = parse_geometry_sphere(sphere, operator_string)
1258 
1259     eidors_msg('@@@ called with "%s" starting.', operator_string, 4);
1260 
1261     if (~isstruct(sphere) || isempty(sphere))
1262         error('Parameter sphere must be a valid structure.');        
1263     else
1264         % Get number of spheres.
1265         n_spheres = numel(sphere);
1266         
1267         % Get structure field names.
1268         field_names = fieldnames(sphere);
1269         n_fields = numel(field_names);
1270   
1271         % Assign default values.
1272         radius          = ones(1, n_spheres);
1273         center          = [0;0;0]*ones(1, n_spheres);
1274         complement_flag = false(1, n_spheres);
1275         
1276         % Parse all structure fields.
1277         for i = 1:n_fields
1278             switch (field_names{i})
1279                 case {'center', 'centre'}
1280                     center = assign_point(sphere, n_spheres, 'sphere', field_names{i}, center);
1281                 case 'radius'
1282                     radius = assign_radius(sphere, n_spheres, 'sphere', field_names{i}, radius);     
1283                 case 'complement_flag'
1284                     complement_flag = assign_flag(sphere, n_spheres, 'sphere', field_names{i}, complement_flag);
1285                 otherwise
1286                     error(['Field name ''%s'' is not valid for a sphere.\nAvailable field names for a sphere are: ' ...
1287                            'center, radius, and complement_flag.'], field_names{i});
1288             end
1289         end
1290         
1291         % Start geo code with an opening parenthesis.
1292         geo_code = '(';
1293 
1294         % Add geo code for each sphere.
1295         for i = 1:n_spheres
1296             if (complement_flag(i))
1297                 geo_code = [geo_code 'not '];
1298             end
1299                 
1300             geo_code = [geo_code sprintf('sphere(%g, %g, %g ; %g)', ...
1301                         center(:,i), radius(i))];
1302                     
1303             if (i < n_spheres)
1304                 geo_code = [geo_code operator_string];
1305             else
1306                 geo_code = [geo_code ')'];             
1307             end
1308         end
1309     end
1310     
1311     eidors_msg('@@@ called with "%s" returning.', operator_string, 4);
1312                              
1313 function write_geo_file(geo_fn, body_solid_code, electrode_solid_code, body_extra_code, electrode_extra_code, body_extra_param, electrode_extra_param)
1314     
1315     % Open geo file for writing.
1316     fid = fopen(geo_fn, 'w');
1317     
1318     if (fid == -1)
1319         error('Unable to open file %s for writing.', geo_fn);
1320     end
1321     
1322     % Write header for geo file.
1323     fprintf(fid, '#Automatically generated by ng_mk_geometric_models\n\n');
1324     fprintf(fid, 'algebraic3d\n\n');
1325     
1326     % Assemble a string to represent the union of all bodies.
1327     total_body_solid = '(';
1328    
1329     for i = 1:numel(body_solid_code)
1330         total_body_solid = [total_body_solid body_extra_param{i}.name];
1331 
1332         if (i < numel(body_solid_code))
1333             total_body_solid = [total_body_solid ' or '];
1334         else
1335             total_body_solid = [total_body_solid ')'];             
1336         end
1337     end
1338     
1339     % Assemble a string to represent the union of all electrodes entering the body.
1340     total_electrode_solid = '(';
1341     n_total_electrode_solid = 0;
1342    
1343     for i = 1:numel(electrode_solid_code)
1344         if (electrode_extra_param{i}.enter_body_flag)
1345             if (n_total_electrode_solid > 0)
1346                 total_electrode_solid = [total_electrode_solid ' or '];
1347             end
1348             total_electrode_solid = [total_electrode_solid electrode_extra_param{i}.name];
1349             n_total_electrode_solid = n_total_electrode_solid + 1;
1350         end
1351     end
1352     total_electrode_solid = [total_electrode_solid ')'];   
1353     
1354     % Write body_extra_code and electrode_extra_code in geo file
1355     for i = 1:numel(body_extra_code)
1356         if (~isempty(body_extra_code{i}))
1357             fprintf(fid, body_extra_code{i});
1358         end
1359     end
1360     for i = 1:numel(electrode_extra_code)
1361         if (~isempty(electrode_extra_code{i}))
1362             fprintf(fid, electrode_extra_code{i});
1363         end
1364     end
1365     fprintf(fid, '\n');
1366  
1367     % Write electrode solids that enter the body in geo file.
1368     for i = 1:numel(electrode_solid_code)
1369         if (~isempty(electrode_solid_code{i}) && electrode_extra_param{i}.enter_body_flag)
1370             fprintf(fid, 'solid %s = %s;\n\n', electrode_extra_param{i}.name, electrode_solid_code{i});
1371         end
1372     end
1373     
1374     % Write body solids in geo file.
1375     for i = 1:numel(body_solid_code)
1376         if (n_total_electrode_solid == 0)
1377             fprintf(fid, 'solid %s = %s;\n\n', body_extra_param{i}.name, body_solid_code{i});
1378         else
1379             fprintf(fid, 'solid %s = not %s and %s;\n\n', body_extra_param{i}.name, total_electrode_solid, body_solid_code{i});            
1380         end
1381     end
1382  
1383     % Write electrode solids that do not enter the body in geo file.
1384     for i = 1:numel(electrode_solid_code)
1385         if (~isempty(electrode_solid_code{i}) && ~electrode_extra_param{i}.enter_body_flag)
1386             fprintf(fid, 'solid %s = not %s and %s;\n\n', electrode_extra_param{i}.name, total_body_solid, electrode_solid_code{i});
1387         end
1388     end
1389     
1390     % Write electrode tlos in geo file.
1391     for i = 1:numel(electrode_solid_code)
1392         if (~isempty(electrode_solid_code{i}))
1393             if (isinf(electrode_extra_param{i}.max_edge_length))
1394                 fprintf(fid, 'tlo %s -col=[1,0,0] -material=%s;\n', electrode_extra_param{i}.name, electrode_extra_param{i}.name);
1395             else
1396                 fprintf(fid, 'tlo %s -col=[1,0,0] -material=%s -maxh=%g;\n', electrode_extra_param{i}.name, electrode_extra_param{i}.name, electrode_extra_param{i}.max_edge_length);          
1397             end
1398         end
1399     end
1400     
1401     % Assume the first object is the main one
1402     mainobj = 'MainNetgenObject';
1403     fprintf(fid, 'solid %s = %s ', mainobj, body_extra_param{1}.name);
1404     for i= 2:numel(body_solid_code)
1405        fprintf(fid, 'and (not %s)', body_extra_param{i}.name);
1406     end
1407     fprintf(fid, ';\n');
1408     body_extra_param{1}.name = mainobj; % rename it for following code
1409      
1410     % Write body tlos in geo file.
1411     for i = 1:numel(body_solid_code)
1412         if (isinf(body_extra_param{i}.max_edge_length))
1413             fprintf(fid, 'tlo %s -col=[0,1,0] -material=%s;\n', body_extra_param{i}.name, body_extra_param{i}.name);
1414         else
1415             fprintf(fid, 'tlo %s -col=[0,1,0] -material=%s -maxh=%g;\n', body_extra_param{i}.name, body_extra_param{i}.name, body_extra_param{i}.max_edge_length);            
1416         end
1417     end
1418     
1419     % Close file.
1420     fclose(fid);
1421 
1422 function mat = read_mat_from_file(fid, nrows, ncols)
1423     mat = fscanf(fid, '%g', [ncols, nrows])';
1424 
1425     % Skip to next line.
1426     if (~isempty(fgetl(fid)))
1427         error('Last line was only partialy read.');
1428     end
1429     
1430 function fmdl = read_vol_file(vol_fn, electrode_extra_param)
1431 
1432     % Open file for reading.
1433     fid = fopen(vol_fn, 'r');
1434 
1435     if (fid == -1)
1436         error('Unable to open file %s for reading.', vol_fn);
1437     end
1438     
1439     % Read a first line in vol file.
1440     line = fgetl(fid);
1441    
1442     % While no EOF or "endmesh" keyword is found.
1443     while (ischar(line) && ~strcmp(line, 'endmesh'))
1444         
1445         % Parse every line if not comment or empty line
1446         if (~isempty(line) && line(1) ~= '#') % Supposing '#' is always the first character of a comment line.
1447             switch(line)
1448                 case 'mesh3d'   % Nothing to do.
1449                 case 'dimension'
1450                     dimension = read_mat_from_file(fid, 1, 1);
1451                     if (dimension ~= 3)
1452                         error('unknown dimension %g in vol file.', dimension);
1453                     end
1454                 case 'geomtype'
1455                     geomtype = read_mat_from_file(fid, 1, 1);
1456                     if (geomtype ~= 0)
1457                         error('unknown %g geomtype in vol file.', geomtype);
1458                     end
1459                 case 'surfaceelements'
1460                     % # surfnr    bcnr   domin  domout      np      p1      p2      p3
1461                     n_surface_elements = read_mat_from_file(fid, 1, 1);
1462                     if (n_surface_elements)
1463                         surface_elements   = read_mat_from_file(fid, n_surface_elements, 8);
1464                     else
1465                         error('vol file contains no surface elements. There is probably something wrong with the provided geometry description.');    
1466                     end
1467                 case 'volumeelements'
1468                     % #  matnr      np      p1      p2      p3      p4
1469                     n_volume_elements = read_mat_from_file(fid, 1, 1);
1470                     if (n_volume_elements)
1471                         volume_elements   = read_mat_from_file(fid, n_volume_elements, 6);
1472                     else
1473                         error('vol file contains no volume elements. There is probably something wrong with the provided geometry description.');    
1474                     end
1475                 case 'edgesegmentsgi2'
1476                     % # surfid  0   p1   p2   trignum1    trignum2   domin/surfnr1    domout/surfnr2   ednr1   dist1   ednr2   dist2
1477                     n_edge_segments_sgi2 = read_mat_from_file(fid, 1, 1);
1478                     edge_segments_sgi2   = read_mat_from_file(fid, n_edge_segments_sgi2, 12);
1479                 case 'points'
1480                     % #          X             Y             Z
1481                     n_points = read_mat_from_file(fid, 1, 1);
1482                     if (n_points)
1483                         points   = read_mat_from_file(fid, n_points, 3);
1484                     else
1485                         error('vol file contains no points. There is probably something wrong with the provided geometry description.');                       
1486                     end
1487                 case 'materials'
1488                     n_materials = read_mat_from_file(fid, 1, 1);
1489                     if (n_materials)
1490                         materials   = cell(n_materials, 2);
1491                         % Read and parse each material line.
1492                         for i = 1:n_materials
1493                             material_line = fgetl(fid);
1494                             sscanf_result = sscanf(material_line, '%g%c%s')';
1495                             materials{i, 1} = sscanf_result(1);
1496                             materials{i, 2} = char(sscanf_result(3:end));
1497                         end
1498                     else
1499                         error('vol file contains no materials. There is probably something wrong with the provided geometry description.');                             
1500                     end
1501                 case 'face_colours'
1502                     % #   Surfnr     Red     Green     Blue
1503                     n_face_colours = read_mat_from_file(fid, 1, 1);
1504                     face_colours   = read_mat_from_file(fid, n_face_colours, 4);
1505                 otherwise
1506                     error('unknown "%s" line in vol file.', line);
1507             end
1508         end
1509         
1510         % Read next line in vol file.
1511         line = fgetl(fid);
1512     end
1513     
1514     % Close file.
1515     fclose(fid);
1516     
1517     if (~exist('points', 'var'))
1518         error('Point description is missing from vol file.');
1519     end
1520  
1521     if (~exist('volume_elements', 'var'))
1522         error('Volume element description is missing from vol file.');
1523     end
1524     
1525     if (~exist('surface_elements', 'var'))
1526         error('Surface element description is missing from vol file.');
1527     end
1528     
1529     if (~exist('materials', 'var'))
1530         error('Material description is missing from vol file.');
1531     end
1532     
1533     % Find electrode and body material indices.
1534     electrode_material = [];
1535     for i = 1:n_materials
1536         material_name   = materials{i, 2};
1537         material_number = materials{i, 1};
1538         
1539 %         if (strncmp(material_name, 'electrode', 9))
1540 %             % Extract electrode number from material_name
1541 %             electrode_number = str2double(material_name(10:end));
1542 %             electrode_material(electrode_number) = material_number;
1543 %         end
1544         for j = 1:numel(electrode_extra_param)
1545             if (strcmp(material_name, electrode_extra_param{j}.name))
1546                 electrode_material(j) = material_number;
1547             end
1548         end
1549     end
1550    
1551     % Remove electrode material if necessary
1552     original_n_nodes     = size(points, 1);
1553     original_n_elements  = size(volume_elements, 1);
1554     original_n_surfaces  = size(surface_elements, 1);
1555     original_n_materials = size(materials, 1);
1556 
1557     for i = 1:numel(electrode_material)
1558         if (~electrode_extra_param{i}.keep_material_flag)
1559             % Remove unwanted volume elements
1560             volume_elements(volume_elements(:, 1) == electrode_material(i), :) = [];
1561 
1562             % Remove unwanted surface elements
1563             surface_elements(surface_elements(:, 3) == electrode_material(i) & ...
1564                              surface_elements(:, 4) == 0 | ...
1565                              surface_elements(:, 4) == electrode_material(i) & ...
1566                              surface_elements(:, 3) == 0, :) = [];
1567         end
1568     end
1569 
1570     % Find nodes that are now unused.
1571     unused_nodes = true(1, size(points, 1));
1572     unused_nodes(volume_elements(:, 3:6))  = false;
1573     unused_nodes(surface_elements(:, 6:8)) = false;     
1574 
1575     % Remove unused points.
1576     points(unused_nodes, :) = [];
1577 
1578     % Compute new node indices after node removal.
1579     new_node_index = (1:original_n_nodes) - cumsum(unused_nodes);   
1580 
1581     % Update node indices for surface and volume elements.
1582     surface_elements(:, 6:8) = new_node_index(surface_elements(:, 6:8));
1583     volume_elements(:, 3:6)  = new_node_index(volume_elements(:, 3:6));
1584 
1585     % Find materials that are now unused.
1586     unused_materials = true(1, size(materials, 1));
1587     unused_materials(volume_elements(:, 1)) = false;
1588 
1589     % Remove unused materials.
1590     materials(unused_materials, :) = [];
1591 
1592     % Compute new material indices after material removal.
1593     new_material_index = (1:original_n_materials) - cumsum(unused_materials);   
1594 
1595     % Update material indices for volume elements.
1596     volume_elements(:, 1)  = new_material_index(volume_elements(:, 1));
1597 
1598     eidors_msg('@@@ Removed %d nodes, %d elements, %d surfaces and %d materials', ...
1599         original_n_nodes     - size(points, 1), ...
1600         original_n_elements  - size(volume_elements, 1), ...
1601         original_n_surfaces  - size(surface_elements, 1), ...
1602         original_n_materials - size(materials, 1), 3);
1603     
1604     % Assign mesh data to fmdl structures.
1605     fmdl.nodes            = points;
1606     fmdl.elems            = volume_elements(:, 3:6);
1607     fmdl.boundary         = surface_elements(:, 6:8);
1608     fmdl.boundary_numbers = surface_elements(:, 2);
1609     for i=1:max(volume_elements(:,1))
1610        fmdl.mat_idx{i}    = find( volume_elements(:, 1) == i);
1611     end
1612     fmdl.mat_name         = materials(:, 2);
1613     
1614     % Find electrode surfaces and nodes.
1615     for i = 1:numel(electrode_material)
1616         % Find surfaces that are part of the electrodes.
1617         if (electrode_extra_param{i}.keep_material_flag)
1618             electrode_boundary = ...
1619                sort(find(surface_elements(:, 3) == 0 & ...
1620                          surface_elements(:, 4) == electrode_material(i) | ...
1621                          surface_elements(:, 4) == 0 & ...
1622                          surface_elements(:, 3) == electrode_material(i)))';
1623         else
1624             electrode_boundary = ...
1625                 sort(find(surface_elements(:, 3) == electrode_material(i) | ...
1626                           surface_elements(:, 4) == electrode_material(i)))';
1627         end
1628      
1629         if (isempty(electrode_boundary))
1630             eidors_msg('WARNING: Electrode #%04d has been removed since it does not contact any body.', i, 2);
1631         else
1632             fmdl.electrode(i).boundary = electrode_boundary;
1633             
1634             % Find nodes that are part of the electrodes.
1635             fmdl.electrode(i).nodes = ...
1636                 unique(fmdl.boundary(fmdl.electrode(i).boundary(:), :))';                          
1637 
1638             % Assign default contact impedance.
1639             fmdl.electrode(i).z_contact = 0.01;
1640             
1641             if (~isempty(electrode_extra_param{i}.name))
1642                 fmdl.electrode(i).name = electrode_extra_param{i}.name;
1643             end
1644         end
1645     end
1646 
1647 function fmdl = complete_fmdl(fmdl, electrode_extra_param)
1648  
1649     % Find center point of domain.
1650     domain_center  = (max(fmdl.nodes)-min(fmdl.nodes))/2 + min(fmdl.nodes);
1651     domain_centers = ones(size(fmdl.nodes, 1), 1)*domain_center;
1652     
1653     % Find node closest to center for ground node.
1654     [unused, min_idx] = min(sum((fmdl.nodes - domain_centers).^2, 2));
1655     fmdl.gnd_node     = min_idx(1);
1656 
1657     fmdl.np_fwd_solve.perm_sym = '{n}';
1658 
1659     fmdl.name = 'ng_mk_geometric_models';
1660 
1661     fmdl.solve=      'eidors_default';
1662     fmdl.jacobian=   'eidors_default';
1663     fmdl.system_mat= 'eidors_default';
1664 
1665     fmdl.normalize_measurements = 0;
1666     
1667     for i = 1:numel(electrode_extra_param)
1668         if (isfield(electrode_extra_param{i}, 'point'))
1669             % Find center point of domain.
1670             electrode_points = ones(size(fmdl.nodes, 1), 1)*electrode_extra_param{i}.point';
1671 
1672             % Find node closest to the electrode point.
1673             [unused, min_idx]       = min(sum((fmdl.nodes - electrode_points).^2, 2));
1674             fmdl.electrode(i).nodes = min_idx(1);
1675             fmdl.electrode(i).boundary = [];
1676 
1677             % Assign default contact impedance.
1678             fmdl.electrode(i).z_contact = 0.01;
1679         end
1680     end
1681 
1682     fmdl = eidors_obj('fwd_model', fmdl);
1683 
1684 function do_unit_test
1685     for tn = 1:do_test_number(0)
1686         eidors_msg('ng_mk_geometric_models: unit_test %02d', tn, 1);
1687         [fmdl, opts] = do_test_number(tn);
1688         subplot(1,3,1)
1689         show_fem(fmdl);
1690         title('show_fem', 'Interpreter', 'none', 'FontWeight', 'bold');
1691         subplot(1,3,2);
1692         show_fem_enhanced(fmdl, opts);
1693         title({'show_fem_enhanced'; '(default options)'}, 'Interpreter', 'none', 'FontWeight', 'bold');
1694         subplot(1,3,3);
1695         opts.edge.color = [0 0 1];
1696         opts.edge.width = 0;
1697         opts.edge.significant.color = [1 0 0];
1698         opts.edge.significant.width = 1.5;
1699         opts.edge.significant.viewpoint_dependent.color = [0 1 0];
1700         opts.edge.significant.viewpoint_dependent.width = 1.5;
1701         show_fem_enhanced(fmdl, opts);
1702         title({'show_fem_enhanced'; '(with some options)'}, 'Interpreter', 'none', 'FontWeight', 'bold');
1703         drawnow;
1704         %pause
1705     end
1706 
1707 function [fmdl, opts] = do_test_number(tn)
1708     opts = struct;
1709     switch tn
1710         % Simple 3D cylinder. Radius = 1 with no electrodes
1711         case 1;
1712             body_geometry.cylinder = struct;
1713             fmdl = ng_mk_geometric_models(body_geometry);
1714         % Simple 3D cylinder. Radius = 1 with 16 spherical electrodes.
1715         case 2;
1716             body_geometry.cylinder = struct;
1717             n_elect = 16;
1718             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1719             for i = 1:n_elect
1720                 electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
1721                 electrode_geometry{i}.sphere.radius = 0.1;
1722             end
1723             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1724         % Simple 3D cylinder. Radius = 1 with 16 cylindrical electrodes.
1725         case 3;
1726             body_geometry.cylinder = struct;
1727             n_elect = 16;
1728             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1729             for i = 1:n_elect
1730                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1731                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1732                 electrode_geometry{i}.cylinder.radius = 0.1;
1733             end
1734             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1735         case 4;
1736             body_geometry.cylinder = struct;
1737             n_elect = 16;
1738             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1739             for i = 1:n_elect
1740                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1741                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1742                 electrode_geometry{i}.cylinder.radius = 0.1;
1743                 electrode_geometry{i}.keep_material_flag = 1;
1744             end
1745             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1746         case 5;
1747             body_geometry.cylinder = struct;
1748             n_elect = 16;
1749             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1750             for i = 1:n_elect
1751                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1752                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1753                 electrode_geometry{i}.cylinder.radius = 0.1;
1754                 electrode_geometry{i}.enter_body_flag = 1;
1755             end
1756             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1757         case 6;
1758             body_geometry.cylinder = struct;
1759             n_elect = 16;
1760             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1761             for i = 1:n_elect
1762                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1763                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1764                 electrode_geometry{i}.cylinder.radius = 0.1;
1765                 electrode_geometry{i}.keep_material_flag = 1;
1766                 electrode_geometry{i}.enter_body_flag = 1;                
1767             end
1768             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1769         case 7;
1770             body_geometry.cylinder = struct;
1771             body_geometry.sphere.center = [0 0 1];
1772             n_elect = 16;
1773             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1774             for i = 1:n_elect
1775                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1776                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1777                 electrode_geometry{i}.cylinder.radius = 0.1;
1778             end
1779             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1780         case 8;
1781             body_geometry.cylinder  = struct;
1782             body_geometry.sphere(1) = struct;  
1783             body_geometry.sphere(2).center = [0 0 1];         
1784             n_elect = 16;
1785             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1786             for i = 1:n_elect
1787                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1788                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1789                 electrode_geometry{i}.cylinder.radius = 0.1;
1790             end
1791             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);   
1792         case 9;
1793             body_geometry.intersection.cylinder(1) = struct;
1794             body_geometry.intersection.cylinder(2).radius     = 0.5;
1795             body_geometry.intersection.cylinder(2).complement_flag = 1;   
1796             n_elect = 16;
1797             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1798             for i = 1:n_elect
1799                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1800                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1801                 electrode_geometry{i}.cylinder.radius = 0.1;
1802             end
1803             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1804         case 10;
1805             body_geometry.intersection(1).sphere(1).radius     = 0.5;
1806             body_geometry.intersection(1).sphere(1).center     = [0 0 2];
1807             body_geometry.intersection(1).sphere(1).complement_flag = 1;
1808             body_geometry.intersection(1).sphere(2).center     = [0 0 2];
1809             body_geometry.intersection(2).cylinder(1).top_center = [0 0 2];
1810             body_geometry.intersection(2).cylinder(2).radius     = 0.5;
1811             body_geometry.intersection(2).cylinder(2).top_center = [0 0 2];
1812             body_geometry.intersection(2).cylinder(2).complement_flag = 1;   
1813             n_elect = 16;
1814             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1815             for i = 1:n_elect
1816                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1817                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1818                 electrode_geometry{i}.cylinder.radius = 0.1;
1819             end
1820             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1821         case 11;
1822             body_geometry.intersection.union(1).sphere.radius = 0.5;
1823             body_geometry.intersection.union(1).sphere.center = [0 0 2];
1824             body_geometry.intersection.union(1).cylinder.radius = 0.5;
1825             body_geometry.intersection.union(1).cylinder.top_center = [0 0 2];
1826             body_geometry.intersection.union(1).complement_flag = 1;
1827             body_geometry.intersection.union(2).sphere.center = [0 0 2];
1828             body_geometry.intersection.union(2).cylinder.top_center = [0 0 2]; 
1829             n_elect = 16;
1830             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1831             for i = 1:n_elect
1832                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1833                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1834                 electrode_geometry{i}.cylinder.radius = 0.1;
1835             end
1836             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1837         case 12;
1838             body_geometry.cone = struct; 
1839             n_elect = 16;
1840             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1841             for i = 1:n_elect
1842                 electrode_geometry{i}.cylinder.top_center    = [0.85*cos(theta(i)) 0.85*sin(theta(i)) 0.5];
1843                 electrode_geometry{i}.cylinder.bottom_center = [0.65*cos(theta(i)) 0.65*sin(theta(i)) 0.5];
1844                 electrode_geometry{i}.cylinder.radius = 0.1;
1845             end
1846             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1847         case 13;
1848             body_geometry.cone = struct; 
1849             n_elect = 16;
1850             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1851             for i = 1:n_elect
1852                 electrode_geometry{i}.sphere.center = [0.75*cos(theta(i)) 0.75*sin(theta(i)) 0.5];
1853                 electrode_geometry{i}.sphere.radius = 0.1;
1854             end
1855             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1856         case 14;
1857             body_geometry.cone(1).top_center = [0 0 1.5];
1858             body_geometry.cone(1).bottom_center = [0 0 0.5];
1859             body_geometry.cone(2).top_center = [0 0 -1.5];
1860             body_geometry.cone(2).bottom_center = [0 0 -0.5];
1861             body_geometry.cylinder.top_center    = [0, 0, 0.5];
1862             body_geometry.cylinder.bottom_center = [0, 0, -0.5];
1863             n_elect = 16;
1864             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1865             for i = 1:n_elect
1866                 electrode_geometry{i}.sphere.center = [0.75*cos(theta(i)) 0.75*sin(theta(i)) 1.0];
1867                 electrode_geometry{i}.sphere.radius = 0.1;
1868                 electrode_geometry{i + n_elect}.sphere.center = [cos(theta(i)) sin(theta(i)) 0];
1869                 electrode_geometry{i + n_elect}.sphere.radius = 0.15;
1870                 electrode_geometry{i + 2*n_elect}.sphere.center = [0.75*cos(theta(i)) 0.75*sin(theta(i)) -1.0];
1871                 electrode_geometry{i + 2*n_elect}.sphere.radius = 0.1;
1872             end
1873             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1874             opts.edge.significant.angle = 15;
1875         case 15
1876             body_geometry.ortho_brick = struct;
1877             fmdl = ng_mk_geometric_models(body_geometry);
1878         case 16
1879             body_geometry.intersection.ortho_brick.opposite_corner_a = [0 0 0];
1880             body_geometry.intersection.ortho_brick.opposite_corner_b = [5 5 4];
1881             for i = 1:4; 
1882                 for j = 1:4; 
1883                     body_geometry.intersection.cylinder(i,j).radius = 0.15;
1884                     body_geometry.intersection.cylinder(i,j).top_center = [i, j, 4];
1885                     body_geometry.intersection.cylinder(i,j).bottom_center = [i, j, 2];
1886                     body_geometry.intersection.cylinder(i,j).complement_flag = 1;
1887                 end; 
1888             end;
1889             fmdl = ng_mk_geometric_models(body_geometry);    
1890         case 17
1891             body_geometry.intersection.ortho_brick.opposite_corner_a = [0 0 0];
1892             body_geometry.intersection.ortho_brick.opposite_corner_b = [5 5 4];
1893             for i = 1:4; 
1894                 for j = 1:4; 
1895                     body_geometry.intersection.cylinder(i, j).radius = 0.15;
1896                     body_geometry.intersection.cylinder(i, j).top_center    = [i, j, 4];
1897                     body_geometry.intersection.cylinder(i, j).bottom_center = [i, j, 2];
1898                     body_geometry.intersection.cylinder(i, j).complement_flag = 1;
1899                     electrode_geometry{i, j, 1}.cylinder.radius        = 0.2;
1900                     electrode_geometry{i, j, 1}.cylinder.top_center    = [i, j, 3.1];
1901                     electrode_geometry{i, j, 1}.cylinder.bottom_center = [i, j, 2.9];
1902                     electrode_geometry{i, j, 2}.cylinder.radius        = 0.2;
1903                     electrode_geometry{i, j, 2}.cylinder.top_center    = [i, j, 2.2];
1904                     electrode_geometry{i, j, 2}.cylinder.bottom_center = [i, j, 2.0];
1905                 end; 
1906             end;
1907             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1908         case 18
1909             body_geometry.parallelepiped  = struct;
1910             body_geometry.max_edge_length = 0.15;
1911             fmdl = ng_mk_geometric_models(body_geometry);
1912         case 19
1913             body_geometry.parallelepiped.vertex   = [ 0;  0;  0];
1914             body_geometry.parallelepiped.vector_a = [ 1;  1;  0];
1915             body_geometry.parallelepiped.vector_b = [ 0;  1;  1];
1916             body_geometry.parallelepiped.vector_c = [ 1;  0;  1];
1917             body_geometry.max_edge_length = 0.15;
1918             fmdl = ng_mk_geometric_models(body_geometry);
1919         case 20
1920             body_geometry.intersection.ortho_brick.opposite_corner_a = [-15, -15, 0];
1921             body_geometry.intersection.ortho_brick.opposite_corner_b = [15, 15, 5];
1922             body_geometry.intersection.half_space.point = [0, 0, 5];
1923             body_geometry.intersection.half_space.outward_normal_vector = [-1, -1, 5];
1924             
1925             fmdl = ng_mk_geometric_models(body_geometry);
1926         case 21
1927             body_geometry.ellipsoid.axis_a = [1 0 0];
1928             body_geometry.ellipsoid.axis_b = [0 2 0];
1929             body_geometry.ellipsoid.axis_c = [0 0 3];
1930             fmdl = ng_mk_geometric_models(body_geometry);   
1931         case 22
1932             body_geometry.ellipsoid.axis_a = [1 0 0];
1933             body_geometry.ellipsoid.axis_b = [0 1 1];
1934             body_geometry.ellipsoid.axis_c = [0 -2 2];
1935             fmdl = ng_mk_geometric_models(body_geometry);   
1936         case 23
1937             body_geometry.elliptic_cylinder.top_center = [0, 0, 10];
1938             body_geometry.elliptic_cylinder.bottom_center = [0, 0, 0];           
1939             body_geometry.elliptic_cylinder.axis_a = [1 0 0];
1940             body_geometry.elliptic_cylinder.axis_b = [0 2 0];  
1941             fmdl = ng_mk_geometric_models(body_geometry);
1942         case 24
1943             body_geometry.elliptic_cylinder.top_center = [0, 5, 5];
1944             body_geometry.elliptic_cylinder.bottom_center = [0, 0, 0];           
1945             body_geometry.elliptic_cylinder.axis_a = [1 0 0];
1946             body_geometry.elliptic_cylinder.axis_b = [0 -2 2];  
1947             fmdl = ng_mk_geometric_models(body_geometry);
1948         case 25
1949             body_geometry.body_of_revolution = struct;
1950             body_geometry.max_edge_length = 0.15;
1951             fmdl = ng_mk_geometric_models(body_geometry);
1952         case 26
1953             body_geometry.body_of_revolution.points   = [1 1; 1 2; 2 1.5; 2 1];
1954             body_geometry.body_of_revolution.segments = [1 2; 2 3; 3 4; 4 1];
1955             body_geometry.max_edge_length = 0.15;
1956             fmdl = ng_mk_geometric_models(body_geometry);
1957         case 27
1958             n_points = 24;
1959             theta = linspace(0, 2*pi, n_points+1)'; theta(end) = [];
1960             body_geometry.body_of_revolution.points   = 2 + [sin(theta) cos(theta)];
1961             body_geometry.body_of_revolution.segments = [(1:n_points)' [(2:n_points) 1]'];
1962             body_geometry.max_edge_length = 0.15;
1963             fmdl = ng_mk_geometric_models(body_geometry);
1964         case 28
1965             n_points = 24;
1966             theta = linspace(0, 2*pi, n_points+1)'; theta(end) = [];
1967             body_geometry.body_of_revolution.points   = 2 + [sin(theta) cos(theta)];
1968             body_geometry.body_of_revolution.segments = [(1:2:n_points)' (2:2:n_points)' [(3:2:n_points) 1]'];
1969             body_geometry.max_edge_length = 0.15;
1970             fmdl = ng_mk_geometric_models(body_geometry);
1971         case 29
1972             body_geometry{1}.cylinder(1).radius        = 0.5;
1973             body_geometry{1}.cylinder(1).top_center    = [0 0 0.75];
1974             body_geometry{1}.cylinder(1).bottom_center = [0 0 0.25];
1975             body_geometry{1}.name                      = 'Object';           
1976             body_geometry{2}.cylinder(2).radius        = 1;
1977             body_geometry{2}.name                      = 'Tank';
1978             n_elect = 16;
1979             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1980             for i = 1:n_elect
1981                 electrode_geometry{i}.cylinder.top_center    = [1.03*cos(theta(i)) 1.03*sin(theta(i)) 0.5];
1982                 electrode_geometry{i}.cylinder.bottom_center = [0.97*cos(theta(i)) 0.97*sin(theta(i)) 0.5];
1983                 electrode_geometry{i}.cylinder.radius = 0.1;
1984             end
1985             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1986         case 30
1987             body_geometry{1}.sphere.radius     = 0.25;
1988             body_geometry{1}.sphere.center     = [0 0 0.5];
1989             body_geometry{1}.name              = 'Sphere';
1990             body_geometry{2}.cylinder.radius   = 1;
1991             body_geometry{2}.name              = 'Tank';           
1992             n_elect = 16;
1993             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
1994             for i = 1:n_elect
1995                 electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
1996                 electrode_geometry{i}.sphere.radius = 0.1;
1997             end
1998             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
1999        case 31
2000             n_sphere = 8;
2001             theta = linspace(0, 2*pi, n_sphere+1); theta(end) = [];   
2002             for i = 1:n_sphere
2003                 body_geometry{i}.sphere.radius   = 0.2;
2004                 body_geometry{i}.sphere.center   = [0.65*cos(theta(i)) 0.65*sin(theta(i)) 0.5];  
2005                 body_geometry{i}.max_edge_length = 0.025*(1 + rem(i,2));
2006                 body_geometry{i}.name            = sprintf('Sphere%d', i);  
2007             end        
2008             body_geometry{n_sphere+1}.cylinder.radius = 1;
2009             body_geometry{n_sphere+1}.name            = 'Tank';  
2010             n_elect = 16;
2011             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
2012             for i = 1:n_elect
2013                 electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
2014                 electrode_geometry{i}.sphere.radius = 0.1;
2015                 electrode_geometry{i}.max_edge_length = 0.025*(1 + rem(i,2));
2016             end
2017             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
2018        case 32
2019             body_geometry.cylinder = struct;
2020             n_elect = 16;
2021             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
2022             for i = 1:n_elect
2023                 electrode_geometry{i}.point = [cos(theta(i)) sin(theta(i)) 0.5];
2024             end
2025             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
2026        case 33     
2027             body_geometry.cylinder = struct;
2028             n_elect = 16;
2029             theta = linspace(0, 2*pi, n_elect+1); theta(end) = [];
2030             for i = 1:n_elect
2031                 if (rem(i,2))
2032                     electrode_geometry{i}.point = [cos(theta(i)) sin(theta(i)) 0.5];
2033                     electrode_geometry{i}.name  = sprintf('Point_Electrode%d', ceil(i/2));
2034                 else
2035                     electrode_geometry{i}.sphere.center = [cos(theta(i)) sin(theta(i)) 0.5];
2036                     electrode_geometry{i}.sphere.radius = 0.1;
2037                     electrode_geometry{i}.name          = sprintf('Circular_Electrode%d', floor(i/2));
2038                 end
2039             end
2040             fmdl = ng_mk_geometric_models(body_geometry, electrode_geometry);
2041        case 34
2042             body_geometry.body_of_extrusion = struct;
2043             body_geometry.max_edge_length = 0.15;
2044             fmdl = ng_mk_geometric_models(body_geometry);
2045         case 35
2046             body_geometry.body_of_extrusion.path_points   = [0 0 0; 0.25 0 1; 0.25 0 2; 0.25 0 3; 0 0 4];
2047             body_geometry.body_of_extrusion.path_segments = [1 2; 2 3; 3 4; 4 5];
2048             body_geometry.max_edge_length = 0.15;
2049             fmdl = ng_mk_geometric_models(body_geometry);
2050        case 36
2051             n_points = 16;
2052             theta = linspace(0, 2*pi, n_points+1)'; theta(end) = [];
2053             body_geometry.body_of_extrusion.profile_points   = 0.2*(2 + [0.75*sin(theta) cos(theta)]);
2054             body_geometry.body_of_extrusion.profile_segments = [(1:n_points)' [(2:n_points) 1]'];
2055             n_points = 32;
2056             theta = linspace(0, 2*pi, n_points+1)'; theta(end) = [];          
2057             body_geometry.body_of_extrusion.path_points   = 1*(2 + [sin(theta) 1.5*cos(theta) zeros(n_points, 1)]);
2058             body_geometry.body_of_extrusion.path_segments = [(1:n_points)' [(2:n_points) 1]'];
2059             body_geometry.body_of_extrusion.vector_d      = [0; 0; 1];
2060             body_geometry.max_edge_length = 0.15;
2061             fmdl = ng_mk_geometric_models(body_geometry); 
2062         case 0; fmdl = 36; % Return maximum number of tests.
2063         otherwise;
2064             error('Invalid test number.')
2065     end