## Copyright (C) 2013 fotios
## 
## 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
## (at your option) 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 Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## buildmesh

## Author: fotios <address@hidden>
## Created: 2013-05-30

function [p, e, t, hmax] = buildmesh (fix, varargin)
  
  for cnt = 1 : ( nargin - 1 )
    if ( !ischar (varargin{cnt}) )
      error ("arguments expected to be strings");
    endif
  endfor
  
  fid = fopen ("ff2O.edp", "w");
    if ( !exist ("ff2O_bound.idp") )
      error ("boundary file is not found");
    else
      fprintf (fid, "include \"ff2O_bound.idp\";\n");
    endif
    
    bnd = [];
    for cnt = 1 : ( nargin - 2 )
      bnd = [bnd, varargin{cnt}, " + "];
    endfor
    bnd = [bnd, varargin{end}];
    fprintf (fid, ["mesh Th = buildmesh(", bnd, ");\n"]);

    if ( fix == 1 )
      fprintf (fid, "include \"ff2O_fix.idp\";\n");
    endif
    
    fprintf (fid, "savemesh(Th, \"ff2O_mesh.msh\");");
  fclose (fid);

  system ("FreeFem++ ff2O.edp");
  
  [p, t, e] = readmesh ();
  
  not = numel (t(:, 1));
  hmax = zeros (1, not);
  
  for cnt = 1 : not
    tp = p(t(cnt, 1:3), 1:2);
    tp = [tp; tp(1,:)];
    h = max (sqrt ( diff (tp(:,1)).^2 + diff (tp(:,2)).^2 ));
    hmax = max ([hmax, h]);
  endfor
  
  p = p(:, 1:2);
  t = t(:, 1:3);
  e = e(:, 1:3);
  
  system ("rm ff2O*.*");

endfunction