Finding the nodes coordinates after displacement.

6 months ago by

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,

# 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

6 months ago by
To see mesh coordnates you can use these commands 

# Mesh related commands

Thanks for your answer. Do you know how to get the same information for the deformed mesh?
written 5 months ago by Nima  
6 months ago by
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.
6 months ago by
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  
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
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
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.

Similar posts:
Search »