fem.cpp

//Echo 2015-08-04
//FEM Prototype

#include <iostream>
#include <cmath>
#include <gsl/gsl_linalg.h>
#include <vector>

//number of triangles in lattice
#define Nt 17

//constant in poisson equation
#define r 4

//access to coordinates in points of a triangle
#define xi (e.p1.x)
#define yi (e.p1.y)

#define xj (e.p2.x)
#define yj (e.p2.y)

#define xk (e.p3.x)
#define yk (e.p3.y)

//values for L functions (below)
#define ai (xj*yk-xk*yj)
#define bi (yj-yk)
#define ci (xk-xj)

#define aj (xk*yi-xi*yk)
#define bj (yk-yi)
#define cj (xi-xk)

#define ak (xi*yj-xj*yi)
#define bk (yi-yj)
#define ck (xj-xi)

//general area of a triangle
#define a abs((xj*yk+xi*yj+xk*yi-(xj*yi+xi*yk+xk*yj))/2)

//area for equilateral triangles
#define ae sqrt(3)*s*s/2

//function of value 1 at p1, p2, p3 resp.
#define Li (ai + bi*x + ci*y)/(2a)
#define Lj (aj + bj*x + cj*y)/(2a)
#define Lk (ak + bk*x + ck*y)/(2a)

using namespace std;

int pointIndex = 0;

struct point{
	double x;
	double y;
	int index;

	point(){}
	point(double x, double y): x(x), y(y){
		index = pointIndex++;
	};
};

struct triangle{
	point p1;
	point p2;
	point p3;

	triangle(){}
	triangle(point p1, point p2, point p3): p1(p1), p2(p2), p3(p3){};

};

//array of triangles
triangle triangleList[Nt];

//finds the final number of vertices given the 
//final number of triangles to make the final array
int findNumVertices(){
	int i = 0;
	int num = 4;
	while(num<Nt){
		num*=4;
		i++;
	}
	return 2*i*(i+1)+1;
}

//current number of triangles
static int c = 0;

//loops through all triangles in the triangle array
//makes 3 new triangles using the midpoints between the vertices
//puts one in the location of the old triangle
//the others at the end
void explode(){
	for(int i=0; i<c; i++){
		triangle e = triangleList[i];
		
		point m1((xi+xj)/2, (yi+yj)/2);
		point m2((xi+xk)/2, (yi+yk)/2);
		point m3((xj+xk)/2, (yj+yk)/2);

		triangleList[i] = triangle(m1, m2, m3);
		triangleList[c+3*i] = triangle(m1, e.p1, e.p2);
		triangleList[c+3*i+1] = triangle(m2, e.p1, e.p3);
		triangleList[c+3*i+2] = triangle(m3, e.p2, e.p3);
	}
	c*=4;
}

//initializes triangle array given the two corners of a rectangle
//NOTE: no rotated rectangles allowed b.c only 2 pt given
//c2 c3 are the other corners, o is the midpt of the diagonal
void initializeRectangleTriangleList(point c0, point c1){

	point c2(c0.x, c1.y);
	point c3(c1.x, c0.y);
	point o((c0.x+c1.x)/2, (c0.y+c1.y)/2);

	triangleList[0] = triangle(c0, c2, o);
	triangleList[1] = triangle(c0, c3, o);
	triangleList[2] = triangle(c1, c2, o);
	triangleList[3] = triangle(c1, c3, o);

	c=4;
	while(c*4<Nt){
		explode();
	}
}

//solves equation Ax=b for an nxn matrix, 
//where A=aVect, b=bVect, and n=dim (from MLib)
  double * LU(double * a_data, double * b_data, int dim){

	gsl_matrix_view m = gsl_matrix_view_array (a_data, dim, dim);
	gsl_vector_view b = gsl_vector_view_array (b_data, dim);

	gsl_vector * x = gsl_vector_alloc (dim);
  
	int s;

	gsl_permutation * p = gsl_permutation_alloc (dim);

	gsl_linalg_LU_decomp (&m.matrix, p, &s);
	gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x); 

    return x->data;
  }

//sets matrix for each triangle
void setTriangleMatrix(triangle e, double * matrix){
	double m[] = 	{bi*bi+ci*ci, bi*bj+ci*cj, bi*bk+ci*ck,
			 bi*bj+ci*cj, bj*bj+cj*cj, bj*bk+cj*ck,
			 bi*bk+ci*ck, bj*bk+cj*ck, bk*bk+ck*ck,
			 -4*r*a*a/3,   -4*r*a*a/3,  -4*r*a*a/3 };

	for(int i=0; i<12; i++){
		*matrix = m[i];
		matrix++;
	}
}

//adds all singular triangle matrices onto a final matrix to be solved
void setFinalMatrix(double *pFinalArray, int v){
	double triangleMatrix[12];
	triangle e;

	int i1;
	int i2;
	int i3;
	
	for(int i=0; i<c; i++){
		e = triangleList[i];

		i1 = e.p1.index;
		i2 = e.p2.index;
		i3 = e.p3.index;

		setTriangleMatrix(e, triangleMatrix);

		/* if anyone knows a better way to do this...*/ 
		*(pFinalArray + i1*v+i1) += triangleMatrix[0];
		*(pFinalArray + i1*v+i2) += triangleMatrix[1];
		*(pFinalArray + i1*v+i3) += triangleMatrix[2];
		*(pFinalArray + i2*v+i1) += triangleMatrix[3];
		*(pFinalArray + i2*v+i2) += triangleMatrix[4];
		*(pFinalArray + i2*v+i3) += triangleMatrix[5];
		*(pFinalArray + i3*v+i1) += triangleMatrix[6];
		*(pFinalArray + i3*v+i2) += triangleMatrix[7];
		*(pFinalArray + i3*v+i3) += triangleMatrix[8];
		*(pFinalArray + (i1+v*v)) += triangleMatrix[9];
		*(pFinalArray + (i2+v*v)) += triangleMatrix[10];
		*(pFinalArray + (i3+v*v)) += triangleMatrix[11];	
	}
}

//given finalArray after solving setFinalMatrix and number of vertices v
//splits the augmented matrix into the coefficient matrix and b vector
//uses gsl method to solve final matrix
void solveFinalArray(double * finalArray, int v){
	LU(finalArray, &finalArray[v*v], v);
}

//prints list of triangle points
//in form: (x1_coor, y1_coor) (x2_coor, y2_coor) (x3_coor, y3_coor)
void printTriangle(triangle e){
	cout << "(" << xi << ", " << yi << ")" << "\t";
	cout << "(" << xj << ", " << yj << ")" << "\t";
	cout << "(" << xk << ", " << yk << ")" << "\n";
}

//prints an indexed list of triangles and their associated vertices
//in form: index (x1_coor, y1_coor), (x2_coor, y2_coor), (x3_coor, y3_coor)
void printTriangleList(){
	for(int i=0; i<c; i++){
		cout << i << "\t";
		printTriangle(triangleList[i]);
	}
}

int main(){

	/* testing triangle initialization */
	initializeRectangleTriangleList(point(0,0), point(1,1));
	printTriangleList();

	/* test number vertices */
	cout << "num vertices: " << findNumVertices() << endl;

	/* test triangle area */
	triangle e(point(0, 0), point(0, 1), point(1, 0));
	cout << "test area: " << a << endl;
	e = triangle(point(0, 0), point (0, 2), point(1, 0));
	cout << "test area: " << a << endl;


	/* test ai, aj, ... etc, compared to textbk excercise */
	point p1(.5, .5);
	point p2(0, 0);
	point p3(1, 0);

	e = triangle(p1, p2, p3);
	
	cout << "test abcs: " << endl;
	cout << ai << "\t" << aj << "\t" << ak << endl;
	cout << bi << "\t" << bj << "\t" << bk << endl;
	cout << ci << "\t" << cj << "\t" << ck << endl;


	//second test matrix
	p1 = point(.5, .5);
	p2 = point(1, 0);
	p3 = point(1, 1);

	e = triangle(p1, p2, p3);

	cout << "test abcs: " << endl;
	cout << ai << "\t" << aj << "\t" << ak << endl;
	cout << bi << "\t" << bj << "\t" << bk << endl;
	cout << ci << "\t" << cj << "\t" << ck << endl;

	/* test setTriangleMatrix, compared to textbook exercise */
	double m[12];
	double * p = &m[0];
	setTriangleMatrix(e, m);

	cout << "test setTriangleMatrix:" << endl;
	for(int j=0; j<4; j++){
		for(int i=0; i<3; i++){
			cout << *p << "\t";
			p++;
		}
		cout << endl;
	}
	
	/* test setFinalMatrix */
	int v = findNumVertices();
	double *finalArray = new double[v*(v+1)];
	for(int j=0; j<v+1; j++){
		for(int i=0; i<v; i++){
			finalArray[j*v+i]=0;
		}
	}

	for(int j=0; j<v+1; j++){
		for(int i=0; i<v; i++){
			cout << finalArray[j*v+i] << "\t";
		}
		cout << endl;
	}

	setFinalMatrix(finalArray, v);
	
	for(int j=0; j<v+1; j++){
		for(int i=0; i<v; i++){
			cout << finalArray[j*v+i] << "\t";
		}
		cout << endl;
	}
	cout << endl << endl;

	/* test solveFinalArray */
	solveFinalArray(finalArray, v);
 	for(int j=0; j<v+1; j++){
		for(int i=0; i<v; i++){
			cout << finalArray[j*v+i] << "\t";
		}
		cout << endl;
	}	

	delete finalArray;
	return 0;
}