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