I have been trying to use FEniCS to solve moving mesh problems and have run into an error while trying to evaluate the function at mesh points after moving the mesh. As an example I have worked with the Poisson Equation provided in the documented demo. What I have tried to do is to solve the Poisson equation after moving the mesh using ALE.move command. Here is the code that I am working with
from dolfin import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x, on_boundary):
return on_boundary
u0 = Constant(1.0)
bc = DirichletBC(V, u0, boundary)
class f_exp(UserExpression):
def eval(self, values, x):
values[0] = 10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)
def value_shape(self):
return ()
class g_exp(UserExpression):
def eval(self, values, x):
values[0] = sin(5*x[0])
def value_shape(self):
return ()
u = TrialFunction(V)
v = TestFunction(V)
f = f_exp(degree = 2)
g = g_exp(degree = 2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx
u = Function(V)
solve(a == L, u, bc)
plt.figure()
plot(u)
*Now the mesh is modified by using the ALE.move command
boundary = BoundaryMesh(mesh, "exterior")
for x in boundary.coordinates():
x[0] *= 2
x[1] *= 2
ALE.move(mesh, boundary)
u = Function(V)
solve(a == L, u, bc)
plt.figure()
plot(u)
The function values can be calculated at the mesh points which is obtained from
mesh.coordinates()
array([[ 0. , 0. ],
[ 0.0625, 0. ],
[ 0.125 , 0. ],
...,
[ 1.875 , 2. ],
[ 1.9375, 2. ],
[ 2. , 2. ]])
u (2, 2)
gives the output as
1.0
Now the mesh is being modified and the problem is being solved on the mesh again
boundary = BoundaryMesh(mesh, "exterior")
for x in boundary.coordinates():
x[0] *= 2
x[1] *= 2
ALE.move(mesh, boundary)
u = Function(V)
solve(a == L, u, bc)
plt.figure()
plot(u)
The mesh points in the new mesh configuration are obtained using the following command
mesh.coordinates()
array([[ 0. , 0. ],
[ 0.125, 0. ],
[ 0.25 , 0. ],
...,
[ 3.75 , 4. ],
[ 3.875, 4. ],
[ 4. , 4. ]])
Then trying to calculate the Function at the mesh points results in the the following error
u (3.75 , 4 )
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-34-feee2d8d526e> in <module>()
----> 1 u (3.75 , 4. )
/usr/lib/python3/dist-packages/dolfin/function/function.py in __call__(self, *args, **kwargs)
335
336 # The actual evaluation
--> 337 self._cpp_object.eval(values, x)
338
339 # If scalar return statement, return scalar value.
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 evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
The point is clearly inside the domain but still it results in an error explaining that the point lies outside the domain. Can anyone kindly explain what I might have done wrong in the code that might be causing the error?
Thank You