ng_mk_geometric_models

PURPOSE ^

NG_MK_GEOMETRIC_MODELS: create geometric mesh models using Netgen. Body

SYNOPSIS ^

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

DESCRIPTION ^

 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);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005