Extracting the values at quadrature points

General
1

Is there any way to extract/store the values of solution/a function of solution at quadrature points?
I want to print/plot the values of q(u) at quadrature points in the following code. How can I achieve that? I have projected q(u) to a quadrature function space. After that, I am not able to print/plot those values.

"""
FEniCS tutorial demo program:
Nonlinear Poisson equation with Dirichlet conditions
in x-direction and homogeneous Neumann (symmetry) conditions
in all other directions. The domain is the unit hypercube in
of a given dimension.

-div(q(u)*nabla_grad(u)) = 0,
u = 0 at x=0, u=1 at x=1, du/dn=0 at all other boundaries.
q(u) = (1+u)^m

Solution method: Newton iteration at the algebraic level.
"""

from dolfin import *
from fenics import *
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
import numpy as np, sys
import matplotlib.pyplot as plt

# Create mesh and define function space
degree = int(2)
divisions = 10

mesh = UnitSquareMesh(10, 10)
Vq = FiniteElement('Lagrange', mesh.ufl_cell(), degree, quad_scheme='default')
V = FunctionSpace(mesh, Vq)

# Quadrature function space
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)

fq = Function(W0)

metadata = {"quadrature_degree": degree, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

coordinates = mesh.coordinates()
print("node number", len(coordinates))
dofs = V.dofmap().dofs()
print("dof_u", dofs)

# Define boundary conditions for initial guess
tol = 1E-14
def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol

Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
bcs = [Gamma_0, Gamma_1]




def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dxm
    b_proj = inner(v, v_)*dxm
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return
  




# Define variational problem for initial guess (q(u)=1, i.e., m=0)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(nabla_grad(u), nabla_grad(v))*dx
f = Constant(0.0)
L = f*v*dx
A, b = assemble_system(a, L, bcs)
u_k = Function(V)
solve(A, u_k.vector(), b, 'lu')
    

# Note that all Dirichlet conditions must be zero for
# the correction function in a Newton-type method
Gamma_0_du = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1_du = DirichletBC(V, Constant(0.0), right_boundary)
bcs_du = [Gamma_0_du, Gamma_1_du]

# Choice of nonlinear coefficient
m = 2

def q(u):
    return (1+u)**m

def Dq(u):
    return m*(1+u)**(m-1)

# Define variational problem for the matrix and vector
# in a Newton iteration
du = TrialFunction(V) # u = u_k + omega*du
a = inner(q(u_k)*nabla_grad(du), nabla_grad(v))*dx + \
    inner(Dq(u_k)*du*nabla_grad(u_k), nabla_grad(v))*dx
L = -inner(q(u_k)*nabla_grad(u_k), nabla_grad(v))*dx

# Newton iteration at the algebraic level
du = Function(V)
u  = Function(V)  # u = u_k + omega*du
omega = 1.0       # relaxation parameter
eps = 1.0
tol = 1.0E-5
iter = 0
maxiter = 25
# u_k must have right boundary conditions
while eps > tol and iter < maxiter:
    iter += 1
    print (iter, 'iteration')
    A, b = assemble_system(a, L, bcs_du)
    print ("size of jacobian",(A.array()).shape)
    print("condition number", np.linalg.cond(A.array()))
    solve(A, du.vector(), b)
    eps = b.norm("l2")
    print ('Norm:', eps)
    u.vector()[:] = u_k.vector() + omega*du.vector()
    u_k.assign(u)


convergence = 'convergence after %d Newton iterations at the algebraic level' % iter
if iter >= maxiter:
    convergence = 'no ' + convergence

print ("""
Solution of the nonlinear Poisson problem div(q(u)*nabla_grad(u)) = f,
with f=0, q(u) = (1+u)^m, u=0 at x=0 and u=1 at x=1.
%s
%s
""" % (mesh, convergence))

plt.figure()
p = plot(u, title="Solution u")
plt.colorbar(p)

plt.show()


local_project(q(u), W0, fq)
2

If you have projected it into a quadrature space, you can get the points by call V.tabulate_dof_coordinates()
as done in the following post How to determine gauss point coordinates in a mesh - #4 by dokken
Here the ith row corresponds to the ith entry of fq.vector().get_local()

2 Likes