How to translate a 2D mesh by a vector?

1

I need to translate the mesh nodes in a 2D geometry by a vector \textbf{u} of size (2,1). Since mesh coordinates are geometry-type objects and not a function object, I tried to project the mesh_node coordinates to a vector function space and then tried adding \textbf{u}. Consider the below MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

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 *
from mshr import *
import numpy as np

        
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), 5)

mesh = generate_mesh(C, 50)

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

mesh_nodes = mesh.coordinates()
print(mesh_nodes)
u = Function(V)
coords = project(mesh_nodes, V)
phi = coords + u

But this raises the error:

coords = project(mesh_nodes, 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: [[-4.99899902  0.10004404]
 [-4.99899902 -0.10004404]
 [-4.99099357 -0.29997192]
 ...
 [ 1.81895092  4.25595418]
 [ 1.60783276 -2.54317265]
 [ 1.12755794 -2.61691696]] can not be converted to any UFL type.

What does the error mean and what can be done to overcome that?

2

You can displace the vertices of a mesh by a Function in a linear continuous VectorFunctionSpace by using the ALE.move() function, e.g.,

V = VectorFunctionSpace(mesh,"CG",1)
u = Function(V)

# (Set `u` to some desired nonzero vector field)

ALE.move(mesh,u)

where the vector field u has the interpretation of a displacement added to the initial mesh vertex positions.

3 Likes
3

Hi @kamensky , thanks for the suggestion. Will V = VectorFunctionSpace(mesh,"Lagrange",1) work?

4

Yes, "CG" (standing for “continuous Galerkin”) and "Lagrange" are equivalent function spaces.

5

The code seems to work for the MWE I created. But in my actual code, it raises an error. Here is my original code followed by the gmsh script:

Script:

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
from ufl import cofac, sqrt
import numpy
import matplotlib.pyplot as plt

# 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
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

import meshio
msh = meshio.read("annulusgeometry.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)

# Create mesh and define function space
V = VectorFunctionSpace(mesh, "Lagrange", 1) 
# Definining the normal to each of the faces
n=FacetNormal(mesh)
# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, 0.0))  # Body force per unit volume
T = Constant(-5.0)
c = Constant((0.0,0.0))

uout=DirichletBC(V, c, mf, 1)

bcs=[]

# Kinematics
d = mesh.geometry().dim() 
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor
M = cofac(F)


# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

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

phi = ALE.move(mesh,u)

Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds_custom(2)
F = derivative(Pi, u, v)
J = derivative(F, u, du)

solve(F == 0, u, bcs, J=J,
     form_compiler_parameters=ffc_options)
file = File("test.pvd");
file << u;

gmsh script:

//+
SetFactory("OpenCASCADE");
Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Line Loop(1) = {1};
//+
Line Loop(2) = {2};
//+
Plane Surface(1) = {1, 2};
//+
Physical Surface("S1", 3) = {1};
//+
Physical Line("C1", 1) = {1};
//+
Physical Line("C2", 2) = {2};

Error message:

Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds_custom(2) 
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 173, in dot
    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: None can not be converted to any UFL type.

But the same code works if I replace your part of the code phi = ALE.move(mesh,u) with the following:

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

ALE.move implicitly moves your mesh, and does not return any field.
If you want to integrate over the new coordinates of the mesh, I suggest you use ALE.move to move your mesh, and then use

x=SpatialCoordinate(mesh)
phi = as_vector((x[0], x[1], x[2]))
1 Like
7

Hi @dokken , that worked. However, I have a question just to clarify myself.

When you say “does not return any field”, do you mean doesn’t return coordinates of the moved mesh?
Also the mesh argument in x=SpatialCoordinate(mesh) is the already moved mesh, right?

8

ALE.move changes the coordinates of the input mesh (and thus does not return an explicit array).

SpatialCoordinate holds a reference to the mesh, and is only evaluated at assembly, and thus contains the coordinates of the moved mesh if used after ALE.move.

9

Got it. Thanks a lot!