augustjohansson / customquad Goto Github PK

View Code? Open in Web Editor NEW

4.0 4.0 2.0 241 KB

Library for using custom quadrature rules at runtime in FEniCSx

License: GNU Lesser General Public License v3.0

Python 90.66% C++ 6.07% Dockerfile 1.19% Shell 2.08%

customquad's People

Contributors

Stargazers

Watchers

Forkers

dbaroli daveb-dev

customquad's Issues

Add tests for algoim qr

Eg qr_pts length match cut_cells, but xyz not

Only double is supported

At least complex should be supported

Running demo_elasticity.py doesn't work

RuntimeError: Function element interpolation points has different shape to Expression interpolation points

We only handle first order derivatives in the basix call

This should be deduced from the form

How to use a c++ compiler when compiling call_basix?

It also looks like we need to set the -fpermissive flag for the cffi to pass.

Support simplices

What is perm in the kernel call?

It's currently ignored

Multiple new calls in call_basix without delete

Small memory leak

Replace FIAT quadrature calls with basix

cell_type = basix.CellType.quadrilateral
points, w = basix.make_quadrature(basix.QuadratureType.Default, cell_type, degree)

Locking the inactive dofs is slow

We need to allocate the matrix properly

Parallel

Make sure the code works in parallel

Create dockerfile with algoim

Make sure vector elements work

Only linear lagrange elements on quads and hexes are supported

Why does the quadrilateral quadrature points need to be mirrored in y=x?

Test code for integrating f = 2*x+y using (0.5, 0.1) as quadrature rule. The exact integral should be 2*0.5+0.1 = 1.1, which is obtained if we mirror the quadrature points.

NB: we do this mirroring in call_basix.

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl
import customquad as cq

cell_type = dolfinx.mesh.CellType.quadrilateral
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, cell_type)
cells = [0]

x = ufl.SpatialCoordinate(mesh)
integrand = 2 * x[0] + x[1]
form = dolfinx.fem.form(
    integrand * ufl.dx(metadata={"quadrature_rule": "runtime"}, domain=mesh)
)

qr_pts = np.array([[0.5, 0.1]])
qr_w = np.array([[1.0]])
m_exact = 1.1 # 2*0.5 + 0.1

qr_data = [(cells, qr_pts, qr_w)]
m = cq.assemble_scalar(form, qr_data)

# Mirror qr in y=x
qr_pts2 = qr_pts.copy()
qr_pts2[:, [1, 0]] = qr_pts2[:, [0, 1]]
qr_data2 = [(cells, qr_pts2, qr_w)]
m2 = cq.assemble_scalar(form, qr_data2)

print(qr_pts, m, m2, m_exact)

This gives [[0.5 0.1]] 0.6999999999999997 1.0999999999999996 1.1.