const double fac_lin[] = {0.5,0.5,0.5,1.0/24.0,1.0/6.0};
const double ms1_lin[] = {1,-1,0,-1,1,0,0,0,0};
const double ms2_lin[] = {2,-1,-1,-1,0,1,-1,1,0};
const double ms3_lin[] = {1,0,-1,0,0,0,-1,0,1};
const double ms4_lin[] = {2,1,1,1,2,1,1,1,2};
const double vs1_lin[] = {1,1,1};

const double fac_quad[] = {1.0/6.0,1.0/6.0,1.0/6.0,1.0/360.0,1.0/6.0};
const double ms1_quad[] = {3,1,0,-4,0,0,1,3,0,-4,0,0,0,0,0,0,0,0,-4,-4,0,8,0,0,0,0,0,0,8,-8,0,0,0,0,-8,8};
const double ms2_quad[] = {6,1,1,-4,0,-4,1,0,-1,-4,4,0,1,-1,0,0,4,-4,-4,-4,0,8,-8,8,0,4,4,-8,8,-8,-4,0,-4,8,-8,8};
const double ms3_quad[] = {3,0,1,0,0,-4,0,0,0,0,0,0,1,0,3,0,0,-4,0,0,0,8,-8,0,0,0,0,-8,8,0,-4,0,-4,0,0,8};
const double ms4_quad[] = {6,-1,-1,0,-4,0,-1,6,-1,0,0,-4,-1,-1,6,-4,0,0,0,0,-4,32,16,16,-4,0,0,16,32,16,0,-4,0,16,16,32};
const double vs1_quad[] = {0,0,0,1,1,1};

    for(int y=0; y<MSize; y++) {
        for(int x=0; x<MSize; x++) {
            gsl_matrix_set(S1,x,y,ms1[y*MSize+x]);
            gsl_matrix_set(S2,x,y,ms2[y*MSize+x]);
            gsl_matrix_set(S3,x,y,ms3[y*MSize+x]);
            gsl_matrix_set(S4,x,y,ms4[y*MSize+x]);
        }
        gsl_vector_set(s1,y,vs1[y]);
    }

    gsl_matrix_scale(S1,fac[0]);
    gsl_matrix_scale(S2,fac[1]);
    gsl_matrix_scale(S3,fac[2]);
    gsl_matrix_scale(S4,fac[3]);
    gsl_vector_scale(s1,fac[4]);

    gsl_matrix_memcpy(Me,S4);
    Stot = gsl_matrix_calloc(numMeshVertices,numMeshVertices);
    Mtot = gsl_matrix_calloc(numMeshVertices,numMeshVertices);

    gsl_matrix_set(S5l,0,0,1.0/3.0); gsl_matrix_set(S5l,0,1,1.0/6.0);
    gsl_matrix_set(S5l,1,1,1.0/3.0); gsl_matrix_set(S5l,1,0,1.0/6.0);

    gsl_matrix_set(S5q,0,0,2.0);  gsl_matrix_set(S5q,0,1,1.0); gsl_matrix_set(S5q,0,2,-0.5);
    gsl_matrix_set(S5q,1,0,1.0);  gsl_matrix_set(S5q,1,1,8.0); gsl_matrix_set(S5q,1,2,1.0);
    gsl_matrix_set(S5q,2,0,-0.5); gsl_matrix_set(S5q,2,1,1.0); gsl_matrix_set(S5q,2,2,2.0);

    gsl_matrix_scale(S5q,1.0/15.0);