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:
which finally form a system of
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:
# 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,
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?