FemSolver/GradientMatrix.cpp at main · cad420/FemSolver · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
//
//  GradientMatrix.cpp
//
//
//

#include "GradientMatrix.h"

GradientMatrix::GradientMatrix()
{
    //create();
}

GradientMatrix::GradientMatrix(const double a_size, const double b_size, const double c_size, const int n_num) : a(a_size), b(b_size), c(c_size), n(n_num)
{
    create();
}

GradientMatrix::~GradientMatrix()
{
    //destroy();
}

PetscErrorCode GradientMatrix::destroy()
{
    ierr = MatDestroy(&matrix);  CHKERRQ(ierr);
    ierr = MatDestroy(&transpose); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode GradientMatrix::create()
{
    double temp;

    double tempA = a, tempB = b, tempC = c;

    // now the gradient (or B) matrices
    ierr = MatCreateSeqDense(PETSC_COMM_SELF, 6, 24, PETSC_NULL, &matrix);  CHKERRQ(ierr);

	double x[8],y[8],z[8];
    double a, b, c, d, e, f, g, h, l;
    double difxs, difys, difzs, difxt, difyt, difzt, difxn, difyn, difzn;
    for (int i = 0; i < 8; i++)
    {
        x[i] = (1 + SIGN(integMtx(i, 0))) / 2.0 * tempA;
        y[i] = (1 + SIGN(integMtx(i, 1))) / 2.0 * tempB;
        z[i] = (1 + SIGN(integMtx(i, 2))) / 2.0 * tempC;
    }
    difxs = difxt = difxn = difys = difyt = difyn = difzs = difzt = difzn = 0;
    for (int i = 0; i < 8; i++)
    {
        difxs += diffMtx(n, 0, i) * x[i];
        difxt += diffMtx(n, 1, i) * x[i];
        difxn += diffMtx(n, 2, i) * x[i];
        difys += diffMtx(n, 0, i) * y[i];
        difyt += diffMtx(n, 1, i) * y[i];
        difyn += diffMtx(n, 2, i) * y[i];
        difzs += diffMtx(n, 0, i) * z[i];
        difzt += diffMtx(n, 1, i) * z[i];
        difzn += diffMtx(n, 2, i) * z[i];
    }

    Jdet = difxs * difyt * difzn + difxt * difyn * difzs + difxn * difys * difzt - difxn * difyt * difzs - difxt * difys * difzn - difxs * difyn * difzt;

    a = difyt * difzn - difzt * difyn;
    b = difys * difzn - difzs * difyn;
    c = difys * difzt - difzs * difyt;
    d = difxt * difzn - difzt * difxn;
    e = difxs * difzn - difzs * difxn;
    f = difxs * difzt - difzs * difxt;
    g = difxt * difyn - difyt * difxn;
    h = difxs * difyn - difys * difxn;
    l = difxs * difyt - difys * difxt;

    for (int i=0; i< NODES_PER_ELEMENT; i++)
    {// 8 integration points

        /*tempA = 1.0 / tempA * SIGN(integMtx(i, 0)) *
            (0.5 + integMtx(n, 1) * SIGN(integMtx(i, 1))) *
            (0.5 + integMtx(n, 2) * SIGN(integMtx(i, 2)));

        tempB = 1.0 / tempB * SIGN(integMtx(i, 1)) *
            (0.5 + integMtx(n, 0) * SIGN(integMtx(i, 0))) *
            (0.5 + integMtx(n, 2) * SIGN(integMtx(i, 2)));

        tempC = 1.0 / tempC * SIGN(integMtx(i, 2)) *
            (0.5 + integMtx(n, 0) * SIGN(integMtx(i, 0))) *
            (0.5 + integMtx(n, 1) * SIGN(integMtx(i, 1)));*/


        temp = a * diffMtx(n, 0, i) - b * diffMtx(n, 1, i) + c * diffMtx(n, 2, i);
        temp = temp / Jdet;

        ierr = MatSetValue(matrix, 0, 0 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 3, 1 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 5, 2 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);

        temp = -d * diffMtx(n, 0, i) + e * diffMtx(n, 1, i) - f * diffMtx(n, 2, i);
        temp = temp / Jdet;


        ierr = MatSetValue(matrix, 1, 1 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 3, 0 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 4, 2 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);

        temp = g * diffMtx(n, 0, i) - h * diffMtx(n, 1, i) + l * diffMtx(n, 2, i);
        temp = temp / Jdet;


        ierr = MatSetValue(matrix, 2, 2 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 4, 1 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
        ierr = MatSetValue(matrix, 5, 0 + 3 * i, temp, INSERT_VALUES); CHKERRQ(ierr);
    }

    ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

    // create transpose
    MatTranspose(matrix, MAT_INITIAL_MATRIX, &transpose);

    return 0;
}