#!/usr/bin/perl

## Copyright 2008 by James Rath and the University of Texas at Austin
## This code is available under a Creative Commons Attribution-ShareAlike license.  You
## are free to copy and share this as you see fit; however, you must attribute any copy
## to me.  You are also free to create derivative works; however, you must acknowledge
## the use of this original, and you must share your work under the same terms you
## received this copy.  For details, see http://creativecommons.org/licenses/by-sa/3.0/
## For permission for any other use, please contact me at ratjamm@alum.mit.edu


## one of these days i'll write usage instructions...
@usage = ("blah\n",
		  "blah\n",
		  "blah\n");


## name on the command line tells us what file(s) to look in
if ($#ARGV+1) {
	($ARGV[0] eq "--help")  &&  die @usage;
	$title = shift @ARGV;
}


## slurp up the nodes
open(NODES,"< $title.node");
# First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>
# skip possible comment lines
$headerline = <NODES>;
while ($headerline =~ /^\s*#/ || $headerline =~ /^\s*$/) {
	$headerline = <NODES>;
	die "all comments in nodes file?!" if (!$headerline);
}
# read and check the information
($nnodes,$dim,$nattr,$bdymarkp) = ($headerline =~ /\s*(\d+)\s*(\d+)\s*(\d+)\s*(\d+)\s*/);
die "naughty nodes file\n" if (!$nnodes || $dim!=2 || $nattr!=0 || $bdymarkp!=1);
# preallocate space for data
$#indices = $nnodes - 1;
$#xs = $nnodes - 1;
$#ys = $nnodes - 1;
$#bdymarks = $nnodes - 1;
# Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]
$nLines = 0;
while (<NODES>) {
	next if ((/^\s*#/) || (/^\s*$/));
	($indices[$nLines],$xs[$nLines],$ys[$nLines],$bdymarks[$nLines]) = split(' ', $_);
	$nLines++;
}
die "didn't read all the nodes?!\n" if ($nLines!=$nnodes);
close(NODES);

# sort the data just in case the nodes weren't written in order
@revIndices = sort { $indicies[$a] <=> $indicies[$b] } 0..$#indices;
@xs = @xs[@revIndices];
@ys = @ys[@revIndices];
@bdymarks = @bdymarks[@revIndices];

## open a file for the mass matrix data
open(MM,"> $title.mm");


## slurp up the elements and write 'em out
open(ELTS,"< $title.ele");
# for checking, read the number of elements from the first line
# First line: <# of triangles> <nodes per triangle> <# of attributes>
$headerline = <ELTS>;
while ($headerline =~ /^\s*#/ || $headerline =~ /^\s*$/) {
	$headerline = <ELTS>;
	die "all comments in elements file?!" if (!$headerline);
}
($ntri,$nodepertri,$nattr) = ($headerline =~ /\s*(\d+)\s*(\d+)\s*(\d+)\s*/);
die "naughty elements file\n" if (!$ntri || $nodepertri!=3 || $nattr!=0);
# read the elements, processing one-by-one
$nLines = 0;
while (<ELTS>) {
	next if ((/^\s*#/) || (/^\s*$/));
# Remaining lines: <triangle #> <node> <node> <node> ... [attributes]
	($tri,$v1,$v2,$v3) = split(' ', $_);
#	darn zero indexing in perl and one indexing in matlab
	$l1sq = ($xs[$v2-1]-$xs[$v3-1])*($xs[$v2-1]-$xs[$v3-1]) + ($ys[$v2-1]-$ys[$v3-1])*($ys[$v2-1]-$ys[$v3-1]);
	$l2sq = ($xs[$v3-1]-$xs[$v1-1])*($xs[$v3-1]-$xs[$v1-1]) + ($ys[$v3-1]-$ys[$v1-1])*($ys[$v3-1]-$ys[$v1-1]);
	$l3sq = ($xs[$v1-1]-$xs[$v2-1])*($xs[$v1-1]-$xs[$v2-1]) + ($ys[$v1-1]-$ys[$v2-1])*($ys[$v1-1]-$ys[$v2-1]);
	$area = &heron(sqrt($l1sq),sqrt($l2sq),sqrt($l3sq));
	printf MM "%d %d %.5g\n", $v1,$v1, 0.25*$l1sq/$area  if  ($bdymarks[$v1-1]==0);
	printf MM "%d %d %.5g\n", $v2,$v2, 0.25*$l2sq/$area  if  ($bdymarks[$v2-1]==0);
	printf MM "%d %d %.5g\n", $v3,$v3, 0.25*$l3sq/$area  if  ($bdymarks[$v3-1]==0);
	printf MM "%d %d %.5g\n", $v1,$v2, 0.125*($l3sq-$l1sq-$l2sq)/$area  if  ($bdymarks[$v1-1]==0 && $bdymarks[$v2-1]==0);
	printf MM "%d %d %.5g\n", $v2,$v1, 0.125*($l3sq-$l1sq-$l2sq)/$area  if  ($bdymarks[$v2-1]==0 && $bdymarks[$v1-1]==0);
	printf MM "%d %d %.5g\n", $v2,$v3, 0.125*($l1sq-$l2sq-$l3sq)/$area  if  ($bdymarks[$v2-1]==0 && $bdymarks[$v3-1]==0);
	printf MM "%d %d %.5g\n", $v3,$v2, 0.125*($l1sq-$l2sq-$l3sq)/$area  if  ($bdymarks[$v3-1]==0 && $bdymarks[$v2-1]==0);
	printf MM "%d %d %.5g\n", $v3,$v1, 0.125*($l2sq-$l1sq-$l3sq)/$area  if  ($bdymarks[$v3-1]==0 && $bdymarks[$v1-1]==0);
	printf MM "%d %d %.5g\n", $v1,$v3, 0.125*($l2sq-$l1sq-$l3sq)/$area  if  ($bdymarks[$v1-1]==0 && $bdymarks[$v3-1]==0);
	$nLines++;
}
# check that we got 'em all
die "didn't read all the triangles?!\n" if ($nLines!=$ntri);
close(ELTS);

## close up shop
close(MM);



### helper functions

sub heron {
	# sort the lengths in descending order
	($a,$b,$c) = sort {$b <=> $a} @_;
	# NB: parens and order are important here for stability on thin triangles
	0.25 * sqrt( ($a+($b+$c)) * ($c-($a-$b)) * ($c+($a-$b)) * ($a+($b-$c)) );
}