Projection onto P1 CG space defined on P2 geometry

3 months ago by

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))

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
I thought that high order geometries were a thing to be implemented in fenicsX. Does high order mean very high order? I thought it was just higher than 1.
written 3 months ago by Miguel  

1 Answer

3 months ago by
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),

# 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})

Absolutely wonderful! Thanks!
written 3 months ago by Nate  
Yes, this is known issue:
written 3 months ago by Jan Blechta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »