Hi @hernan , here is the MWE. The problem is the solution “u” is a n1 column vector, where n is the number of nodes of the mesh. But the vector “translation” is n2. This records the X and the Y components of the translation for each of the nodes. Please suggest what should be done here. Need help with this.
MWE:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import cofac
from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
# 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}
import numpy as np
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
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;
