Transform dolfinx assemble matrix to numpy array

dokken, I’ll show you the picture after open “tag_triangle.xdmf”. The error should happened when transforming mesh file to xdmf file.

This is clearly a issue with either:
a) Pygmsh
b) Gmsh
c) Meshio
d) Paraview
as it works for me with:

  • pygmsh 6.0.8
  • gmsh 4.5.1
  • meshio 4.0.9
  • paraview version 5.7.0

dokken, it’s the following code questionable, but I cannot figure it out.

import pygmsh
import meshio
import numpy as np

filename = 'strip.msh'
mesh = meshio.read(
    filename,  # string, os.PathLike, or a buffer/open file
)
points, cells, cell_data, field_data = mesh.points, mesh.cells, mesh.cell_data, mesh.field_data

triangle_data = None
line_data = None
tri_cell_idx = []
line_cell_idx = []

for i, cell in enumerate(cells):
    if cell.type == 'triangle':
        print(i, 'triangle')
        tri_cell_idx.append(i)
        if triangle_data is None:
            triangle_data = cell.data
        else:
            triangle_data = np.vstack([triangle_data, cell.data])
    else:
        print(i, 'line')
        line_cell_idx.append(i)
        if line_data is None:
            line_data = cell.data
        else:
            line_data = np.vstack([line_data, cell.data])


meshio.xdmf.write('tag_triangle.xdmf', meshio.Mesh(
    points=points,
    cells={'triangle': triangle_data},
    cell_data={
        'triangle': [np.concatenate([cell_data['gmsh:physical'][i] for i in tri_cell_idx])]
    },
    field_data=field_data
))
meshio.xdmf.write("tag_all.xdmf", meshio.Mesh(
    points=points,
    cells={"line": line_data},
    cell_data={
        "line": [np.concatenate([cell_data['gmsh:physical'][i] for i in line_cell_idx])]
    }
))

dokken, can you check the physical boundary of the mesh is correct for your software version? You can run the below code. The first value should be 3.680109307299481.

import numpy as np
import pygmsh
import meshio
from dolfinx.io import XDMFFile
from dolfinx import (MPI, FacetNormal, Function, FunctionSpace, 
                     has_petsc_complex, DirichletBC, cpp)
from dolfinx.fem.assemble import assemble_scalar, assemble_matrix
from dolfinx.fem import locate_dofs_topological
from ufl import (TestFunctions, TrialFunctions, grad, inner, curl, split,
                  Measure, FiniteElement, dx, Dx, dot)
import dolfinx
from scipy.constants import c
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs
from slepc4py import SLEPc
from petsc4py import PETSc

with XDMFFile(MPI.comm_world, "tag_triangle.xdmf") as xdmf_infile:
    mesh = xdmf_infile.read_mesh(cpp.mesh.GhostMode.none)
    mvc_subdomain = xdmf_infile.read_mvc_size_t(mesh)
    # tag_info = xdmf_infile.read_information_int()

with XDMFFile(MPI.comm_world, "tag_all.xdmf") as xdmf_infile:
    mvc_boundaries = xdmf_infile.read_mvc_size_t(mesh)

mf_triangle = cpp.mesh.MeshFunctionSizet(mesh, mvc_subdomain, 0)
mf_line = cpp.mesh.MeshFunctionSizet(mesh, mvc_boundaries, 0)

V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)

S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()

for i in range(x.shape[0]):
    if mf_triangle.values[i] == 2:
        e_r.vector.setValueLocal(i, 8.875)
    else:
        e_r.vector.setValueLocal(i, 1)

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

f = 25e9
k0_sqr = (2*np.pi*f/c)**2

electric_wall = Function(W)

with electric_wall.vector.localForm() as bc_local:
    bc_local.set(0)

bndry_facets = np.where(mf_line.values == 1)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc= DirichletBC(electric_wall, bdofs)

theta_sqr = k0_sqr * 8.875

a = 1/theta_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx 
b = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
c = b + a 

S = assemble_matrix(b, [bc])
T = assemble_matrix(c, [bc])

S.assemble()
S = S.realPart()

T.assemble()
T = T.realPart()

es = SLEPc.EPS(which='LR').create(comm=MPI.comm_world)
es.setDimensions(100)
es.setOperators(S, T)
es.setFromOptions()
es.solve()

vr, vi = S.getVecs()
lam = es.getEigenpair(0, vr, vi)
print(lam)

The script gives me:

(35.7681974852247-1.2044409474322886e-12j)

But as I said. This is clearly a meshing issue. I suggest either to raise an issue at pygmsh or meshio’s github. But then you need to figure out what is wrong with your mesh.

As I said in the previous post, to do so, write the working and failing mesh to file using: meshio.xdmf.write("filename.xdmf", mesh, data_format="XML") . Then you can compare to similar meshes in plaintext (use coarse meshes).

ok, I’ll try. Thank you.