/* integrate-demo.c
 *
 * A demo file to illustrate integration over triangulated domains.
 *
 * Author: Rouben Rostamian <rostamian@umbc.edu>
 * November 2005
 * Revised for Math 625, April 2006
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "read-mesh-data.h"
#include "twb-quad.h"
#include "integrate.h"
#include "fem.h"

double testfunc1(double x, double y)
{
	return 1.0;	/* answer should be the domain's area */
}

double testfunc2(double x, double y)
{
	return x*y*y;	/* exact integral over [0,1]x[0,1] is 1/6 */
}

double testfunc3(double x, double y)
{
	return sin(x*y);  /* exact integral over [0,1]x[0,1] is 0.2398117420 */
}

static void show_usage(const char *prog)
{
	fprintf(stderr, "Usage:  %s quad_degree\n", prog);
	fprintf(stderr, "        where quad_degree is the desired "
			"quarature degree\n");
}

int main(int argc, char **argv)
{
	Fem *fem;
	int quad_degree;
	char *endptr;

	if (argc != 2) {
		show_usage(argv[0]);
		return EXIT_FAILURE;
	}

	quad_degree = strtol(argv[1], &endptr, 10);
	if (*endptr != '\0') {
		show_usage(argv[0]);
		return EXIT_FAILURE;
	}

	fem = malloc(sizeof *fem);
	if (fem==NULL) {
		fprintf(stderr, "%s: out of memory!\n", argv[0]);
		return EXIT_FAILURE;
	}

	fem->qdat = twb_qdat(quad_degree);

	/*
	fem->mesh = read_mesh_data("Square/square.1");
	fem->mesh = read_mesh_data("Annulus/annulus.1");
	fem->mesh = read_mesh_data("Ell/ell.1");
	*/
	fem->mesh = read_mesh_data("Square/square.2");

	if (fem->mesh==NULL)
		return EXIT_FAILURE;	/* diagnostics printed
					   in read_mesh_data */
	fem->f = testfunc3;

	printf("integral = %g\n", integrate(fem));

	free_mesh(fem->mesh);
	return EXIT_SUCCESS;
}