How to move a solution by a translation?

1

I have an annulus with inner radius 1 and outer radius 2. I am applying a traction of 5 units on the inner boundary. Since there are no essential boundary conditions, there is a translation of the solution. I want to get rid of the effect of translation on the solution. I calculated the amount of translation and formed a translation matrix of dimension (n,2) where n is the number of nodes. So each row of the translation matrix records the amount of translation on the coordinate of the nodes

MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import cofac
from dolfin import *
import numpy as np
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File

# 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}

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C = C2 - C1

mesh = generate_mesh(C, 100)

class inner(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class outer(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class annulus(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)
outer().mark(boundary_markers, 2)
annulus().mark(surface_markers, 3)

#V = VectorFunctionSpace(mesh, "Lagrange", 1)
#mesh11 = project(mesh1,V)

plot(mesh)
plt.show()

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


V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)

du = TrialFunction(V)
v  = TestFunction(V)             
u  = Function(V)                 
T = Constant(5.0)

bcs=[]

# Kinematics
d = mesh.geometry().dim()
I = Identity(d)             
F = I + grad(u)             
C = F.T*F                   
M = cofac(F)
Finvt = (inv(F)).T
Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.499   
mu, lmbda = 27.9, Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) + ((nu)*(mu)/(1-2*nu))*(J-1)**2 - mu*(ln(J))

Coords = Expression(('x[0]','x[1]'),degree=0)
coords = project(Coords, V)
phi = coords+u

Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds(1)
ddPI = derivative(Pi, u, v)
solve(ddPI == 0, u, bcs,
     form_compiler_parameters=ffc_options)
     
mesh1 = Mesh(mesh) 
X = mesh1.coordinates()
X += np.vstack(map(u, X))

mesh_moved = Mesh(mesh1)
outer_nodes = BoundaryMesh(mesh_moved, "exterior", True)
translated_center = np.mean(outer_nodes.coordinates(), axis=0)

n = np.shape(mesh_moved.coordinates())[0]
translation = np.zeros([n, 2])
translation[:, 0] = translated_center[0]; translation[:, 1] = translated_center[1]
mesh_moved.coordinates()[:] = mesh_moved.coordinates() - translation

trans = project(translation, V)
u1 = Function(V)
u1.vector()[:] = u.vector() - trans.vector()

file = File("MWE.pvd");
file << u1;

But there is a problem with the project function and the error is:

trans = project(translation, V)
  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: [[-3.91119991 -4.4667194 ]
 [-3.91119991 -4.4667194 ]
 [-3.91119991 -4.4667194 ]
 ...
 [-3.91119991 -4.4667194 ]
 [-3.91119991 -4.4667194 ]
 [-3.91119991 -4.4667194 ]] can not be converted to any UFL type.

Any suggestions would be extremely helpful.

2

You should make a Vector cg1 function space, and a corresponding function. With this function use the dolfin.vertex_to_dof_map to assign the data to u.vector().

1 Like