# include <math.h>
# include <stdio.h>
# include <time.h>

# include "ppc_array.h"
# include "ppc_mesh.h"
# include "ppc_problem_spec.h"
# include "ppc_triangle.h"
# include "ppc_xmalloc.h"

/******************************************************************************/

static struct triangulateio *problem_spec_to_triangle ( struct problem_spec *spec )
{
	int i;
	struct triangulateio *in = xmalloc(sizeof *in);

	/* process points */
	in->numberofpoints = spec->npoints;
	make_vector(in->pointlist, 2 * in->numberofpoints);
	for (i = 0; i < in->numberofpoints; i++) {
		in->pointlist[2*i]   = spec->points[i].x;
		in->pointlist[2*i+1] = spec->points[i].y;
	}

	make_vector(in->pointmarkerlist, in->numberofpoints);
	for (i = 0; i < in->numberofpoints; i++)
		in->pointmarkerlist[i] = spec->points[i].bc;

	in->numberofpointattributes = 0;

	/* process segments */
	in->numberofsegments = spec->nsegments;
	make_vector(in->segmentlist, 2 * in->numberofsegments);
	for (i = 0; i < in->numberofsegments; i++) {
		in->segmentlist[2*i]   = spec->segments[i].point_no_1;
		in->segmentlist[2*i+1] = spec->segments[i].point_no_2;
	}

	make_vector(in->segmentmarkerlist, in->numberofsegments);
	for (i = 0; i < in->numberofsegments; i++)
		in->segmentmarkerlist[i] = spec->segments[i].bc;

	/* process holes */
	in->numberofholes = spec->nholes;
	if (in->numberofholes != 0) {
		make_vector(in->holelist, 2 * in->numberofholes);
		for (i = 0; i < in->numberofholes; i++) {
			in->holelist[2*i]   = spec->holes[i].x;
			in->holelist[2*i+1] = spec->holes[i].y;
		}
	} else
		in->holelist = NULL;

	in->numberofregions = 0;

	return in;
}
/******************************************************************************/

static struct triangulateio *do_triangulate ( struct triangulateio *in, double a )
{
	char *opts_string = "Qzpeq30a%f";
	char opts[64];
	struct triangulateio *out = xmalloc(sizeof *out);

	out->pointlist = NULL;
	out->pointmarkerlist = NULL;
	out->edgelist = NULL;
        out->edgemarkerlist = NULL;
	out->trianglelist = NULL;
	out->segmentlist = NULL;
        out->segmentmarkerlist = NULL;
	sprintf(opts, opts_string, a);
	triangulate ( opts, in, out, NULL );

	return out;
}

/******************************************************************************/

static void assign_elem_edges(
		struct elem *elems, int nelems,
		struct edge *edges, int nedges)
{
	for (int r = 0; r < nelems; r++) {
		for (int i = 0; i < 3; i++) {	/* i: vertex index */
			int j = (i+1)%3;
			int k = (i+2)%3;
			int n1 = elems[r].n[j]->nodeno;
			int n2 = elems[r].n[k]->nodeno;
			for (int s = 0; s < nedges; s++) {
				int m1 = edges[s].n[0]->nodeno;
				int m2 = edges[s].n[1]->nodeno;
				if ((m1 == n1 && m2 == n2)
						|| (m1 == n2 && m2 == n1)) {
					elems[r].e[i] = &edges[s];
					break;
				}
			}
		}
	}
}
/******************************************************************************/

static void set_element_edge_vectors(struct elem *ep)
{
	for (int i = 0; i < 3; i++) {	/* i: vertex index */
		int j = (i+1)%3;
		int k = (i+2)%3;
		ep->ex[i] = ep->n[k]->x - ep->n[j]->x;
		ep->ey[i] = ep->n[k]->y - ep->n[j]->y;
	}
}
/******************************************************************************/

static void set_element_area(struct elem *ep)
{
	ep->area = (ep->ex[0]*ep->ey[1] - ep->ex[1]*ep->ey[0])/2.0;
}
/******************************************************************************/

static void set_edge_vectors_and_areas(struct elem *elems, int nelems)
{
	for (int i = 0; i < nelems; i++) {
		struct elem *ep = &elems[i];
		set_element_edge_vectors(ep); 
		set_element_area(ep); 
	}
}
/******************************************************************************/

static struct mesh *triangle_to_mesh(struct triangulateio *out)
{
	struct node *nodes;
	struct edge *edges;
	struct elem *elems;
	int i, nnodes, nedges, nelems;
	struct mesh *mesh = xmalloc(sizeof *mesh);

	nnodes = out->numberofpoints;
	make_vector(nodes, nnodes);
	for (i = 0; i < out->numberofpoints; i++) {
		nodes[i].nodeno = i;
		nodes[i].x = out->pointlist[2*i];
		nodes[i].y = out->pointlist[2*i+1];
		nodes[i].z = 0.0;
		nodes[i].bc = out->pointmarkerlist[i];
	}

	nedges = out->numberofedges;
	make_vector ( edges, nedges );
	for (i = 0; i < nedges; i++) {
		edges[i].edgeno = i;
		edges[i].n[0] = &nodes[out->edgelist[2*i]];
		edges[i].n[1] = &nodes[out->edgelist[2*i+1]];
		edges[i].bc = out->edgemarkerlist[i];
	}

	nelems = out->numberoftriangles;
	make_vector(elems, nelems);
	for (i = 0; i < nelems; i++) {
		elems[i].elemno = i;
		elems[i].n[0] = &nodes[out->trianglelist[3*i]];
		elems[i].n[1] = &nodes[out->trianglelist[3*i+1]];
		elems[i].n[2] = &nodes[out->trianglelist[3*i+2]];
	}

	assign_elem_edges(elems, nelems, edges, nedges);
	set_edge_vectors_and_areas(elems, nelems);

	mesh->nnodes = nnodes;
	mesh->nedges = nedges;
	mesh->nelems = nelems;
	mesh->nodes = nodes;
	mesh->edges = edges;
	mesh->elems = elems;

	return mesh;
}
/******************************************************************************/

static void free_triangle_in_structure(struct triangulateio *in)
{
	free_vector(in->pointlist);
	free_vector(in->pointmarkerlist);
	free_vector(in->segmentlist);
	free_vector(in->segmentmarkerlist);
	free_vector(in->holelist);
	free(in);
}
/******************************************************************************/

static void free_triangle_out_structure(struct triangulateio *out)
{

	free(out->pointlist);
	free(out->pointmarkerlist);
	free(out->edgelist);
	free(out->edgemarkerlist);
	free(out->trianglelist);
	free(out->segmentlist);
	free(out->segmentmarkerlist);
	free(out);
}
/******************************************************************************/

void free_mesh ( struct mesh *mesh )
{
  if ( mesh == NULL )
  {
    return;
  }

  free_vector ( mesh->nodes );
  free_vector ( mesh->edges );
  free_vector ( mesh->elems );
  free ( mesh );

  return;
}
/******************************************************************************/

struct mesh *make_mesh ( struct problem_spec *spec, double a )
{
  struct triangulateio *in, *out;
  struct mesh *mesh;

  in = problem_spec_to_triangle(spec);
  out = do_triangulate ( in, a );
  mesh = triangle_to_mesh(out);
  free_triangle_in_structure(in);
  free_triangle_out_structure(out);
  return mesh;
}
/******************************************************************************/

void mesh_to_eps ( struct mesh *mesh, char *outfile )
/*
The "physical" coordinates range over xmin..xmax, and ymin..ymax.
The postscript bounding box is (0,0) to (W,H), where max(W,H) = D is
define below.

The physical domain is scaled to the postscript domain by a factor of s.
Thus, the width of the image area is s*w.  We add margins of p*s*w
from the left and right, so that 2*p*s*w + s*w = W.  Similarly,
2*p*s*h + s*h = H.  Therefore: W = (1+2*p)*s*w, and H = (1+2*p)*s*h, and
max(W,H) = (1+2*p)*s*max(w,h).  Consequently, D = (1+2*p)*s*d, where
d=max(w,h).  From here we get:
   s = D/((1+2*p)*d).
In particular, (1+2*p)s = D/d, therefore W = D/d*w and H = D/d*h.

The linear mapping from the physical domain (x,y) to the postscript
domain (X,Y) takes xmin to the images left margin, that is p*s*w,
therefore X = p*s*w + s*(x - xmin).  Similarly Y = p*s*h + s*(y - ymin).

These are all incorporated into the function mesh_to_eps().

2012-12-29
*/
{
#define D 400	/* the larger of the W & H of the bounding box */

	FILE *fp;
	struct node *nodes = mesh->nodes;
	struct edge *edges = mesh->edges;
	struct elem *elems = mesh->elems;
	time_t now;
	double xmin, xmax, ymin, ymax, w, h, W, H, d, s;
	double p = 0.01;
	int i, k;

	if ((fp = fopen(outfile, "w")) == NULL) {
		fprintf(stderr, "cannot open file %s for writing\n", outfile);
		return;
	}

	xmin = xmax = nodes[0].x;
	ymin = ymax = nodes[0].y;
	for (i = 1; i < mesh->nnodes; i++) {
		if (nodes[i].x < xmin)
			xmin = nodes[i].x;
		else if(nodes[i].x > xmax)
			xmax = nodes[i].x;
		if (nodes[i].y < ymin)
			ymin = nodes[i].y;
		else if(nodes[i].y > ymax)
			ymax = nodes[i].y;
	}
	w = xmax - xmin;
	h = ymax - ymin;
	d = fmax ( w, h );
	W = D/d*w;		/* bounding box width */
	H = D/d*h;		/* bounding box height */
	s = D/((1+2*p)*d);	/* scale */

	fputs("%!PS-Adobe-3.0 EPSF-3.0\n", fp);
	fprintf(fp, "%%%%BoundingBox: 0 0 %g %g\n", W, H);
	fprintf(fp, "%%%%Title: (%s)\n", outfile);
	fprintf(fp, "%%%%Creator: C Projects, %s\n", __FILE__);
	now = time(NULL);
	fprintf(fp, "%%%%CreationDate: %s", ctime(&now));
	fputs("%%EndComments\n", fp);

	fputs("gsave\n", fp);
	fputs("1 1 0 setrgbcolor\n", fp);
	for (i = 0; i < mesh->nelems; i++) {
		double x[3], y[3];
		for (k = 0; k < 3; k++) {
			x[k] = p*s*w + s*(elems[i].n[k]->x - xmin);
			y[k] = p*s*h + s*(elems[i].n[k]->y - ymin);
		}
		fprintf(fp, "%g %g moveto %g %g lineto %g %g lineto "
				"closepath fill\n",
				x[0], y[0],
				x[1], y[1],
				x[2], y[2]);
	}

	fputs("0 setgray\n", fp);
	for (i = 0; i < mesh->nedges; i++) {
		double x1 = edges[i].n[0]->x;
		double y1 = edges[i].n[0]->y;
		double x2 = edges[i].n[1]->x;
		double y2 = edges[i].n[1]->y;
		double X1 = p*s*w + s*(x1 - xmin);
		double Y1 = p*s*h + s*(y1 - ymin);
		double X2 = p*s*w + s*(x2 - xmin);
		double Y2 = p*s*h + s*(y2 - ymin);
		fprintf(fp, "%g %g moveto %g %g lineto stroke\n",
			X1, Y1, X2, Y2);
	}
	fputs("grestore\n", fp);
	fputs("showpage\n", fp);
	fputs("%%EOF\n", fp);
	fclose(fp);
	fprintf(stderr, "postscript output written to file %s\n", outfile);

  return;
}