Construct M matrix from given discrete velocity set order for MRT D3Q19 lattice model in lattice Boltzmann method

getM.c

/*
 *construct M from c
 *refer to D'Humières, Dominique Ginzburg, Irina Krafczyk, Manfred Lallemand, Pierre Luo, Li-Shi
  2012
 */

#include <stdio.h>


#define Q 19

int e[Q][3] = {
    { 0,  0,  0},

    { 1,  0,  0},
    {-1,  0,  0},
    { 0,  1,  0},
    { 0, -1,  0},
    { 0,  0,  1},
    { 0,  0, -1},

    { 0,  1,  1},
    { 0, -1, -1},
    { 0, -1,  1},
    { 0,  1, -1},

    {-1,  0, -1},
    { 1,  0,  1},
    {-1,  0,  1},
    { 1,  0, -1},

    {-1,  1,  0},
    { 1, -1,  0},
    {-1, -1,  0},
    { 1,  1,  0}
};

int main()
{

    int i;
    int easq;
    //# 0
    for(i=0; i<Q; i++){
        printf("% 4d", 1);
    }
    printf("\n");

    //# 1
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(19*(e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2]) - 30));
    }
    printf("\n");

    //# 2
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)((21*easq*easq - 53*easq + 24)/2));
    }
    printf("\n");

    //# 3
    for(i=0; i<Q; i++){
        printf("% 4d", (int)e[i][0]);
    }
    printf("\n");

    //# 4
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)((5*easq - 9)*e[i][0]));
    }
    printf("\n");

    //# 5
    for(i=0; i<Q; i++){
        printf("% 4d", (int)e[i][1]);
    }
    printf("\n");

    //# 6
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)((5*easq - 9)*e[i][1]));
    }

    printf("\n");
    //# 7
    for(i=0; i<Q; i++){
        printf("% 4d", (int)e[i][2]);
    }
    printf("\n");

    //# 8
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)((5*easq - 9)*e[i][2]));
    }
    printf("\n");

    //# 9
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)(3*e[i][0]*e[i][0] - easq));
    }
    printf("\n");

    //#10
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)(3*easq -5)*(3*e[i][0]*e[i][0] - easq));
    }
    printf("\n");

    //#11
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][1]*e[i][1] - e[i][2]*e[i][2]));
    }
    printf("\n");

    //#12
    for(i=0; i<Q; i++){
        easq = e[i][0]*e[i][0] + e[i][1]*e[i][1] + e[i][2]*e[i][2];
        printf("% 4d", (int)(3*easq -5)*(e[i][1]*e[i][1]  - e[i][2]*e[i][2]));
    }
    printf("\n");

    //#13
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][0]*e[i][1]));
    }
    printf("\n");


    //#14
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][1]*e[i][2]));
    }
    printf("\n");


    //#15
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][0]*e[i][2]));
    }
    printf("\n");

    //#16
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][1]*e[i][1] - e[i][2]*e[i][2])*e[i][0]);
    }
    printf("\n");

    //#17
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][2]*e[i][2] - e[i][0]*e[i][0])*e[i][1]);
    }
    printf("\n");

    //#18
    for(i=0; i<Q; i++){
        printf("% 4d", (int)(e[i][0]*e[i][0] - e[i][1]*e[i][1])*e[i][2]);
    }
    printf("\n");



    return 0;
}