Import MATLAB matrices to Fenics using vector functions


113
views
0
3 months ago by

I am trying to import a series of MATLAB matrices into Fenics using vector functions. I have succeeded doing so using scalar functions (code attached below), but I am having trouble with vector functions. The problem arises from the ordering of the nodes in Fenics functions and vector functions. The below code finds the coordinates in the mesh corresponding to the coordinates in the function space, and sets the value for each node (please note that our MATLAB matrices are ordered according to their coordinates starting from (0,0,0) --> (0,2,0)-->...--> (126,126,38)). I would like to do the same for vector functions, but I don't know how to use the one vertex index that I find to attribute the 3 values of the vector for that node. I would appreciate any help.

Please note that the below code works for the scalar ur (or any other scalar) in the V space, but I would like to make it work for the vector dr = (ur, vr, wr) and a second vector di = (ui, vi, wi) in the Vu space.

import scipy.io as sio
from dolfin import *
import numpy as np

mat_vars = sio.loadmat("our_matrices.mat")
disp = [mat_vars['ux_r'], mat_vars['ux_i'], mat_vars['uy_r'], mat_vars['uy_i'], mat_vars['uz_r'], mat_vars['uz_i']]
mesh = BoxMesh(Point(0.,0.,0.),Point(126.,126.,38.), 63, 63, 19)
V  = FunctionSpace(mesh, 'CG', 1)
Vu = VectorFunctionSpace(mesh, 'CG', 1)
ur = Function(V)

dofmap = V.dofmap()
nvertices = mesh.ufl_cell().num_vertices()

# Set up a vertex_2_dof list
indices = [dofmap.tabulate_entity_dofs(0, i)[0] for i in range(nvertices)]

vertex_2_dof = dict()
[vertex_2_dof.update(dict(vd for vd in zip(cell.entities(0), dofmap.cell_dofs(cell.index())[indices]))) for cell in cells(mesh)]

# Get the vertex coordinates
X = mesh.coordinates()

dof_idx = vertex_2_dof[0]
ur.vector()[dof_idx] = disp[0][0][0]
for k in range(20):
    print("slice %d"%(k))
    for i in range(64):
        for j in range(64):
            vertex_idx = np.where((X == (i*2,j*2,k*2)).all(axis = 1))[0] 
            if not vertex_idx:
                print 'No matching vertex!'
            else:
                vertex_idx = vertex_idx[0]
                dof_idx = vertex_2_dof[vertex_idx]
                ur.vector()[dof_idx] = disp[0][(i) + (j*64) + (k*64*64)][0]​
Community: FEniCS Project

1 Answer


2
3 months ago by
I found the answer for anyone that is interested. First, one has to convert their matrices into single columns where the order is as: 1st node-x, 1st node-y, 1st node-z, 2nd node-x, etc. Then one has to save the column as ascii in a text file. The rest is done in python:

from dolfin import *
import numpy as np
# your mesh goes here, mesh = ...
Vu = VectorFunctionSpace(mesh, 'CG', 1)
var = Function(Vu)
var_np = np.loadtxt(#your txt file goes here)
dofs = dof_to_vertex_map(Vu)
var_np = var_np[dofs]
var.vector().set_local(var_np)​
Please login to add an answer/comment or follow this question.

Similar posts:
Search »