Hi group!
This is a crosspost from physics forum, so anything between [tex] tags
is latex stuff. Hope it's not too cluttered :)
Problem is solving the following in software:
[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 thingie is about iterative methods, but that's not
a problem so far. I'm stuck on transforming this problem to a matrix-
vector form.
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.
My own attempt follows:
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:
[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:
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 looks 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!
