TypeError: BoundingBoxTree.init() takes 2 positional arguments but 3 were given

Karlosjusus · November 22, 2023, 12:04pm

Hi everbody!
I’m trying to load a .msh mesh into dolfinxa and call BoundingBoxTree

I get the error

BoundingBoxTree.__init__() takes 2 positional arguments but 3 were given

my code is:

import dolfinx
import gmsh
from mpi4py import MPI

from dolfinx.io import gmshio


domain, cell_tags, ft = gmshio.read_from_msh("msh/venturi_beta08_40.msh", MPI.COMM_WORLD, 0, gdim=2)
ft.name = "Facet markers"

tdim = domain.topology.dim
tree = dolfinx.geometry.BoundingBoxTree(domain, tdim)

Have anyone a idea?
thank you so much

dparsons · November 22, 2023, 12:09pm

dolfinx 0.7 got a change in API.
cf. Fix bounding boxes and more by jorgensd · Pull Request #135 · jorgensd/dolfinx-tutorial · GitHub

Karlosjusus · November 22, 2023, 1:48pm

Thank you @dparsons!
But I have several problems now. My code worked before the changes. But know I have problem with

-write.meshtag

TypeError: XDMFFile.write_meshtags() missing 1 required positional argument: 'x'

-write.function

Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?

Do you know what I have to change?

Thank you again!

here is the full code:

import dolfinx
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
import pyvista
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, create_matrix, set_bc)
import matplotlib.pyplot as plt

from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)

from dolfinx.mesh import create_unit_square
from petsc4py.PETSc import ScalarType

from dolfinx import geometry
from dolfinx.io import gmshio
import meshio
import dolfinx.cpp as _cpp
import dolfinx.fem as _fem

###############################################################################

    #Mesh with GMsh
    
################################################################################
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0

gmsh.initialize()
     
    #Points
n_point = 50
str_point = str(n_point)
Re = 500
str_re = str(Re)

##################################################################################

    # Coordinates
    
##################################################################################


cx = [0,0.225,0.252,0.292,0.3737,0.9450,0.9450,0.3737,0.292,0.252,0.225,0]
cy = [0.025, 0.025,0.02,0.02,0.025,0.025,-0.025, -0.025,-0.02,-0.02,-0.025,-0.025]


bump = 0.05

Toni = False

a = 2.0*(bump-1)/(n_point*1.3)+1 

n_bump = int(math.log(bump,a))
c = 1e-1;
gmsh.model.geo.add_point(0, 0.025, 0, c,1)
gmsh.model.geo.add_point(0.225, 0.025, 0, c,2)
gmsh.model.geo.add_point(0.252, 0.02, 0, c,3)
gmsh.model.geo.add_point(0.292, 0.02, 0, c, 4)
gmsh.model.geo.add_point(0.3737, 0.025, 0, c, 5)
gmsh.model.geo.add_point(0.9450, 0.025, 0, c, 6)
gmsh.model.geo.add_point(0.9450, -0.025, 0, c, 7)
gmsh.model.geo.add_point(0.3737, -0.025, 0, c,8)
gmsh.model.geo.add_point(0.292, -0.02, 0, c, 9)
gmsh.model.geo.add_point(0.252, -0.02, 0, c, 10)
gmsh.model.geo.add_point(0.225, -0.025, 0, c, 11)
gmsh.model.geo.add_point(0, -0.025, 0, c, 12)
gmsh.model.occ.synchronize()
    #Lines
    
gmsh.model.geo.addLine(1, 2 ,1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4,3)
gmsh.model.geo.addLine(4, 5,4)
gmsh.model.geo.addLine(5, 6,5)
gmsh.model.geo.addLine(6, 7,6)
gmsh.model.geo.addLine(7, 8,7)
gmsh.model.geo.addLine(8, 9,8)
gmsh.model.geo.addLine(9, 10,9)
gmsh.model.geo.addLine(10, 11,10)
gmsh.model.geo.addLine(11, 12,11)
gmsh.model.geo.addLine(12, 1,12)
gmsh.model.occ.synchronize()

gmsh.model.geo.add_curve_loop([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],1)
gmsh.model.geo.add_plane_surface([1],1)
gmsh.model.geo.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)


if Toni == False:
    gmsh.model.addPhysicalGroup(1, [12], 1)
    gmsh.model.addPhysicalGroup(1, [6], name = 'inlet')
    gmsh.model.addPhysicalGroup(1, [1,2,3,4,5,7,8,9,10,11], name= 'walls')
    gmsh.model.addPhysicalGroup(2, [1], 1)
    gmsh.model.occ.synchronize()


gmsh.model.geo.mesh.setTransfiniteCurve(6, n_bump, "Bump", bump)
gmsh.model.geo.mesh.setTransfiniteCurve(12, n_bump, "Bump", bump)
gmsh.model.geo.mesh.setTransfiniteCurve(1, int(5.2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(2, int(0.8*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(3, int(1*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(4, int(2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(5, int(10*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(7, int(10*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(8, int(2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(9, int(1*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(10, int(0.8*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(11, int(5.2*n_point))

gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1,6,7,12])


gmsh.model.geo.mesh.setRecombine(2, 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)



domain, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
ft.name = "Facet markers"
gmsh.finalize()

##########################################################################################################################################

    # Define function 

def closest_point(mesh, point,tol):
    points = point
    while True:
        try:
            entity = dolfinx.geometry.compute_closest_entity(tree, mid_tree, mesh, points)
            break
        except RuntimeError:
            print(points)
            for j in range(len(mesh.geometry.x)):
                if (np.abs(mesh_geom[j][0]-points[0])< tol and np.abs(mesh_geom[j][1]-points[1])< tol):
                     return points, j              
    
    geom_dof = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object,tdim,entity, False)
    #geom_dof = dolfinx.cpp.mesh.compute_incident_entities(mesh,[entity], 2,1)
    #print('geom',geom_dof)
    mesh_nodes = mesh_geom[geom_dof][0]
    #print('mesh_nodes', mesh_nodes)
    dis = dolfinx.geometry.compute_distance_gjk(points, mesh_nodes)
    index = -100
    for i in range(len(mesh_nodes)):
        #print(mesh_nodes[i])
        if (np.abs(mesh_nodes[i][0]-points[0]+ dis[0]) < tol and np.abs(mesh_nodes[i][1]-points[1]+ dis[1]) < tol):
            index = i
            
    if (index == -100):
        return point, index
    else:
        return points[0] - dis[0], points[1] - dis[1] , geom_dof[0][index]
    
def facet_normal_approximation(V, mt, mt_id, tangent=False): 
                               #jit_options: dict = {}, form_compiler_options: dict = {}):

    comm = V.mesh.comm
    n = ufl.FacetNormal(V.mesh)
    nh = _fem.Function(V)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    ds = ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
    if tangent:
        if V.mesh.geometry.dim == 1:
            raise ValueError("Tangent not defined for 1D problem")
        elif V.mesh.geometry.dim == 2:
            a = ufl.inner(u, v) * ds
            L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
        else:
            def tangential_proj(u, n):
                """
                See for instance:
                https://link.springer.com/content/pdf/10.1023/A:1022235512626.pdf
                """
                return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
            c = fem.Constant(V.mesh, [1, 1, 1])
            a = ufl.inner(u, v) * ds
            L = ufl.inner(tangential_proj(c, n), v) * ds
    else:
        a = (ufl.inner(u, v) * ds)
        L = ufl.inner(n, v) * ds

    # Find all dofs that are not boundary dofs
    imap = V.dofmap.index_map
    all_blocks = np.arange(imap.size_local, dtype=np.int32)
    top_blocks = _fem.locate_dofs_topological(V, V.mesh.topology.dim - 1, mt.find(mt_id))
    deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]

    # Note there should be a better way to do this
    # Create sparsity pattern only for constraint + bc
    bilinear_form = _fem.form(a)
    pattern = _fem.create_sparsity_pattern(bilinear_form)
    pattern.insert_diagonal(deac_blocks)
    pattern.finalize()
    u_0 = _fem.Function(V)
    u_0.vector.set(0)

    bc_deac = _fem.dirichletbc(u_0, deac_blocks)
    A = _cpp.la.petsc.create_matrix(comm, pattern)
    A.zeroEntries()

    # Assemble the matrix with all entries
    form_coeffs = dolfinx.cpp.fem.pack_coefficients(bilinear_form._cpp_object)
    form_consts = _cpp.fem.pack_constants(bilinear_form._cpp_object)
    _cpp.fem.petsc.assemble_matrix(A, bilinear_form._cpp_object,
                                   form_consts, form_coeffs, [bc_deac._cpp_object])
    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore
        _cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac._cpp_object], 1.0)
    A.assemble()
    linear_form = _fem.form(L)
    b = _fem.petsc.assemble_vector(linear_form)

    _fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)  # type: ignore
    _fem.petsc.set_bc(b, [bc_deac])

    # Solve Linear problem
    solver = PETSc.KSP().create(V.mesh.comm)  # type: ignore
    solver.setType("cg")
    solver.rtol = 1e-8
    solver.setOperators(A)
    solver.solve(b, nh.vector)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)  # type: ignore
    return nh


########################################################################################################

    # Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
v_cg2 = VectorElement("CG", domain.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg2, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

fdim = domain.topology.dim -1

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
def u_exact(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 1
    return values

def p_exact(x):
    values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 0
    return values

V1, _ = W.sub(0).collapse()
Q1, _ = W.sub(1).collapse()
u_icondition = Function(V1)
u_wallcondition = Function(V1)
p_icondition = Function(Q1)
u_icondition.interpolate(u_exact)
p_icondition.interpolate(p_exact)

up_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological((W.sub(0),V1), fdim, wall_facets)
down_f = dolfinx.fem.locate_dofs_topological((W.sub(1),Q1), fdim, down_facets)

mu = Constant(domain, PETSc.ScalarType(1.0/Re))
print(mu)
n = FacetNormal(domain)   
nh = facet_normal_approximation(V1, facet_tag, 3, tangent=False)

with dolfinx.io.XDMFFile(domain.comm, "Plot/Functions/normal_wall.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(nh)

dokken · November 22, 2023, 2:29pm

Meshtags now require the geometry as second input

The decision of having to explicitly interpolate into a supported XDMF-file function space (will be a continuous lagrange space of same order as the mesh geometry) was made to highlight that this makes a loss of accuracy and can be an expensive operation.

If you are using higher order Lagrange/DG elements, you should use VTKFile or VTXWriter, as they support this kind of output. If you are using other spaces, such as RT, Nedelec etc, they should be interpolated into an appropriate DG space and use VTXWrtier or VTKFile for outputting.

Topic		Replies	Views
Bounding box tree error Errors	1	20	October 26, 2024
TypeError: FiniteElement.__init__() missing 3 required positional arguments: 'reference_value_shape', 'pullback', and 'sobolev_space' General	3	133	May 27, 2024
Meshtags_from_entities(), incompatible function arguments dolfinx mesh , dolfinx	2	53	July 4, 2024
Error after gmshio.model_to_mesh dolfinx	8	954	September 25, 2023
TypeError: __init__(): incompatible function arguments. [...] dolfinx mesh	7	71	September 25, 2024

TypeError: BoundingBoxTree.__init__() takes 2 positional arguments but 3 were given

Related topics

TypeError: BoundingBoxTree.init() takes 2 positional arguments but 3 were given