writing the solution in to the txt

Asked by Bahram on 2013-05-06

Hello,
I am very new to Fenics, so sorry if this is question is very simple.
I want to save the solution in a format that I can read automatically by matlab. the problem is that I can use u_nodal_values =

u.vector()
u_array = u_nodal_values.array()
coor = mesh.coordinates()
numpy.savetxt('PointCoord.txt', u_nodal_values, delimiter=',')
numpy.savetxt('MeshCoord.txt', coor, delimiter=',')

to save the coordinate and the displacement in to two different file but since I working on 3D domain, I have 3 displacement for each node and I am not sure how these are connected.
I convert the PVD file to the CSV format in paraview manually but I am not sure how can I save the nodes and corresponding displacment to each node in to a txt file.
any suggestion?

Thanks

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Nick Davies
Solved:
2013-05-08
Last query:
2013-05-08
Last reply:
2013-05-08
Best Nick Davies (ntd14) said : #1

I will try to help, but I dont use matlab much, and not 100 % sure what you are asking, but I think it is how can I link my displacements to my coordinates. The problem is that the u_vec and the coordinates dont have the same ordering. I have a script which does a similar thing to what you are wanting and you are welcome to use and abuse, it is modified from a script a different member (Johan Hake) on here posted. For reference what this code actually does is take the coordinates of the mesh, updates them with the displacements and returns them as new coordinates.

#these are probably already in your code as the problem setup
V = VectorFunctionSpace(mesh1, "Lagrange", 1)
VV = FunctionSpace(mesh1, "Lagrange", 1)
.........
solve()
####

coor_int = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), V).vector().array()
    #turing displacement into numpy array
    u_vec = u.vector().array()
    #calc new coordinates by subtaracting the displacement at each interploated
    #node from the oringonal mesh
    new_coor = coor_int+u_vec

    #initalising matries
    dofs_to_vert = np.zeros(mesh1.num_vertices(), dtype=np.uintp)
    vectordofs_to_vert = np.zeros((mesh1.num_vertices()*3), dtype=np.uintp)
    vectordofs_to_subvert = np.zeros((mesh1.num_vertices()*3),dtype=np.uintp)
    cellinds = np.zeros((mesh1.num_cells(),4), dtype=np.uintp)
    #creating degrees of freedom maps
    dm = VV.dofmap()
    dms = [V.sub(i).dofmap() for i in range(3)]

    for cell in cells(mesh1):
        cell_ind = cell.index()
        vert_inds = cell.entities(0)
        cellinds[cell_ind,:] = vert_inds #retuns the local indacies for each cell
        dofs_to_vert[dm.cell_dofs(cell_ind)] = vert_inds #dofs to vertcies map
        for i, (dms_i, dmcs_i) in enumerate(zip(dms, dms)):
            vectordofs_to_vert[dms_i.cell_dofs(cell_ind)] = vert_inds #gives map for coords to cells, not sencitive to xyz position
            vectordofs_to_subvert[dms_i.cell_dofs(cell_ind)] = i #gives map for each x,y,z to to particular cells for map above
            #initalise more arrays
    map_mat = np.zeros((len(new_coor),3), dtype=float)
    coorxyz = np.zeros((len(new_coor)/3,3), dtype=float)
    #make matrix of maps above
    map_mat[:,0] = vectordofs_to_vert
    map_mat[:,1] = vectordofs_to_subvert
    map_mat[:,2] = new_coor
    #sort entries using maps matrix above
    for ij in range(0,len(map_mat)):
            if map_mat[ij,1] == 0:
                coorxyz[map_mat[ij,0],0] = map_mat[ij,2]
            if map_mat[ij,1] == 1:
                coorxyz[map_mat[ij,0],1] = map_mat[ij,2]
            if map_mat[ij,1] == 2:
                coorxyz[map_mat[ij,0],2] = map_mat[ij,2]

Note this only works for 3D and linear elements.
coorxyz is a matrix with all of the coordinates in x, y and z columns after the displacements have been added.
You could probably modify this to deal with the displacments and coordinates separately but keep them all in the right order so that the two matrices match when imported into matlab if you want them separate.

here is a link to the original code which might help too
https://answers.launchpad.net/fenics/+question/215827

Good luck

Bahram (behinoo) said : #2

Thanks Nick Davies, that solved my question.