Source file meshF.ml

# 1 "src/meshFC.ml"
(* Functions for layout fortran_layout.
 ***********************************************************************)

open Printf
open Bigarray
open Mesh_common

type mesh = fortran_layout t
type vec = fortran_layout Mesh_common.vec
type mat = fortran_layout Mesh_common.mat
type int_mat = fortran_layout Mesh_common.int_mat
type int_vec = fortran_layout Mesh_common.int_vec

let layout = fortran_layout;;

let empty_vec = Array1.create int layout 0
let empty_mat2 = Array2.create (float64) fortran_layout (2) (0)
let empty_mat4 = Array2.create (float64) fortran_layout (4) (0)
let empty_int_mat2 = Array2.create (int) fortran_layout (2) (0)
let empty_int_mat3 = Array2.create (int) fortran_layout (3) (0)

let check_point name point =
  if Array2.dim2(point) = 0 then
    invalid_arg(name ^ ": points cannot be empty");
  if Array2.dim1(point) <> 2 then
    invalid_arg(name ^ ": dim1 points must be 2")

let check_point_marker name point = function
  | None -> empty_vec
  | Some m ->
     let n = Array1.dim m in
     if 0 < n && n < Array2.dim2(point) then
       invalid_arg(name ^ ": point_marker too small");
     m

let check_segment_marker name segment = function
  | None -> empty_vec
  | Some m ->
     let n = Array1.dim m in
     if 0 < n && n < Array2.dim2(segment) then
       invalid_arg(name ^ ": segment_marker too small");
     m

let check_hole name = function
  | None -> empty_mat2
  | Some h ->
     if Array2.dim2(h) > 0 && Array2.dim1(h) <> 2 then
       invalid_arg(name ^ ": dim1 hole must be 2");
     h

let check_region name = function
  | None -> empty_mat4
  | Some r ->
     if Array2.dim2(r) > 0 && Array2.dim1(r) <> 4 then
       invalid_arg(name ^ ": dim1 region must be 4");
     r

let pslg ~hole ~region ~point_marker ~point ~segment_marker ~segment =
  check_point "Mesh.pslg" point;
  let point_marker = check_point_marker "Mesh.pslg" point point_marker in
  let segment_marker =
    check_segment_marker "Mesh.pslg" segment segment_marker in
  let hole = check_hole "Mesh.pslg" hole in
  let region = check_region "Mesh.pslg" region in
  (object
      method point = point
      method point_marker = point_marker
      method segment = segment
      method segment_marker = segment_marker
      method hole = hole
      method region = region
    end : fortran_layout pslg)

(* Similar to [make_mesh] but with some elementary checks. *)
let create ~hole ~region ~point_marker ~point ~segment_marker ~segment
      ~neighbor ~edge ~edge_marker ~triangle =
  check_point "Mesh.create" point;
  let point_marker = check_point_marker "Mesh.create" point point_marker in
  let segment = match segment with
    | None -> empty_int_mat2
    | Some s ->
       if Array2.dim2(s) > 0 && Array2.dim1(s) <> 2 then
         invalid_arg "Mesh.create: dim1 segment must be 2";
       s in
  let segment_marker =
    check_segment_marker "Mesh.create" segment segment_marker in
  let hole = check_hole "Mesh.create" hole in
  let region = check_region "Mesh.create" region in
  if Array2.dim2(triangle) = 0 then
    invalid_arg "Mesh.create: triangle cannot be empty";
  if Array2.dim1(triangle) < 3 then
    invalid_arg "Mesh.create: dim1 triangle must be at least 3";
  let neighbor = match neighbor with
    | None -> empty_int_mat3
    | Some nbh ->
       if Array2.dim2(nbh) > 0 then (
         if Array2.dim2(nbh) <> Array2.dim2(triangle) then
           invalid_arg "Mesh.create: dim2 neighbor <> dim2 triangle";
         if Array2.dim1(nbh) <> 3 then
           invalid_arg "Mesh.create: dim1 neighbor <> 3";
       );
       nbh in
  let edge = match edge with
    | None -> empty_int_mat2
    | Some e ->
       if Array2.dim2(e) > 0 && Array2.dim1(e) <> 2 then
         invalid_arg "Mesh.create: dim1 edge <> 2";
       e in
  let edge_marker = match edge_marker with
    | None -> empty_vec
    | Some e ->
       if Array1.dim e > 0 && Array1.dim e <> Array2.dim2(edge) then
         invalid_arg "Mesh.create: dim2 edge_marker <> dim2 edge";
       e in
  (object
     method point = point
     method point_marker = point_marker
     method segment = segment
     method segment_marker = segment_marker
     method hole = hole
     method region = region
     method triangle = triangle
     method neighbor = neighbor
     method edge = edge
     method edge_marker = edge_marker
   end : fortran_layout t)


(** Return the smaller box (xmin, xmax, ymin, ymax) containing the [mesh]. *)
let bounding_box (mesh: mesh) =
  let xmin = ref infinity
  and xmax = ref neg_infinity
  and ymin = ref infinity
  and ymax = ref neg_infinity in
  let point = mesh#point in
  for i = 1 to Array2.dim2(point) do
    let x = point.{1,i}
    and y = point.{2,i} in
    if x > !xmax then xmax := x;
    if x < !xmin then xmin := x;
    if y > !ymax then ymax := y;
    if y < !ymin then ymin := y;
  done;
  (!xmin, !xmax, !ymin, !ymax)

let latex_write ?edge:(edge_color=fun _ -> Some black) (mesh: mesh) fh =
  let edge = mesh#edge in
  let pt = mesh#point in
  if Array2.dim2(edge) = 0 then invalid_arg "Mesh.latex: mesh#edge must be nonempty";
  if Array2.dim1(edge) <> 2 then
    invalid_arg "Mesh.latex: mesh#edge must have 2 rows (fortran)";
  if Array2.dim2(pt) = 0 then invalid_arg "Mesh.latex: mesh#point must be nonempty";
  if Array2.dim1(pt) <> 2 then
    invalid_arg "Mesh.latex: mesh#point must have 2 rows (fortran)";
  let xmin, xmax, ymin, ymax = bounding_box mesh in
  latex_begin fh (xmax -. xmin) (ymax -. ymin) xmin ymin;
  (* Write lines *)
  fprintf fh "  %% %i triangles\n" (Array2.dim2(mesh#triangle));
  for e = 1 to Array2.dim2(edge) do
    match edge_color e with
    | None -> ()
    | Some color ->
      let i1 = edge.{1,e}
      and i2 = edge.{2,e} in
      let p1 = { x = pt.{1,i1};  y = pt.{2,i1} }
      and p2 = { x = pt.{1,i2};  y = pt.{2,i2} } in
      line fh color p1 p2
  done;
  (* Write points *)
  fprintf fh "  %% %i points\n" (Array2.dim2(pt));
  for i = 1 to Array2.dim2(pt) do
    point_xy fh i (pt.{1,i}) (pt.{2,i});
  done;
  latex_end fh

let latex ?edge mesh filename =
  let fh = open_out filename in
  try latex_write ?edge mesh fh;
      close_out fh
  with e -> close_out fh;
           raise e

let scilab (mesh: mesh) ?(longitude=70.) ?(azimuth=60.)
      ?(mode=`Triangles) ?(box=`Full) ?edgecolor
      (z: vec) fname =
  let triangle = mesh#triangle in
  let pt = mesh#point in
  if Array2.dim2(triangle) = 0 then
    invalid_arg "Mesh.scilab: mesh#triangle must be nonempty";
  if Array2.dim1(triangle) < 3 then
    invalid_arg "Mesh.scilab: mesh#triangle must have at least \
	         3 rows (fortran)";
  if Array2.dim2(pt) = 0 then invalid_arg "Mesh.scilab: mesh#point must be nonempty";
  if Array2.dim1(pt) <> 2 then
    invalid_arg "Mesh.scilab: mesh#point must have 2 rows (fortran)";
  if Array1.dim z < Array2.dim2(pt) then
    invalid_arg "Mesh.scilab: vector too small";
  let fname =
    if Filename.check_suffix fname ".sci" then Filename.chop_extension fname
    else fname in
  let sci = fname ^ ".sci"
  and xf = fname ^ "_x.dat"
  and yf = fname ^ "_y.dat"
  and zf = fname ^ "_z.dat" in
  let mode = match mode with
    | `Triangles -> 1
    | `Triangles_only -> 0
    | `No_triangles -> -1 in
  let box = match box with
    | `None -> 0
    | `Behind -> 2
    | `Box_only -> 3
    | `Full -> 4 in
  let edgecolor, er, eg, eb = match edgecolor with
    | None -> false, 0., 0., 0.
    | Some(`Color c) ->
       (true,
        float((c lsr 16) land 0xFF) /. 255.,
        float((c lsr 8) land 0xFF) /. 255.,
        float(c land 0xFF) /. 255.)
    | Some(`Grey c) ->
       if c <= 0. || c > 1. then (false, 0., 0., 0.)
       else (true, c, c, c) in
  let fh = open_out sci in
  (* Put the edge color at the bottom of the colormap so it is usually
     hidden.  Moreover, put enough color in the map so the edge color
     is seldom drawn. *)
  fprintf fh "mode(0);\n\
              // Run in Scilab with: exec('%s')\n\
              // Written by the OCaml Mesh module (version 0.9.5).\n\
              // mesh: %i triangles, %i points.\n\
              ocaml = struct('f', scf(), 'e', null, \
                             'x', fscanfMat('%s'), 'y', fscanfMat('%s'), \
                             'z', fscanfMat('%s'));\n\
              clf();\n\
              ocaml.e = gce();\n\
              ocaml.e.hiddencolor = -1;\n\
              ocaml.f.color_map = jetcolormap(100);\n"
    sci (Array2.dim2(triangle)) (Array2.dim2(pt))
    (Filename.basename xf) (Filename.basename yf) (Filename.basename zf);
  if edgecolor && mode >= 0 then
    fprintf fh "ocaml.f.color_map(1,:) = [%g, %g, %g];\n\
                xset('color', 1);\n"
      er eg eb;
  fprintf fh "plot3d1(ocaml.x, ocaml.y, ocaml.z, theta=%g, alpha=%g, \
                flag=[%d,2,%d]);\n\
              disp('Save: xs2pdf(ocaml.f, ''%s.pdf'')');\n"
    longitude azimuth mode box (Filename.basename fname);
  close_out fh;
  let save_mat fname coord =
    let fh = open_out fname in
    (* We traverse several times the triangles but Scilab will not
       have to transpose the matrices. *)
    for point = 1 to 3 do
      for t = 1 to Array2.dim2(triangle) do
        fprintf fh "%.16e " (coord (triangle.{point,t}))
      done;
      fprintf fh "\n";
    done;
    close_out fh in
  save_mat xf (fun i -> pt.{1,i});
  save_mat yf (fun i -> pt.{2,i});
  save_mat zf (fun i -> z.{i})

let is_allowed_matlab c =
  ('0' <= c && c <= '9') || ('a' <= c && c <= 'z') || ('A' <= c && c <= 'Z')
  || c = '_'

let matlab (mesh: mesh) ?(edgecolor=`Color 0) ?(linestyle="-") ?(facealpha=1.)
           (z: vec) fname =
  let tr = mesh#triangle in
  let pt = mesh#point in
  if Array2.dim2(tr) = 0 then
    invalid_arg "Mesh.matlab: mesh#triangle must be nonempty";
  if Array2.dim1(tr) < 3 then
    invalid_arg "Mesh.matlab: mesh#triangle must have at least \
	         3 rows (fortran)";
  if Array2.dim2(pt) = 0 then invalid_arg "Mesh.matlab: mesh#point must be nonempty";
  if Array2.dim1(pt) <> 2 then
    invalid_arg "Mesh.matlab: mesh#point must have 2 rows (fortran)";
  let dir = Filename.dirname fname
  and base = Filename.basename fname in
  let base =
    if Filename.check_suffix base ".m" then
      Bytes.unsafe_of_string(String.sub base 0 (String.length base - 2))
    else Bytes.of_string base in
  (* Matlab filenames can contain only alphanumeric characters and
     underscores.  Convert all other characters to underscore *)
  for i = 0 to Bytes.length base - 1 do
    if not(is_allowed_matlab (Bytes.get base i)) then Bytes.set base i '_'
  done;
  let base = Bytes.unsafe_to_string base in
  let mat = Filename.concat dir (base ^ ".m") in
  let save_xy fh coord =
    for p = 1 to Array2.dim2(pt) do fprintf fh "%.13g " (pt.{coord,p}) done;
    fprintf fh "\n" in
  let fh = open_out mat in
  fprintf fh "%% Run in Matlab with: run %s\n\
              %% Created by the OCaml Mesh module (version 0.9.5).\n\
              %% print -painters -dpdf -r600 %s.pdf\n" mat base;
  fprintf fh "mesh_x = [" ;
  save_xy fh 1;
  fprintf fh "];\nmesh_y = [";
  save_xy fh 2;
  fprintf fh "];\nmesh_z = [";
  for i = 1 to Array1.dim(z) do fprintf fh "%.13f " z.{i} done;
  fprintf fh "];\nmesh_triangles = [";
  for t = 1 to Array2.dim2(tr) do
    fprintf fh "%i %i %i; " (tr.{1,t}) (tr.{2,t}) (tr.{3,t})
  done;
  let edgecolor = match edgecolor with
    | `None -> "'none'"
    | `Flat -> "'flat'"
    | `Interp -> "'interp'"
    | `Color c ->
       if c < 0 then "'none'"
       else let b = float(c land 0xFF) /. 255.
            and g = float((c lsr 8) land 0xFF) /. 255.
            and r = float((c lsr 16) land 0xFF) /. 255. in
            sprintf "[%g,%g,%g]" r g b in
  let facealpha = if facealpha < 0. then 0.
                  else if facealpha > 1. then 1.
                  else facealpha in
  (* FIXME: protect against strings containing "'". *)
  fprintf fh "];\ntrisurf(mesh_triangles, mesh_x, mesh_y, mesh_z, \
              'FaceAlpha', %f, 'EdgeColor', %s, 'LineStyle', '%s');\n"
          facealpha edgecolor linestyle;
  close_out fh
;;


(* Sort the vertices at node [n0] by increasing (counterclockwise)
   angle w.r.t. the base vertex [i0].  [TriangularSurfacePlot] (not
   [PlanarGraphPlot] it seems) requires the vertices to be ordered. *)
let sort_counterclockwise (pt: mat) n0 = function
  | ([] | [_]) as adj -> adj
  | n1 :: tl ->
    let x0 = pt.{1, n0} and y0 = pt.{2, n0} in
    let dx1 = pt.{1, n1} -. x0 and dy1 = pt.{2, n1} -. y0 in
    (* Since [atan2] returns an angle in ]-pi, pi], the angle of
       (dx1,dy1) will be set to pi so that the order given by the
       angles is correct.  Also there is no need to norm the vectors
       [(dx1,dy1)] and [(dx,dy)] because that will only dilate
       [(e1,e2)] which does not change the value of [atan2]. *)
    let angle n =
      let dx = pt.{1, n} -. x0 and dy = pt.{2, n} -. y0 in
      let e1 = -. dx *. dx1 -. dy *. dy1
      and e2 = dx *. dy1 -. dy *. dx1 in
      atan2 e2 e1 in
    (* Add angles *)
    let tl = List.map (fun n -> (n, angle n)) tl in
    let tl = List.fast_sort (fun (_,a1) (_,a2) -> compare a1 a2) tl in
    n1 :: List.map (fun (n,_) -> n) tl
;;

(* Return an array [adj] such that [adj.(i)] is the list of the
   adjacent nodes to [i]. *)
let adjacency (mesh: mesh) =
  let pt = mesh#point in
  let n = Array2.dim2(pt) in
  let adj = Array.make (n + 1) [] in
  let edge = mesh#edge in
  for e = 1 to Array2.dim2(edge) do
    let i1 = edge.{1,e}
    and i2 = edge.{2,e} in
    adj.(i1) <- i2 :: adj.(i1);
    adj.(i2) <- i1 :: adj.(i2);
  done;
  (* This is important for TriangularSurfacePlot (that uses the order
     for orientation?).  *)
  Array.mapi (fun n0 adj -> sort_counterclockwise pt n0 adj) adj

let is_allowed_mathematica c =
  ('0' <= c && c <= '9') || ('a' <= c && c <= 'z') || ('A' <= c && c <= 'Z')

let count_mathematica_allowed base =
  let n = ref 0 in
  for i = 0 to String.length base - 1 do
    if is_allowed_mathematica (String.unsafe_get base i) then incr n
  done;
  !n

(* Remove all chars that are not alphanumeric. *)
let mathematica_safe base =
  let len = count_mathematica_allowed base in
  if len = String.length base then base
  else (
    let base' = Bytes.create len in
    let j = ref 0 in
    for i = 0 to String.length base - 1 do
      let c = String.unsafe_get base i in
      if is_allowed_mathematica c then (
        Bytes.set base' !j c;
        incr j;
    )
    done;
    Bytes.unsafe_to_string base'
  )

let mathematica_print_float fh f =
  let s = Bytes.unsafe_of_string(sprintf "%.16g" f) in
  try
    let e = Bytes.index s 'e' in
    output fh s 0 e;  output_string fh "*^";
    output fh s (e + 1) (Bytes.length s - e - 1)
  with Not_found ->
    output fh s 0 (Bytes.length s)

let mathematica (mesh: mesh) (z: vec) fname =
  let pt = mesh#point in
  if Array2.dim2(pt) = 0 then
    invalid_arg "Mesh.mathematica: mesh#point must be nonempty";
  if Array2.dim1(pt) <> 2 then
    invalid_arg "Mesh.mathematica: mesh#point must have 2 rows (fortran)";
  if Array2.dim2(mesh#edge) = 0 then
    invalid_arg "Mesh.mathematica: mesh#edge must be nonempty";
  if Array2.dim1(mesh#edge) <> 2 then
    invalid_arg "Mesh.mathematica: mesh#edge must have 2 rows (fortran)";
  let base = Filename.basename fname in
  let pkg, fname =
    if Filename.check_suffix base ".m" then
      mathematica_safe(String.sub base 0 (String.length base - 2)), fname
    else mathematica_safe base, fname ^ ".m" in
  let pkg = String.capitalize_ascii pkg in
  let fh = open_out fname in
  fprintf fh "(* Created by the OCaml Mesh module (version 0.9.5)) \
              *)\n";
  fprintf fh "%s`xyz = {" pkg;
  output_string fh "{";
  mathematica_print_float fh pt.{1, 1};  output_string fh ", ";
  mathematica_print_float fh pt.{2, 1};  output_string fh ", ";
  mathematica_print_float fh z.{1};        output_string fh "}";
  for i = 1 + 1 to Array2.dim2(pt) do
    output_string fh ", {";
    mathematica_print_float fh pt.{1, i};  output_string fh ", ";
    mathematica_print_float fh pt.{2, i};  output_string fh ", ";
    mathematica_print_float fh z.{i};        output_string fh "}"
  done;
  fprintf fh "};\n\n";
  let adj = adjacency mesh in
  let output_adj i =
    (* mathematica indices start at 1 *)
    match adj.(i) with
    | [] -> fprintf fh "{%i, {}}" (i)
    | n :: tl ->
      fprintf fh "{%i, {%i" (i) (n);
      List.iter (fun n -> fprintf fh ", %i" (n)) tl;
      fprintf fh"}}" in
  fprintf fh "%s`adj = {" pkg;
  output_adj 1;
  for i = 1 + 1 to Array.length adj - 1 do
    output_string fh ", "; output_adj i
  done;
  fprintf fh "};\n\n";
  fprintf fh "Needs[\"ComputationalGeometry`\"];\n";
  fprintf fh "TriangularSurfacePlot[%s`xyz, %s`adj, Axes -> True]\n" pkg pkg;
  close_out fh
;;

(************************************************************************)
(* mesh_level_curvesF.ml included by "make_FC_code.ml" with Mesh = "Mesh". *)
(* Generic code to draw level cuves.  To be included in a file that
   defines the drawing primitives. *)

module M = Map.Make(struct
                      type t = int
                      let compare x y = compare (x:int) y
                    end)

(* Module to build a structure helping to determine when the segment
   joining 2 points are on the boundary. *)
module Edge =
struct
  let make() = ref M.empty

  let add_edge t i1 i2 =
    assert(i1 < i2);
    try
      let v = M.find i1 !t in
      v := i2 :: !v
    with Not_found ->
      t := M.add i1 (ref [i2]) !t

  (* Declare the segment joining the points of indexes [i1] and [i2]
     as being part of the boundary.   It is auusmed that [i1 <> i2]. *)
  let add t i1 i2 =
    if i1 < i2 then add_edge t i1 i2 else add_edge t i2 i1

  let on_boundary t i1 i2 =
    assert(i1 < i2);
    try
      let v = M.find i1 !t in List.mem i2 !v
    with Not_found -> false

  (* Tells whether the segment (if any) joining the points of indices
     [i1] and [i2] is on the boundary (according to the information in
     [t]).  It is assumed that [i1 <> i2]. *)
  let on t i1 i2 =
    if i1 < i2 then on_boundary t i1 i2 else on_boundary t i2 i1
end;;

let default_level_eq l1 l2 =
  abs_float(l1 -. l2) <= 1E-8 *. (abs_float l1 +. abs_float l2)

let mid p q = {x = 0.5 *. (p.x +. q.x);  y = 0.5 *. (p.y +. q.y) }

(* Intersection of the curve et level [l] and the line passing through
   (x1,y1) and (x2,y2).  [z1 <> z2] assumed. *)
let intercept {x=x1; y=y1} z1 {x=x2; y=y2} z2 l =
  let d = z1 -. z2 and a = l -. z2 and b = z1 -. l in
  {x = (a *. x1 +. b *. x2) /. d;  y = (a *. y1 +. b *. y2) /. d }

let draw_levels ~boundary (mesh: mesh) (z: vec)
    ?(level_eq=default_level_eq) levels surf =
  let edge = mesh#edge in
  let marker = mesh#edge_marker in
  let pt = mesh#point in
  if Array2.dim2(edge) = 0 then
    invalid_arg("Mesh.level_curves: mesh#edge must be nonempty");
  if Array2.dim1(edge) <> 2 then
    invalid_arg("Mesh.level_curves: mesh#edge must have 2 rows (fortran)");
  if Array1.dim marker < Array2.dim2(edge) then
    invalid_arg("Mesh.level_curves: dim mesh#edge_marker < number edges");
  if Array2.dim2(pt) = 0 then
    invalid_arg("Mesh.level_curves: mesh#point must be nonempty");
  if Array2.dim1(pt) <> 2 then
    invalid_arg("Mesh.level_curves: mesh#point must have 2 rows (fortran)");
  let bd = Edge.make() in
  (* Draw the boundary edges *)
  for e = 1 to Array2.dim2(edge) do
    let m = marker.{e} in
    if m <> 0 (* not an interior point *) then begin
      let i1 = edge.{1,e}
      and i2 = edge.{2,e} in
      Edge.add bd i1 i2; (* collect boundary points *)
      match boundary m with
      | None -> ()
      | Some color ->
          let p1 = { x = pt.{1,i1};  y = pt.{2,i1} }
          and p2 = { x = pt.{1,i2};  y = pt.{2,i2} } in
          line surf color p1 p2
    end
  done;
  let tr = mesh#triangle in
  if Array2.dim2(tr) = 0 then
    invalid_arg("Mesh.level_curves: mesh#triangle must be nonempty");
  if Array2.dim1(tr) < 3 then
    invalid_arg("Mesh.level_curves: mesh#triangle must have at least 3 \
      rows (fortran) or 3 columns (C)");
  let marker = mesh#point_marker in
  for t = 1 to Array2.dim2(tr) do
    let i1 = tr.{1,t}
    and i2 = tr.{2,t}
    and i3 = tr.{3,t} in
    let p1 = { x = pt.{1,i1};  y = pt.{2,i1} }
    and z1 = z.{i1} in
    let p2 = { x = pt.{1,i2};  y = pt.{2,i2} }
    and z2 = z.{i2} in
    let p3 = { x = pt.{1,i3};  y = pt.{2,i3} }
    and z3 = z.{i3} in
    List.iter
      (fun (l, color) ->
         (* Draw the level curve [l] on the triangle [t] except if
            that curve is on the boundary. *)
         if level_eq l z1 then (
           if level_eq l z2 then (
             if level_eq l z3 then
               (* The entire triangle is at the same level.  Try to
                  remove boundary edges. *)
               if Edge.on bd i1 i2 then
                 if Edge.on bd i1 i3 || Edge.on bd i2 i3 then
                   triangle surf color p1 p2 p3 (* Full triangle ! *)
                 else line surf color p3 (mid p1 p2)
               else (* i1-i2 not on boundary *)
                 if Edge.on bd i1 i3 then
                   if Edge.on bd i2 i3 then triangle surf color p1 p2 p3
                   else line surf color p2 (mid p1 p3)
                 else (* i1-i3 not on boundary *)
                   if Edge.on bd i2 i3 then line surf color p1 (mid p2 p3)
                   else triangle surf color p1 p2 p3 (* Full triangle ! *)
             else (* l = z1 = z2 <> z3 *)
               if not(Edge.on bd i1 i2) then line surf color p1 p2
           )
           else (* l = z1 <> z2 *)
             if level_eq l z3 then (* l = z1 = z3 <> z2 *)
               (if not(Edge.on bd i1 i3) then line surf color p1 p3)
             else
               if (z2 < l && l < z3) || (z3 < l && l < z2) then
                 line surf color p1 (intercept p2 z2 p3 z3 l)
         )
         else if l < z1 then (
           if level_eq l z2 then
             if level_eq l z3 then
               (if not(Edge.on bd i2 i3) then line surf color p2 p3)
             else if l > z3 then (* z3 < l = z2 < z1 *)
               line surf color p2 (intercept p1 z1 p3 z3 l)
             else (* corner point, inside the domain.  Ususally this
                     happens because the level line passes through a
                     triangle corner. *)
               (if marker.{i2} = 0 then point surf i2 p2)
           else if l < z2 then (
             if level_eq l z3 then
               (if marker.{i3} = 0 then point surf i3 p3)
             else if l > z3 then
               line surf color (intercept p1 z1 p3 z3 l)
                 (intercept p2 z2 p3 z3 l)
           )
           else (* z2 < l < z1 *)
             line surf color (intercept p1 z1 p2 z2 l)
               (if level_eq l z3 then p3
                else if l < z3 then intercept p2 z2 p3 z3 l
                else (* l > z3 *)   intercept p1 z1 p3 z3 l)
         )
         else (* l > z1 *) (
           (* Symmetric of [l < z1] with all inequalities reversed *)
           if level_eq l z2 then
             if level_eq l z3 then
               (if not(Edge.on bd i2 i3) then line surf color p2 p3)
             else if l < z3 then (* z1 < l = z2 < z3 *)
               line surf color p2 (intercept p1 z1 p3 z3 l)
             else (* corner point, inside the domain *)
               (if marker.{i2} = 0 then point surf i2 p2)
           else if l > z2 then (
             if level_eq l z3 then
               (if marker.{i3} = 0 then point surf i3 p3)
             else if l < z3 then
               line surf color (intercept p1 z1 p3 z3 l)
                 (intercept p2 z2 p3 z3 l)
           )
           else (* z1 < l < z2 *)
             line surf color (intercept p1 z1 p2 z2 l)
               (if level_eq l z3 then p3
                else if l > z3 then intercept p2 z2 p3 z3 l
                else (* l < z3 *)   intercept p1 z1 p3 z3 l)
         )
      ) levels
  done
;;

type polygon_fill =
| Tri123 (* triangle with edge 1 and cut in edges 2, 3 *)
| Tri231 | Tri312
| Quad123 (* Quadrilateral with edges 1-2 and 1-3 of the triangle cut *)
| Quad231 | Quad312
| Whole | Empty;;

(* base 3: c1 + 1 + 3(c2 + 1) + 9(c3 + 1).  The [c1], [c2] and [c3]
   are the comparisons of the 3 corners with the desired level. *)
let index c1 c2 c3 = c1 + 3 * c2 + 9 * c3 + 13

let super =
  let d = Array.make 27 Empty in
  d.(index( 1) ( 1) ( 1)) <- Whole;
  d.(index( 1) ( 1) ( 0)) <- Whole;
  d.(index( 1) ( 1) (-1)) <- Quad312;
  d.(index( 1) ( 0) ( 1)) <- Whole;
  d.(index( 1) ( 0) ( 0)) <- Whole;
  d.(index( 1) ( 0) (-1)) <- Tri123;
  d.(index( 1) (-1) ( 1)) <- Quad231;
  d.(index( 1) (-1) ( 0)) <- Tri123;
  d.(index( 1) (-1) (-1)) <- Tri123;
  d.(index( 0) ( 1) ( 1)) <- Whole;
  d.(index( 0) ( 1) ( 0)) <- Whole;
  d.(index( 0) ( 1) (-1)) <- Tri231;
  d.(index( 0) ( 0) ( 1)) <- Whole;
  d.(index( 0) ( 0) ( 0)) <- Empty; (* > 0 required *)
  d.(index( 0) ( 0) (-1)) <- Empty;
  d.(index( 0) (-1) ( 1)) <- Tri312;
  d.(index( 0) (-1) ( 0)) <- Empty;
  d.(index( 0) (-1) (-1)) <- Empty;
  d.(index(-1) ( 1) ( 1)) <- Quad123;
  d.(index(-1) ( 1) ( 0)) <- Tri231;
  d.(index(-1) ( 1) (-1)) <- Tri231;
  d.(index(-1) ( 0) ( 1)) <- Tri312;
  d.(index(-1) ( 0) ( 0)) <- Empty;
  d.(index(-1) ( 0) (-1)) <- Empty;
  d.(index(-1) (-1) ( 1)) <- Tri312;
  d.(index(-1) (-1) ( 0)) <- Empty;
  d.(index(-1) (-1) (-1)) <- Empty;
  d

let sub =
  let d = Array.make 27 Empty in
  d.(index( 1) ( 1) ( 1)) <- Empty;
  d.(index( 1) ( 1) ( 0)) <- Empty;
  d.(index( 1) ( 1) (-1)) <- Tri312;
  d.(index( 1) ( 0) ( 1)) <- Empty;
  d.(index( 1) ( 0) ( 0)) <- Empty;
  d.(index( 1) ( 0) (-1)) <- Tri312;
  d.(index( 1) (-1) ( 1)) <- Tri231;
  d.(index( 1) (-1) ( 0)) <- Tri231;
  d.(index( 1) (-1) (-1)) <- Quad123;
  d.(index( 0) ( 1) ( 1)) <- Empty;
  d.(index( 0) ( 1) ( 0)) <- Empty;
  d.(index( 0) ( 1) (-1)) <- Tri312;
  d.(index( 0) ( 0) ( 1)) <- Empty;
  d.(index( 0) ( 0) ( 0)) <- Empty; (* < 0 required *)
  d.(index( 0) ( 0) (-1)) <- Whole;
  d.(index( 0) (-1) ( 1)) <- Tri231;
  d.(index( 0) (-1) ( 0)) <- Whole;
  d.(index( 0) (-1) (-1)) <- Whole;
  d.(index(-1) ( 1) ( 1)) <- Tri123;
  d.(index(-1) ( 1) ( 0)) <- Tri123;
  d.(index(-1) ( 1) (-1)) <- Quad231;
  d.(index(-1) ( 0) ( 1)) <- Tri123;
  d.(index(-1) ( 0) ( 0)) <- Whole;
  d.(index(-1) ( 0) (-1)) <- Whole;
  d.(index(-1) (-1) ( 1)) <- Quad312;
  d.(index(-1) (-1) ( 0)) <- Whole;
  d.(index(-1) (-1) (-1)) <- Whole;
  d

let draw_xxx_level decision name ?(boundary=(fun _ -> Some black))
    (mesh: mesh) (z: vec) l color surf =
  let edge = mesh#edge in
  let edge_marker = mesh#edge_marker in
  let pt = mesh#point in
  if Array2.dim2(edge) = 0 then
    invalid_arg("Mesh" ^ name ^ ": mesh#edge must be nonempty");
  if Array2.dim1(edge) <> 2 then
    invalid_arg("Mesh" ^ name ^ ": mesh#edge must have 2 rows (fortran)");
  if Array1.dim edge_marker < Array2.dim2(edge) then
    invalid_arg("Mesh" ^ name ^ ": dim mesh#edge_marker < number edges");
  if Array2.dim2(pt) = 0 then
    invalid_arg("Mesh" ^ name ^ ": mesh#point must be nonempty");
  if Array2.dim1(pt) <> 2 then
    invalid_arg("Mesh" ^ name ^ ": mesh#point must have 2 rows (fortran)");
  let tr = mesh#triangle in
  if Array2.dim2(tr) = 0 then
    invalid_arg("Mesh" ^ name ^ ": mesh#triangle must be nonempty");
  if Array2.dim1(tr) < 3 then
    invalid_arg("Mesh" ^ name ^ ": mesh#triangle must have at least 3 \
      rows (fortran) or 3 columns (C)");
  for t = 1 to Array2.dim2(tr) do
    let i1 = tr.{1,t}
    and i2 = tr.{2,t}
    and i3 = tr.{3,t} in
    let p1 = { x = pt.{1,i1};  y = pt.{2,i1} }
    and z1 = z.{i1} in
    let p2 = { x = pt.{1,i2};  y = pt.{2,i2} }
    and z2 = z.{i2} in
    let p3 = { x = pt.{1,i3};  y = pt.{2,i3} }
    and z3 = z.{i3} in
    match decision.(index (compare z1 l) (compare z2 l) (compare z3 l)) with
    | Tri123 -> fill_triangle surf color p1 (intercept p1 z1 p2 z2 l)
                                           (intercept p1 z1 p3 z3 l)
    | Tri231 -> fill_triangle surf color p2 (intercept p2 z2 p3 z3 l)
                                           (intercept p2 z2 p1 z1 l)
    | Tri312 -> fill_triangle surf color p3 (intercept p3 z3 p1 z1 l)
                                           (intercept p3 z3 p2 z2 l)
    | Quad123 -> fill_quadrilateral surf color (intercept p1 z1 p2 z2 l)
                                              (intercept p1 z1 p3 z3 l) p3 p2
    | Quad231 -> fill_quadrilateral surf color (intercept p2 z2 p3 z3 l)
                                              (intercept p2 z2 p1 z1 l) p1  p3
    | Quad312 -> fill_quadrilateral surf color (intercept p3 z3 p1 z1 l)
                                              (intercept p3 z3 p2 z2 l) p2 p1
    | Whole -> fill_triangle surf color p1 p2 p3
    | Empty -> ()
  done;
  (* Draw the boundary edges (over the filled area) *)
  for e = 1 to Array2.dim2(edge) do
    let m = edge_marker.{e} in
    if m <> 0 (* not an interior point *) then begin
      match boundary m with
      | None -> ()
      | Some color ->
        let i1 = edge.{1,e}
        and i2 = edge.{2,e} in
        let p1 = { x = pt.{1,i1};  y = pt.{2,i1} }
        and p2 = { x = pt.{1,i2};  y = pt.{2,i2} } in
        line surf color p1 p2
    end
  done

let draw_super_level ?boundary mesh z level color surf =
  draw_xxx_level super ".super_level" ?boundary mesh z level color surf

let draw_sub_level ?boundary mesh z level color surf =
  draw_xxx_level sub ".sub_level" ?boundary mesh z level color surf
;;
(************************************************************************)

let level_curves ?(boundary=(fun _ -> Some black)) (mesh: mesh) (z: vec)
    ?level_eq levels fname =
  let fh = open_out fname in
  let xmin, xmax, ymin, ymax = bounding_box mesh in
  latex_begin fh (xmax -. xmin) (ymax -. ymin) xmin ymin;
  draw_levels ~boundary mesh z ?level_eq levels fh;
  latex_end fh;
  close_out fh

let super_level ?boundary (mesh: mesh) (z: vec) level color fname =
  let fh = open_out fname in
  let xmin, xmax, ymin, ymax = bounding_box mesh in
  latex_begin fh (xmax -. xmin) (ymax -. ymin) xmin ymin;
  draw_super_level ?boundary mesh z level color fh;
  latex_end fh;
  close_out fh

let sub_level ?boundary (mesh: mesh) (z: vec) level color fname =
  let fh = open_out fname in
  let xmin, xmax, ymin, ymax = bounding_box mesh in
  latex_begin fh (xmax -. xmin) (ymax -. ymin) xmin ymin;
  draw_sub_level ?boundary mesh z level color fh;
  latex_end fh;
  close_out fh


(* Determine the number of superdiagonals + 1 main diagonal *)
let band_height_P1 filter (mesh: mesh) =
  let tr = mesh#triangle in
  let kd = ref 0 in
  match filter with
  | None ->
    for t = 1 to Array2.dim2(tr) do
      let i1 = tr.{1,t}
      and i2 = tr.{2,t}
      and i3 = tr.{3,t} in
      kd := max4 !kd (abs(i1 - i2)) (abs(i2 -i3)) (abs(i3 - i1))
    done;
    !kd + 1
  | Some cond ->
    for t = 1 to Array2.dim2(tr) do
      let i1 = tr.{1,t}
      and i2 = tr.{2,t}
      and i3 = tr.{3,t} in
      if cond i1 then (
        if cond i2 then
          if cond i3 then
            kd := max4 !kd (abs(i1 - i2)) (abs(i2 -i3)) (abs(i3 - i1))
          else (* exlude i3 *)
            kd := max2 !kd (abs(i2 - i1))
        else (* exclude i2 *) if cond i3 then
          kd := max2 !kd (abs(i3 - i1))
      )
      else (* exclude i1 *) if cond i2 && cond i3 then
        kd := max2 !kd (abs(i3 - i2))
    done;
    !kd + 1

(* Return the index with the lowest nonnegative [deg] (negative
   degrees are ignored).  Return [-1] if all degrees are < 0. *)
let min_deg (deg: int array) =
  let i = ref(-1) in
  let degi = ref(max_int) in
  for j = 1 to Array.length deg - 1 do
    if deg.(j) >= 0 && deg.(j) < !degi then (i := j;  degi := deg.(j))
  done;
  !i


(* sub
 ***********************************************************************)

(* Iterator with indices adapted to the current layout. *)
let rec iteri f i = function
  | [] -> ()
  | x :: tl -> f i x;  iteri f (succ i) tl

let iteri f l = iteri f 1 l


let filter_columns_shift (m: int_mat) select shift =
  let cols = ref [] in
  let nselected = ref 0 in (* length of [cols] *)
  for c = 1 to Array2.dim2(m) do
    if select m c then (cols := c :: !cols;  incr nselected)
  done;
  let cols = List.rev !cols in
  let m' = Array2.create (int) fortran_layout (Array2.dim1(m)) (!nselected) in
  iteri (fun i pi ->
         for j = 1 to Array2.dim1(m') do
           m'.{j,i} <- m.{j,pi} - shift
         done
        ) cols;
  m', !nselected, cols

let sub_markers (v: int_vec) n cols =
  if Array1.dim v = 0 then v (* no markers *)
  else (
    let v' = Array1.create (int) fortran_layout (n) in
    iteri (fun i pi ->  v'.{i} <- v.{pi}) cols;
    v'
  )

let internal_sub (mesh: fortran_layout #t) ?pos len =
  let pos = match pos with
    | None -> 1
    | Some pos ->
       if pos < 1 then invalid_arg "Mesh.sub: pos < 1";
       pos in
  if len <= 0 then invalid_arg "Mesh.sub: len <= 0";
  if pos + len > Array2.dim2(mesh#point) then
    invalid_arg "Mesh.sub: len too large";
  let shift = pos - 1 in
  let max_point_idx = pos + len - 1 in
  let sub_point i = pos <= i && i <= max_point_idx in
  (* Points *)
  let point = Array2.sub_right mesh#point pos len in
  let point_marker = Array1.sub mesh#point_marker pos len in
  (* Segments *)
  let select2 (m: int_mat) i = sub_point m.{1,i}
                               && sub_point m.{2,i} in
  let new_seg, n, cols = filter_columns_shift mesh#segment select2 shift in
  let new_seg_marker = sub_markers mesh#segment_marker n cols in
  (* Triangles *)
  let select3 (m: int_mat) t = sub_point m.{1,t} && sub_point m.{2,t}
                               && sub_point m.{3,t} in
  let new_tr, n_tr, cols_tr = filter_columns_shift mesh#triangle
                                                   select3 shift in
  (* Neighbors corresponding to the selected triangles. *)
  let new_neighbor =
    let old_nbh = mesh#neighbor in
    if Array2.dim2(old_nbh) = 0 then old_nbh
    else (
      let nbh = Array2.create (int) fortran_layout (3) (n_tr) in (* new neighbor *)
      let trans = Array1.create (int) fortran_layout (Array2.dim2(mesh#triangle)) in (* old idx → new *)
      Array1.fill trans (-1); (* default: no corresponding index *)
      iteri (fun i pi -> trans.{pi} <- i) cols_tr;
      iteri (fun i pi ->
             nbh.{1,i} <- trans.{old_nbh.{1,pi}};
             nbh.{2,i} <- trans.{old_nbh.{2,pi}};
             nbh.{3,i} <- trans.{old_nbh.{3,pi}};
            ) cols_tr;
      nbh
    ) in
  (* Edges *)
  let new_edge, n, cols = filter_columns_shift mesh#edge select2 shift in
  let new_edge_marker = sub_markers mesh#edge_marker n cols in
  (make_mesh
     ~point: point
     ~point_marker: point_marker
     ~segment: new_seg
     ~segment_marker: new_seg_marker
     ~hole: mesh#hole (* keep *)
     ~region: mesh#region (* keep *)
     ~triangle: new_tr
     ~neighbor: new_neighbor
     ~edge: new_edge
     ~edge_marker: new_edge_marker,
   n_tr, cols_tr)


let sub (mesh: mesh) ?pos len =
  let m, _, _ = internal_sub mesh ?pos len in
  m


(* Permutations
 ***********************************************************************)

(** Apply the permutation [perm] to the [mesh]. *)
let do_permute_points name (mesh: mesh) (perm: int_vec) (inv_perm: int_vec)
    : mesh =
  (* Build the new mesh *)
  let old_pt = mesh#point in
  let n = Array2.dim2(old_pt) in
  if n <> Array1.dim perm then
    invalid_arg(sprintf "%s: dim2 #point = %i <> dim perm = %i"
                        name n (Array1.dim perm));
  let pt = Array2.create (float64) fortran_layout (2) (n) in
  let last_pt_idx = Array2.dim2(pt) in
  for i = 1 to last_pt_idx do
    let old_i = perm.{i} in
    pt.{1,i} <- old_pt.{1,old_i};
    pt.{2,i} <- old_pt.{2,old_i};
  done;
  let old_ptm = mesh#point_marker in
  let ptm = Array1.create int layout (Array1.dim old_ptm) in
  for i = 1 to Array1.dim(ptm) do ptm.{i} <- old_ptm.{perm.{i}} done;
  let old_seg = mesh#segment in
  let seg = Array2.create (int) fortran_layout (2) (Array2.dim2(old_seg)) in
  for s = 1 to Array2.dim2(seg) do
    let i1 = old_seg.{1,s} in
    if i1 < 1 || i1 > last_pt_idx then
      failwith(sprintf "%s: mesh#segment.{%i} = %i not in [%i..%i]"
                       name s i1 1 last_pt_idx);
    seg.{1,s} <- inv_perm.{i1};
    let i2 = old_seg.{2,s} in
    if i2 < 1 || i2 > last_pt_idx then
      failwith(sprintf "%s: mesh#segment.{%i} = %i not in [%i..%i]"
                       name s i2 1 last_pt_idx);
    seg.{2,s} <- inv_perm.{i2};
  done;
  let old_tr = mesh#triangle in
  let tr = Array2.create (int) fortran_layout (Array2.dim1(old_tr)) (Array2.dim2(old_tr)) in
  for t = 1 to Array2.dim2(tr) do
    for c = 1 to Array2.dim1(tr) do
      tr.{c,t} <- inv_perm.{old_tr.{c,t}}
    done;
  done;
  let old_edge = mesh#edge in
  let edge = Array2.create (int) fortran_layout (2) (Array2.dim2(old_edge)) in
  for e = 1 to Array2.dim2(edge) do
    edge.{1,e} <- inv_perm.{old_edge.{1,e}};
    edge.{2,e} <- inv_perm.{old_edge.{2,e}};
  done;
  make_mesh
    ~point: pt
    ~point_marker: ptm
    ~segment: seg
    ~segment_marker: mesh#segment_marker
    ~hole: mesh#hole
    ~region: mesh#region
    ~triangle: tr
    ~neighbor: mesh#neighbor
    ~edge: edge
    ~edge_marker: mesh#edge_marker


let permute_points_name = "Mesh.permute_points"

let permute_points_unsafe mesh perm =
  let n = Array2.dim2(mesh#point) in
  (* Inverse perm *)
  let inv_perm = Array1.create int layout n in
  for i = 1 to Array1.dim(perm) do inv_perm.{perm.{i}} <- i done;
  do_permute_points permute_points_name mesh perm inv_perm

let inverse_perm name (perm: int_vec) =
  (* Inverse perm and check that [perm] is indeed a permuation. *)
  let inv_perm = Array1.create int layout (Array1.dim perm) in
  Array1.fill inv_perm (-1); (* never an index *)
  let last_el = Array1.dim(perm) in
  for i = 1 to last_el do
    let pi = perm.{i} in
    if pi < 1 || pi > last_el then
      invalid_arg(sprintf "%s: perm.{%i} = %i not in [%i..%i]"
                          name i pi 1 last_el)
    else if inv_perm.{pi} < 0 then inv_perm.{pi} <- i
    else invalid_arg(sprintf "%s: not a permutation (perm.{%i} = %i = \
                              perm.{%i})" name inv_perm.{pi} pi i)
  done;
  inv_perm

let permute_points (mesh: mesh) ~inv perm =
  let inv_perm = inverse_perm permute_points_name perm in
  if inv then do_permute_points permute_points_name mesh inv_perm perm
  else do_permute_points permute_points_name mesh perm inv_perm


let do_permute_triangles name (mesh: mesh) (perm: int_vec) =
  let old_tr = mesh#triangle in
  let n = Array2.dim2(old_tr) in
  if n <> Array1.dim perm then
    invalid_arg(sprintf "%s: dim2 #triangle = %i <> dim perm = %i"
                        name n (Array1.dim perm));
  let tr = Array2.create (int) fortran_layout (Array2.dim1(old_tr)) (n) in
  let last_tr_idx = Array2.dim2(tr) in
  for i = 1 to last_tr_idx do
    for j = 1 to Array2.dim1(tr) do
      tr.{j,i} <- old_tr.{j,perm.{i}}
    done
  done;
  let old_nbh = mesh#neighbor in
  let nbh =
    if Array2.dim2(old_nbh) = 0 then old_nbh
    else (
      if Array2.dim1(old_nbh) <> 3 then
        invalid_arg(sprintf "%s: invalid mesh: ROW #neighbor <> 3" name);
      if n <> Array2.dim2(old_nbh) then
        invalid_arg(sprintf "%s: invalid mesh: COL #neighbor = %i <> \
                             COL #triangle = %i" name (Array2.dim2(old_nbh)) n);
      let nbh = Array2.create (int) fortran_layout (3) (n) in
      for i = 1 to last_tr_idx do
        let old_i = perm.{i} in
        nbh.{1,i} <- old_nbh.{1,old_i};
        nbh.{2,i} <- old_nbh.{2,old_i};
        nbh.{3,i} <- old_nbh.{3,old_i};
      done;
      nbh
    ) in
  make_mesh
    ~point: mesh#point
    ~point_marker: mesh#point_marker
    ~segment: mesh#segment
    ~segment_marker: mesh#segment_marker
    ~hole: mesh#hole
    ~region: mesh#region
    ~triangle: tr
    ~neighbor: nbh
    ~edge: mesh#edge
    ~edge_marker: mesh#edge_marker


let permute_triangles_name = "Mesh.permute_triangles"

let permute_triangles (mesh: mesh) ~inv perm =
  let inv_perm = inverse_perm permute_triangles_name perm in
  if inv then do_permute_triangles permute_triangles_name mesh inv_perm
  else do_permute_triangles permute_triangles_name mesh perm


(* Band
 ***********************************************************************)

(* http://ciprian-zavoianu.blogspot.com/2009/01/project-bandwidth-reduction.html
*)
let cuthill_mckee ~rev perm (mesh: mesh) : mesh =
  let n = Array2.dim2(mesh#point) in
  let perm = match perm with
    | None -> Array1.create int layout n
    | Some p ->
      if Array1.dim p <> n then
        invalid_arg "Mesh.cuthill_mckee: dim perm <> number of points";
      p in
  let deg = Array.make (n + 1) 0 in (* degree of adjacency of each node *)
  let nbh = Array.make (n + 1) [] in (* list of adjacent nodes *)
  let edge = mesh#edge in
  for e = 1 to Array2.dim2(edge) do
    let i1 = edge.{1,e}
    and i2 = edge.{2,e} in
    nbh.(i1) <- i2 :: nbh.(i1);
    deg.(i1) <- deg.(i1) + 1;
    nbh.(i2) <- i1 :: nbh.(i2);
    deg.(i2) <- deg.(i2) + 1;
  done;
  let free = ref(1) in (* first free position in [perm] *)
  let q = Queue.create () in
  let add node =
    perm.{!free} <- node;
    incr free;
    deg.(node) <- -1; (* [i] put in the final vec. *)
    let nbhs = List.filter (fun i -> deg.(i) >= 0) nbh.(node) in
    let nbhs = List.fast_sort (fun i1 i2 -> compare deg.(i1) deg.(i2)) nbhs in
    List.iter (fun i -> Queue.add i q) nbhs
  in
  let last_pt = Array1.dim(perm) in
  while !free <= last_pt do
    add (min_deg deg);
    while not(Queue.is_empty q) do
      let c = Queue.take q in
      if deg.(c) >= 0 then add c
    done
  done;
  if rev then (
    let s = if 1 = 0 then n-1 else n+1 in (* FIXME: cond known at compil. *)
    for i = 1 to n/2 -1 + 1 do
      let t = perm.{i} in
      perm.{i} <- perm.{s-i};
      perm.{s-i} <- t;
    done
  );
  permute_points_unsafe mesh perm

(* A Generalized GPS Algorithm For Reducing The Bandwidth And Profile
   Of A Sparse Matrix, Q. Wang, Y. C. Guo, and X. W. Shi
   http://www.jpier.org/PIER/pier90/09.09010512.pdf *)
let ggps (mesh: mesh) perm : mesh =
  let n = Array2.dim2(mesh#point) in
  let perm = match perm with
    | None -> Array1.create int layout n
    | Some p ->
      if Array1.dim p <> n then
        invalid_arg "Mesh.ggps: dim perm <> number of points";
      p in
  let deg = Array.make (n + 1) 0 in (* degree of adjacency of each node *)
  let edge = mesh#edge in
  for e = 1 to Array2.dim2(edge) do
    let i1 = edge.{1,e}
    and i2 = edge.{2,e} in
    deg.(i1) <- deg.(i1) + 1;
    deg.(i2) <- deg.(i2) + 1;
  done;
  let _v = min_deg deg in
  (* FIXME *)
  permute_points_unsafe mesh perm

(* Local Variables: *)
(* compile-command: "make -k -C .." *)
(* End: *)