dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 21, 2005 at 10:51:21AM -0500, Andy Ray Terrel wrote:
> Hello Anders,
> 
> Okay so after getting everything to compile and running I am getting an 
> error whenever I assemble my second equation.  It is as follows:
> 
> *** Assertion (test_element.spacedim() == trial_element.spacedim()) 
> failed [FEM.cpp:30: assemble()]
> Deliberately raising a segmentation fault.

Aha! This is interesting. I have deliberately put in this assertion in
the assembly code since there are some assumptions that the trial and
test spaces are of the same dimension. This is not the case here,
since you want to assemble a matrix A that is not square (and then use
this rectangular matrix as a building block for a square matrix).

You need to modify the assembly code and remove this assumption. The
code is pretty clean and I think you should be able to see what's
going on. Take a look at the first function in the file
src/kernel/fem/FEM.cpp. I think the only thing you need to do is to
add a pair of variables m, M in addition to n, N. Something like

uint m = test_element.spacedim();
uint n = trial_element.spacedim();
real* block = new real[m*n];
int* test_dofs = new int[m];
int* trial_dofs = new int[n];

uint M = size(mesh, test_element);
uint N = size(mesh, test_element);
A.init(M, N, 1);
A = 0.0;

and so on. When it works, you can either send me a patch or just
commit your changes to CVS. Remember to add a line

// Modified by Andy Terrel, 2005.

at the top of the file.

> Here is the form:
> 
> vector = FiniteElement("Vector Lagrange", "triangle", 1)
> scalar = FiniteElement("Lagrange", "triangle", 1)
> 
> q = BasisFunction(scalar)
> u = BasisFunction(vector)
> 
> a = q * u[i].dx(i) * dx
> 
> Also, how do I call the solver for these guys?  I don't have a right 
> hand side to give to the GMRES solver.   Should is there a way to 
> generate the appropriate right hand side from the size of the matrix.

You need to define a form for the right-hand side as well, something
like

v = BasisFunction(appropriate element)
f = Function(appropriate element)

L = f*v*dx

You then get a Function f in DOLFIN that you can set to whatever you
want (0, 1, sin(x) etc).

> Something like MatGetSize(A1.mat(),&m,&n), and then build a the RHS 
> vector from m and n?
> 
> Finally, how do I get more control over the PETSc solver? Say I want to 
> use the Jacobi Preconditioner as apposed to the ILU default. 

You should be able to set the preconditioner by getting the PETSc KSP
pointer by doing something like

GMRES gmres;
KSP ksp = gmres.solver();

Then you can do whatever you want with the ksp pointer:

PC pc;
KSPGetPC(ksp, &pc);
PCSetType(pc, PCILU); // Or something else

When you then call gmres.solve(A, x, b), the solver will use the
preconditioner you specified. We should probably add a function to
GMRES that allows a user to specify the preconditioner without going
through the PETSc calls. (Add one if you like.)

/Anders