gmsh_mk_2d_model

GMSH_MK_2D_MODEL create a 2D mesh with GMSH
 mdl = gmsh_mk_2d_model(shape)

 SHAPE can be:
  - xy (Nx2)             : a counter- clockwise list of points in 2D 
                           defining the outer contour
  - {xy, xy1, xy2, ...}  : allows specifying additional counter-clockwise 
                           loops  xy1, xy2, etc, which represent holes in  
                           the bigger contour xy contour
  - {..., maxsz}         : specifies the characteristic length at each
                           point (default: 0.1);

 For sufficently large MAXSZ, no new points are created on the boundary.

  See also NG_MK_2D_MODEL

0001 function mdl = gmsh_mk_2d_model(varargin)
0002 %GMSH_MK_2D_MODEL create a 2D mesh with GMSH
0003 % mdl = gmsh_mk_2d_model(shape)
0004 %
0005 % SHAPE can be:
0006 %  - xy (Nx2)             : a counter- clockwise list of points in 2D
0007 %                           defining the outer contour
0008 %  - {xy, xy1, xy2, ...}  : allows specifying additional counter-clockwise
0009 %                           loops  xy1, xy2, etc, which represent holes in
0010 %                           the bigger contour xy contour
0011 %  - {..., maxsz}         : specifies the characteristic length at each
0012 %                           point (default: 0.1);
0013 %
0014 % For sufficently large MAXSZ, no new points are created on the boundary.
0015 %
0016 %  See also NG_MK_2D_MODEL
0017 
0018 %% (C) 2021 Bartlomiej Grychtol, License: GPL version 2 or version 3
0019 % $Id: gmsh_mk_2d_model.m 6420 2022-11-29 15:08:09Z aadler $
0020 
0021 if nargin==1 && ischar(varargin{1}) && strcmp(varargin{1},'UNIT_TEST'); do_unit_test; return; end
0022 
0023 [shape, cl] = process_input(varargin{:});
0024 
0025 namestem = tempname;
0026 write_geo_file(shape, cl, namestem);
0027 call_gmsh([namestem '.geo']);
0028 [srf,vtx,fc,bc,simp,edg,mat_ind] = gmsh_read_mesh([namestem '.msh']);
0029 mdl.type     = 'fwd_model';
0030 mdl.name = 'gmsh_2d_model';
0031 mdl.nodes    = vtx;
0032 mdl.elems    = simp;
0033 mdl= eidors_obj('fwd_model', mdl);
0034 
0035 delete([namestem '.geo']);
0036 delete([namestem '.msh']);
0037 
0038 
0039 function fname = write_geo_file(shape, cl, namestem)
0040 fname = [namestem '.geo'];
0041 n_points = 0;
0042 n_lines = 0;
0043 loop_idx = [];
0044 fid= fopen(fname,'w');
0045 for i = 1:length(shape)
0046     n_sh_pts = size(shape{i},1);
0047     for p = 1:n_sh_pts
0048         fprintf(fid,'Point(%d) = {%f, %f, 0, %f};\n',...
0049             n_points + p, shape{i}(p,1), shape{i}(p,2), cl);
0050     end
0051     
0052     for p = 1:n_sh_pts-1
0053         fprintf(fid,'Line(%d) = {%d, %d};\n',...
0054             n_lines + p, n_points + p, n_points + p + 1);
0055     end
0056     fprintf(fid,'Line(%d) = {%d, %d};\n',...
0057         n_lines + n_sh_pts, n_points + n_sh_pts, n_points + 1);
0058     
0059     fprintf(fid,'Line Loop (%d) = {%d%s};\n',...
0060         n_lines + n_sh_pts + 1, n_lines + 1,...
0061         sprintf(', %d', n_lines + (2:n_sh_pts)));
0062     
0063     loop_idx = [loop_idx, n_lines + n_sh_pts + 1];
0064     n_points = n_points + n_sh_pts;
0065     n_lines = n_lines + n_sh_pts + 1;
0066 end
0067 str = '';
0068 if numel(loop_idx) > 1
0069     str = sprintf(', %d', loop_idx(2:end));
0070 end
0071 
0072 fprintf(fid,'Plane Surface (%d) = {%d%s};\n',...
0073         n_lines + 1, loop_idx(1), str);
0074         
0075 fclose(fid);
0076 
0077 
0078 function [shape, cl] = process_input(shape)
0079 cl = 0.1;
0080 if ~iscell(shape)
0081    shape = {shape};
0082 end
0083 
0084 if numel(shape) > 1 && numel(shape{end}) == 1
0085     cl = shape{end};
0086     shape(end) = [];
0087 end
0088 
0089 function do_unit_test
0090   th=linspace(0,2*pi,9);
0091   th(end) = [];
0092   xy = [cos(th);sin(th)]';
0093   fmdl=gmsh_mk_2d_model(xy);
0094   subplot(331); show_fem(fmdl) 
0095   title 'standard'
0096 
0097   unit_test_cmp('stop sign 01', max(fmdl.nodes), [1,1]);
0098   unit_test_cmp('stop sign 02', min(fmdl.nodes),-[1,1]);
0099 
0100   fmdl=gmsh_mk_2d_model({xy,0.5});
0101   subplot(332); show_fem(fmdl) 
0102   title 'coarse'
0103 
0104   fmdl=gmsh_mk_2d_model({xy,xy*0.5,0.2});
0105   subplot(333); show_fem(fmdl) 
0106   title 'hole'