### Finding the nodes coordinates after displacement.

318

views

0

Hi,

For the hyperelasticity problem below:

For the hyperelasticity problem below:

```
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True}
# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"))
r = Expression(("scale*0.0",
"scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
"scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3)
bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(F, u, du)
# Solve variational problem
solve(F == 0, u, bcs, J=J,
form_compiler_parameters=ffc_options)
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
# Plot and hold solution
plot(u, mode = "displacement", interactive = True)
```

How can I see the coordinates of the original mesh and the deformed mesh (possibly with the same order)?
Community: FEniCS Project

### 3 Answers

2

To see mesh coordnates you can use these commands

# Mesh related commands

print(mesh.coordinates())

print(mesh.num_cells())

print(mesh.num_vertices())

print(str(mesh))

Thanks for your answer. Do you know how to get the same information for the deformed mesh?

# Mesh related commands

print(mesh.coordinates())

print(mesh.num_cells())

print(mesh.num_vertices())

print(str(mesh))

written
5 months ago by
Nima

0

also see fenics book and tutorials .. there is mention of

**mesh editor**and**dofmap**which can clarify your concept about the way mesh nodes and coordinates are linked and ordered.0

If you just want to see the deformed shape, the most convenient way would be to open your displacement pvd-file in Paraview and use the "Warp By Vector" filter on the displacement field.

Thank you so much. But the problem is that I need to see the coordinates for the deformed mesh as a list.

written
5 months ago by
Nima

1

How about taking your mesh coordinates and simply adding your displacement field?

Can't try it out myself right now but something along the lines

Can't try it out myself right now but something along the lines

```
x = SpatialCoordinate(mesh)
deformed_coords = project(x+u, V)
print deformed_coords.vector().array()[:]
```

may work, where `u`

is your displacement and `V`

the function space `u`

comes from.

written
5 months ago by
klunkean

Thanks again. I tried your code. However, there is something that concerns me. For the original mesh I used

As for your code I get an array of 21675 numbers and clearly this is not consistent. Because if resize into 3-tuples (which I don't know how) it will give us much less 3D points than dof_coordinates.

```
N = V.dim()
d = mesh.geometry().dim()
dof_coordinates = V.tabulate_dof_coordinates()
dof_coordinates.resize( ( N , d ) )
```

which first gives me an array of 65025 numbers for dof_coordinates and then resizes it to 3-tuples as points in 3D, therefore a 21675x3 matrix. As for your code I get an array of 21675 numbers and clearly this is not consistent. Because if resize into 3-tuples (which I don't know how) it will give us much less 3D points than dof_coordinates.

written
5 months ago by
Nima

Please login to add an answer/comment or follow this question.