use strict;
use lib '.';
use generate_cmgui_mesh qw(generate_cmgui_mesh);
use Math::Trig qw(pi);

my $input_node_filename = "refine.exnode";
my $input_element_filename = "refine.exelem";

my $output_node_filename = "refine_volume.exnode";
my $output_element_filename = "refine_volume.exelem";

# These will be calculated from the mesh read in
my $xi1_discretization;
my $xi2_discretization;
my $xi3_discretization;

my $pi = pi;

#If $generate_cubics is non zero then the resulting mesh will be cubic Hermite, otherwise linear.
my $generate_cubics = 1;

my $group_name = "medrect";

my $start_node = 193;

# The node file is assumed to have 
# 1) coordinates, coordinate, rectangular cartesian, #Components=3
#   x.  Value index= 1, #Derivatives= 3 (d/ds1,d/ds2,d2/ds1ds2)
#   y.  Value index= 5, #Derivatives= 3 (d/ds1,d/ds2,d2/ds1ds2)
#   z.  Value index= 9, #Derivatives= 3 (d/ds1,d/ds2,d2/ds1ds2)
my @nodes = ();
my $line;

open INPUT, "<$input_node_filename";

while (defined ($line = <INPUT>))
  {
	 if ($line =~ m/Node:\s*(\d+)/)
		{
		  my @node_array = $1;

		  #Read the next three lines assuming these have four values each.
		  for (my $i = 0 ; $i < 3 ; $i++)
			 {
				$line = <INPUT>;
				if ($line =~ m/([\-\dE\+\.]+)\s+([\-\dE\+\.]+)\s+([\-\dE\+\.]+)\s+([\-\dE\+\.]+)/)
				  {
					 push @node_array, $1, $2, $3, $4;
				  }
				else
				  {
					 die "Unable to read coordinates and derivatives for node $node_array[0]\n";
				  }
			 }
		  push @nodes, \@node_array;
		}
  }

close INPUT;

print "Read " . scalar @nodes . " nodes.\n";


my @elements = ();

open INPUT, "<$input_element_filename";

while (defined ($line = <INPUT>))
  {
	 # Ignore everything but the node indicies.
	 if ($line =~ m/Nodes:/)
		{
		  my @node_array;

		  $line = <INPUT>;
		  if ($line =~ m/(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/)
			 {
				push @node_array, $1, $2, $3, $4;
			 }
		  else
			 {
				die "Unable to read node indices : \n$_\n";
			 }
		  push @elements, \@node_array;
		}
  }

close INPUT;

print "Read " . scalar @elements . " elements.\n";

my @node_spread;

my $search_node_xi1 = 0;
my $search_node_xi2 = $start_node;
my $xi1_dimension = 0;
my $xi2_dimension = 0;

my $element;

$node_spread[$xi1_dimension][$xi2_dimension] = $start_node;

while ($search_node_xi2)
  {
	 #Find the element with this node in the first position
	 undef $element;
	 my $i = 0;
	 while ($elements[$i] && (${$elements[$i]}[0] != $search_node_xi2))
		{
		  $i++;
		}
	 $search_node_xi2 = 0;
	 if ($element = $elements[$i])
		{
		  #Remove this element from the list
		  splice @elements, $i, 1;

		  $node_spread[$xi1_dimension + 1][$xi2_dimension] = 
			 $$element[1];
		  $search_node_xi1 = $$element[1];
		  $node_spread[$xi1_dimension][$xi2_dimension + 1] = 
			 $$element[2];
		  $search_node_xi2 = $$element[2];

		  # Do the xi1 search
		  while ($search_node_xi1)
			 {
				undef $element;
				my $i = 0;
				while ($elements[$i] && (${$elements[$i]}[0] != $search_node_xi1))
				  {
					 $i++;
				  }
				$search_node_xi1 = 0;
				if ($element = $elements[$i])
				  {
					 $xi1_dimension++;
					 
					 #Remove this element from the list
					 splice @elements, $i, 1;
					 $node_spread[$xi1_dimension + 1][$xi2_dimension] = 
						$$element[1];
					 $search_node_xi1 = $$element[1];
				  }
			 }
		  $xi1_dimension = 0;
		  $xi2_dimension++;
		}
  }

if (@elements)
  {
	 die scalar @elements . " elements remaining after node spread.";
  }

$xi1_discretization = scalar @node_spread - 1;
$xi3_discretization = 1;
if ($xi3_discretization > 1)
  {
	 die "The location function code will need to be expanded to interpolate nodes and more.";
  }
$xi2_discretization = (scalar @{$node_spread[0]} + 1) / 2 - $xi3_discretization;

print ("xi1 $xi1_discretization\nxi2 $xi2_discretization\nxi3 $xi3_discretization\n");

sub location_function
{
  my ($xi1_step, $xi2_step, $xi3_step) = @_;

  my @coordinates;
  my @derivatives;

  my $phi = $pi * $xi3_step / $xi3_discretization - $pi / 2.0;

  my $xi2 = $xi2_step / $xi2_discretization;
  my $psi = $pi * $xi2;

  my $xi1 = $xi1_step / $xi1_discretization;
  my $theta = $pi * $xi1;

  my $xi1_spread = $xi1_step;
  my $xi2_spread;
  my $xi2_spread_opposite;
  my $xi2_xi30_spread = $xi2_step;
  my $xi2_xi31_spread = (2 * $xi2_discretization - $xi2_step) % (2 * $xi2_discretization);
  if ($xi3_step == 0)
	 {
		$xi2_spread = $xi2_xi30_spread;
		$xi2_spread_opposite = $xi2_xi31_spread;
	 }
  elsif ($xi3_step == 1)
	 {
		$xi2_spread = $xi2_xi31_spread;
		$xi2_spread_opposite = $xi2_xi30_spread;
	 }
  else
	 {
		die "xi3_step is $xi3_step";
	 }
  my $node_number = $node_spread[$xi1_spread][$xi2_spread];
  #print ("  $xi1_spread $xi2_spread $node_number\n");
  if (!$node_number)
	 {
		die "No node was stored for position $xi1_spread $xi2_spread";
	 }

  my $node;
  my $i = 0;
  while ($nodes[$i] && (${$nodes[$i]}[0] != $node_number))
  {
	 $i++;
  }

  if (!($node = $nodes[$i]))
  {
	 die "Unable to find a node entry for node $node_number";
  }

  #We don't need to calculate this all the time but it is easier to start
  #this way and will be helpful if we end up adding inbetween layers in xi3.
  my $node_number_opposite = $node_spread[$xi1_spread][$xi2_spread_opposite];
  #print ("  $xi1_spread $xi2_spread_opposite $node_number_opposite\n");
  if (!$node_number_opposite)
	 {
		die "No node was stored for opposite position $xi1_spread $xi2_spread_opposite";
	 }

  my $node_opposite;
  $i = 0;
  while ($nodes[$i] && (${$nodes[$i]}[0] != $node_number_opposite))
  {
	 $i++;
  }

  if (!($node_opposite = $nodes[$i]))
  {
	 die "Unable to find a node entry for opposite node $node_number_opposite";
  }

  for (my $i = 0 ; $i < 3 ; $i++)
	 {
		$coordinates[$i] = ${$node}[$i * 4 + 1];

		# d/ds1
		$derivatives[$i * 7] = ${$node}[$i * 4 + 2];

		if ($xi3_step == 0)
		  {
			 # d/ds2
			 $derivatives[$i * 7 + 1] =
				${$node}[$i * 4 + 3];

			 # d2/ds1ds2
			 $derivatives[$i * 7 + 2] =
				${$node}[$i * 4 + 4];
		  }
		else
		  {
			 # d/ds2
			 $derivatives[$i * 7 + 1] =
				-${$node}[$i * 4 + 3];

			 # d2/ds1ds2
			 $derivatives[$i * 7 + 2] =
				-${$node}[$i * 4 + 4];				
		  }

		if ($xi2_step == 0)
		  {
			 # d/ds3
			 $derivatives[$i * 7 + 3] =
				0.0;

			 # d2/ds1ds3
			 $derivatives[$i * 7 + 4] =
				0.0;
		  }
		elsif ($xi2_step == $xi2_discretization)
		  {
			 # d/ds3
			 $derivatives[$i * 7 + 3] =
				0.0;

			 # d2/ds1ds3
			 $derivatives[$i * 7 + 4] =
				0.0;
		  }
		else
		  {
			 #Want a linear difference across the element
			 if ($xi3_step == 0)
				{
				  # d/ds3
				  $derivatives[$i * 7 + 3] =
					 - ${$node}[$i * 4 + 1]
						+ ${$node_opposite}[$i * 4 + 1];
				}
			 else
				{
				  # d/ds3
				  $derivatives[$i * 7 + 3] =
					 ${$node}[$i * 4 + 1]
						- ${$node_opposite}[$i * 4 + 1];
				}


			 # d2/ds1ds3
			 $derivatives[$i * 7 + 4] =
				0.0;
		  }


		# d2/ds2ds3
		$derivatives[$i * 7 + 5] =
		  0.0;

		# d3/ds1ds2ds3
		$derivatives[$i * 7 + 6] =
		  0.0;

	 }

  return(\@coordinates, \@derivatives);
}


generate_cmgui_mesh( $output_node_filename, $output_element_filename,
     $xi1_discretization, $xi2_discretization, $xi3_discretization,
     \&location_function, $generate_cubics, $group_name );