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 6913 2024-05-29 20:07:53Z bgrychtol $
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 fprintf(fid,'Mesh 2;\n');
0076         
0077 fclose(fid);
0078 
0079 
0080 function [shape, cl] = process_input(shape)
0081 cl = 0.1;
0082 if ~iscell(shape)
0083    shape = {shape};
0084 end
0085 
0086 if numel(shape) > 1 && numel(shape{end}) == 1
0087     cl = shape{end};
0088     shape(end) = [];
0089 end
0090 
0091 function do_unit_test
0092   th=linspace(0,2*pi,9);
0093   th(end) = [];
0094   xy = [cos(th);sin(th)]';
0095   fmdl=gmsh_mk_2d_model(xy);
0096   subplot(331); show_fem(fmdl) 
0097   title 'standard'
0098 
0099   unit_test_cmp('stop sign 01', max(fmdl.nodes), [1,1]);
0100   unit_test_cmp('stop sign 02', min(fmdl.nodes),-[1,1]);
0101 
0102   fmdl=gmsh_mk_2d_model({xy,0.5});
0103   subplot(332); show_fem(fmdl) 
0104   title 'coarse'
0105 
0106   fmdl=gmsh_mk_2d_model({xy,xy*0.5,0.2});
0107   subplot(333); show_fem(fmdl) 
0108   title 'hole'