I’m experimenting on stablization techniques of the advection-diffusion-reaction equation, on which I found some code by Evan Michael Cummings. He uses a Bubble Enriched method for stabilization:
from fenics import *
mesh = IntervalMesh(50 , 0, 1)
Q = FunctionSpace(mesh , 'CG', 1)
B = FunctionSpace(mesh , 'B', 2)
M = Q + B
def left(x, on_boundary ):
return on_boundary and x[0] == 0
uD = project ( Constant(0.0), M)
leftBC = DirichletBC(Q, 0.0, left )
leftBC_b = DirichletBC(Q, uD , left )
kappa = Constant(1.0/100.0)
s = Constant(5.0)
d = Constant(10.0)
f = Function(Q)
f.vector()[25] = 1000 # this is about the middle for a 50 element mesh
# bubble - enriched solution :
u = TrialFunction(M)
v = TestFunction(M)
us = Function(M)
uf2 = Function(Q)
a = + kappa * u.dx(0) * v.dx(0) * dx \
+ d * u.dx(0) * v * dx \
+ s * u * v * dx
L = f * v * dx
solve(a == L, us , leftBC_b )
uf2.interpolate (us)
However, this formulation doesn’t seem to work (anymore) in FEniCS, because I get the following error:
TypeError: unsupported operand type(s) for +: ‘FunctionSpace’ and ‘FunctionSpace’
Does anyone know how to adapt the formulation to obtain the Bubble Enriched functionspace? I found some info on NodalEnrichedElements here, but I can’t figure out how to implement this in this case because they do not apply it to a FunctionSpace in the example.
Thanks, this works indeed! My next step is to apply this to a mixed functionspace, for the concentrations of two different substances:
QE = FiniteElement('CG',mesh.ufl_cell(),1)
BE = FiniteElement('B',mesh.ufl_cell(),2)
QME = MixedElement([QE, QE])
BME = MixedElement([BE, BE])
M = FunctionSpace(mesh,QME+BME)
Q = FunctionSpace(mesh,QME)
However the form above doesn’t work, because in the defenition of M it gives the error: AttributeError: ‘MixedElement’ object has no attribute ‘degree’. Any idea how to solve this?
Try reversing the order of the steps: enrich the scalar elements first, then create the MixedElement:
QE = FiniteElement('CG',mesh.ufl_cell(),1)
BE = FiniteElement('B',mesh.ufl_cell(),2)
M = FunctionSpace(mesh,MixedElement(2*[QE+BE,]))
Q = FunctionSpace(mesh,MixedElement(2*[QE,]))
I was wondering if you could walk me through the case that instead of the scalar element (CG) we have vector element (NED or RT). MWE of the Unstable form:
mesh = UnitSquareMesh(10, 10)**
CGE = FiniteElement(“CG”, mesh.ufl_cell(), degreeCG)
CurlE = FiniteElement(“N2curl”, mesh.ufl_cell(), degreeCurl) # N1curl and N2curl
DivE = FiniteElement(“BDM”, mesh.ufl_cell(), degreeDiv) # RT and BDM
Z = FunctionSpace(mesh, MixedElement([CGE, CGE, CurlE, CurlE, DivE, DivE]))
But I couldn’t understand how define bubble function for vector elements.
The first few things I tried gave me errors, so I don’t have any immediate answer. However, if nothing else works, you could try to just include the enrichment as a separate component in a mixed space, then add it to the vector field in UFL, and use the resulting UFL Sum in variational forms. For example, the following runs
from dolfin import *
N = 10
mesh = UnitSquareMesh(N,N)
DivE = FiniteElement("BDM", mesh.ufl_cell(), degree=1)
BVE = VectorElement("B", mesh.ufl_cell(), degree=3)
V = FunctionSpace(mesh,MixedElement((DivE,BVE)))
# BDM and bubble space are separate components of mixed space:
w = Function(V)
# Use UFL to manually form enriched vector:
bdm_part = as_vector([w[0],w[1]])
bubble_part = as_vector([w[2],w[3]])
u_enriched = bdm_part + bubble_part
although I can’t vouch for whether this is a useful/stable element for any purpose.
I could not understand the last part (Use UFL to manually form enriched vector).
could you please let me know what would be this u_enriched?
And why we are not considering our MixedElement like a normal two fields space?
The goal of the original post was to obtain a bubble-enriched element for a single vector field, rather than two different fields. However, the desired base element type did not support bubble enrichment. This is just a mathematically-equivalent workaround.
Thanks for clarification.
I understood fenics does not support bubble enrichment for two or more vector fields. right?
So If this is the case, could you please give me a hint on how can I implement that manually?
Following my example above, the idea is just to take your original variational form, but replace the solution field with u_enriched, and also replace the test function with a similarly-defined v_enriched, which is constructed from a TestFunction in the mixed space, in an analogous way to how u_enriched is constructed from the Function (or, in a linear problem, TrialFunction) w.
The issue here is that the element enrichment is canonically done on the reference element, when defining the shape functions. However, the vector-valued shape functions for RT-type elements are pushed forward to the physical element using the Piola transform (i.e., the transformation law for vector proxies of 2-forms), whereas the bubble element uses the same scalar shape functions for every component, each of which is pushed forward by simple composition with the mapping (i.e., each component is treated like a 0-form). Thus, the two elements are in some sense “incompatible”, in that the enrichment cannot be done on the reference element, then pushed forward. The workaround I proposed above essentially pushes each set of basis functions forward, then does the enrichment in physical space. As I mentioned, this is not obviously a good idea to me. (If what you really want is to enrich, then push forward, you could also try to apply a Piola transform yourself, on bubble_part, using ufl.Jacobian(mesh), but, for non-degenerated constant-Jacobian elements, I believe the transformed bubble_part will just span the same space, since the components will just be scaled and rotated by element-wise constant amounts.)
It seems that maybe you’re unclear on my workaround, so I put together a simple example, using a scalar problem where the enrichment can be done in both ways:
from dolfin import *
# Set up mesh and elements:
mesh = UnitIntervalMesh(3)
CGE = FiniteElement("CG",mesh.ufl_cell(),1)
BE = FiniteElement("Bubble",mesh.ufl_cell(),2)
# Simple variational problem, for demonstration purposes:
def a(u,v):
return inner(grad(u),grad(v))*dx + u*v*dx
def L(v):
x = SpatialCoordinate(mesh)
return sin(x[0])*v*dx
# Solve doing the enrichment "manually":
V_mixed = FunctionSpace(mesh,MixedElement([CGE,BE]))
w = TrialFunction(V_mixed)
dw = TestFunction(V_mixed)
u_enriched = w[0] + w[1]
v_enriched = dw[0] + dw[1]
wh = Function(V_mixed)
solve(a(u_enriched,v_enriched)==L(v_enriched),wh)
uh_enriched = wh[0] + wh[1]
# Solve doing the enrichment the standard way:
V = FunctionSpace(mesh,CGE+BE)
uh = Function(V)
v = TestFunction(V)
solve(a(TrialFunction(V),v)==L(v),uh)
# Solutions are exactly identical:
print(assemble((uh-uh_enriched)**2*dx))
Sorry for bumping this thread, I have a question about this method. While searching for stabilization methods for drift-diffusion equation that can be easily implemented in FEniCS, I found this one as a possible candidate. I tried to adapt the code for 2D case, but it crashes with an error (MWE and error are given below). So, is it possible to use Bubble enriched stabilization in 2D and how to implement it?
MWE is given below:
from dolfin import *
# mesh = UnitIntervalMesh(3) #It works in 1D
mesh = RectangleMesh(Point(0.0, 0.0),Point(1.25e-2, 1.35e-2), 15, 15, "crossed") #It crashes when using 2D mesh
CGE = FiniteElement("CG",mesh.ufl_cell(),1)
BE = FiniteElement("Bubble",mesh.ufl_cell(),2)
def a(u,v):
return inner(grad(u),grad(v))*dx + u*v*dx
def L(v):
x = SpatialCoordinate(mesh)
return sin(x[0])*v*dx
V = FunctionSpace(mesh,CGE+BE)
uh = Function(V)
v = TestFunction(V)
solve(a(TrialFunction(V),v)==L(v),uh)
and the error is:
RuntimeError: Bubble element of degree 2 and codimension 0 has no dofs