Output of conditional statement

I/O
1

Hello all,

I wanted to know how I can get the output of a conditional UFL as an integer or float in an example like below:
Check = conditional( lt(dot(d,d), 0.01), 0, 1)

In the above script, d is a function of VectorFunctionSpace in a 2D mesh and has different values at different points. For getting its integral over the domain, I know that I can use assemble(Check * dx) . But I couldn’t find how I can find its output at each point of space.

Thank you in advance for your help!

2

I believe you can project your functional into a proper function space. Try something like:

V = FunctionSpace(mesh,'CG',1) # or whatever function space
f = project(Check, V)
File('f.pvd') << f
3

Thank you so much for your very helpful response. I am trying to run a code like this:

import numpy as np
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 20, 20)
V = VectorFunctionSpace(mesh, "CG", 1)

left =  CompiledSubDomain("near(x[0], side)", side = -1.0)
y_BC = Expression(("1.5*x[0]+ x[1]", "0*x[0]+1*x[1]"), degree=1)
bcs = [DirichletBC(V, y_BC, left)]

# Define functions
dy = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
y  = Function(V)                 # Displacement from previous iteration
d = Function(V)
d_initial= Expression(("x[0]","x[1]"), degree=1)
d=project(d_initial,V)
y_reference = Expression(("x[0]","x[1]"),degree=1)
y.assign(project(y_BC, V)) 
V1 = FunctionSpace(mesh,'CG',1)

def W1(y,d): 
    F = grad(y)             
    C = F.T*F
    Check = conditional(lt(dot(d,d), 0.4), 0, 1)
    f = project(Check, V1)
    print(f.vector().get_local())
    if f == 0: #Evaluate Check the value at each point
        return  0.5*tr(F.T*F)      
    else:  
        return tr(F.T*F)
    
# Total potential energy
Pi = (dot(d,d))*W1(y, d)*dx 
residual = derivative(Pi, y, v) # Compute first variation of Pi
jacobian = derivative(residual, y, dy) # Compute Jacobian of residual
solve(residual == 0, y, bcs, J=jacobian) # Solve variational problem

And now my problem is that when I print f at each point using print(f.vector().get_local()) , I don’t get 0 -1 values, which I was expecting as outputs of Check.

4

I think the issue here is that you are trying to project a discontinuous functional into a continuous function space (see this discussion for instance [my bad at choosing a CG space of order 1]).
Depending on what you want, you can define V1 as a Discontinuous space and then project or try the following as proposed in the above link I shared:

v = TestFunction(V1)
lhs = assemble(inner(Constant(1.0),v)*dx)
rhs = assemble(inner(Check,v)*dx)
f2 = Function(V1)
as_backend_type(f2.vector())\
    .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                           as_backend_type(lhs).vec())
print(f2.vector().get_local())
5

Thank you so much for your great help!