main.cpp

#include <iostream>
#include <stdlib.h>
#include <fstream>
#include <cstdlib>
using namespace std;
#include <string.h>
#include <math.h>
#include <omp.h>
#include <chrono>
#include <vector>
#include <set>
#include <algorithm>

#include "topol.h"
#include "gmres.h"
#include "utilities.h"


// Main---------------------------------------------------------------------------------------
int main(int argc, const char* argv[]){

    //Argument check using the terminal
    if (argc < 6) {
        printf("Too few arguments. \n Usage: [tetra file] [coord file] [np] [maxit] [restart] \n");
	exit(1);
    }

    int np = atoi(argv[3]);        // number of threads
    int maxit = atoi(argv[4]);     // max iteration GMRES
    int restart = atoi(argv[5]);   // Restart parameter GMRES
    // Read tetrahedrons
    ifstream tetra_file(argv[1]);
    // Read header
    int ntet, nodes_per_elem, base;
    tetra_file >> ntet >> nodes_per_elem >> base;         // number of tetra elements
    int *tetra_buf = (int*) malloc(4*ntet*sizeof(int));   // Continue allocation of tetra elements
    int **tetra = (int**) malloc(ntet*sizeof(int*));      // Puntatori ai nodi di ciascun tetraedro
    int k = 0;
    int junk;

    for (int i=0; i<ntet; i++){
        tetra[i] = &(tetra_buf[k]);
	k += 4;
	tetra_file >> junk;
        for (int j=0; j < 4; j++){
	    tetra_file >> tetra[i][j];
	}
	tetra_file >> junk;
    }

    // Close the input file
    tetra_file.close();

    //Set C-style of tet: needed to convert file in base 1 (Fortran style) in base 0
    for (int i = 0; i < ntet; i++)
        for (int j = 0; j < 4; j++) tetra[i][j]--;
        
    //Number of equations
    int nn = 0;
    for (int i = 0; i < ntet; i++)
        for (int j = 0; j < 4; j++)
	    if (tetra[i][j] > nn) nn = tetra[i][j];  
    nn++;


    printf("Number of DOF: %d \n", nn);
    
   

    // Topology -----------------------------------------------------------------
    int nterm;
    int n1 = 600;
    int *iat = nullptr;
    int *ja = nullptr;
    topol(nn, ntet, n1, tetra, nterm, iat, ja);
    //topol2(nn, ntet, tetra, nterm, iat, ja);
    
    printf("Topology created \n");

    // Now we have iat and ja ready to be used
    // -----------------------------------------
    
    ifstream coord_file(argv[2]);
    
    // Read header 
    int nnodes, dim, junk1, junk2; 
    coord_file >> nnodes >> dim >> junk1 >> junk2;

    double *coord_buff = (double*) malloc(3*nnodes*sizeof(double)); // we are in 3D
    double **coord = (double**) malloc(nnodes*sizeof(double*));
    
    for (int i = 0; i < nn; i++) coord[i] = nullptr;
    
    for (int i = 0; i < nnodes; i++) {
        int id;
        coord_file >> id;
        id--;
        
        coord[id] = &coord_buff[3 * id];
        
        coord_file >> coord[id][0]
                   >> coord[id][1]
                   >> coord[id][2];
    }

    
    coord_file.close();
    
    for (int k = 0; k < ntet; k++) {
        for (int i = 0; i < 4; i++) {
            if (tetra[k][i] >= nnodes) {
                printf("Tetra %d use nodo %d ma coor has only %d nodi\n", k, tetra[k][i], nnodes);
                exit(1);
            }
        }
    }
    
    for (int i = 0; i < nn; i++) {
        if (coord[i] == nullptr) {
             printf("Node FEM %d NO coordinates\n", i);
             exit(1);
        }
    }

    //------------------------------------------

    // Diffusion and flow velocity
    double* D = (double*) malloc(3*sizeof(double));
    D[0] = 0.1; D[1] = 0.1; D[2] = 0.6;
    double* v = (double*) malloc(3*sizeof(double));
    v[0]=1.0; v[1] = 3.0; v[2] = 1.0;

    printf("D: ");
    printf("%f %f %f\n", D[0],D[1],D[2]);
    printf("v: ");
    printf("%f %f %f\n", v[0],v[1],v[2]);

    //-----------------------------------------
    // Assembly: we have to build H B P csr matrix -> coefH, coefB, coefP
    double start;
    double end;

    start = omp_get_wtime();

    double* coefH = (double*) malloc (nterm *sizeof(double));
    double* coefB = (double*) malloc (nterm *sizeof(double)); 
    double* coefP = (double*) malloc (nterm *sizeof(double));

    for (int i = 0; i<nterm; i++){
        coefH[i]= 0;
	coefB[i] = 0;
	coefP[i] = 0;
    }

    # pragma omp parallel num_threads(np)
    {
        double* a = (double*) malloc(4*sizeof(double));
	double* b = (double*) malloc(4*sizeof(double));
	double* c = (double*) malloc(4*sizeof(double));
	double* d = (double*) malloc(4*sizeof(double));
	double *Hloc_buf = (double*) malloc(4*4*sizeof(double));
	double **Hloc = (double**) malloc(4*sizeof(double*));
	double *Ploc_buf = (double*) malloc(4*4*sizeof(double));
	double **Ploc =(double**) malloc(4*sizeof(double*));
	double *Bloc_buf = (double*) malloc(4*4*sizeof(double));
	double **Bloc = (double**) malloc(4*sizeof(double*));
	int kk =0;

	for (int i = 0; i<4; i++){
            Hloc[i] = &(Hloc_buf[kk]);
	    Ploc[i] = &(Ploc_buf[kk]);
	    Bloc[i] = &(Bloc_buf[kk]);
	    kk += 4;
	}
        # pragma omp for
	for (int jj=0; jj < ntet; jj++){
	    int* tet = tetra[jj];
	    
	    double* n0 = coord[tet[0]]; 
	    double* n1 = coord[tet[1]]; 
	    double* n2 = coord[tet[2]]; 
	    double* n3 = coord[tet[3]];

	    double vol = (det3(n1,n2,n3) - 
			    det3(n0,n2,n3) +
			    det3(n0,n1,n3) -
			    det3(n0,n1,n2))/6;

	    for (int ii=0; ii < 4; ii++){
	        int idx[3]; int cnt = 0;
	        for (int jj = 0; jj < 4; jj++){
		    if (ii == jj){
		        continue;
	            }
                    idx[cnt] = jj;
                    cnt += 1;
                }
                double segno = pow(-1, ii);
                double* nj = coord[tet[idx[0]]];
		double* nk = coord[tet[idx[1]]];
		double* nm = coord[tet[idx[2]]];

		a[ii] = segno * det3(nj,nk,nm);
		b[ii] = segno * (-1) * (nk[1]*nm[2] -nm[1]*nk[2] - (nj[1]*nm[2]-nj[2]*nm[1]) + nj[1]*nk[2]-nj[2]*nk[1]);
		c[ii] = segno * (nk[0]*nm[2]-nm[0]*nk[2] - (nj[0]*nm[2]-nj[2]*nm[0])+nj[0]*nk[2]-nj[2]*nk[0]);
		d[ii] = segno * (nk[1]*nm[0]-nm[1]*nk[0] - (nj[1]*nm[0]-nj[0]*nm[1])+nj[1]*nk[0]-nj[0]*nk[1]); // I swap col 1 with col 2 so the det change sign


            }
            
            for (int i = 0; i<4; i++){
	        double* Hi = Hloc[i];
		double* Pi = Ploc[i];
		double* Bi = Bloc[i];
		
		for (int j=0; j<4; j++){
		    Hi[j] = (D[0]*b[i]*b[j] + D[1]*c[i]*c[j] + D[2]*d[i]*d[j])/ (36 * abs(vol));
	            Pi[j] = abs(vol)/20;
	            Bi[j] = sign(vol)/24 * (v[0]*b[j]+v[1]*c[j]+v[2]*d[j]);	    
		}
	        Pi[i] *=2;	
            }
            
            // locate them inside coefH, coefP, coefB

            for (int i=0; i<4; i++) {
	        // tet[i] gives the row index
		double* Hi = Hloc[i];
                double* Pi = Ploc[i];
                double* Bi = Bloc[i];
                for (int j=0; j<4; j++){
		    // in ja from iat[tet[i]] look for tet[j]
		    int curr_node = tet[i];
		    for (int ii=iat[curr_node]; ii < iat[curr_node+1]; ii++){
			int jjj = tet[j];
			if (ja[ii] == jjj){
			    #pragma omp atomic
		            coefH[ii] += Hi[j];
                            #pragma omp atomic
			    coefP[ii] += Pi[j];
                            #pragma omp atomic
			    coefB[ii] += Bi[j];
			
			}

		    }
	    
		}

	    }
	}
	
	
	free(a);
	free(b);
	free(c);
	free(d);
	free(Hloc);
	free(Bloc);
	free(Ploc);
	free(Hloc_buf);
	free(Ploc_buf);
	free(Bloc_buf);
    }

    end = omp_get_wtime();
    printf("Assembly time: %f\n", end - start);

    // GMRES-----------------------------------------------------------------------------------
    
    double* coefA = (double*) malloc (nterm * sizeof(double));
    double* coefAmod = (double*) malloc(nterm *sizeof(double)); //matrix to modify in order to impose the DC bc.
    
    double dt = 1;
    for (int i = 0; i < nterm; i++){
	coefA[i] = coefB[i] + coefH[i] + coefP[i]/dt; //
	coefAmod[i] = coefA[i];

    }

    double* bdval = (double*) malloc(nnodes*sizeof(double));
    for (int i = 0; i < nnodes; i++){
        if ((coord[i][0] < 1e-10)&&(coord[i][1] < .5 + 1e-10)&&(coord[i][2]< 1e-10)){
	    bdval[i] = 1.0;
	}
	else {
	    bdval[i] = 0.0;
	}	    
    }

    double* Abdval = (double*) calloc(nnodes, sizeof(double));

    matcsrvecprod(nn, iat, ja, coefA, bdval, Abdval, np);

    // Impose BC----------------------------------------------------------------------------
    for (int iii =0; iii < nn; iii++){
	if (bdval[iii]){ // Set the row equal to zero except for the diagonal element
	    for (int j = iat[iii]; j < iat[iii+1]; j++){
	        if (ja[j] == iii){
	    	    coefAmod[j] = 1.0;
		}
		else {
		    coefAmod[j] = 0.0;
		}
	    }
	}
	else { // Set the columns where q[j] = 1 equal to 0
	    for (int j = iat[iii]; j < iat[iii+1]; j++){
	        if (bdval[ja[j]]){
	    	    coefAmod[j] = 0.0;
		}
	    }
	}
    }
    

    // RHS assembly-----------------------------------------------------------------------------
    double* rhs = (double*) calloc(nnodes, sizeof(double));

    matcsrvecprod(nn, iat, ja, coefP, bdval, rhs, np);
    for (int i = 0; i < nnodes; i++){
	if (bdval[i]){
	    rhs[i] = bdval[i];
	}
	else {
	    rhs[i] = rhs[i]/dt - Abdval[i]; // (P/dt)u_prev - Abdval
	}
    }
    

    // Solution vector---------------------------------------------------------------------------
    double* S = (double*) calloc(nnodes, sizeof(double));

    double tol = 1e-9;
    
    start = omp_get_wtime();
    gmres(nnodes, iat, ja, coefAmod, rhs, tol, maxit, restart, np, S);
    end = omp_get_wtime();
    
    printf("GMRES time: %f\n", end - start);
    
    
    // Relative error---------------------------------------------------------------------------
    double* AS = (double*) malloc(nnodes*sizeof(double));
    double* rtrue = (double*) calloc(nnodes, sizeof(double));
    matcsrvecprod(nn, iat, ja, coefAmod, S, AS, np);
    
    double num = 0.0, den = 0.0;
    for (int i=0; i<nnodes; i++) rtrue[i] = rhs[i] - AS[i];
    printf("Final residual ||b - A*x|| = %1.4e\n", norm(rtrue, nnodes, np));
    

    // free memory
    free(coord);
    free(coord_buff);
    free(tetra);
    free(tetra_buf);
    free(iat);
    free(ja);
    free(D);
    free(v);
    
    free(coefH);
    free(coefP);
    free(coefB);
    free(coefA);
    free(coefAmod);
    free(bdval);
    free(Abdval);
    free(S);
    free(rhs);
    free(AS);
    free(rtrue);

}