FEniCS / dolfinx

Next generation FEniCS problem solving environment
https://fenicsproject.org
GNU Lesser General Public License v3.0
731 stars 177 forks source link

[BUG]: can't integrate the inner product of the gradient of a 3D vector and itself over the domain using dolfinx.fem.assemble #2659

Closed ruf10 closed 1 year ago

ruf10 commented 1 year ago

How to reproduce the bug

I try to calculate the volume of my domain/mesh, and fem.assemble_scalar is not working, also I try to calculate the integral of nabla_grad of a vector inner product with itself over the domain, its also returns me the error: Can anyone help me with that quickly?

Thank you!

Best, RUi

Minimal Example (Python)

F = (1./dt)*inner(u,v)*dx + eps*p*q*dx +b(un,u,v)*dx +2.*nu*a_sym(u,v)*dx
        F += - c(p,v)*dx+ c(q,u)*dx +nu_t*a_sym(u,v)*dx
        F -= (1./dt)*(inner(un,v))*dx
        unPlus1, pnPlus = solve_u(F, t)

        def grad_sym(u):
            return 0.5 * (nabla_grad(u) + nabla_grad(u).T)

        def frobenius_norm_squared(T):
            return ufl.inner(T, T)

        u_Sym = grad_sym(unPlus1)
        frobenius_norm_squared_u_Sym = frobenius_norm_squared(u_Sym)

        # Define measure for integration over the domain
        dx = ufl.dx(domain=msh)

        # Define form for integral
        M_form = frobenius_norm_squared_u_Sym * dx

        # Assemble integral over the domain
        M = dolfinx.fem.assemble.assemble_scalar(M_form)

Output (Python)

Traceback (most recent call last):
  File "/Users/ruifang/ct.py", line 481, in <module>
    M = dolfinx.fem.assemble.assemble_scalar(M_form)
  File "/Users/ruifang/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/assemble.py", line 125, in assemble_scalar
    constants = constants or _pack_constants(M)
TypeError: pack_constants(): incompatible function arguments. The following argument types are supported:
    1. (form: dolfinx::fem::Form<float>) -> numpy.ndarray[numpy.float32]
    2. (e: dolfinx::fem::Expression<float>) -> numpy.ndarray[numpy.float32]
    3. (form: dolfinx::fem::Form<double>) -> numpy.ndarray[numpy.float64]
    4. (e: dolfinx::fem::Expression<double>) -> numpy.ndarray[numpy.float64]
    5. (form: dolfinx::fem::Form<std::__1::complex<float> >) -> numpy.ndarray[numpy.complex64]
    6. (e: dolfinx::fem::Expression<std::__1::complex<float> >) -> numpy.ndarray[numpy.complex64]
    7. (form: dolfinx::fem::Form<std::__1::complex<double> >) -> numpy.ndarray[numpy.complex128]
    8. (e: dolfinx::fem::Expression<std::__1::complex<double> >) -> numpy.ndarray[numpy.complex128]

Invoked with: Form([Integral(Inner(ComponentTensor(Product(FloatValue(0.5), Indexed(Sum(NablaGrad(Coefficient(FunctionSpace(Mesh(VectorElement(Basix element (P, tetrahedron, 1, equispaced, unset, False), 3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3)), 115)), Transposed(NablaGrad(Coefficient(FunctionSpace(Mesh(VectorElement(Basix element (P, tetrahedron, 1, equispaced, unset, False), 3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3)), 115)))), MultiIndex((Index(3700), Index(3701))))), MultiIndex((Index(3700), Index(3701)))), ComponentTensor(Product(FloatValue(0.5), Indexed(Sum(NablaGrad(Coefficient(FunctionSpace(Mesh(VectorElement(Basix element (P, tetrahedron, 1, equispaced, unset, False), 3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3)), 115)), Transposed(NablaGrad(Coefficient(FunctionSpace(Mesh(VectorElement(Basix element (P, tetrahedron, 1, equispaced, unset, False), 3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3)), 115)))), MultiIndex((Index(3700), Index(3701))))), MultiIndex((Index(3700), Index(3701))))), 'cell', Mesh(VectorElement(Basix element (P, tetrahedron, 1, equispaced, unset, False), 3), 0), 'everywhere', {}, None)])

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

IgorBaratta commented 1 year ago

Could you try calling form first, eg:

from dolfinx.fem import form
....
M_form = form(M_form) # JIT
M = dolfinx.fem.assemble.assemble_scalar(M_form)
ruf10 commented 1 year ago

thank you so much!!! it works!