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)

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

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!

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