Hi Community,
I am solving a nonlinear elasticity problem and I’m going to show stress on the deformed configuration.
For this aim I am getting a VTK file for my deformation as follows:
#generating 2D mesh
domain_vertices = [Point(0.0, 0.0),
Point(48.0, 44.0),
Point(48.0, 60.0),
Point(0.0, 44.0)]
domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,5)
CGG = FiniteElement(“CG”, mesh.ufl_cell(), 1)
Curll = FiniteElement(“N2curl”, mesh.ufl_cell(), 1) # N1curl and N2curl
Divv = FiniteElement(“BDM”, mesh.ufl_cell(), 1) # RT and BDM
ZZ= FunctionSpace(mesh, MixedElement([CGG, CGG, Curll, Curll, Divv, Divv]))
Build space for vector functions
vFE = VectorElement(‘Lagrange’, mesh.ufl_cell(),1)
W = FunctionSpace(mesh, vFE)
parameters[‘allow_extrapolation’] = True
UU_array = as_vector([U1, U2])
UU = project(UU_array, W)
vtkfile = File(‘deformation.pvd’)
vtkfile << UU
But I am confused on how can I get my stress as a VTK file to be able to show both of them on one figure in paraview.
I have tried to get scalar values of stress as follows:
PP = np.sqrt(norm(P1, norm_type=“L2”)**2
+ norm(P2, norm_type=“L2”)**2)
which P1 and P2 are row vector components of the stress , FYI; P1 = (P1[0], P1[1]) and P2 = (P2[0], P2[1])
and they try to project them into the FunvtionSpace:
Stress = project(PP, W)
But It is not working.
The Error:
ufl.log.UFLException: Shapes do not match: and FloatValue d=139818763805864>.
I think the problem with my projection arises from FunctionSpace.
I appreciate any hint or comment to resolved this issue.
Thanks.
Hi,
W is a vectorial function space where PP is scalar, also don’t use Numpy sqrt but UFL
Best
Hi,
beside what bleyerj pointed out: Regarding your original question I would suggest using xdmf files instead of VTK. Although with ParaView filters it is possible to visualize on the deformed configuration, it is a whole lot easier with xdmf.
Check out the following example code where I made up a displacement field and calculate the linear elastic stress response (notice that I define the Frobenius norm of the stress in directly UFL):
from dolfin import *
mesh = UnitSquareMesh(16,16)
dFE = FiniteElement("DG", mesh.ufl_cell(), 0)
FE = FiniteElement("Lagrange", mesh.ufl_cell() ,1)
vFE = VectorElement(FE)
tFE = TensorElement(dFE)
U = FunctionSpace(mesh, dFE)
V = FunctionSpace(mesh, vFE)
W = FunctionSpace(mesh, tFE)
u = Function(V) # displacement
uVals = Expression(('x[0]','x[0]*pow(x[1],2)'), element=vFE) # Just an example displacement field
u.interpolate(uVals)
eps = sym(grad(u)) # UFL definition of strains
mu = 80e+9
lam = 120e+9
sig = 2.*mu*eps + lam*tr(eps)*Identity(2) # UFL definition of stress
fileResults = XDMFFile("output.xdmf")
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True
disp = Function(V, name='Displacement')
disp.assign(u)
fileResults.write(disp,0.)
stress = Function(W, name='Stress')
stress.assign(project(sig, W))
fileResults.write(stress,0.)
stressNorm = Function(U, name="StressNorm")
stressNorm.assign(project(sqrt(inner(sig,sig)),U))
fileResults.write(stressNorm,0.)
Opened in ParaView, you can access all fields you saved to the same output file as seen on the top left in the following screenshot.
Then you can warp by vector and color by scalars.
4 Likes
