0001 function [fmdl, mat_idx] = ng_mk_geometric_models(body_geometry, electrode_geometry)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132 if (ischar(body_geometry) && strcmp(body_geometry, 'UNIT_TEST'))
0133 do_unit_test;
0134 return;
0135 end
0136
0137
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
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
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
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
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
0172 n_body = numel(body_geometry);
0173 n_electrode = numel(electrode_geometry);
0174
0175
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
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
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
0206 fn_prefix = tempname;
0207 geo_fn = [fn_prefix, '.geo'];
0208 vol_fn = [fn_prefix, '.vol'];
0209
0210
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
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
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
0228 call_netgen(geo_fn, vol_fn);
0229
0230
0231 fmdl_mat_idx{1} = read_vol_file(vol_fn, electrode_extra_param);
0232
0233
0234 if ~eidors_debug('query','ng_mk_geometric_models:keep_temp_files')
0235 delete(geo_fn);
0236 delete(vol_fn);
0237 end
0238
0239
0240 fmdl_mat_idx{1} = complete_fmdl(fmdl_mat_idx{1}, electrode_extra_param);
0241
0242
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
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
0362 if (isfield(geometry, 'point'))
0363
0364 if (numel(geometry) ~= 1)
0365 error('Field name "point" must define only a single point.');
0366 end
0367
0368
0369 field_names = fieldnames(geometry);
0370 n_fields = numel(field_names);
0371
0372
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
0382 if (numel(geometry.point) ~= 3)
0383 error('geometry.point value is not valid.');
0384 end
0385
0386
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
0395 extra_code = '';
0396
0397
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
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
0415 if (~isstruct(geometry) || isempty(geometry))
0416 error('Parameter geometry must be a valid structure.');
0417 else
0418
0419 n_geometries = numel(geometry);
0420
0421
0422 field_names = fieldnames(geometry);
0423 n_fields = numel(field_names);
0424
0425
0426 geo_code = '(';
0427 for i = 1:n_geometries
0428
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
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
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
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
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
0562 n_body_of_extrusions = numel(body_of_extrusion);
0563
0564
0565 field_names = fieldnames(body_of_extrusion);
0566 n_fields = numel(field_names);
0567
0568
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
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
0598 geo_code = '(';
0599 extra_code = '';
0600
0601
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
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
0653
0654
0655
0656
0657
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
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
0685 n_body_of_revolutions = numel(body_of_revolution);
0686
0687
0688 field_names = fieldnames(body_of_revolution);
0689 n_fields = numel(field_names);
0690
0691
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
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
0718 geo_code = '(';
0719 extra_code = '';
0720
0721
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
0766 n_cones = numel(cone);
0767
0768
0769 field_names = fieldnames(cone);
0770 n_fields = numel(field_names);
0771
0772
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
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
0799 geo_code = '(';
0800
0801
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
0830 n_cylinders = numel(cylinder);
0831
0832
0833 field_names = fieldnames(cylinder);
0834 n_fields = numel(field_names);
0835
0836
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
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
0860 geo_code = '(';
0861
0862
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
0891 n_ellipsoids = numel(ellipsoid);
0892
0893
0894 field_names = fieldnames(ellipsoid);
0895 n_fields = numel(field_names);
0896
0897
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
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
0924 geo_code = '(';
0925
0926
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
0961 n_elliptic_cylinders = numel(elliptic_cylinder);
0962
0963
0964 field_names = fieldnames(elliptic_cylinder);
0965 n_fields = numel(field_names);
0966
0967
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
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
0994 geo_code = '(';
0995
0996
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
1033 n_half_spaces = numel(half_space);
1034
1035
1036 field_names = fieldnames(half_space);
1037 n_fields = numel(field_names);
1038
1039
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
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
1060 geo_code = '(';
1061
1062
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
1089 n_ortho_bricks = numel(ortho_brick);
1090
1091
1092 field_names = fieldnames(ortho_brick);
1093 n_fields = numel(field_names);
1094
1095
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
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
1116 geo_code = '(';
1117
1118
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
1146 n_parallelepipeds = numel(parallelepiped);
1147
1148
1149 field_names = fieldnames(parallelepiped);
1150 n_fields = numel(field_names);
1151
1152
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
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
1179 geo_code = '(';
1180
1181
1182 for i = 1:n_parallelepipeds
1183 if (complement_flag(i))
1184 geo_code = [geo_code 'not '];
1185 end
1186
1187
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
1193 opposite_vertex = vertex(:,i) + vector_a(:,i) + vector_b(:,i) + vector_c(:,i);
1194
1195
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
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
1242 n_spheres = numel(sphere);
1243
1244
1245 field_names = fieldnames(sphere);
1246 n_fields = numel(field_names);
1247
1248
1249 radius = ones(1, n_spheres);
1250 center = [0;0;0]*ones(1, n_spheres);
1251 complement_flag = false(1, n_spheres);
1252
1253
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
1269 geo_code = '(';
1270
1271
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
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
1300 fprintf(fid, '#Automatically generated by ng_mk_geometric_models\n\n');
1301 fprintf(fid, 'algebraic3d\n\n');
1302
1303
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
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
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
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
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
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
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
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;
1386
1387
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
1397 fclose(fid);
1398
1399 function mat = read_mat_from_file(fid, nrows, ncols)
1400 mat = fscanf(fid, '%g', [ncols, nrows])';
1401
1402
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
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
1417 line = fgetl(fid);
1418
1419
1420 while (ischar(line) && ~strcmp(line, 'endmesh'))
1421
1422
1423 if (~isempty(line) && line(1) ~= '#')
1424 switch(line)
1425 case 'mesh3d'
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
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
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
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
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
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
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
1488 line = fgetl(fid);
1489 end
1490
1491
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
1511 electrode_material = [];
1512 for i = 1:n_materials
1513 material_name = materials{i, 2};
1514 material_number = materials{i, 1};
1515
1516
1517
1518
1519
1520
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
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
1537 volume_elements(volume_elements(:, 1) == electrode_material(i), :) = [];
1538
1539
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
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
1553 points(unused_nodes, :) = [];
1554
1555
1556 new_node_index = (1:original_n_nodes) - cumsum(unused_nodes);
1557
1558
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
1563 unused_materials = true(1, size(materials, 1));
1564 unused_materials(volume_elements(:, 1)) = false;
1565
1566
1567 materials(unused_materials, :) = [];
1568
1569
1570 new_material_index = (1:original_n_materials) - cumsum(unused_materials);
1571
1572
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
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
1592 for i = 1:numel(electrode_material)
1593
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
1612 fmdl.electrode(i).nodes = ...
1613 unique(fmdl.boundary(fmdl.electrode(i).boundary(:), :))';
1614
1615
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
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
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
1647 electrode_points = ones(size(fmdl.nodes, 1), 1)*electrode_extra_param{i}.point';
1648
1649
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
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
1682 end
1683
1684 function [fmdl, opts] = do_test_number(tn)
1685 opts = struct;
1686 switch tn
1687
1688 case 1;
1689 body_geometry.cylinder = struct;
1690 fmdl = ng_mk_geometric_models(body_geometry);
1691
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
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;
2040 otherwise;
2041 error('Invalid test number.')
2042 end