Retrieve vector size from c file

dolfinx
1

Hi,

can I retrieve the vector size (block size) from a generated c-file?

For example:

from pathlib import Path

import ffcx.main

cwd = Path.cwd()
infile = cwd / "elasticity.py"

ffcx.main.main(["-o", str(cwd), "--visualise", str(infile)])

with elasticity.py

from dolfinx import fem, mesh
from mpi4py import MPI
from ufl import Identity, TestFunction, TrialFunction, dx, inner, nabla_grad, tr

domain = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.quadrilateral)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = TrialFunction(V)
v = TestFunction(V)
E = 7e9
nu = 0.33
lam = E*nu/((1+nu)*(1-2*nu))
mu = 0.5*E/(1+nu)
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(v):
    return 2.0 * mu * epsilon(v) + lam * tr(epsilon(v)) * Identity(len(v))
a = inner(sigma(u), epsilon(v)) * dx

From elasticity.c I would like to retrieve the vector size (=2).
(Fenics 0.9)

Thanks in advance!

2

I don’t think you can do that.
Why would you need to do this?

3

It would have been helpful to me for testing element matrices with input data of the right size. I found some information in ufcx_form and I thought I might have missed it.

4

What exactly is your workflow with FFCx/DOLFINx?

In the example above, you mix the DOLFINx python interface with the “pure” ufl/basix/ffcx interface.

5

In the example above, you mix the DOLFINx python interface with the “pure” ufl/basix/ffcx interface.

With the pure ufl formulation I ran into the following problem:

MWE:

import basix.ufl
from ufl import FunctionSpace, Identity, Mesh, TestFunction, TrialFunction, dx, grad, inner, tr

element = basix.ufl.element("Lagrange", "quadrilateral", 1)
domain = Mesh(basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(2,)))
V = FunctionSpace(domain, element)

u = TrialFunction(V)
v = TestFunction(V)
E = 7e9
nu = 0.33
lam = E*nu/((1+nu)*(1-2*nu))
mu = 0.5*E/(1+nu)
def epsilon(u):
    return 0.5*(grad(u) + grad(u).T)
def sigma(v):
    return 2.0 * mu * epsilon(v) + lam * tr(epsilon(v)) * Identity(len(v))
a = inner(sigma(u), epsilon(v)) * dx

error message:

Traceback (most recent call last):
  File "/home/bjoern/ffcx/demo/elasticity.py", line 18, in <module>
    a = inner(sigma(u), epsilon(v)) * dx
              ^^^^^^^^
  File "/home/bjoern/ffcx/demo/elasticity.py", line 17, in sigma
    return 2.0 * mu * epsilon(v) + lam * tr(epsilon(v)) * Identity(len(v))
                      ^^^^^^^^^^
  File "/home/bjoern/ffcx/demo/elasticity.py", line 15, in epsilon
    return 0.5*(grad(u) + grad(u).T)
                          ^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 365, in _transpose
    return Transposed(self)
           ^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 106, in __init__
    raise ValueError("Transposed is only defined for rank 2 tensors.")
ValueError: Transposed is only defined for rank 2 tensors.

But would the original example not be the right way to formulate the problem using dolfinx?

What exactly is your workflow with FFCx/DOLFINx?

I generate c-files of different order/ element type and then I call he tabulate_tensor_integral for different inputs (to test vs. comsol). However, since I only test a few problems, I can hard code the vector size as well.

6

You get your error message due to the fact that you want a vector element, i.e.

element = basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(2, ))

Usually users do not manually generate the C files to call tabulate_tensor_integral. It is usually achieved through dolfinx.fem.form or dolfinx.fem.create_form, as shown in:
Create domain independent forms and call dolfinx.fem.create_form

**Create problem specific forms and call dolfinx.fem.form:

1 Like
7

Thank you for the clarification! The links are very helpful and I can print the matrices. For the domain independent form, I was wondering how I reduce the block size to 1. Or is this not possible in that case?

MWE (the commented line would raise a segmentation fault)

import basix.ufl
import dolfinx
import ufl
from mpi4py import MPI

element = basix.ufl.element("Lagrange", "quadrilateral", 1, shape=(2,))
domain = ufl.Mesh(element)
V = ufl.FunctionSpace(domain, element)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, cell_type=dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.functionspace(mesh,element)

# V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
a_compiled = dolfinx.fem.compile_form(MPI.COMM_WORLD, a)


bilinear_form = dolfinx.fem.create_form(a_compiled, [V, V], mesh, {}, {}, {})

A = dolfinx.fem.create_matrix(bilinear_form)
A_scipy = A.to_scipy()
print(f"{A_scipy=}")
dolfinx.fem.assemble_matrix(A, bilinear_form, bcs=[])
print(A_scipy.todense())
8

Then you need to change this to a scalar element, ie the element here should be a blocked element.

1 Like
9

Got it, thanks again!