Using the output instead of saving it in .pvd

mesh
1

Hi,
For exemple in the hyperelasticity exemple from the doc the result of a simulation is a file .PVD containing the initial mesh ‘coord’ and a displacement vector ‘f_9’
Instead of saving this result into a .pvd file I’d like tu use directly : coord + coef * f_9
with a choosen coefficient. In a format like a simple matrix containing [x,y,z] of the result point.
How can it be done in python?
thanks

2

Hi manu,

without a MWE it is hard to see what you need, but if you have a displacement vector u, you can define a new function which interpolates a position function like this:

V=VectorFunctionSpace(mesh, ‘CG’, 1) # With a given 2D mesh
coords = Expression((“x[0]”, “x[1]”), degree=1)
coord_fun = Function(V)
coord_fun.interpolate(V)
u_new = Function(V)
u_new.vector()[:] += coord_fun.vector() + coef * u.vector()

To make it more efficient, you can use the axpy functionality from petsc instead of the last line:

u_new.vector().vec().axpy(1, coord_un.vector().vec())
u_new.vector().vec().axpy(coef, u.vector().vec())

Best regards

4

Hi nabarnaf,
I think I want the same thing but without vector at the end but an array [[x,y,z], …].
Here is the code I took From https://fenicsproject.org/docs/dolfin/2019.1.0/python/demos/hyperelasticity/demo_hyperelasticity.py.html

import matplotlib.pyplot as plt
from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

c = Constant((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, degree=2)
bcl = DirichletBC(V, c, left)
bcs = [bcl, bcr]

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

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

Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
F = derivative(Pi, u, v)
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 solution
plot(u)
plt.show()

And from this I want to do something like this to extract in the array data my result :

arr = u.vector().get_local()
coor = mesh.coordinates()
vtd = dof_to_vertex_map(V)

values = []
for i, dum in enumerate(coor):
    values.append([arr[vtd[3*i]],arr[vtd[3*i+1]],arr[vtd[3*i+2]]])

coeff=5
data=[]
for i, dum in enumerate(coor):
    data.append([coor[i][0]+coeff*values[i][0],coor[i][1]+coeff*values[i][1],coor[i][2]+coeff*values[i][2]])

but I still don’t have the good result. Is my use of dof_to_vertex_map ok?

5

Well it does work, the error was else where.