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) + 2*muepsilon(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)