Source file mesh_triangleC.ml

# 1 "triangle/mesh_triangleFC.ml"
(* Binding to Triangle for layout c_layout. *)

open Printf
open Bigarray
open Mesh_triangle_common

type layout = Bigarray.c_layout
type mesh = layout t
type mat = layout Mesh.mat
type vec = layout Mesh.vec
type int_mat = layout Mesh.int_mat
type int_vec = layout Mesh.int_vec

let layout = Bigarray.c_layout
let default_switches = "z"

let empty_vec = Array1.create int layout 0
let empty_mat0 = Array2.create (float64) c_layout (0) (0)
let empty_mat2 = Array2.create (float64) c_layout (0) (2)
let empty_mat4 = Array2.create (float64) c_layout (0) (4)

let pslg ~hole ~region ~point_attribute ~point_marker ~point
         ~segment_marker ~segment =
  let point_marker = match point_marker with
    | None -> empty_vec
    | Some m ->
       let n = Array1.dim m in
       if 0 < n && n < Array2.dim1(point) then
         invalid_arg "Mesh_triangle.pslg: point_marker too small";
       m in
  let point_attribute = match point_attribute with
    | None -> empty_mat0
    | Some a ->
       if Array2.dim2(a) > 0 && Array2.dim1(a) <> Array2.dim1(point) then
         invalid_arg "Mesh_triangle.pslg: dim1 point_attribute <> dim1 point";
       a in
  let segment_marker = match segment_marker with
    | None -> empty_vec
    | Some m ->
       let n = Array1.dim m in
       if 0 < n && n < Array2.dim1(segment) then
         invalid_arg "Mesh_triangle.pslg: segment_marker too small";
       m in
  let hole = match hole with
    | None -> empty_mat2
    | Some h ->
       if Array2.dim1(h) > 0 && Array2.dim2(h) <> 2 then
         invalid_arg "Mesh_triangle.pslg: dim2 hole must be 2";
       h in
  let region = match region with
    | None -> empty_mat4
    | Some r ->
       if Array2.dim1(r) > 0 && Array2.dim2(r) <> 4 then
         invalid_arg "Mesh_triangle.pslg: dim2 region must be 4";
       r in
  (object
      method point = point
      method point_marker = point_marker
      method point_attribute = point_attribute
      method segment = segment
      method segment_marker = segment_marker
      method hole = hole
      method region = region
    end : c_layout pslg)


external triangle :
  string ->                        (* options *)
  layout t ->
  vec                             (* trianglearea *)
  -> mat * mat * int_vec * int_mat * mat * int_mat * int_mat * int_vec
    * (* edge *) int_mat * int_vec
    * (* voronoi *) mat * mat * int_mat * mat
  = "triangulate_c_layout"


let empty_vec = Array1.create float64 layout 0 (* not used => global *)

(* check that all C "triexit" have been avoided. *)

let triangulate ?(delaunay=true) ?min_angle ?max_area ?(region_area=false)
    ?max_steiner ?(voronoi=false) ?(edge=true) ?(neighbor=false)
    ?(subparam=false) ?triangle_area ?triunsuitable
    ?(check_finite=true) ?(debug=true) ?verbose
    ~pslg ~refine (mesh: layout t) =
  (* Check points *)
  let point = mesh#point in
  if Array2.dim2(point) <> 2 then invalid_arg("dim2 mesh#point <> 2");
  if Array2.dim1(mesh#point_attribute) > 0
    && Array2.dim1(mesh#point_attribute) < Array2.dim1(point) then
    invalid_arg("dim1 mesh#point_attribute < dim1 mesh#point");
  if Array1.dim mesh#point_marker > 0
    && Array1.dim mesh#point_marker < Array2.dim1(point) then
    invalid_arg("dim mesh#point_marker < dim1 mesh#point");
  if check_finite then (
    (* Check that no point contains NaN (or infinities).  Triangle
       seems to go into an infinite loop with these which can easily
       be confused with other difficulties. *)
    for i = 0 to Array2.dim1(point) do
      if not(is_finite(point.{i,0})) then
        invalid_arg(sprintf "mesh#point.{%i, %i} is not finite"
                            (i)(0));
      if not(is_finite(point.{i,1})) then
        invalid_arg(sprintf "mesh#point.{%i, %i} is not finite"
                            (i)(1));
    done;
  );
  let switches = Buffer.create 20 in
  Buffer.add_string switches default_switches;
  (* Check for PSLG *)
  if pslg then (
    if Array2.dim1(mesh#segment) > 0 then begin
      if Array2.dim2(mesh#segment) <> 2 then invalid_arg("dim2 segment <> 2");
      if Array1.dim mesh#segment_marker > 0
        && Array1.dim mesh#segment_marker < Array2.dim1(mesh#segment) then
        invalid_arg("dim mesh#segment_marker < dim1 mesh#segment");
    end;
    if not refine then (
      let hole = mesh#hole in
      if Array2.dim1(hole) > 0 && Array2.dim2(hole) <> 2 then
        invalid_arg("dim2 hole <> 2");
      let region = mesh#region in
      if Array2.dim1(region) > 0 then (
        if Array2.dim2(region) <> 4 then invalid_arg("dim2 region <> 4");
        Buffer.add_char switches 'A'; (* regional attributes *)
        if region_area then Buffer.add_char switches 'a'; (* area constraint *)
      );
      if check_finite then (
        for i = 0 to Array2.dim1(hole) do
          if not(is_finite(hole.{i,0})) then
            invalid_arg(sprintf "mesh#hole.{%i, %i} is not finite"
                                (i)(0));
          if not(is_finite(hole.{i,1})) then
            invalid_arg(sprintf "mesh#hole.{%i, %i} is not finite"
                                (i)(1));
        done;
        for i = 0 to Array2.dim1(region) do
          for j = 0 to Array2.dim2(region) do
            if not(is_finite(region.{i,j})) then
              invalid_arg(sprintf "mesh#region.{%i, %i} is not finite"
                                  (i)(j));
          done
        done
      )
    );
    Buffer.add_char switches 'p';
    if Array2.dim2(mesh#segment) = 0 || Array2.dim1(mesh#segment) = 0 then
      Buffer.add_char switches 'c';
  );
  (* Check for refinement -- triangles *)
  if refine then (
    if Array2.dim1(mesh#triangle) > 0 then begin
      if Array2.dim2(mesh#triangle) < 3 then
        invalid_arg("dim2 mesh#triangle < 3");
      if Array2.dim2(mesh#triangle_attribute) > 0
        && Array2.dim1(mesh#triangle_attribute) < Array2.dim1(mesh#triangle) then
        invalid_arg("dim1 mesh#triangle_attribute < dim1 mesh#triangle");
    end;
    Buffer.add_char switches 'r';
    (* Check triangle_area *)
    (match triangle_area with
     | Some a ->
        if Array1.dim a < Array2.dim1(mesh#triangle) then
          invalid_arg("dim triangle_area < dim1 mesh#triangle");
        Buffer.add_char switches 'a';
     | None -> ());
  );
  (* Area constraints *)
  (match max_area with
   | None -> ()
   | Some a -> bprintf switches "a%f" a);
  let triangle_area = match triangle_area with
    | None -> empty_vec
    | Some a -> a (* for refinement only *) in
  (* Check for a triunsuitable function *)
  (match triunsuitable with
  | None -> ()
  | Some f -> register_triunsuitable f;  Buffer.add_char switches 'u');
  (* Other switches *)
  if delaunay then Buffer.add_char switches 'D';
  (match min_angle with
  | None -> ()
  | Some a ->
    if a < 0. || a > 60. then (* required: 3 min_algle <= 180 *)
      Buffer.add_char switches 'q'
    else
      (* Angle may include a decimal point, but not exponential notation. *)
      bprintf switches "d%f" a);
  (match max_steiner with
   | None -> ()
   | Some a -> bprintf switches "S%i" a);
  if voronoi then Buffer.add_char switches 'v';
  if edge then Buffer.add_char switches 'e';
  if neighbor then Buffer.add_char switches 'n';
  if subparam then Buffer.add_string switches "o2";
  if not debug then Buffer.add_char switches 'Q';
  (match verbose with
   | Some `V -> Buffer.add_string switches "V";
   | Some `VV -> Buffer.add_string switches "V";
   | Some `VVV -> Buffer.add_string switches "VVV";
   | None -> ());
  (* Call triangle and build the resulting objects *)
  let point, point_attribute, point_marker, triangle, triangle_attribute,
      neighbor, segment, segment_marker, edge, edge_marker,
      vor_point, vor_point_attribute, vor_edge, vor_normal =
    triangle (Buffer.contents switches) mesh triangle_area in
  let mesh_out : layout t =
    (make_mesh
      ~point:              point
      ~point_attribute:    point_attribute
      ~point_marker:       point_marker
      ~triangle:           triangle
      ~triangle_attribute: triangle_attribute
      ~neighbor:           neighbor
      ~segment:            segment
      ~segment_marker:     segment_marker
      ~edge:               edge
      ~edge_marker:        edge_marker
      ~hole: mesh#hole
      ~region: mesh#region)
  and vor : layout voronoi =
    (object
      method point               = vor_point
      method point_attribute     = vor_point_attribute
      method edge                = vor_edge
      method normal              = vor_normal
     end) in
  (mesh_out, vor)


(* Sub
 ***********************************************************************)

let sub (mesh: mesh) ?(pos=0) len =
  let m, n_tr, cols_tr = Mesh__MeshC.internal_sub (mesh :> Mesh__MeshC.mesh)
                                             ~pos len in
  let point_attribute =
    if Array2.dim2(mesh#point_attribute) = 0 || Array2.dim1(mesh#point_attribute) = 0 then
      mesh#point_attribute
    else
      Array2.sub_left mesh#point_attribute pos len in
  let triangle_attribute =
    let old_att = mesh#triangle_attribute in
    if Array2.dim2(old_att) = 0 || Array2.dim1(old_att) = 0 then old_att
    else (
      let att = Array2.create (float64) c_layout (n_tr) (Array2.dim2(old_att)) in
      Mesh__MeshC.iteri (fun i pi ->
                    for j = 0 to Array2.dim2(att) - 1 do
                      att.{i,j} <- old_att.{pi,j};
                    done
                   ) cols_tr;
      att
    ) in
  extend_mesh m
              ~point_attribute: point_attribute
              ~triangle_attribute: triangle_attribute


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

let permute_points_name = "Mesh_triangle.permute_points"

let do_permute_points (old_mesh: mesh) (perm: int_vec) inv_perm : mesh =
  let mesh = Mesh__MeshC.do_permute_points permute_points_name
                                      (old_mesh :> Mesh__MeshC.mesh)
                                      perm inv_perm in
  (* Permute the attributes *)
  let old_attr : mat = old_mesh#point_attribute in
  let attr = Array2.create (float64) c_layout (Array2.dim1(old_attr)) (Array2.dim2(old_attr)) in
  for i = 0 to Array2.dim1(old_attr) - 1 do
    let old_i = perm.{i} in
    for a = 0 to Array2.dim2(old_attr) - 1 do
      attr.{i,a} <- old_attr.{old_i,a}
    done
  done;
  extend_mesh mesh
              ~point_attribute: attr
              ~triangle_attribute: old_mesh#triangle_attribute


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


let permute_triangles_name = "Mesh_triangle.permute_triangles"

let do_permute_triangles (old_mesh: mesh) (perm: int_vec) : mesh =
  let mesh = Mesh__MeshC.do_permute_triangles permute_triangles_name
                                         (old_mesh :> Mesh__MeshC.mesh) perm in
  (* Permute attributes *)
  let old_attr : mat = old_mesh#triangle_attribute in
  let attr = Array2.create (float64) c_layout (Array2.dim1(old_attr)) (Array2.dim2(old_attr)) in
  for i = 0 to Array2.dim1(old_attr) - 1 do
    let old_i = perm.{i} in
    for a = 0 to Array2.dim2(old_attr) - 1 do
      attr.{i,a} <- old_attr.{old_i,a}
    done
  done;
  extend_mesh mesh
              ~point_attribute: (old_mesh#point_attribute)
              ~triangle_attribute: attr

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