I ran the real version of Dolfin. The newest docker images can be found at:
https://hub.docker.com/u/dolfinx
I get the same problem when trying to run your code with dolfinx/complex.
@wwshunan I’ve been able to reduce your example to a way smaller code:
import numpy as np
import pygmsh
import meshio
from dolfinx.io import XDMFFile
from dolfinx import (MPI, FunctionSpace,
cpp, UnitSquareMesh)
from dolfinx.fem.assemble import assemble_matrix
from dolfinx.fem import locate_dofs_topological
from ufl import (TestFunctions, TrialFunctions, inner,
FiniteElement, dx, TrialFunction, TestFunction)
import dolfinx
mesh = UnitSquareMesh(MPI.comm_world, 35, 35)
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
u,p = TrialFunctions(W)
v,q = TestFunctions(W)
a = inner(u, v)*dx + inner(p, q)*dx
A = assemble_matrix(a)
A.assemble()
print(W.dim())
C = A.convert('dense')
print(C.getDenseArray())
This works code runs for a 30 x 30 UnitSquareMesh, but not a 35x35 UnitSquareMesh, indicating that there are memory issues when converting large systems (The matrix will be ~40 000x40 000). Which is a crazy large matrix to store in python.
Thank you, dokken. Is it then impossible to transform the matrix to numpy array in case dolfinx now don’t provide SLEPcEigenSolver?
The PETScEigen solver is still in dolfinx. It is just not exposed to the Python layer through pybind. This is something that you could do yourself, or make an issue on the dolfinx issue tracker and hope that someone does it for you.
dokken, thank you. It’s a bit hard for me to bind to python and you said dolfinx/real have no problem. I’ll try it.
It’s probably the dense mode that’s killing you. You can convert to a sparse CSR matrix using something like
A = assemble_matrix(a, [bc])
A.assemble()
assert isinstance(A, PETSc.Mat)
ai, aj, av = A.getValuesCSR()
Asp = csr_matrix((av, aj, ai))I also found the silimar solution on the internet, but not aware of what it means. Thank you.
CSR means compressed-sparse-row. So it’s just storing the column index, row index and the value. I’m sure you can find some scipy examples online that will help with your understanding if required.
OK, I see. 
Hi everyone, I generated the mesh using the gmsh software tool and used meshio to create the xdmf files. We can see the generated physical boundaries are not correct. What’s the problem in the transform file? Thank you.
You Need to supply us with a Minimal Working Example (MWE) for anyone to answer this question.
Dear dokken, sorry for not pointing out the example was at the top of this post. I’ll repeat the code with problem for convenience.
gmsh geo file:
//+
d_w = DefineNumber[ 0.0127, Name "Parameters/d_w" ];
//+
d_h = DefineNumber[ 1.27e-3, Name "Parameters/d_h" ];
//+
s_w = DefineNumber[ 1.27e-3, Name "Parameters/s_w" ];
//+
s_h = DefineNumber[ 1.27e-4, Name "Parameters/s_h" ];
//+
Point(1) = {0, 0, 0, 0.0005};
//+
Point(2) = {d_w, 0, 0, 0.0005};
//+
Point(3) = {d_w, d_h, 0, 0.0005};
//+
Point(4) = {d_w, d_w, 0, 0.0005};
//+
Point(5) = {0, d_w, 0, 0.0005};
//+
Point(6) = {0, d_h, 0, 0.0005};
//+
Point(7) = {d_w/2-s_w/2, d_h, 0, 0.00005};
//+
Point(8) = {d_w/2+s_w/2, d_h, 0, 0.00005};
//+
Point(9) = {d_w/2+s_w/2, d_h+s_h, 0, 0.00005};
//+
Point(10) = {d_w/2-s_w/2, d_h+s_h, 0, 0.00005};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {6, 1};
//+
Line(7) = {6, 7};
//+
Line(8) = {7, 8};
//+
Line(9) = {8, 3};
//+
Line(10) = {8, 9};
//+
Line(11) = {9, 10};
//+
Line(12) = {10, 7};
//+
Curve Loop(1) = {5, 7, -12, -11, -10, 9, 3, 4};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {7, 8, 9, -2, -1, -6};
//+
Plane Surface(2) = {2};
//+
Physical Curve("PEC") = {5, 6, 1, 2, 3, 4, 11, 12, 8, 10};
//+
Physical Surface("Air") = {1};
//+
Physical Surface("Dielectric") = {2};
geo mesh to xdmf:
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])]
}
))@wwshunan I do not understand what you mean by saying the physical boundaries are not correct? I’ve opened the mesh in its geo format, and the xdmf format matches the output of the geo format. The triangular mesh consists of two markers, indicating each domain, while you only have one physical tag for the whole boundary.
dokken, did you use the option “XDMF Reader” of paraview 5.8.0 to open the tag_triangle.xdmf file? I used the generated mesh to solve the eigenvalue problem and got the incorrect eigenvalues. Then I created the geometry with the pygmsh module point by point and finally generated the mesh using its generate_mesh function. Everything work. However it’s not convenient to create the geometry by code.
I am using Paraview 5.7 and can use either “XDMF Reader” or “XDMF3ReaderS” to visualize the mesh.
It is also not clear to me what your problem when visualizing the tag_triangle.xdmf is, as you have not explained that in detail.
You need to be able to make a minimal example illustrating that something is wrong. As you have a point by point geometry, you should be able to compare the two meshes and look for differences.
This can be done by using meshio.xdmf.write("filename.xdmf", mesh, data_format="XML"). Then you can compare to similar meshes in plaintext (use coarse meshes).
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).