Setting minimal values to a function object

dolfinx
1

Hi all, is there a way to set fixed values to a function object? My hope is to monitor if the solution drops below a certain eps, say eps = 1e-16, and if the function value is smaller than eps at any point, to return the value eps for that point.

In this way, I try to get rid of negative and/or very low values of the function that often mess up a solution. Alternatives to this “marking” would also be quite welcome.

2

I assume you want this as some kind of post-processing step. Consider the following.

import dolfinx
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 8)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: np.sin(np.pi*x[0]))

print("Before manipulation")
print(u.x.array)

u.x.array[:] = np.where(u.x.array < 0.5, 0.0, u.x.array)
u.x.scatter_forward()

print("After manipulation")
print(u.x.array)

which yields

Before manipulation
[0.00000000e+00 3.82683432e-01 7.07106781e-01 9.23879533e-01
 1.00000000e+00 9.23879533e-01 7.07106781e-01 3.82683432e-01
 1.22464680e-16]
After manipulation
[0.         0.         0.70710678 0.92387953 1.         0.92387953
 0.70710678 0.         0.        ]
1 Like
3

@nate Thnak you for the suggestion. It works on its own, but I guess I am using another type of the Function object, as when I try to use the example in my code it throws the error that AttributeError: 'Function' object has no attribute 'x'

Instead, I found a workaround to do something like this:

from fenics import *

mesh = Mesh("mesh.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 2)
u = interpolate(Expression('sin(pi * x[0])', degree=2, domain=mesh), V)
u.assign(interpolate(Expression('f > err ? f : err', degree=2, f=u, err=0.5), V))

Is that a valid thing to do or is there a better way to tackle this?

Also, it is a somewhat intermittent step, as I am solving a diffusion equation so this is a check that I do on every step to ensure strictly positive numbers. Without it, there are quite a few numerical issues in the solution going haywire and this workaround at least seems to mitigate it a bit, but not perfectly.

4

You are using legacy Dolfin, while Nate suggests a solution for DOLFINx (as you didn’t specify the version you were using in your original post).

You can do

u.vector()[:] = np.where(u.vector()[:] < 0.5, 0.0, u.vector()[:])
u.vector().apply('insert')
2 Likes
5

@dokken Many thanks, that works!
It seems like it works even better than my “workaround”, as it gets rid of those odd numerical issues that I mentioned.