# Unexpected dimensions after matrix assembly

**URL:** <https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090>\
**Category:** General\
**Created:** [June 27, 2024, 4:57pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090 "2024-06-27T16:57:49Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![arash\_kazemi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/arash_kazemi/32/5457_2.png) [@arash\_kazemi](https://fenicsproject.discourse.group/u/arash_kazemi)\
**Post date:** [June 27, 2024, 4:57pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/1 "2024-06-27T16:57:49Z")

</div>

Hey everyone,

I’m trying to delve into my weak form in more detail, but I’m having trouble understanding how assembly works in FEniCSx. Specifically, I want to assemble two matrices, `Kau` and `Kua` , which are defined as follows:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/71e260f78097682560ad7f86e5847a124061ca59.png)  
which finally form a system of  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/c584635e10c371f13db45d258f8c31fd2826fcac.png)

`u` is a vector, and `alpha` is a tensor element (3 by 3). Since I have 8 nodes in my element, I expect `K_{au}` to be an 81 by 72 matrix and `K_{ua}` to be a 72 by 81 matrix. However, this is not the case, and I don’t know why. Here is my code to check the dimensions:

```auto
# Local assembly of cell integral or exterior facet integrals
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT
import cffi
# import numpy as np
# import ufl
# from mpi4py import MPI
# from petsc4py.PETSc import ScalarType
# from ufl import dx, grad, inner, ds
from dolfinx import fem, mesh
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
                         locate_dofs_topological, assemble, locate_dofs_topological, Expression)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type

_type_to_offset_index = {fem.IntegralType.cell: 0, fem.IntegralType.exterior_facet: 1}

class LocalAssembler():

    def __init__ (self, form, integral_type: fem.IntegralType =fem.IntegralType.cell):
        self.form = fem.form(form)
        self.integral_type = integral_type
        self.update_coefficients()
        self.update_constants()

        subdomain_ids = self.form._cpp_object.integral_ids(integral_type)
        assert(len(subdomain_ids) == 1)
        assert(subdomain_ids[0] == -1)
        is_complex = np.issubdtype(ScalarType, np.complexfloating)
        nptype = "complex128" if is_complex else "float64"
        o_s = self.form.ufcx_form.form_integral_offsets[_type_to_offset_index[integral_type]]
        o_e = self.form.ufcx_form.form_integral_offsets[_type_to_offset_index[integral_type]+1]
        assert o_e - o_s == 1

        self.kernel = getattr(self.form.ufcx_form.form_integrals[o_s], f"tabulate_tensor_{nptype}")
        self.active_cells = self.form._cpp_object.domains(integral_type, -1)
        assert len(self.form.function_spaces) == 2
        self.local_shape = [0,0]
        for i, V in enumerate(self.form.function_spaces):
            self.local_shape[i] = V.dofmap.dof_layout.block_size * V.dofmap.dof_layout.num_dofs

        e0 = self.form.function_spaces[0].element
        e1 = self.form.function_spaces[1].element
        needs_transformation_data = e0.needs_dof_transformations or e1.needs_dof_transformations or \
            self.form._cpp_object.needs_facet_permutations
        if needs_transformation_data:
            raise NotImplementedError("Dof transformations not implemented")

        self.ffi = cffi.FFI()
        V = self.form.function_spaces[0]
        self.x_dofs = V.mesh.geometry.dofmap

    def update_coefficients(self):
        self.coeffs = fem.assemble.pack_coefficients(self.form._cpp_object)[(self.integral_type, -1)]

    def update_constants(self):
        self.consts = fem.assemble.pack_constants(self.form._cpp_object)

    def update(self):
        self.update_coefficients()
        self.update_constants()

    def assemble_matrix(self, i:int, local_index:int = 0):

        x = self.form.function_spaces[0].mesh.geometry.x
        x_dofs = self.x_dofs[i]
        geometry = np.zeros((len(x_dofs), 3), dtype=np.float64)
        geometry[:, :] = x[x_dofs]

        A_local = np.zeros((self.local_shape[0], self.local_shape[0]), dtype=ScalarType)
        facet_index = np.array([local_index], dtype=np.intc)
        facet_perm = np.zeros(0, dtype=np.uint8)
        if self.coeffs.shape == (0, 0):
            coeffs = np.zeros(0, dtype=ScalarType)
        else:
            coeffs = self.coeffs[i,:]
        ffi_fb = self.ffi.from_buffer
        self.kernel(ffi_fb(A_local), ffi_fb(coeffs), ffi_fb(self.consts), ffi_fb(geometry),
               ffi_fb(facet_index), ffi_fb(facet_perm))
        return A_local

#####################################################################
# ......... Define parameteres
Fl = 1e-8
h = 0.76e-3 / Fl
H = h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
                  [2, 2, 2], cell_type=CellType.hexahedron)
# .... Define Elements
u_el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  
lm_el = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  
# .... Define Function Spaces
V = functionspace(mesh, u_el)
AL = functionspace(mesh, lm_el)
# .... Define Trail and Test Functions
u, lamu = TrialFunction(V), TrialFunction(AL)
del_u, del_lamu = TestFunction(V), TestFunction(AL)
# ..... Create forms
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)

kua = -lamu[j, k] * del_u[k].dx(j) * dx
kau = del_lamu[j, k] * u[k].dx(j) * dx
#...................
assembler = LocalAssembler(kua)
assemblerB = LocalAssembler(kau)
for cell in range(mesh.topology.index_map(mesh.topology.dim).size_local):
    if cell == 7:
        A = assembler.assemble_matrix(cell)
        B = assemblerB.assemble_matrix(cell)
        print(A.shape)
        print(B.shape)

```

Here is the output I get,  
 ![Screenshot from 2024-06-27 10-38-18](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8d8bb67da0002868c67fbc97ddb2247b2db6a409.png)

The first part of the code (Assembly) I found in another post, and the second part is a simple variational form. Can anyone explain how FEniCSx is computing these matrices?

---

<div class="post-metadata">

**Author:** ![arash\_kazemi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/arash_kazemi/32/5457_2.png) [@arash\_kazemi](https://fenicsproject.discourse.group/u/arash_kazemi)\
**Post date:** [June 28, 2024, 5:17pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/2 "2024-06-28T17:17:00Z")

</div>

Could it be because it should be a square matrix, and we cannot have a matrix with a different number of rows and columns?

---

<div class="post-metadata">

**Author:** ![hherlyng](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hherlyng/32/5843_2.png) [@hherlyng](https://fenicsproject.discourse.group/u/hherlyng)\
**Post date:** [June 28, 2024, 5:39pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/3 "2024-06-28T17:39:01Z")

</div>

Hi Arash,

Just out of curiosity, is this a code you have written yourself, or is the whole code written by Jørgen? (I see he is appropriately credited at the top of the script)

Anyways, in the line  
`A_local = np.zeros((self.local_shape[0], self.local_shape[0]), dtype=ScalarType)`,  
both the row and column dimensions of the matrix are set to the first element of `local_shape`, which is why you get a square matrix.

I’m not really sure about the context here, and specifically what exactly the operator \delta in your expressions mean, so it’s hard to attempt explaining why without speculation.

Cheers,  
Halvor

---

<div class="post-metadata">

**Author:** ![arash\_kazemi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/arash_kazemi/32/5457_2.png) [@arash\_kazemi](https://fenicsproject.discourse.group/u/arash_kazemi)\
**Post date:** [June 28, 2024, 5:56pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/4 "2024-06-28T17:56:22Z")

</div>

I didn’t write the code for the assembly; as I mentioned, I borrowed it from another post. I just used it to check some portions of the variational form in my problem. In my expression, δ means the test function. In Kau, δu is the test finction and in Kua, δα is the test function.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 28, 2024, 7:11pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/5 "2024-06-28T19:11:45Z")

</div>

I probably made an assumption on the test and trial space, and this should be adjusted. I’ll have a look tomorrow or over the weekend.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 1, 2024, 2:50pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/6 "2024-07-01T14:50:04Z")

</div>

> [@arash\_kazemi](#):
>
> ` A_local = np.zeros((self.local_shape[0], self.local_shape[0]), dtype=ScalarType)`

Modifying this to:

```python
        A_local = np.zeros((self.local_shape[0], self.local_shape[1]), dtype=ScalarType)

```

should make it work.

---

<div class="post-metadata">

**Author:** ![arash\_kazemi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/arash_kazemi/32/5457_2.png) [@arash\_kazemi](https://fenicsproject.discourse.group/u/arash_kazemi)\
**Post date:** [July 1, 2024, 5:09pm UTC](https://fenicsproject.discourse.group/t/unexpected-dimensions-after-matrix-assembly/15090/7 "2024-07-01T17:09:20Z")

</div>

It works now, Thanks.