How to implement a "project" function correctly?

1

In the following MWE, I am solving the Poisson equation (with the right-hand side multiplied by T) on a disc of radius Sqrt(3). T is a scalar that ranges from 0 to 5. Then in order to apply the solution to every mesh coordinate, I applied the following commands:

Y = mesh.coordinates()
X = project(Y,W1)
X += np.vstack(map(u,X))

MWE:

from fenics import *
from mshr import Circle, generate_mesh
import numpy as np

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary        


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

C = Circle(Point(0, 0), sqrt(3))

mesh = generate_mesh(C, 100)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(1*dx(1)), assemble(1*ds(2)))

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

n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D = Constant(0.0)
d = mesh.geometry().dim()
I = Identity(d)            
T = Constant(0.0)
c = Constant((1.0,0.0))
bc = DirichletBC(W1, u_D, boundary_markers, 2)
print(assemble(1*ds(2)))
for k in range(0,5,1):
# Define variational problem
  u = TrialFunction(W1)
  v = TestFunction(W1)
  f = Constant(-G/mu) # f=-G/mu
  a = dot(grad(u), grad(v))*dx
  L = T*f*v*dx 
  T = T + k
# Compute solution
  u = Function(W1)
  solve(a == L, u, bc)
  Y = mesh.coordinates()
  X = project(Y,W1)
  X += np.vstack(map(u,X))

Error:

 X = project(Y,W1)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
    L = ufl.inner(w, v) * dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 155, in inner
    b = as_ufl(b)
  File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
    raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: [[-1.73197184  0.01653895]
 [-1.73197184 -0.01653895]
 [-1.73134016 -0.04961083]
 ...
 [ 0.97718622  1.04771416]
 [ 0.09944235  0.30137168]
 [ 0.06897023  0.31762106]] can not be converted to any UFL type.

What’s wrong here with my implementation of the “project” function or is there another problem with my code? Any help or suggestion will be helpful.

2

The project function takes in an Expression, Function, or ufl.Form g, and solves the problem:
\int_{\Omega} u \cdot v \mathrm{d}x = \int_{\Omega} g \cdot v \mathrm{d}x where u is the TrialFunction, v the TestFunction.

What do you mean when you say you want to apply the solution to every mesh coordinate? The field of your solution u is a scalar CG 1 field, while the mesh coordinates is a vector field.

3

Sorry for the late response.

Yes, you are right. I made a mistake. There should be a Vector function space instead of a function space for the solution u. I have made the changes to my MWE

By “apply the solution to every mesh coordinate”, I meant to say that if u=(u_1, u_2) is the solution and suppose (x_1,y_1) is one my mesh coordinate then the modified mesh coordinate will be X +=(x_1 + u_1,x_2 + u_2).

Here is the corrected MWE:

from fenics import *
from mshr import Circle, generate_mesh
import numpy as np

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary        


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

C = Circle(Point(0, 0), sqrt(3))

mesh = generate_mesh(C, 100)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(1*dx(1)), assemble(1*ds(2)))

W1 = VectorFunctionSpace(mesh, "Lagrange", 1)

n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D = Constant((0.0,0.0))
d = mesh.geometry().dim()
I = Identity(d)            
T = Constant(0.0)
bc = DirichletBC(W1, u_D, boundary_markers, 2)

print(assemble(1*ds(2)))

for k in range(0,5,1):

# Define variational problem
  u = TrialFunction(W1)
  v = TestFunction(W1)
  f = Constant(-G/mu) # f=-G/mu
  a = inner(grad(u), grad(v))*dx
  L = T*f*inner(Constant((1.0,1.0)),v)*dx 
  T = T + k

# Compute solution
  u = Function(W1)
  solve(a == L, u, bc)
  Y = mesh.coordinates()
  Z = project(Y, W1)
  Y += np.vstack(map(u,Z))

Still the error reads:

Z = project(Y, W1)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
    L = ufl.inner(w, v) * dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 155, in inner
    b = as_ufl(b)
  File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
    raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: [[-1.73197184  0.01653895]
 [-1.73197184 -0.01653895]
 [-1.73134016 -0.04961083]
 ...
 [ 0.97718622  1.04771416]
 [ 0.09944235  0.30137168]
 [ 0.06897023  0.31762106]] can not be converted to any UFL type.
4

Use ALE.move(mesh, u) as explained in How to translate a 2D mesh by a vector? - #6 by dokken