1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Stiffness matrix / load vector calculation

  1. Jun 20, 2008 #1
    1. The problem statement, all variables and given/known data
    Problem is (to be solved in software, this is actually a programming assignment):

    [tex]-\triangle u = f[/tex] on [tex]\Omega = (0,1)^2[/tex]
    [tex]u = 0[/tex] on [tex]\delta\Omega[/tex]
    And
    [tex]f = 2x(1-x)+2y(1-y)[/tex]
    And the exact solution known to be [tex]u(x,y) = x(1-x)y(1-y)[/tex]
    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.

    2. Relevant 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.


    3. 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:
    Code (Text):

    a   b   1
    a+h b   1
    a   b+h 1
     
    And with the point down:
    Code (Text):

    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:

    [tex]a(\phi(i),\phi(j)) = \int_{triangle}\nabla\phi(i) . \nabla\phi(j)[/tex]

    Leading to the following two element stiffness matrices:

    Point up:
    Code (Text):

    1   -.5 -.5
    -.5  .5   0
    -.5   0  .5
     
    Point down:
    Code (Text):

    .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.

    Code (Text):

        // 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:
    Code (Text):

     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):
    Code (Text):

    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!
     
    Last edited: Jun 20, 2008
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you help with the solution or looking for help too?



Similar Discussions: Stiffness matrix / load vector calculation
  1. Matrix Question (Replies: 0)

  2. Vandermonde matrix (Replies: 0)

  3. Vector Potential (Replies: 0)

Loading...