### Projection onto P1 CG space defined on P2 geometry

123

views

0

Consider the following:

```
from dolfin import *
geom_degree = 2
mesh = UnitDiscMesh.create(MPI.comm_world, 1, geom_degree, 2)
f = project(Constant(1.0), FunctionSpace(mesh, "CG", 1))
print(f.vector()[:])
```

My numerical analysis is rather poor. Naively I expect `f.vector()[:]`

to contain entirely `1.0`

. What I see is

```
[0.9141741 1.18315029 0.9141741 0.9141741 0.9141741 0.9141741
0.9141741 ]
```

Can anyone provide some insight into what I should expect?

Community: FEniCS Project

written
3 months ago by
Miguel

### 1 Answer

5

Here is what I think is going on. The default quadrature degree for the RHS vector is 1, because it is a linear shape function multiplied by a constant, but that is not sufficient to exactly integrate the non-constant Jacobian determinant from using a geometrical degree of 2 in the unit disc mesh. Take a look at the following code:

```
from dolfin import *
# Produces expected result with default quadrature:
#geom_degree = 1
# Requires more quadrature for expected result:
geom_degree = 2
# (Changed for 2017.2)
#mesh = UnitDiscMesh.create(MPI.comm_world, 1, geom_degree, 2)
mesh = UnitDiscMesh.create(mpi_comm_world(), 1, geom_degree, 2)
f_default = project(Constant(1.0), FunctionSpace(mesh, "CG", 1))
f = project(Constant(1.0), FunctionSpace(mesh, "CG", 1),
form_compiler_parameters={'quadrature_degree':2})
# This seriously under-integrates the projection LHS, and produces really
# incorrect results:
#f = project(Constant(1.0), FunctionSpace(mesh, "CG", 1),
# form_compiler_parameters={'quadrature_degree':1})
print(f_default.vector().get_local())
print(f.vector().get_local())
```

Absolutely wonderful! Thanks!

written
3 months ago by
Nate

Yes, this is known issue: https://bitbucket.org/fenics-project/ufl/issues/80/degree-estimation-doesnt-cater-for-non

written
3 months ago by
Jan Blechta

Please login to add an answer/comment or follow this question.