function strFile = data2geo(data,lc,varargin)
%% fileStr = data2geo(data,lc)
% Takes data assuming certain order and builds a .geo file compatible with gmsh.
%
%!demo
%! points = [0 0 0; 0.1 0 0; 0.1 .3 0; 0 0.3 0];
%! strFile = data2geo(points,0.009);
%! disp(strFile)
% -------------------
% Copyright 2010 Juan Pablo Carbajal
% carbajal@ifi.uzh.ch
% Monday, November 01 2010
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.


% This variable can then be used in the definition of Gmsh's simplest
% `elementary entity', a `Point'. A Point is defined by a list of
% four numbers: three coordinates (X, Y and Z), and a characteristic
% length (lc) that sets the target element size at the point:
% The distribution of the mesh element sizes is then obtained by
% interpolation of these characteristic lengths throughout the
% geometry. Another method to specify characteristic lengths is to
% use a background mesh (see `t7.geo' and `bgmesh.pos').

pointStr = @(n,xyz,l)sprintf('Point(%d) = {%f,%f,%f,%f};\n',n,xyz,l);

% We can then define some additional points as well as our first
% curve.  Curves are Gmsh's second type of elementery entities, and,
% amongst curves, straight lines are the simplest. A straight line is
% defined by a list of point numbers. In the commands below, for
% example, the line 1 starts at point 1 and ends at point 2:

lineStr = @(n,i,j) sprintf('Line(%d) = {%d,%d};\n',n,i,j);
 
% The third elementary entity is the surface. In order to define a
% simple rectangular surface from the four lines defined above, a
% line loop has first to be defined. A line loop is a list of
% connected lines, a sign being associated with each line (depending
% on the orientation of the line):

lineLoopStr = @(nl,lns) sprintf(['Line Loop(%d) = {%d' ...
                                        repmat(',%d',1,nl-2) '};\n'],nl,lns);
nl = @()sprintf('\n');

[n dim] = size(data);
if dim == 2
    data(:,3) = zeros(n,1);
end

header = ' // File created with Octave';
strFile = [];
strFile = [strFile header nl()];

% Points
strFile = [strFile '// Points' nl()];

for i=1:n
    strFile = [strFile pointStr(i,data(i,:),lc)];
end

% Lines
strFile = [strFile '// Lines' nl()];
for i=1:n-1
    strFile = [strFile lineStr(i,i,i+1)];
end
%strFile = [strFile lineStr(n,n,1)];

% Loop
strFile = [strFile lineLoopStr(n,1:n-1)];

if nargin > 2
    outfile = varargin{1};
    fid = fopen(outfile,'w');
    fprintf(fid,'%s',strFile);
    fclose(fid);
end    

