Point-wise Conditional

4 months ago by
Hello everyone,

I'm trying to use conditional for an function that contains information for each node. Let's say we have 50 nodes in our mesh and function "T" has scalar information for each node (So, we can say T as (50,1) array). By using the "T" function I want to use conditional. Is it possible to check each of the node value in that function by conditional ? I could not get any solution for it. It seems like FEniCS does not allow it.Do you have any suggestions ? Thanks.

Community: FEniCS Project

2 Answers

4 months ago by
The Conditional in FEniCS (more precisely in UFL and FFC) is treated as any other ufl expression in a form - i.e. evaluated at quadrature point in order to assemble a tensor. Since currently implemented quadrature rules produce points that are not at nodes you cannot simply implement a conditional checking values at nodes (=at degrees of freedom).

Solution is not simple nor straightforward. Perhaps the easiest way is to interpolate (or project) the function "T" to DG0 space. Then asking conditional(T > 1.0, ..., ...) has apparent meaning, since the value of that function at quadrature points coincide with value at nodes. But you are loosing a bit information if "T" was originally in richer space.

One can also experiment with the special Quadrature finite element family where degrees of freedom coincide with quadrature points. But here, the only advantage is that if you print vector of degrees of freedom, T.vector()[:] you can see the values that will be checked by conditional. However, "T" will be evaluated inside mesh cell's, not at vertices of mesh.
4 months ago by
Hello Michal,

Thank you for your answer. However I still could not use the conditional operator in a right way. I used also DG0 space but it did not work. I am sharing a sample code here, if you can look into that I will be appreciate.
from __future__ import print_function
from fenics import *
%matplotlib inline

# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
g = 0.4*delta**2
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 1,1,1)
V = VectorFunctionSpace(mesh, 'P', 1)
Va = FunctionSpace(mesh, 'DG', 0)

tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

# Define variational problem
u = Function(V)
d = 3  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
h = Constant((0, 0, 0))


Pi=Psi*dx-dot(f, u)*dx - dot(h, u)*ds



# Plot solution
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, Va)

the error that I have got :

RuntimeError                              Traceback (most recent call last)
<ipython-input-6-24a059050e86> in <module>()
     46 u_magnitude = sqrt(dot(u, u))
     47 u_magnitude = project(u_magnitude, Va)
---> 48 T=conditional(le(u_magnitude,0.007),project((0.0),Va),project(u_magnitude/(2*mu+(2/3)),Va))

/usr/lib/python3/dist-packages/dolfin/fem/projection.py in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
    140     # Assemble linear system
    141     A, b = assemble_system(a, L, bcs=bcs,
--> 142                            form_compiler_parameters=form_compiler_parameters)
    144     # Solve linear system for projection

/usr/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble_system(A_form, b_form, bcs, x0, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend)
    364     # Call C++ assemble function
--> 365     assembler = cpp.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
    366     assembler.add_values = add_values
    367     assembler.finalize_tensor = finalize_tensor

/usr/lib/python3/dist-packages/dolfin/cpp/fem.py in __init__(self, a, L, bcs)
   2695         """
-> 2696         _fem.SystemAssembler_swiginit(self, _fem.new_SystemAssembler(a, L, bcs))
   2698     def assemble(self, *args) -> "void":


*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***     fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to assemble system.
*** Reason:  expected a linear form for L.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0
*** DOLFIN version: 2017.1.0
*** Git changeset:  7215d3707edbf33a51e7b4bc3dfc8d67922651a2
*** -------------------------------------------------------------------------

Additionally, is there any example about how can I use "Quadrature" finite element family for a such application ?

@Christian:  there are a few symbols that are predefined in UFL specification.  Looking at the error, 'L' is typically a linear form, and you have defined it as a constant.  Try changing this to 'Lval' or something else.
See the ufl specification for details on keyword symbols:
written 4 months ago by jwinkle  
Thanks for the advice ! I got the solution.
written 4 months ago by Christian  
It is not a problem with conditional. The problem is in projection. For a technical reason just replace project((0.0), Va) with project(Constant(0.0), Va) inside the conditional.

The reason: FEniCS creates a simple bilinear and unilinear forms for `L^2` projection putting the argument of projection into the right-hand-side of the problem. See line https://bitbucket.org/fenics-project/dolfin/src/5c18533c286f81a177e6f2a3a27ca4ca57e07f8e/python/dolfin/fem/projection.py?at=master&fileviewer=file-view-default#projection.py-129
written 4 months ago by Michal Habera  
You are totally right, thanks for your response! I got what you mean from the link that you shared. Lastly, is there any method to print the vector comes from conditional by using DG0 space ? T.vector().array() does not work.
written 4 months ago by Christian  
Of course you can't access vector of values, because what conditional produces is just a symbolic UFL representation that can only be used in a finite element form. You can as well project the conditional again to Va space, but that would be a bit dirty...

I see now what you wanted to do. You can apply such conditional with the following piece of code
u_magnitude = project(sqrt(dot(u, u)), Va)
T = Function(Va)

for cell in cells(mesh):
    dofindex = Va.dofmap().cell_dofs(cell.index())
    if u_magnitude.vector()[dofindex] <= 0.007:
        T.vector()[dofindex] = 0.0
        T.vector()[dofindex] =\
            u_magnitude.vector()[dofindex] / (2 * mu + (2.0 / 3))​

I iterate through each cell, get a dof index for that cell and ask for value of u_magnitude in that cell. Depending on that the function T is prepared. Now you can save T to a file (xdmf, vtk) and visualize...
written 4 months ago by Michal Habera  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »