Failed to interpolate Gaussian distribution

Hi all,

I have a question regarding interpolation. I wanted to interpolate a Gaussian distribution for a boundary condition. However, I could not obtain a correct interpolation of the expression defining the distribution. I came across a similar issue when interpolating a function defined by exponential function. Do you have any ideas on what is going on here? Note that the version of Dolphin is 2019.2.0.

import dolfin as dl
import numpy as np
import matplotlib.pyplot as plt
N = 100 # number of elements
T = 10
mesh = dl.IntervalMesh(N, 0, T)
degree = 1
V = dl.FunctionSpace(mesh, “Lagrange”, degree)
mu = 5
sigma = 1.0
u_ini = dl.Expression(‘1/(sigma* sqrt(2* pi))* exp(-1/2* (pow((x[0] - mu)/sigma, 2)))’, mu = mu, sigma = sigma, degree = 1)
u_n = dl.interpolate(u_ini, V)

time = np.linspace(0, T, 101)

u_true = 1/(sigma* np.sqrt(2* 3.14))* np.exp(-1/2* ((time - mu)/sigma)**2)

dl.plot(u_n)
plt.plot(time, u_true)
plt.ylim(-0.1, 1.1)

I appreciate your help.

The problem is in your definition of u_ini. You should be careful of order of evaluation. Whichever of * or / lies on the left of the other gets evaluated first, and always before + and -. For this simple case, you should be able to see by evaluating the below expression and check that

(pow((x[0] - mu)/sigma, 2.0))*1/2

evaluates to your desired result, whereas

1/2*(pow((x[0] - mu)/sigma, 2.0))

evaluates to zero (0). And exp(0) always evaluates to 1 no matter what your x[0] is. You can explicitly specify floating point division as a default in your expression (1./2 vs 1/2) by

dl.Expression('1/(sigma*sqrt(2*pi))*exp(-1./2*(pow((x[0] - mu)/sigma, 2.0)))', mu = mu, sigma = sigma, degree = 1)

Also, you can enclose your code in triple backticks ``` to ensure proper readability and execution.

3 Likes

Hi Bhavesh,

Thank you for your help! I did not notice the problem is in the floating division since I am not familiar with C++ syntax (so used to Python easy one…). Thank you so much!

1 Like

thanks for the awesome information.

thanks my issue has been fixed.

Hello @dokken

How can I make a graph of the distribution calculated through the expression?

Project the expression into an appropriate function space, and then plot it with either matplotlib or paraview

This is my code. My mesh is in 2D vertical.

from dolfin import *
import numpy as np

L = 100
mesh = UnitSquareMesh(L, L)
    
V = FunctionSpace(mesh, "Lagrange", 1)
    
L = 11000
mu_phi = 0.5
x0_phi = 0.5 * L
phi0 = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi)
phia = interpolate(phi0, V)

file = File("phi/phi0.pvd")
file << phia

I am not sure what you expect. Just plotting this with wolfram alpha gives you 0 everywhere. Due to your large L (are you sure L shouldn’t be 1?): plot e^(-0.5*(x-5500)**2) for x in 0, 1 - Wolfram|Alpha

Sorry, I had a 11000 long mesh but for the example I used a fenics mesh. I just want to know if I did the distribution right by looking at the graph. But I can’t do it with the fenics plot nor with the paraview.

from dolfin import *
import numpy as np

mesh = RectangleMesh(Point(0.0, 0.0), Point(100.0, 50.0), 10, 10)

V = FunctionSpace(mesh, "Lagrange", 1)
mu_phi = 0.5
L = 100

V = FunctionSpace(mesh, "Lagrange", 1)

mu_phi = 0.5
x0_phi = 0.5 * L
phi0 = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi)
phia = interpolate(phi0, V)

file = File("phi/phi0.pvd")
file << phia

Hello @dokken

I have an issue with the implementation of a Gaussian bell. My goal is to create degradation with respect to the z-axis (my domain is x, z). After a certain height, I want to assign a specific value. The first ‘if’ statement works, but the ‘else’ part does not give me the desired value. Thank you for the help.

from dolfin import *
import matplotlib.pyplot as plt
Nx = Ny = 20
mesh = UnitSquareMesh(Nx, Ny)
Q = FunctionSpace(mesh, "CG", 1)
mu = 0.5
H = 0.5
A = 0.4
sigma = 0.1
class porosity(UserExpression):
    def __init__(self, mu, sigma, H, A, **kwargs):
        super().__init__(**kwargs)
        self.mu = mu
        self.sigma = sigma
        self.H = H
        self.A = A

    def eval(self, values, x):
        phi0 = 1 - A * exp(-1./2 * (pow((x[0] - mu)/sigma, 2)))
        #phi0 = 0.6
        if x[1] <= 0.5:
            values[0] = (1.0 - phi0) * sqrt(x[1] / H) + phi0
        else:
            values[0] = 1.0

    def value_shape(self):
        return ()
phi = porosity(mu=mu, sigma=sigma, H=H, A=A)
cp = plot(interpolate(phi, Q))
plt.colorbar(cp)
plt.show()

update

I visualized it in ParaView, and it does show what I expected. I would like to know if pyplot can help me achieve that visualization similar to ParaView. Also, I would like to know how I can view the value of a constant that depends on phi (from the previous code) like this.

K = pow(10, 2) * pow(phi, 3) / (150 * pow(1 - phi, 2) + 1e-2)

In legacy dolfin this is tricky. It is easier in dolfinx with dolfinx.fem.Expression

Where can I find information on how to do the installation and any tutorials?

Installation notes for dolfinx can be found at

with demos at
https://docs.fenicsproject.org/dolfinx/main/python/demos.html
and corresponding paper

Ive written several tutorials at
https://jsdokken.com/FEniCS23-tutorial/README.html
and
https://jsdokken.com/dolfinx-tutorial/

@dokken

Hello, I’m trying to install Fenicsx, and I’m encountering these issues. Can you help me?

mkdir build
cd build
cmake ..
make install
mkdir: cannot create directory ‘build’: File exists
CMake Error: The source directory "/home/gerard" does not appear to contain CMakeLists.txt.
Specify --help for usage, or press the help button on the CMake GUI.
make: *** No rule to make target 'install'.  Stop.

Where is the build directory compared to the dolfinx source dir?

Note that there are many other ways to install dolfinx than from source, for instance:

  • conda
  • docker
  • spack
  • apt-get