wells-group / asimov-custom-assemblers Goto Github PK

Assembly of matrices using dof transformations instead of permutations is expensive, and should be avoided at all costs.

This is an optimization after #12 is merged

Suggestion: optional input parameter quadrature_degree for function create_distance_map. If quadrature degree unchanged, check if _qp_ref_facet and _phi_ref_facets already computed (these are the same independent of origin_meshtag).

Add run ffcx to cmake of demos

Wrapping the basix::make_quadrature in a nice format.

All kernels should then take in the quadrature rule, instead of the quadrature degree.
Avoids alot of recomputation and mixed inputs.

Python tests started failing on 15/05/2021

Support either pulling back dolfin-x functions, or pulling back lambda functions

Adapt to array and use xtensor allclose

I am getting a compilation error for the Python module with the latest DolfinX. This worked without problems last week so I think it is due to the refactoring of the pack coefficients that was introduced a few days ago:

FEniCS/dolfinx#1966

This is with gcc 9.3.0 on CentOS 8. The relevant error is below. Any suggestions would be greatly appreciated!

dolfinx_cuas/include/dolfinx_cuas/utils.hpp:212:71: error: no matching function for call to 'pack_coefficient_entity<double, std::tuple<int, int, int, int> >(tcb::span, int, __gnu_cxx::__alloc_traits<std::allocator<tcb::span >, tcb::span >::value_type&, tcb::span&, const dolfinx::fem::DofMap&, tcb::span<const std::tuple<int, int, int, int> >&, dolfinx_cuas::pack_coefficients(std::vector<std::shared_ptr<const dolfinx::fem::Function > >, std::variant<tcb::span<const int, -1>, tcb::span<const std::pair<int, int>, -1>, tcb::span<const std::tuple<int, int, int, int>, -1> >) [with T = double]::<lambda(auto:164&&)> [with auto:164 = tcb::span<const std::tuple<int, int, int, int> >&]::<lambda(const std::tuple<int, int, int, int>&)>&, __gnu_cxx::__alloc_traits<std::allocator, int>::value_type, int, const std::function<void(const tcb::span&, const tcb::span&, int, int)>&)'
211 | dolfinx::fem::impl::pack_coefficient_entity<
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212 | T, std::tuple<std::int32_t, int, std::int32_t, int>>(
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
213 | xtl::span(c), 2 * cstride, v[coeff], cell_info, *dofmaps[coeff], entities,
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
214 | fetch_cell1, coeffs_offsets[coeff] + coeffs_offsets[coeff + 1],
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
215 | elements[coeff]->space_dimension(), transform);

These functions are duplicates of
https://github.com/FEniCS/dolfinx/blob/3c6c1ce9a1956dcb6b7066d438ec5013e47fb733/cpp/dolfinx/fem/CoordinateElement.h#L99-L117
and can be used with

dolfinx::fem::CoordinateElement::compute_....

which would require less maintenance from our end.

Any objections? @IgorBaratta @SarahRo

• Scalar Mass Matrix on Quads
• Scalar Mass Matrix on Hexes
• Vector Mass Matrix on Hexes
• Scalar Stiffness Matrix on Quads
• Scalar Stiffness Matrix on Hexes
• Vector Stiffness Matrix on Hexes

Tensor product Gauss-Legendre (GL) with q=p+2, p is the polynomial order of the element.
Use nodal basis with p+1 Gauss-Legendre-Lobatto (GLL) points in each spatial dimension.

Currently depends on a form to get alot of data. As these are custom kernels, it should only depend on function space + coeffs and constants.

Interface should thus take in function space, and potential coeffs and constants to be used in the kernel.

Required for gradients of surface integrals

Currently we time the whole assembly, including creating the matrix (and sparistypattern).
This is not optimal, as the sparsity pattern might be alot smaller if assembling over set of sub entities (for instance exterior facets).

If we want to have accurate comparisons, we need to change the signature of the assembly function to have the matrix as input (dolfin-x already has this option).

To generalize the assembler, we should allow it to integrate over a set of sub entities (similar to having a ufl measure with meshtags in dolfinx). This is currently implemented for the surface assembler, and should be easy to generalize to the cell assemblers.

• Mass Matrix - Lagrange - a = inner(u, v)*dx
• Mass Matrix - Vector Lagrange - a = inner(u, v)*dx
• Mass Matrix - Nédélec - a = inner(u, v)*dx
• Stiffness Matrix - Lagrange - a = inner(grad(u), grad(v))*dx
• Stiffness Matrix - Vector Lagrange
• Stiffness Matrix - Nédélec - a = inner(curl(u), curl(v))*dx
• Surface Integration - Vector Lagrange - a = inner(u, v*n)*dx

Currently we can use dolfinx with complex petsc, as the python installation of the header only library fails to detect that the build is complex (Similar problem in the dolfin-x CI): https://github.com/FEniCS/dolfinx/runs/3929272622
In Build Python interface:

...
     -- Looking for PETSC_USE_COMPLEX
    -- Looking for PETSC_USE_COMPLEX - not found

while
Configure C++ yields

-- Looking for PETSC_USE_COMPLEX
-- Looking for PETSC_USE_COMPLEX - found

The implementation for non-templated code should be in a .cpp file rather than header files.

Would avoid copying of basis functions, coordinates etc, as we could use xt::view (not allowed by xt::linalg::dot).

Challenging task that requires a few steps.

• Find local facet index
• Get quadrature points on reference facet
• Push quadrature points from reference facet to reference cell
• Pull back data from physical element to these quadrature points