Volume after deformation

1

Hi,
I am currently stuck at the following problem. I modified a demo file where rectangular mesh is deformed under its own weight:

mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(L, W, H), 20, 8, 4)
V = VectorFunctionSpace(mesh, ‘P’, 1)

with fixed boundaries:

left = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 0.0)
right = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 10.0)

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]

The solution is saved:
u = Function(V)
solve(a == L, u, bcs)

File(‘elasticity/displacement.pvd’) << u

Now i would like to know the local volume change, lets say from the four mesh-volumes in the center of the rectangle. So i would need to calculate the inital volume (which is ok, one can do that just using the dimension and resolution) and the volume after displacement.

See the picture below, i guess then it gets more clear.
Thank you

hmmm
hmmm.png914×442 149 KB

2

Hi, consider the following

from dolfin import *
import numpy as np

# Undeformed configuration
mesh0 = UnitSquareMesh(32, 32)
File('before.pvd') << mesh0
# Displacement
V = VectorFunctionSpace(mesh0, 'CG', 2)
u = interpolate(Expression(('x[0]*x[1]', '0'), degree=2), V)

# Deformation takes the form x = F(X) = X + u(X) so then
# volume element, dv, in deformed configuration is dv = det(F)*dV,
# with dV the volume element of reference configuration. One way of 
# getting the deformed volume is therefore to integrate det(F)*dV over
# undeformed domain. 
volume_after0 = assemble(det(Identity(2) + grad(u))*dx)

# Alternatively actually deform the mesh by displacement 
mesh1 = Mesh(mesh0)
X = mesh1.coordinates()
X += np.vstack(map(u, X))

File('after.pvd') << mesh1
# and compute the volume of the deformed mesh
volume_after = assemble(Constant(1)*dx(domain=mesh1))

# Should be the same
print(volume_after0, volume_after)

# For doing this over just part of the domains ...
cell_f0 = MeshFunction('size_t', mesh0, mesh0.topology().dim(), 0)
# Left half
CompiledSubDomain('x[0] < 0.5 + DOLFIN_EPS').mark(cell_f0, 1)
dV = Measure('dx', domain=mesh0, subdomain_data=cell_f0)
volume_after0 = assemble(det(Identity(2) + grad(u))*dV(1))

# Assume here that no cells disappeared/got entangled
cell_f1 = MeshFunction('size_t', mesh1, mesh1.topology().dim(), 0)
cell_f1.array()[:] = cell_f0.array()

dV = Measure('dx', domain=mesh1, subdomain_data=cell_f1)
# and compute the volume of the deformed mesh
volume_after = assemble(Constant(1)*dV(1))

print(volume_after0, volume_after)
2 Likes
3

Hi,
thank you very much. I am starting to understand more :slight_smile:

And it works just fine, I just had to change to Identity(3) (for a box) and adjust the parameters.
Thank you very much :slight_smile:

4

Sorry, me again. The results seem to be wrong.
Using:

CompiledSubDomain(‘x[0] < 11.0 + DOLFIN_EPS’).mark(cell_f0, 1)
dV = Measure(‘dx’, domain=mesh0, subdomain_data=cell_f0)
volume_before1 = assemble(det(Identity(3) + grad(u))*dV(1))
print(“Volume before of left half”, volume_before1)

I find a volume of 92.8. But that should be the volume before deformation, right? So 80.

mesh0 = BoxMesh(Point(0.0, 0.0, 0.0), Point(10.0, 4.0, 2.0), 20, 8, 4)
V = VectorFunctionSpace(mesh, “Lagrange”, 1)

Thanks for your help again :slight_smile:

Edit:
Replacing

volume_before1 = assemble(det(Identity(3) + grad(u))*dV(1))

with:

volume_before1 = assemble(Constant(1)*dV(1))

seems to work.
Best

5

Can you post full (and formated) code? Anyways, both volumes are volumes of the deformed object. The difference is that one is computed without actually deforming it, while in the other case you deform the mesh and integrate.

6

Aa ok, then i misunderstood.
Here the Code:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np

Scaled variables

Len = 10.0; W = 4.0; H=2.0;
mu = 1
rho = 1
delta = W/Len
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, “Lagrange”, 1)

Define boundary condition

tol = 1E-14

Mark boundary subdomians

left = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 0.0)
right = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 10.0)

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]

Define strain and stress

def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))

def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(u)

Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

Compute solution

u = Function(V)
solve(a == L, u, bcs)

Save solution to file in VTK format

File(‘elasticity/displacement.pvd’) << u

#volume of a cell
#c1=Cell(mesh, 7)
#print(c1.volume())

######Volume before and after deformation###############
volume_before0 = assemble(Constant(1)*dx(domain=mesh))

Deformation takes the form x = F(x) = X + u(X) so then

volume element, dv, in deformed configuration is dv = det(F)*dV,

with dV the volume element of reference configuration. One way of

getting the deformed volume is therefore to integrate det(F)*dV over

undeformed domain.

volume_after0 = assemble(det(Identity(3) + grad(u))*dx)
print("Volume before ", volume_before0)
print("Volume after ", volume_after0)
###########
###########
##another way is to make new mesh, and then actually deform it:
mesh1 = Mesh(mesh)
X = mesh1.coordinates()
X += np.vstack(map(u, X))
volume_after = assemble(Constant(1)*dx(domain=mesh1))
print(“Volume after 2nd way”, volume_after)
##########################################################

For doing this over just part of the domains …

##########################################################
cell_f0 = MeshFunction(‘size_t’, mesh, mesh.topology().dim(), 0)

Right half

CompiledSubDomain(‘x[0] < 4.0 + DOLFIN_EPS’).mark(cell_f0, 1)
dV = Measure(‘dx’, domain=mesh, subdomain_data=cell_f0)
#volume_before1 = assemble(det(Identity(3) + grad(u))*dV(1))
volume_before1 = assemble(Constant(1)*dV(1))
print(“Volume before small half”, volume_before1)

Assume here that no cells disappeared/got entangled

cell_f1 = MeshFunction(‘size_t’, mesh1, mesh1.topology().dim(), 0)
cell_f1.array()[:] = cell_f0.array()
dV = Measure(‘dx’, domain=mesh1, subdomain_data=cell_f1)

and compute the volume of the deformed mesh

volume_after1 = assemble(Constant(1)*dV(1))
print(“Volume after of small half:”, volume_after1)

and compute the volume of the deformed mesh

cell_f0 = MeshFunction(‘size_t’, mesh, mesh.topology().dim(), 0)

Right half

CompiledSubDomain(‘x[0] < 6.0 + DOLFIN_EPS’).mark(cell_f0, 1)
dV = Measure(‘dx’, domain=mesh, subdomain_data=cell_f0)
#volume_before1 = assemble(det(Identity(3) + grad(u))*dV(1))
volume_before2 = assemble(Constant(1)*dV(1))
print(“Volume before large half”, volume_before2)
cell_f1 = MeshFunction(‘size_t’, mesh1, mesh1.topology().dim(), 0)
cell_f1.array()[:] = cell_f0.array()
dV = Measure(‘dx’, domain=mesh1, subdomain_data=cell_f1)

and compute the volume of the deformed mesh

volume_after2 = assemble(Constant(1)*dV(1))
print(“Volume after of large half:”, volume_after2)
print("--------------")
print(“Volume before center:”, volume_before2-volume_before1)
print(“Volume after center:”, volume_after2-volume_after1)

7 
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np
# Scaled variables

Len = 10.0; W = 4.0; H=2.0;
mu = 1
rho = 1
delta = W/Len
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
# Define boundary condition

tol = 1E-14
# Mark boundary subdomians

left = CompiledSubDomain('near(x[0], side) && on_boundary', side = 0.0)
right = CompiledSubDomain('near(x[0], side) && on_boundary', side = 10.0)

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]
# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
# Compute solution

u = Function(V)
solve(a == L, u, bcs)


######Volume before and after deformation###############
volume_before0 = assemble(Constant(1)*dx(domain=mesh))
volume_after0 = assemble(det(Identity(3) + grad(u))*dx)
print("Volume before ", volume_before0)
print("Volume after ", volume_after0)

##another way is to make new mesh, and then actually deform it:
mesh1 = Mesh(mesh)
X = mesh1.coordinates()
X += np.vstack(map(u, X))
volume_after = assemble(Constant(1)*dx(domain=mesh1))
print('Volume after 2nd way', volume_after)


cell_f0 = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)

CompiledSubDomain('x[0] < 4.0 + DOLFIN_EPS').mark(cell_f0, 1)
dV = Measure('dx', domain=mesh, subdomain_data=cell_f0)

volume_before = assemble(Constant(1)*dV(1))
print('Volume before small half', volume_before)

volume_after0 = assemble(det(Identity(3) + grad(u))*dV(1))
print('Volume after (0) small half', volume_after0)

cell_f1 = MeshFunction('size_t', mesh1, mesh1.topology().dim(), 0)
cell_f1.array()[:] = cell_f0.array()
dv = Measure('dx', domain=mesh1, subdomain_data=cell_f1)

volume_after1 = assemble(Constant(1)*dv(1))
print('Volume after (1) small half', volume_after1)
8

Nice, thanks alot :slight_smile:

9

Hello,
may i ask another question.
How can I save the mesh of the deformed subdomain only?

I have tried the following:

mesh1 = Mesh(mesh)
X = mesh1.coordinates()
X += np.vstack(map(u, X))
cell_f5 = MeshFunction(‘size_t’, mesh1, mesh.topology().dim(), 0)
CompiledSubDomain(‘x[0] < 4.0 + DOLFIN_EPS’).mark(cell_f5, 1)
File(‘liner_def/subdomain1.pvd’) << cell_f5

But that saves the whole mesh of the deformation, with value 1 for the subdomain and value 0 for the rest. I would need only the deformed submesh.
Also consider that i would need to save more infomation that just the coordinates, e.g. the stress, like in the demo:

s = sigma(u) - (1./3)*tr(sigma(u))Identity(d) # deviatoric stress
von_Mises = sqrt(3./2inner(s, s))
V = FunctionSpace(mesh, ‘P’, 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title=‘Stress intensity’)

If only 0 or 1 is written therem, it would not be os much use ^^.

Thank you very much for your help again :slight_smile:

10

There is another Problem: With

cell_f0 = MeshFunction(‘size_t’, mesh, mesh.topology().dim(), 0)
CompiledSubDomain(‘x[2] > 1.0 + DOLFIN_EPS’).mark(cell_f0, 1)

one gets the layer of which the conditions applies after deformation. Actually the condition was meant for the inital setup. For example look at this picture:

test1
test1888×420 18.2 KB

The fist upper layer of the mesh is marked, but only unti x[2]>1.0. But what I would need is that this condition is applied to the undeformed mesh, i.e. that the complete upper layer is part of the subdomain.
Thanks again ^^

Edit: Would be possible to define such subdomains from the very beginning? Because in the end i want to have exactly this, different layers/domains with different material properties…

11

Me again.
So I can use the subdomains defined on the inital Mesh to calculate the Volume:

dV = Measure(‘dx’, domain=mesh, subdomain_data=subdomains)
volume_before0 = assemble(Constant(1)*dV(0))
print(“Volume of upper Layer before deformation”, volume_before0)
volume_after1 = assemble(det(Identity(3) + grad(u))*dV(0))
print(“Volume of upper Layer after deformation”, volume_after1)
volume_before0 = assemble(Constant(1)*dV(1))
print(“Volume of bottom Layer before deformation”, volume_before0)
volume_after1 = assemble(det(Identity(3) + grad(u))*dV(1))
print(“Volume of bottom Layer after deformation”, volume_after1)

and save the displacement in a file:

File(‘liner_def/domains_after.pvd’) << u

Or show the marked Subdomains of the inital mesh:

File(‘liner_def/domains_before.pvd’) << subdomains

Now i would like to save the marked subdomains of the deformed mesh, i.e. the solution u with the subdomains defined before.
Is there a way?
Thanks again :slight_smile:

ping