How to define trial and test functions for mixed problems

General
1

I want to define the trial and test functions for a mixed problem:

Blockquote
from ufl import TrialFunction, TestFunction, dot, div, inner ,dx
from dolfinx import mesh, fem
from mpi4py import MPI
import numpy as np
from ufl import (VectorElement, FiniteElement, MixedElement )
Lx,Ly = 1,1.1
xL=[0,0]
xR=[Lx,Ly]
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[4, 4], mesh.CellType.triangle)
RTelement = FiniteElement(“RT”,domain.ufl_cell(), 1)
DGelement = FiniteElement(“DG”,domain.ufl_cell(), 0)
Mixed_element = MixedElement([RTelement, DGelement])
MixedSpace = fem.FunctionSpace(domain, Mixed_element)
‘’‘Trial and test functions’‘’
sigma, u = TrialFunction(MixedSpace)
tau, v = TestFunction(MixedSpace)

I get the following error:
sigma, u = TrialFunction(MixedSpace)
ValueError: too many values to unpack (expected 2)

2

use TrialFunctions and TestFunctions

3

Thank you. It works fine now.

4

One more question please:
How can I turn PETScMatrix to a numpy array?

5

Please consider reading @dokken solution on this post Conversion from array to PETSc and from PETSc to array.

1 Like
6

I used the mentioned code and I get :

‘petsc4py.PETSc.Mat’ object has no attribute ‘shape’. Did you mean: ‘scale’?

7

I have found the solution. For others to benafit (With fenicsx, I don’t know if it works with fenics) if you have a petsc matrix B and you want to convert it to numpy, just use

convertToDense(B))

8

Note that this only works for small problems, see Conversion from array to PETSc and from PETSc to array - #4 by dokken

9

Thank you.
But I encountered now a new problem:
when I change the domain to quads I get the following:
Error:
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cellname “quadrilateral” invalid for “Raviart-Thomas” finite element.

This is the code:

Blockquote
from dolfinx.mesh import create_unit_square, create_box, CellType
from dolfinx import mesh, plot, fem
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
import pyvista # visualizing the mesh using pyvista, an interface to the VTK toolkit.
print(pyvista.global_theme.jupyter_backend)
from ufl import (VectorElement, FiniteElement, MixedElement )
###############################################################
Lx,Ly = 1,1.1
xL=[0,0]
xR=[Lx,Ly]
num_elemnts = 4
RT_order = 1
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.quadrilateral)
‘’’
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.triangle)
#############################################################
print(‘Constructing the mesh …’)
tdim = domain.topology.dim
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()
print(‘done…’)
##############################################################
‘’‘The function spaces’‘’
print(‘Constructing the spaces …’)
RTelement = FiniteElement(“RT”,domain.ufl_cell(), RT_order)
DGelement = FiniteElement(“DG”,domain.ufl_cell(), RT_order-1)
Mixed_element = MixedElement([RTelement, DGelement])
MixedSpace = fem.FunctionSpace(domain, Mixed_element)

10

RT is not defined on quadrilaterals DefElement: Raviart–Thomas
You Need to use RTCF DefElement: Q H(div)

11

There is something that I don’t get or I’m missing.
I have previously defined it on a different problem with traig/quads elements and I get the correct results:
The code runs for both domains:

Blockquote
from dolfinx.mesh import create_unit_square, create_box, CellType
from dolfinx import mesh, plot, fem
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
import pyvista # visualizing the mesh using pyvista, an interface to the VTK toolkit.
print(pyvista.global_theme.jupyter_backend)
from ufl import (VectorElement, FiniteElement, MixedElement )
###############################################################
Lx,Ly = 1,1.1
xL=[0,0]
xR=[Lx,Ly]
num_elemnts = 4
RT_order = 1
##############################################################
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.quadrilateral)
‘’’
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.triangle)
‘’’
#############################################################
print(‘Constructing the mesh …’)
tdim = domain.topology.dim
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()
print(‘done…’)
##############################################################
‘’‘The function spaces’‘’
print(‘Constructing the spaces …’)
space = fem.FunctionSpace(domain, (“RT”, 1))

12

Please don’t use blockquotes. Instead use 3x`, ie

```python 
import dolfinx
#add the rest of your code here 
```

Please also add any error messages that you get with the following encapsulation

```bash
Add error message here
```
13

In the below code, I use the RT element with quads and it works fine. But in the previous code for the mixed space, I had to run it with RTCF with quads. I would like to understand why does it work in the below code and not with the mixed case?

from dolfinx.mesh import create_unit_square, create_box, CellType
from dolfinx import mesh, plot, fem
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
import pyvista # visualizing the mesh using pyvista, an interface to the VTK toolkit.
print(pyvista.global_theme.jupyter_backend)
from ufl import (VectorElement, FiniteElement, MixedElement )
###############################################################
Lx,Ly = 1,1.1
xL=[0,0]
xR=[Lx,Ly]
num_elemnts = 4
RT_order = 1
##############################################################
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.quadrilateral)
‘’’
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array(xL), np.array(xR)],
[num_elemnts, num_elemnts], mesh.CellType.triangle)
‘’’
#############################################################
print(‘Constructing the mesh …’)
tdim = domain.topology.dim
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()
print(‘done…’)
##############################################################
‘’‘The function spaces’‘’
print(‘Constructing the spaces …’)
space = fem.FunctionSpace(domain, (“RT”, 1))
14

That is because when you create the first function space, you are not creating a finite element explicitly, but implicitly with a tuple

This calls basix directly (https://github.com/FEniCS/dolfinx/blob/v0.6.0/python/dolfinx/fem/function.py#L481C21-L484), which allows you to use RT for quads (https://github.com/FEniCS/basix/blob/v0.6.0/python/basix/ufl_wrapper.py#L1251, https://github.com/FEniCS/basix/blob/0fc130f7d790f793ffa57bf8d056232ff6f12be8/python/basix/finite_element.py#L27).

ufl.FiniteElement is stricter in this sense, and requires you to use RTCF (https://github.com/FEniCS/ufl/blob/7b24065a5f39ab4296c4ab8246f2b493d7bcd803/ufl/finiteelement/elementlist.py#L105-L106, https://github.com/FEniCS/ufl/blob/7b24065a5f39ab4296c4ab8246f2b493d7bcd803/ufl/finiteelement/elementlist.py#L182-L183)