Source file mesh_easymeshF.ml

# 1 "easymesh/mesh_easymeshFC.ml"
(* FORTRAN/C functions *)

open Printf
open Scanf
open Bigarray

type pslg = fortran_layout Mesh.pslg
type mesh = fortran_layout Mesh.t

(* Write the [pslg] to the channel [fh]. *)
let output_pslg fh (pslg: pslg) area =
  let pt = pslg#point
  and seg = pslg#segment
  and pt_marker = pslg#point_marker
  and seg_marker = pslg#segment_marker in
  let pt_marker =
    if Array1.dim pt_marker > 0 then (fun i -> pt_marker.{i})
    else (fun _i -> 1) in
  let seg_marker =
    if Array1.dim seg_marker > 0 then (fun i -> seg_marker.{i})
    else (fun _i -> 1) in
  (* Save points coordinates *)
  fprintf fh "%i\n" (Array2.dim2(pt)); (* number of nodes *)
  for i = 1 to Array2.dim2(pt) do
    (* EasyMesh expects indexes in 0 .. nnodes-1 *)
    fprintf fh "%i: %.13g %.13g %.13g %i\n"
      (i - 1) (pt.{1,i}) (pt.{2,i})  area (pt_marker i)
  done;
  (* Save segments *)
  fprintf fh "%i\n" (Array2.dim2(seg)); (* number of segments *)
  for i = 1 to Array2.dim2(seg) do
    fprintf fh "%i: %i %i %i\n"
      (i - 1)  (seg.{1,i} - 1) (seg.{2,i} - 1)
      (seg_marker i)
  done

(* FIXME: comments are possible in the files *)
(* FIXME: if a file does not exists, return empty array? *)
(* [read_fortran fname] reads the collection of filenames [fname].n,
   [fname].e and [fname].s and creates a mesh struture with fortran
   layout.  This function can throw a variety of exceptions depending
   of what goes wrong.  *)
let read (pslg: pslg) fname : mesh =
  (* Read nodes *)
  let fh = open_in (fname ^ ".n") in
  let nnodes = int_of_string(input_line fh) in
  let pt = Array2.create (float64) fortran_layout (2) (nnodes) in
  let pt_marker = Array1.create (int) fortran_layout (nnodes) in
  let sb = Scanning.from_channel fh in
  for _i = 1 to Array2.dim2(pt) do
    bscanf sb " %i %g %g %i" (fun i x y m ->
                                 let i = i + 1 in
                                 pt.{1,i} <- x;
                                 pt.{2,i} <- y;
                                 pt_marker.{i} <- m)
  done;
  close_in fh;
  (* Read triangles *)
  let fh = open_in (fname ^ ".e") in
  let n = int_of_string(input_line fh) in
  let tr = Array2.create (int) fortran_layout (3) (n)
  and tr_nbh = Array2.create (int) fortran_layout (3) (n) in
  let sb = Scanning.from_channel fh in
  for _i = 1 to Array2.dim2(tr) do
    bscanf sb " %i %i %i %i %i %i %i %_i %_i %_i %_f %_f %_i"
      (fun e i j k ei ej ek ->
         let e = e + 1 in
         tr.{1,e} <- i + 1;
         tr.{2,e} <- j + 1;
         tr.{3,e} <- k + 1;
         tr_nbh.{1,e} <- ei + 1;
         tr_nbh.{2,e} <- ej + 1;
         tr_nbh.{3,e} <- ek + 1;
      )
  done;
  close_in fh;
  let edge, edge_marker =
    try
      (* Read edges, if file exists *)
      let fh = open_in (fname ^ ".s") in
      let n = int_of_string(input_line fh) in
      let edge = Array2.create (int) fortran_layout (2) (n)
      and edge_marker = Array1.create (int) fortran_layout (n) in
      let sb = Scanning.from_channel fh in
      for _i = 1 to Array2.dim2(edge) do
        bscanf sb " %i %i %i %_i %_i %i" (fun s c d m ->
            let s = s + 1 in
            edge.{1,s} <- c + 1;
            edge.{2,s} <- d + 1;
            edge_marker.{s} <- m;
          )
      done;
      close_in fh;
      (edge, edge_marker)
    with Sys_error _ ->
      ( (Array2.create (int) fortran_layout (2) (0)), (Array1.create (int) fortran_layout (0)) ) in
  let segment = pslg#segment in
  let segment_marker = pslg#segment_marker in
  let hole = pslg#hole in
  let region = pslg#region in
  (object
     method point = pt
     method point_marker = pt_marker
     method triangle = tr
     method neighbor = tr_nbh
     method edge = edge
     method edge_marker = edge_marker
     (* From PSLG *)
     method segment = segment
     method segment_marker = segment_marker
     method hole = hole
     method region = region
   end : fortran_layout Mesh.t)

let empty_pslg : pslg =
  let empty_mat = Array2.create (float64) fortran_layout (2) (0)
  and empty_int_mat = Array2.create (int) fortran_layout (2) (0)
  and empty_vec = Array1.create (int) fortran_layout (0) in
  (object
     method point = empty_mat
     method point_marker = empty_vec
     method segment = empty_int_mat
     method segment_marker = empty_vec
     method hole = empty_mat
     method region = empty_mat
   end)

(* [sort3 a b c] sort [a], [b], [c] from the larger to the smaller *)
let sort3 (a:int) b c =
  if a > b then
    if a <= c then (c, a, b)
    else if b > c then (a, b, c) else (a, c, b)
  else (* a <= b *)
    if b <= c then (c, b, a)
    else if a > c then (b, a, c) else (b, c, a)

let write (mesh: mesh) file =
  (* Compute the information needed by easymesh *)
  let pt = mesh#point in
  let n = Array2.dim2 pt in
  let tr = mesh#triangle in
  let ed = mesh#edge in
  (* For two nodes (i,j), i > j, return the (at most) two triangles
     sharing that edge. *)
  let pt_tr = Hashtbl.create (Array2.dim2(ed)) in
  for t = 1 to Array2.dim2(tr) do
    let i1, i2, i3 = (* i1 > i2 > i3 *)
      sort3 (tr.{1,t}) (tr.{2,t}) (tr.{3,t}) in
    if i1 = i2 || i1 = i3 || i2 = i3 then
      invalid_arg "Easymesh.write: illegal mesh (triangle corner indices)";
    Hashtbl.add pt_tr (i1,i2) t;
    Hashtbl.add pt_tr (i1,i3) t;
    Hashtbl.add pt_tr (i2,i3) t;
  done;
  (* For two nodes (i,j), i > j, give the index of the edge, if exists. *)
  let pt_ed = Hashtbl.create (Array2.dim2(ed)) in
  for e = 1 to Array2.dim2(ed) do
    let i = ed.{1,e} and j = ed.{2,e} in
    if i = j then invalid_arg "Easymesh.write: illegal mesh (egde enpoints)"
    else if i > j then Hashtbl.add pt_ed (i,j) e
    else Hashtbl.add pt_ed (j,i) e
  done;
  (* Write [file].n *)
  let pt_marker = mesh#point_marker in
  let marker =
    if Array1.dim pt_marker = 0 then (fun _ -> 1) else (fun i -> pt_marker.{i}) in
  let fh = open_out(file ^ ".n") in
  fprintf fh "%i\n" n;
  for i = 1 to Array2.dim2(pt) do
    fprintf fh "%i: %.17g %.17g %i\n" (i - 1)
      (pt.{1,i}) (pt.{2,i}) (marker i)
  done;
  close_out fh;
  (* Write [file].e *)
  let fh = open_out(file ^ ".e") in
  fprintf fh "%i\n" (Array2.dim2(tr));
  let other_triangle i12 t = (* assume i12 = (i1, i2) with i1 > i2 *)
    match Hashtbl.find_all pt_tr i12 with
    | [ _ ] -> -1 (* no other triangle *)
    | [t1; t2] -> if t1 = t then t2 - 1 else t1 - 1
    | _ -> assert false in
  for t = 1 to Array2.dim2(tr) do
    let i1, i2, i3 = (* i1 > i2 > i3 *)
      sort3 (tr.{1,t}) (tr.{2,t}) (tr.{3,t}) in
    let x1 = pt.{1,i1} and y1 = pt.{2,i1}
    and x2 = pt.{1,i2} and y2 = pt.{2,i2}
    and x3 = pt.{1,i3} and y3 = pt.{2,i3} in
    (* Compute the circumcenter.
       See http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html *)
    let x21 = x2 -. x1 and y21 = y2 -. y1
    and x31 = x3 -. x1 and y31 = y3 -. y1 in
    (* FIXME: compute the determinant more accurately. *)
    let det = 2. *. (x21 *. y31 -. x31 *. y21) in
    let d21 = x21 *. x21 +. y21 *. y21
    and d31 = x31 *. x31 +. y31 *. y31 in
    let mx = x1 -. (y21 *. d31 -. y31 *. d21) /. det
    and my = y1 +. (x21 *. d31 -. x31 *. d21) /. det in
    let i1' = (i2,i3) and i2' = (i1,i3) and i3' = (i1,i2) in
    fprintf fh "%i: %i %i %i %i %i %i %i %i %i %g %g 0\n" (t - 1)
      (i1 - 1) (i2 - 1) (i3 - 1)
      (other_triangle i1' t) (* opposite i1 *)
      (other_triangle i2' t)
      (other_triangle i3' t)
      (Hashtbl.find pt_ed i1' - 1)
      (Hashtbl.find pt_ed i2' - 1)
      (Hashtbl.find pt_ed i3' - 1)
      mx my
  done;
  close_out fh;
  (* Write [file].s *)
  let fh = open_out(file ^ ".s") in
  fprintf fh "%i\n" (Array2.dim2(ed));
  let ed_marker = mesh#edge_marker in
  let marker =
    if Array1.dim ed_marker = 0 then (fun _ -> 1) else (fun e -> ed_marker.{e}) in
  for e = 1 to Array2.dim2(ed) do
    let i1 = ed.{1,e} and i2 = ed.{2,e} in
    let i12 = if i1 > i2 then (i1, i2) else (i2, i1) in
    let t1, t2 = match Hashtbl.find_all pt_tr i12 with
      | [t1] -> t1 - 1, -1
      | [t1; t2] -> t1 - 1, t2 - 1
      | _ -> assert false in
    fprintf fh "%i: %i %i %i %i %i\n" (e - 1) (i1 - 1) (i2 - 1)
      t1 t2 (marker e)
  done;
  close_out fh