- #1
Darkwings
- 1
- 0
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:
Code:
a b 1
a+h b 1
a b+h 1And with the point down:
Code:
a+h b 1
a+h b+h 1
a b+h 1In 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:
Code:
1 -.5 -.5
-.5 .5 0
-.5 0 .5Point down:
Code:
.5 -.5 0
-.5 1 -.5
0 -.5 .5Okay, now I'm trying to work the magic to create the global stiffness matrix from this with the following C code.
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:
Code:
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.000000Which 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:
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.016667Now 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: