Stiffness matrix / load vector calculation

[Darkwings]

Homework Statement

Problem is (to be solved in software, this is actually a programming assignment):

$$-\triangle u = f$$ on $$\Omega = (0,1)^2$$
$$u = 0$$ on $$\delta\Omega$$
And
$$f = 2x(1-x)+2y(1-y)$$
And the exact solution known to be $$u(x,y) = x(1-x)y(1-y)$$
Okay, so the actual assignment is about iterative methods, but that's not a problem so far. I'm stuck on transforming this problem to a matrix-vector form.

Homework Equations

See problem description. Furthermore, the triangulation is with isosceles triangles, and the case I will present here is for a triangulation with just 8 triangles. (So that's just the main square splitted in four little squares, and these squares splitted to triangles by drawing a line from the upper left to the bottom right).

The relevant basis function are linears.

The Attempt at a Solution

Okay, first the stiffness matrix.

I started by creating the element stiffness matrices for the two different types of triangles - the ones with the points up and the ones with the points down. I did this by simply creating two linear systems in maple:

For triangles with point up:

a   b   1
a+h b   1
a   b+h 1

And with the point down:

a+h b   1
a+h b+h 1
a   b+h 1

In which (a,b) is the coordinate of the bottom right corner of the smallest square encompassing the triangle. For triangles pointing up this is just the bottom right corner of the triangle itself.
h Is the size of the smallest sides of the triangle.

Anyway, I think it's clear that you can create three different basis function per triangle type, for which you can calculate:

$$a(\phi(i),\phi(j)) = \int_{triangle}\nabla\phi(i) . \nabla\phi(j)$$

Leading to the following two element stiffness matrices:

Point up:

1   -.5 -.5
-.5  .5   0
-.5   0  .5

Point down:

.5  -.5   0
-.5   1 -.5
0   -.5  .5

Okay, now I'm trying to work the magic to create the global stiffness matrix from this with the following C code.

    // set to zero
    for (int i = 0; i < size; i ++) {
        for (int j = 0; j < size; j ++) A[i][j] = 0.0;
    }

    size = sqrt (size); // now we've got the size of the matrix in one direction

    // triangles pointing up
    double E1[9] = { 1.0,  -0.5, -0.5,
                    -0.5,   0.5,  0.0,
                    -0.5,   0.0,  0.5};

    // triangles points down
    double E2[9] = { 0.5,  -0.5,    0.0,
                    -0.5,   1.0,   -0.5,
                     0.0,  -0.5,    0.5};

    // go over all elements
    int ecount = ((size - 1) * 2) * (size - 1);

    int i,j,k;
    double x,y;
    for (int e = 0; e < ecount; e ++) {
        i = e/2 + ((e/2) / (size-1));
        if (e % 2 == 0) { // pointing up
            j = i + 1;
            k = i + size;

            A[i][i] += E1[0];
            A[i][j] += E1[1];
            A[j][i] += E1[1]; // symmetric anyway
            A[i][k] += E1[2];
            A[k][i] += E1[2]; // symmetric anyway
            A[j][j] += E1[4];
            A[j][k] += E1[5];
            A[k][j] += E1[5]; // symmetric
            A[k][k] += E1[8];
        } else { // pointing down
            i = e/2 + ((e/2) / (size-1)) + 1;
            j = i + size;
            k = j - 1;

            A[i][i] += E2[0];
            A[i][j] += E2[1];
            A[j][i] += E2[1]; // symmetric anyway
            A[i][k] += E2[2];
            A[k][i] += E2[2]; // symmetric anyway
            A[j][j] += E2[4];
            A[j][k] += E2[5];
            A[k][j] += E2[5]; // symmetric
            A[k][k] += E2[8];
        }
    }

As you can see, I'm looping over all elements in the mesh. Then I start by calculating to which point in the global matrix I must add the current elements contributions, and finally I add these points.

For reference, this results in:

 1.000000       -0.500000        0.000000       -0.500000        0.000000        0.000000        0.000000        0.000000        0.000000
-0.500000        2.000000       -0.500000        0.000000       -1.000000        0.000000        0.000000        0.000000        0.000000
 0.000000       -0.500000        1.000000        0.000000        0.000000       -0.500000        0.000000        0.000000        0.000000
-0.500000        0.000000        0.000000        2.000000       -1.000000        0.000000       -0.500000        0.000000        0.000000
 0.000000       -1.000000        0.000000       -1.000000        4.000000       -1.000000        0.000000       -1.000000        0.000000
 0.000000        0.000000       -0.500000        0.000000       -1.000000        2.000000        0.000000        0.000000       -0.500000
 0.000000        0.000000        0.000000       -0.500000        0.000000        0.000000        1.000000       -0.500000        0.000000
 0.000000        0.000000        0.000000        0.000000       -1.000000        0.000000       -0.500000        2.000000       -0.500000
 0.000000        0.000000        0.000000        0.000000        0.000000       -0.500000        0.000000       -0.500000        1.000000

Which look at least like a shot in the right direction to me.

Calculation of the load vector is done in a similar way. I took the same six (three for both triangle types) basis functions. Multiplied all six with the right-hand side of the problem. Then took the integral over a triangle. This results in six functions depending on a,b and h, and results in the following load-vector (again after adding triangle contributions in the same way as with the stiffness matrix, but now contributions are calculated run-time):

0:  0.016667
1:  0.083333
2:  0.045833
3:  0.083333
4:  0.208333
5:  0.083333
6:  0.045833
7:  0.083333
8:  0.016667

Now I've got a system that is solvable, but the solution is nowhere near the exact solution... I don't even know if I've got the stiffness matrix or the load vector wrong (or both), so *any* help would be greatly appreciated.

If more details are needed I can of course give these :)

Thanks!

 

[bob012345]

I worked out the local stiffness matrices and get the same matrix for up and down triangles which matched your case of the down triangle. Is it possible you did not follow the counterclockwise convention for nodes {i,j,k} evaluating both local matrices?