Assigning a tensor to each element for a transversely isotropic material

dolfinx
1

Hello,
In the problem I have at hand I need to assign to each cell a rotation matrix, which will be applied to the strain tensor. From previous topics, I understood how a tensor is assigned to each node Fenicsx and anisotropic material properties - #4 by dokken. And also in Create a function from cell data - #2 by dokken it was clarified how to assign a scalar to each cell. But the data I need for creating the rotation matrix (3 angles of rotation) is from a vtu file. I don’t exactly know how to assign a nonscalar data to each cell, and how to import the vtu file such so, that each vector of the rotation angles is assigned to their corresponding element, not any wrong element.
Thank you in advance

2

Cell-wise double (tensor data) from VTU does not have native support in DOLFINx.

You should be able to modify:

To make it work. However, as you haven’t provided any code to make a minimal reproducible example, i cannot help you further.

3

Thank you for your attention. Sorry it took a while till I answer. I think I’ve got the procedure. But it seems to me that in the snippet, you have assigned to each node a vector while I’m looking for a vector for each cell. For the use of vtu files I have changed the read_mesh_data function as bellow. Also I have made changed to the definition of element object. I thought it should be a discontinuous, zero-degree polynomial. But for the rest of it I don’t know how I can get the cell indices (as it is done with the node indices with mesh.geometry.input_global_indices) and assign each cell the corresponding vector.

def read_mesh_data(file_path:os.path, mesh: dolfinx.mesh.Mesh, data_path:str):

    msh = meshio.read(file_path)
    num_cells_global = mesh.topology.index_map(mesh.topology.dim).size_global
    num_nodes_global = mesh.topology.index_map(mesh.topology.dim-3).size_global
    dataset = msh.get_cell_data(data_path, "tetra")
    print(dataset.shape)
    shape = dataset.shape
    assert shape[0] == num_cells_global, f"Got data of shape {shape}, expected {num_cells_global, shape[1]}"
    dtype = dataset.dtype
    # Read data locally on each process
    local_input_range = adios4dolfinx.comm_helpers.compute_local_range(mesh.comm, num_nodes_global)    
    local_input_data = dataset[local_input_range[0]:local_input_range[1]]

    # Create appropriate function space (based on coordinate map)
    assert len(mesh.geometry.cmaps) == 1, "Mixed cell-type meshes not supported"
    element = basix.ufl.element(
        basix.ElementFamily.P,
        basix.CellType.tetrahedron,
        0,
        basix.LagrangeVariant.equispaced,
        # In the snippet mesh.geometry.cmaps[0].variant = 1, was used instead of the line above
        discontinuous=True,
        shape=(shape[1],),
        gdim=mesh.geometry.dim)

    # Assumption: Same doflayout for geometry and function space, cannot test in python    
    V = dolfinx.fem.FunctionSpace(mesh, element)
    uh = dolfinx.fem.Function(V, name=data_path)
    # Assume that mesh is first order for now
    assert mesh.geometry.cmaps[0].degree == 1, "Only linear meshes supported"
    x_dofmap = mesh.geometry.dofmap
    igi = np.array(mesh.geometry.input_global_indices, dtype=np.int64)
    global_geom_input = igi[x_dofmap]
    global_geom_owner = adios4dolfinx.utils.index_owner(mesh.comm, global_geom_input.reshape(-1), num_nodes_global)
    for i in range(shape[1]):
        arr_i = adios4dolfinx.comm_helpers.send_dofs_and_recv_values(global_geom_input.reshape(-1), global_geom_owner, mesh.comm, local_input_data[:,i], local_input_range[0])
        dof_pos = x_dofmap.reshape(-1)*shape[1]+i
        uh.x.array[dof_pos] = arr_i

    return uh