Update node coordinates in paraview

1

Hi all. I’m trying to solve a problem of a hyperelastic material, which is subjected to a displacement load. I want to visualize the calculated displacement and show the deformation of the material by paraview. But I don’t know how to update the node coordinates to show the deformation of the material. The corresponding codes and pictures are attached below. Does anyone have experience dealing with similar issues？Thanks in advance.

import matplotlib.pyplot as plt
from dolfin import *
from fenics import *
from ufl import indices

Optimization options for the form compiler

parameters[“form_compiler”][“cpp_optimize”] = True
parameters[“form_compiler”][“representation”] = “uflacs”

parameters[“form_compiler”][“quadrature_degree”]=4

#mesh=RectangleMesh(Point(0,0),Point(1,1),1,1,“left/right”)
mesh=UnitSquareMesh(20,20)

mesh = UnitSquareMesh.create(2,2, CellType.Type.quadrilateral)

#plot(mesh)

V = VectorFunctionSpace(mesh, “Lagrange”, 1)

Mark boundary subdomians

bot = CompiledSubDomain(“near(x[1], side) && on_boundary”, side = 0.0)
top = CompiledSubDomain(“near(x[1], side) && on_boundary”, side = 1.0)
applied_disp = Expression(“disp”, disp=0.0, degree=2)

bcl = DirichletBC(V, Constant((0,0)), bot)
bcr = DirichletBC(V.sub(1), applied_disp, top)
bcs = [bcl, bcr]

Define functions

du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
u_old = Function(V)

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

def safeSqrt(x):
return sqrt(x+DOLFIN_EPS)

Kinematics

d = len(u)
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = variable(F.T*F) # Right Cauchy-Green tensor

q1=(F[0,1]-F[1,0])**2
q2=(F[0,0]+F[1,1])**2
q3=(F[0,1]+F[1,0])2
q4=(F[0,0]-F[1,1])2
eig1=(safeSqrt(q1+q2)-safeSqrt(q3+q4))/2
eig2=(safeSqrt(q1+q2)+safeSqrt(q3+q4))/2
eig12=eig12
eig22=eig22

def psign(u):
return conditional(lt(u,1),1,u)

def nsign(u):
return conditional(le(1,u),1,u)

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/(3*(1-2*nu)))
eps=1E-8

psi=(mu/2*(eig12+eig22-2-2ln(eig1)-2ln(eig2))+0.5lmbda(ln(J))**2)*dx

psi1=(mu/2*(psign(eig12)+psign(eig22)-2-2ln(psign(eig1))-2ln(psign(eig2)))+0.5lmbda(ln(psign(J)))**2)dx
psi2=(mu/2(nsign(eig12)+nsign(eig22)-2-2ln(nsign(eig1))-2ln(nsign(eig2)))+0.5lmbda(ln(nsign(J)))**2)*dx

F=derivative(psi1,u,v)+derivative(psi2,u,v)
J = derivative(F, u, du)

t=0
ur=1.0
deltaT=0.1
tol=1e-3
maxiteration=25
AM_tolerance=1E-3

conc_f=File("./hyperelastic2/u.pvd")
fname=open(‘F_D_curve1.txt’,‘w’)
problem=NonlinearVariationalProblem(F,u,bcs,J=J)
solver=NonlinearVariationalSolver(problem)

#staggered scheme
while t<=1.0:
t+=deltaT
if t>=0.4:
deltaT=0.05
applied_disp.disp=t*ur
iter=0
err=1
while err>AM_tolerance and iter<maxiteration:
iter +=1

    solver.solve()
    u_err=u.vector()-u_old.vector()
    err=u_err.norm('linf')
    u_old.assign(u)
    if err < tol:
        print('Iterations:', iter, ',Total time', t)
        conc_f << u

微信截图_20210426215415
微信截图_20210426215415771×603 118 KB

2

In Paraview there is a filter called Warp by vector that you can apply to your function.
See for instance Using Paraview for visualization — FEniCS-X tutorial

1 Like
3

Dear dokken, thanks for your response, the problem has been addressed.

4

Dear dokken, I have encountered another question: when I tried to solve a thermoelastic problem, I want to import the obtained temperature and displacement into the vtu file simultaneously, so that I can observe the temperature distribution and material deformation by paraview, the command “conc_f << u,T” (u denotes the displacement, T represents the obtained temperature) seems incorrect for this problem, do you have any suggestion to me?

5

You should use XDMFFile for this:

from dolfin import *

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
u = Function(V)
p = Function(Q)
outfile = XDMFFile("output.xdmf")
t = 0
outfile.write(u, t)
outfile.write(p, t)

and then follow the instructions about opening the file and extracting blocks in my previous post.

1 Like
6

Thanks for your response!