How to integrate per cell?

variational formulation
1

I am writting the navier stokes solution for stationary flow with Galerkin Least Squares Method (GLS). The Least Square terms are weighted by the local reynolds number:

Re_L=\frac{\rho\cdot |\boldsymbol u|\cdot D}{\mu}

Where |\boldsymbol u| is the norm of the velocity, \rho the density, D the CellRadius*2 and \mu the viscosity.

So I need to integrate per cell over the domain

How to do it?

from dolfin import *
from mshr import *
import ufl as ufl
import matplotlib.pyplot as plt
import numpy as np
L=1
Ux=1

channel = Rectangle(Point(0, 0), Point(L, L))
mesh    = generate_mesh(channel , 100)

tol        = 1E-8
mu         = 1
rho        = 1
nu         = mu/rho

# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
V = FunctionSpace(mesh, P2)
Q = FunctionSpace(mesh, P1)
W  = FunctionSpace(mesh, TH)

# Define boundaries
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] > L-tol
class walls(SubDomain):
	def inside(self, x, on_boundary):
		return on_boundary and x[0] > L-tol or on_boundary and x[0] < tol or on_boundary and x[1] < tol
bcu_top    = DirichletBC(W.sub(0), Constant((Ux, 0)), top())
bcu_walls     = DirichletBC(W.sub(0), Constant((0, 0)), walls())

bcs = [bcu_walls, bcu_top]
# Define variational problem
v, q   = TestFunctions(W)

up   = Function(W)
u, p = split(up)
f      = Constant((0, 0)) # Possible forcing term
nu     = Constant(nu)

c = Circumradius(mesh)
c_vector=project(c,FunctionSpace(mesh, "DG", 0)).vector()[:] * 2
ReL=rho*c_vector/mu#*norm(v)

G      = inner(dot(grad(u), u), v)*dx+2*nu*inner(sym(grad(u)), sym(grad(v)))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
continuity=inner(div(u),div(v))*dx
movement=inner((-grad(p)+div(2*nu*sym(grad(u)))),-grad(q)+div(2*nu*sym(grad(v))))*dx
GLS=(continuity+movement)

c = Circumradius(mesh)
c_vector=project(c,FunctionSpace(mesh, "DG", 0)).vector()[:] * 2
ReL=rho*c_vector/mu#*norm(v) #( I still dont know how to extract the norm of the velocity and put into inside local reynolds)

# IF I UNCOMMENT HERE BELOW THE ReL I got the error
F=G+GLS#*ReL

solve(F == 0, up, bcs, solver_parameters = {"newton_solver":{ "linear_solver" : "mumps"}})

u, p = up.split(deepcopy=True)

c=plot(u.sub(0),levels=np.linspace(-0.3,1,14),cmap='rainbow')
plt.colorbar(c)
plt.show()
print(c_vector)
print(GLS)
print(ReL)
2

See for Element stiffness matrix - #2 by bhaveshshrimali

3

I tried, without even weighting :

for i, cell in enumerate(cells(mesh)):
	A = assemble_local(GLS, cell)

F=G+A

Not expecting keyword arguments when solving linear algebra problem.