Quadrature Projection

4 months ago by

I'm tring to project a vector by using quadrature as following :

elementS = VectorElement("Quadrature", mesh.ufl_cell(), degree=QUAD_DEG, quad_scheme="default",dim=6)
VS = FunctionSpace(mesh,elementS)
def sigma_vec(s):
    return as_vector([s[0], s[1], s[2],s[3], s[4], s[5]])

But then I got the following error :

[[ 0.5854102 0.1381966 0.1381966]
[ 0.1381966 0.5854102 0.1381966]
[ 0.1381966 0.1381966 0.5854102]
[ 0.1381966 0.1381966 0.1381966]]

quad points:
[[ 0.25 0.25 0.25]]
Points must be equal to coordinates of quadrature points

I tried to change degree of quadrature but it did not help. What the "points" correspond here ? When I increase the degree of quadrature points also increases but points and quad points are not equal.


Community: FEniCS Project
Could you post the complete code (so that we see how you define your mesh and what you are printing)?
It may also be that you want to use interpolate instead of project.  See:
written 4 months ago by jwinkle  

3 Answers

4 months ago by
If I change the projection to


it works for me in 2017.2 (with some deprecation warnings).  Another way to set things at quadrature points that I've found necessary in some situations is the following:

def projectToQuadPts(toProject,VS):
    u = TrialFunction(VS)
    v = TestFunction(VS)
    lhsForm = inner(u,v)*dx
    rhsForm = inner(toProject,v)*dx
    p = {"representation":"quadrature"}
    A = assemble(lhsForm,form_compiler_parameters=p)
    B = assemble(rhsForm,form_compiler_parameters=p)
    u = Function(VS)
    # Because the matrix is diagonal, it should be inverted exactly by the
    # Jacobi preconditioner on the first iteration.  
    return u

sigma = projectToQuadPts(sigma_vec(s),VS)

I need this, for instance, when toProject is some complicated UFL algebra object that project() will not accept as an argument.  Also, with the above, UFLACS representation can be used in the rest of the code, which is more robust when differentiating very complicated forms.

Hello David,

Thank you a lot, it works fine. But as you mentioned it also gives some warnings and makes the compile slower. By adding "form_compiler_parameters={'quadrature_degree':QUAD_DEG}" what did we change actually on projection ?
written 4 months ago by Christian  
My guess is that the project() function defaults to some quadrature rule other than the one that is used to define the Quadrature-type FunctionSpace, so it ends up trying to evaluate functions in that space at points where they are not defined, unless we manually specify that the projection use a quadrature rule of the same degree as the FunctionSpace.
written 4 months ago by David Kamensky  
4 months ago by
What are you trying to project onto your 3D FunctionSpace?  Your extracted vector has 6 elements but your mesh will be 5x3x3 nodes, and further VS is vector-valued.  What are the 6 elements returned from sigma_vec supposed to represent?

Also, to help users see your problem, see this pinned suggestion:

--e.g., I didn't even realize your error was from the compiler (I thought you printed it yourself), so it's good to post the full error message.
Actually I was trying to test the projection on quadrature for future purpose to calculate sigma on quadratures. Regarding the error, I just copied the beginning of the error. It directly starts like that but you are right I should have added the rest of them. Thanks.

By the way, it worked what David suggested below.
written 4 months ago by Christian  
4 months ago by
Hello Jwinkle,

This is the code that gives error while projection :

from dolfin import *
from mshr import *
import numpy as np
import warnings

parameters["form_compiler"]["representation"] = "quadrature"  

dx = dx(metadata={'quadrature_degree': QUAD_DEG})

L = 2; W = 1
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 4, 2, 2)

V = VectorFunctionSpace(mesh,"Lagrange",2)
elementT = VectorElement("Quadrature", mesh.ufl_cell(), degree=QUAD_DEG, quad_scheme="default",dim=36)
VT = FunctionSpace(mesh,elementT)
elementS = VectorElement("Quadrature", mesh.ufl_cell(), degree=QUAD_DEG, quad_scheme="default",dim=6)
VS = FunctionSpace(mesh,elementS)
du, v = TrialFunction(V), TestFunction(V)
t, s = Function(VT), Function(VS)

def sigma_vec(s):
    return as_vector([s[0], s[1], s[2],s[3], s[4], s[5]])


Thanks !

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

Similar posts:
Search »