How are local and global dofs linked in Mixed spaces?

jt509 · November 22, 2022, 5:22pm

FEniCS version: 2019.1.0 installed via Conda

Consider the following MWE on mixed spaces:

from dolfin import *

## Spaces
mesh = UnitSquareMesh(2, 2)
P1 = FiniteElement('CG', triangle, 1)
P0 = FiniteElement('DG', triangle, 0)
elem = MixedElement([P1, P0, P0])
W  = FunctionSpace(mesh, elem)

## Trial and Test functions
w1, w2, w3 = TrialFunctions(W)
t1, t2, t3 = TestFunction(W)

## global dofs
dofs0 = W.sub(0).dofmap().dofs()
dofs1 = W.sub(1).dofmap().dofs()
dofs2 = W.sub(2).dofmap().dofs()
print(dofs0, dofs1, dofs2)

## Interpolate
f1 = interpolate(Expression('-x[1]', degree = 1), W.sub(1).collapse())
assigner0 = FunctionAssigner(W.sub(1), f1.function_space())
f1V = Function(W)
assigner0.assign(f1V.sub(1), f1)

f3 = interpolate(Constant(1), W.sub(1).collapse())
assigner2 = FunctionAssigner(W.sub(1), f3.function_space())
z0 = Function(W)
assigner2.assign(z0.sub(1), f3)

f4 = interpolate(Constant(-1.7), W.sub(2).collapse())
assigner3 = FunctionAssigner(W.sub(2), f4.function_space())
alp0 = Function(W)
assigner3.assign(alp0.sub(2), f4)

## Working on local spaces
z = z0.vector()[dofs1]
alp = alp0.vector()[dofs2]
f1VV = f1V.vector()[dofs1]
# Condition check
A = alp + (z-f1VV) 
marked_cellsA = []
CellMarker = MeshFunction('size_t', mesh, mesh.topology().dim())
CellMarker.set_all(0)     # 0 -> all elements
cell_indices = CellMarker.where_equal(0)
Wloc = W.sub(1).collapse().dofmap()
for cell in cell_indices:
    cell_dofs = Wloc.cell_dofs(cell)
    values = A[cell_dofs]
    for value in values:
        if value >= 0:
            marked_cellsA.append(cell)
    CellMarker.array()[marked_cellsA] = 1            
print(marked_cellsA)     

## Mark the cells
dAct = Measure('dx', subdomain_data=CellMarker)

## Integrations
Eq1 = t2*w2*dx
Eq2 = t2*w2*dAct(1)

Output

[0, 1, 2, 5, 6, 11, 14, 17, 22] [3, 7, 9, 12, 15, 18, 20, 23] [4, 8, 10, 13, 16, 19, 21, 24]
[0, 3]

Questions:

What does dx represent in the case of mixed spaces? Is it still integration over elements, if yes then how is it related to the global dofs?
How are the marke_cellsA local dofs belonging to P0, i.e., [{\color{blue}0}\; {\color{green}3}] linked with the global dofs [{\color{blue}3}, 7, 9, {\color{green}12}, 15, 18, 20, 23] and [{\color{blue}4}, 8, 10, {\color{green}13}, 16, 19, 21, 24], respectively (Is there a 1-1 mapping between them represented by blue and green colours).
What should one expect from the internal Eq2 = t2w2dAct(1)?

dokken · November 22, 2022, 6:52pm

dx is not related to a function space, it is related to the mesh. It means that you will integrate over every cell in your mesh, per term in your variational forms.

I am not really sure what you are trying to do here.
You are trying to find all cells that are satisfying some condition. over every cell I guess.
You could use a ufl conditional for this, i.e.

condition = ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0)
Eq2 = t2*w2*condition*dx

As far as I can tell you would only integrate over cell 0 and 3.
You could check the integration area with the following:

# Mark the cells
dAct = Measure('dx', subdomain_data=CellMarker, domain=mesh)

# Integrations
condition = ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0)

print(assemble(1*dAct(1)),
      assemble(condition*dx), assemble(1*dx(domain=mesh)))

jt509 · November 23, 2022, 4:49am

Thanks for your suggestions @dokken. I was trying to implement the ufl.conditional in my MWE, so I added the following lines of code to my MWE

from ufl import *

##  towards the end
condition = ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0)
Eq2 = t2*w2*condition*dx

But it gives the following error

Exception has occurred: AttributeError
'FunctionSpace' object has no attribute 'sub'
     File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2glb.py", line 17, in <module>
         dofs0 = W.sub(0).dofmap().dofs()

and if I remove the from\;\;ufl\;\;import\;\;* then the conditional statement gives the following error

Exception has occurred: NameError
name 'ufl' is not defined
  File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2glb.py", line 64, in <module>
    condition = ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0)

jt509 · November 23, 2022, 5:06am

Regarding:

print(assemble(1*dAct(1)), assemble(1*dx(domain=mesh)))

when I print:

print(assemble(1*dx(domain=mesh)))

I get the output as

1.0

Shouldn’t it be a column vector with all entries equal to the area of the element?

Moreover when I print

print(assemble(1*dAct(1)))

I get the following error

Exception has occurred: UFLException
This integral is missing an integration domain.
  File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2glb.py", line 65, in <module>
    print(assemble(1*dAct(1)))

dokken · November 23, 2022, 5:41am

jt509:

from ufl import *

##  towards the end
condition = ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0)
Eq2 = t2*w2*condition*dx

But it gives the following error

Exception has occurred: AttributeError
'FunctionSpace' object has no attribute 'sub'
     File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2glb.py", line 17, in <module>
         dofs0 = W.sub(0).dofmap().dofs()

and if I remove the from;;ufl;;import;;* fromuflimport∗from;;ufl;;import;;* then the conditional statement gives the following error

Exception has occurred: NameError
name 'ufl' is not defined
  File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2gl

Use import ufl or from ufl import conditional. You should never use wildcard import (from xyz import *) on more than one package.

No, then you would have to multiply by a TestFunction.

jt509:

print(assemble(1*dAct(1)))

I get the following error

Exception has occurred: UFLException
This integral is missing an integration domain.
  File "/Users/jaitushar/Documents/myFEniCS/Misc/loc2glb.py", line 65, in <module>
    print(assemble(1*dAct(1)))

Look at my redefinition of dAct

Which includes domain=mesh

jt509 · November 23, 2022, 12:26pm

Hi @dokken, thank you for your help. The ufl.conditional(ufl.gt()) is promising.

However, I do not understand the syntax used inside ufl.gt(). I am guessing that alp0[2] + z0[1] - f1V[1] is telling to take alp0, z0 and f1V from their respective subspaces. However, what does 0 in ufl.gt(alp0[2] + z0[1] - f1V[1], 0) and 1, 0 in ufl.conditional(ufl.gt(alp0[2] + z0[1] - f1V[1], 0), 1, 0) mean?

dokken · November 23, 2022, 12:29pm

Please read the documentation: https://fenics.readthedocs.io/projects/ufl/en/latest/manual/

jt509 · December 10, 2022, 11:37am

This question is in continuation of the previous thread @dokken. I tried to extend this to vector spaces, but the results after using the ufl.conditional do not make sense. Please see the following MWE in this regard:

from dolfin import *
import ufl

# Spaces
mesh = UnitSquareMesh(1, 1)
P1 = FiniteElement('CG', triangle, 1)
P0 = VectorElement('DG', triangle, 0)
elem = MixedElement([P1, P0, P0])
W  = FunctionSpace(mesh, elem)

VR = VectorFunctionSpace(mesh, "R", 0)
v = TestFunction(VR)

# Trial and Test functions
w1, w2, w3 = TrialFunctions(W)
t1, t2, t3 = TestFunctions(W)

# print global dofs
dofs0 = W.sub(0).dofmap().dofs()
dofs1 = W.sub(1).dofmap().dofs()
dofs2 = W.sub(2).dofmap().dofs()
print(dofs0, dofs1, dofs2)

# Interpolate
f1 = interpolate(Expression(('-x[1]','-x[1]'), degree = 1), W.sub(1).collapse())
assigner0 = FunctionAssigner(W.sub(1), f1.function_space())
f1V = Function(W)
assigner0.assign(f1V.sub(1), f1)

f3 = interpolate(Constant((1,1)), W.sub(1).collapse())
assigner2 = FunctionAssigner(W.sub(1), f3.function_space())
z0 = Function(W)
assigner2.assign(z0.sub(1), f3)

f4 = interpolate(Constant((-1.,-1.7)), W.sub(2).collapse())
assigner3 = FunctionAssigner(W.sub(2), f4.function_space())
alp0 = Function(W)
assigner3.assign(alp0.sub(2), f4)

lam = 1

## Working on local spaces
z = z0.vector()[dofs1]
alp = alp0.vector()[dofs2]
f1VV = f1V.vector()[dofs1]

A = alp + lam*(z-f1VV); 
print(A)

conditionA = ufl.conditional(ufl.le(alp0[2] + lam*(z0[1]- f1V[1]), 0), 1, 0)

e0 = assemble(inner(v,w2)*dx)
print(e0.array())

e2 = assemble(inner(v,w2)*conditionA*dx)
print(e2.array())

Output:

[0, 1, 2, 7] [3, 4, 8, 9] [5, 6, 10, 11]
[ 0.66666667 -0.03333333  0.33333333 -0.36666667]
[[0.  0.  0.  0.5 0.  0.  0.  0.  0.5 0.  0.  0. ]
 [0.  0.  0.  0.  0.5 0.  0.  0.  0.  0.5 0.  0. ]]
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

The positive values 0.66667 and 0.333333 correspond to the global Dofs 3 and 8, respectively. Then why is e2 = assemble(inner(v,w2)*conditionA*dx) an array of zeros? Shouldn’t it be zero just on the degrees of freedom where the conditionA is not met?

The documentation for the conditional,

@ufl_type(binop="__le__")
class LE(BinaryCondition):
    __slots__ = ()

    def __init__(self, left, right):
        BinaryCondition.__init__(self, "<=", left, right)

    def evaluate(self, x, mapping, component, index_values):
        a = self.ufl_operands[0].evaluate(x, mapping, component, index_values)
        b = self.ufl_operands[1].evaluate(x, mapping, component, index_values)
        return bool(a <= b)

will this perform properly with the vector input?

Thanks!

dokken · December 12, 2022, 8:26am

As you can see, the dimension of alp0 is 5

In [13]: len(alp0)                                                                                                                                                                                                                                                                                            
Out[13]: 5

where the first entry works on the scalar space, the index [1,2] on the second space, [3,4] on the third space.

Im not sure what you expect this value to be, if you expected the conditional being on a vector expression

as with the current variational form, conditionA has to evaluate to a scalar value, while your conditional (as far as I understood it) is a vector expression.

jt509 · December 12, 2022, 8:54am

I was expecting that the conditional is something that when compiled with FFC will act in a form as lamba x: 1 if conditionA holds else 0 (which is the case when my P0 space is scalar and not a vector.)

while your conditional (as far as I understood it) is a vector expression – Yes, I was hoping that the conditional would evaluate at individual dofs. More specifically, for the output

[0, 1, 2, 7] [3, 4, 8, 9] [5, 6, 10, 11]

[ 0.66666667 -0.03333333  0.33333333 -0.36666667]

[[0.  0.  0.  0.5 0.  0.  0.  0.  0.5 0.  0.  0. ]
 [0.  0.  0.  0.  0.5 0.  0.  0.  0.  0.5 0.  0. ]]

dofs 3 and 8 correspond to the x-component of the vector (both 0.5), while dofs 4 and 9 correspond to the y-component of the vector. And observing that the conditionA is only satisfied for the dofs 4 (which is -0.0333333) and 9 (which is -0.3666667) and applying conditionA, I was trying to get,

[[0.  0.  0.  0.  0.   0.  0.  0.  0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.5 0.  0.  0.  0.  0.5 0.  0. ]]

In other words, my conditionA would work like an indicator function. Can this be done using ufl.condtional for the vector space case?

dokken · December 12, 2022, 10:04am

Multiple things to note regarding conditionals:

conditionals are evaluated at quadrature points, not at the location of the dof.
when you call something like inner(u,v)*conditonalSomething*dx this is expanded as \sum_{q=0}^{N_q} w_q\left(\sum_{i=0}^1 u[i](x_q)v[i](x_q) conditionalSomething(x_q)\right)
where x_q is the qth global quadrature point. As you can see it does not make sense variationally for the conditional to be vector valued here.

jt509 · December 12, 2022, 10:25am

Understood. Thank you @Dokken. A workaround then would to evaluate the conditional component wise. But still keeping it as an object of the W space so that it is scalar valued but still has a compatible shape for inner(u,v)*conditionalSomething*dx.

Can the ufl.conditional statement conditionA = ufl.conditional(ufl.le(alp0[2] + lam*(z0[1]- f1V[1]), 0), 1, 0) be modified component wise?

Topic		Replies	Views
Weird behaviour in mixed spaces with vector elements vs scalar elements General	0	256	December 31, 2022
Understanding the mapping of local and global degrees of freedom in Mixed spaces General	0	373	November 13, 2022
Aligning dofs between subspace of mixed element and new function General	6	331	September 13, 2023
Cell markers in a Mixed Space dolfinx	10	619	March 8, 2023
MixedElement formulation leads to unintuitive behaviour of the API dolfinx	3	769	July 3, 2022

How are local and global dofs linked in Mixed spaces?

Related topics