# Leave boundary dofs out of matrix construction

I'm implementing a Navier--Stokes solver for axisymmetric geometries, and in the weak forumulation, the term

u[0] * v[0] / r^2

appears (cf. <https:/

In Dolfin, the resulting matrix would look like

[I 0 ] [u_r] = [0]

[b A] [u ] [f]

where b corresponds to the terms where u[0] is nonzero at the boundary, and thus b would contain NaNs. They don't play a particular role since u_r is 0 there, but I don't if every backend picks that up. Also, applying some sort of LU solver, in fact any method that does something with the matrix entries (e.g., AMG), is asking for all sorts of trouble.

Question:

Is is possible to restrict the trial space as well to directly get an equation system of the form

[A][u] = [f]

?

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Garth Wells

- Solved:
- 2013-05-08

- Last query:
- 2013-05-08

- Last reply:
- 2013-05-05

Jan Blechta (blechta) said : | #1 |

On Sun, 05 May 2013 02:21:17 -0000

Nico Schlömer <email address hidden> wrote:

> New question #228170 on DOLFIN:

> https:/

>

> I'm implementing a Navier--Stokes solver for axisymmetric geometries,

> and in the weak forumulation, the term

>

> u[0] * v[0] / r^2

You should simplify bad terms by multiplication by volume element

2*pi*r*dx - UFL does not reduce fractions except very special cases.

>

> appears (cf.

> <https:/

Note this is a Laplace form - stress form is more suitable for some

problems.

> If the domain contains r=0 -- which is perfectly legal --, e.g. the

> unit square, the 1/r^2 looks a little frightening at first. This

> isn't critical if the trial space and test space only contain

> functions at least proportional to r, corresponding to homogenous

> Dirichlet conditions for u[0] at r=0. This is quite easy to enforce

> on the test space (namely by applying said Dirichlet boundary

> conditions). In Dolfin, the resulting matrix would look like

>

> [I 0 ] [u_r] = [0]

> [b A] [u ] [f]

>

> where b corresponds to the terms where u[0] is nonzero at the

> boundary, and thus b would contain NaNs. They don't play a particular

I think there are no NaNs as Gauss qaudrature is used therefore 1/r

term in not evaluated at r=0. When Dirichlet conditions for u_r are to

be enforced on r=0 bad terms are not evaluated at r=0, only

corresponding dofs are found and manipulated in linear system's matrix

and vector.

> role since u_r is 0 there, but I don't if every backend picks that

> up. Also, applying some sort of LU solver, in fact any method that

> does something with the matrix entries (e.g., AMG), is asking for all

> sorts of trouble.

>

> Question:

> Is is possible to restrict the trial space as well to directly get an

> equation system of the form

>

> [A][u] = [f]

>

> ?

>

Nico Schlömer (nschloe) said : | #2 |

> You should simplify bad terms by multiplication by volume element

> 2*pi*r*dx - UFL does not reduce fractions except very special cases.

Makes sense. -- After all, looking at the problem from a 3D angle, one would really like to weight the equations (and with it the residual) with the volume of the rotational 3D element.

> I think there are no NaNs as Gauss qaudrature is used therefore 1/r term in not evaluated at r=0.

Ah, okay. In face, the PDE isn't defined in r=0 either, so the only thing one has to take care of is that the functional doesn't grow indefinitely towards r=0. Is there any documentation about the Gauss points?

Garth Wells (garth-wells) said : | #3 |

On 5 May 2013 21:11, Nico Schlömer <email address hidden> wrote:

> Question #228170 on DOLFIN changed:

> https:/

>

> Status: Answered => Open

>

> Nico Schlömer is still having a problem:

>> You should simplify bad terms by multiplication by volume element

>> 2*pi*r*dx - UFL does not reduce fractions except very special cases.

>

> Makes sense. -- After all, looking at the problem from a 3D angle, one

> would really like to weight the equations (and with it the residual)

> with the volume of the rotational 3D element.

>

>> I think there are no NaNs as Gauss qaudrature is used therefore 1/r

> term in not evaluated at r=0.

>

> Ah, okay. In face, the PDE isn't defined in r=0 either, so the only

> thing one has to take care of is that the functional doesn't grow

> indefinitely towards r=0. Is there any documentation about the Gauss

> points?

>

Take a look at the FFC code:

ffc/

Garth

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

Nico Schlömer (nschloe) said : | #4 |

Thanks Garth Wells, that solved my question.