Define the function possitive part

General
1

Hello there,

I am new to fenics and python, moving from Matlab/Octave.

I want to define a function which gives the positive part of the trial function. I have tried this

def g(u):
if u>0:
pnu=u
else
pnu=0
return pnu

And then use it in

w=TrialFunction(V)
v=TestFunction(V)
J=g(w) * v* dx

But it does not work. I have also tried to use

pnu=conditional(gt(u,0),u,0)

instead of the if/else block without succes and wonder if I will have to get the nodal values using

nodal_values_u = u.vector()

and then compare componentwise, and then reconstruct somehow, but probably there is a better way.

Some help would be greatly appreciated.

Cheers.

2

What you’re proposing is not a linear problem. Taking the positive part of a ufl.TrialFunction does not make sense in variational formulations.

However, your problem is well suited for solution by iterative methods. For example Newton’s iterative method for which there are demos.

Based on your syntax I assume you’re using legacy dolfin: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/nonlinear-poisson/demo_nonlinear-poisson.py.html.

3

Thank you for your quick response.

I am aware this is a nonlinearity. In fact, I was following the tutorial “Solving PDEs in Python – The FEniCS Tutorial Volume I”, by HP Langtangen and A. Logg, particularly the section 3.2 (page 46) where they do implement a Nonlinear Poisson equation by using a variational formulation. Their code, which can be downladed from

But their code won’t run and issue the error

AttributeError: ‘dolfin.cpp.la.PETScVector’ object has no attribute ‘array’.

I wonder if this is because that code is fenics legacy and my installation is not.

Thank you in advance.

4

The code from the original fenics tutorial is outdated, and .array() as to be replaced with .get_local().
Updated demos for legacy FEniCS can be found at: Fix almost every outdated tutorial by JSS95 ¡ Pull Request #74 ¡ hplgit/fenics-tutorial ¡ GitHub
However, I would strongly suggest to migrate to DOLFINx, as DOLFIN is no longer maintained.

1 Like
5

Thank you for your quick response and useful link. It does compile now.

I still do not know how to implement my positive part function. I am aware I have a lot to learn.
Can you reference a tutorial of the likes of that of Langtangen and Logg but updated to dolfinx?

Thanks a lot.

6

https://jsdokken.com/dolfinx-tutorial/
http://jsdokken.com/FEniCS23-tutorial/README.html
https://docs.fenicsproject.org/dolfinx/v0.7.2/python/demos.html
to mention a few.

You can use a ufl.conditional, ufl.conditional(ufl.gt(u, 0), u, 0). See for instance: Using Conditional function inside a Project function - #4 by dokken

1 Like
7

Thank you for this useful links and remarks.

8

I am answering myself now. It worked by defining the function positive part

def qabs(z):
return (z+ abs(z))*0.5

Then using it in the variational formulation for a Function, not TrialFunction

u = Function(V) # Note: not TrialFunction!
v = TestFunction(V)
F = dot((grad(u)), (grad(v)))dx+kappa * dot(qabs(u),v)dx - fvdx

Here kappa is a penalty constant. And then solving by using

solve(F == 0, u, bc)