### Point-wise Conditional

232

views

0

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.

Regards,

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.

Regards,

Community: FEniCS Project

### 2 Answers

6

The

Solution is not simple nor straightforward. Perhaps the easiest way is to interpolate (or project) the function "T" to

One can also experiment with the special

`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.
0

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.

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

Regards,

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
lambda_=1.25
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))
Psi=(lambda_/2)*(tr(epsilon(u)))**2+mu*tr(epsilon(u)*epsilon(u))
Pi=Psi*dx-dot(f, u)*dx - dot(h, u)*ds
F=derivative(Pi,u,v)
solve(F==0,u,bc)
# Plot solution
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, Va)
T=conditional(le(u_magnitude,0.007),project((0.0),Va),project(u_magnitude/(2*mu+(2/3)),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)
143
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)
363
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)
2694
2695 """
-> 2696 _fem.SystemAssembler_swiginit(self, _fem.new_SystemAssembler(a, L, bcs))
2697
2698 def assemble(self, *args) -> "void":
RuntimeError:
*** -------------------------------------------------------------------------
*** 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 ?

Regards,

1

@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:

https://fenics.readthedocs.io/projects/ufl/en/latest/manual.html

See the ufl specification for details on keyword symbols:

https://fenics.readthedocs.io/projects/ufl/en/latest/manual.html

written
4 months ago by
jwinkle

1

Thanks for the advice ! I got the solution.

written
4 months ago by
Christian

1

It is not a problem with conditional. The problem is in projection. For a technical reason just replace

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

`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

2

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

I see now what you wanted to do. You can apply such conditional with the following piece of code

I iterate through each cell, get a dof index for that cell and ask for value of

`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
else:
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.