添加链接
link管理
链接快照平台
  • 输入网页链接,自动生成快照
  • 标签化管理网页链接

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

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:
            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)
 Karlosjusus:

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

Meshtags now require the geometry as second input

Karlosjusus:

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

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.