I’m using DOLFINx 0.4.1 with ubuntu 22.4. Is there a way to make a fenicsx mesh out of the arrays containing the points and the elements. As a example the two lists are given below:
import numpy as np
import matplotlib.pyplot as plt
pp = np.array([[0.5,1],[2,1],[3,1.5],[3.5,2.5],[2.2,2],[1,2.2]])
ee = np.array([[0,1,5],[1,4,5],[1,2,4],[2,3,4]])
plt.triplot(pp[:,0],pp[:,1],ee)
plt.show()
Consider:
import dolfinx
import ufl
from mpi4py import MPI
import numpy as np
gdim = 2
shape = "triangle"
degree = 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0.5, 1], [2, 1], [3, 1.5], [3.5, 2.5], [2.2, 2], [1, 2.2]])
cells = np.array([[0, 1, 5], [1, 4, 5], [1, 2, 4], [2, 3, 4]], dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
many thanks. Are the meshes typically saved as xdmf files in fenicsx and how to read these files back into fenicsx?
Meshes are stored with xdmf in DOLFINx, which you can read back in.
Thank you for the snippet @dokken ! I tried to extend it to my use case that involves a line mesh along x=1 (a 1D mesh in 2D geometry) on which to interpolate another vector-valued Function object like so :
from setup import *
from dolfinx.mesh import create_mesh
from mpi4py.MPI import COMM_WORLD as comm
from dolfinx.fem import FunctionSpace, Function
cell = ufl.Cell("interval", geometric_dimension=2)
VE = ufl.VectorElement("Lagrange", cell, 2, 3)
domain = ufl.Mesh(VE)
n=3
x = np.vstack((np.ones(n),np.linspace(0,2,n))).T
L = np.arange(n, dtype=np.int32)
cells = np.vstack((L[:-1],L[1:])).T
mesh = create_mesh(comm, cells, x, domain)
u = Function(FunctionSpace(mesh, VE))
But I keep getting errors
Traceback (most recent call last):
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 179, in create_mesh
Traceback (most recent call last):
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 179, in create_mesh
mesh = _cpp.mesh.create_mesh(comm, cells, cmap, x, partitioner)
TypeError: create_mesh(): incompatible function arguments. The following argument types are supported:
1. (comm: MPICommWrapper, cells: dolfinx::graph::AdjacencyList<long>, element: dolfinx::fem::CoordinateElement, x: numpy.ndarray[numpy.float64], partitioner: Callable[[MPICommWrapper, int, int, dolfinx::graph::AdjacencyList<long>], dolfinx::graph::AdjacencyList<int>]) -> dolfinx::mesh::Mesh
Invoked with: <mpi4py.MPI.Intracomm object at 0x7f7d45d7fc30>, array([[0, 1],
[1, 2]], dtype=int32), <dolfinx.cpp.fem.CoordinateElement object at 0x7f7d3e9d5830>, array([[1., 0.],
[1., 1.],
[1., 2.]]), <built-in method of PyCapsule object at 0x7f7d2f8da3d0>
During handling of the above exception, another exception occurred:
mesh = _cpp.mesh.create_mesh(comm, cells, cmap, x, partitioner)
TypeError: create_mesh(): incompatible function arguments. The following argument types are supported:
1. (comm: MPICommWrapper, cells: dolfinx::graph::AdjacencyList<long>, element: dolfinx::fem::CoordinateElement, x: numpy.ndarray[numpy.float64], partitioner: Callable[[MPICommWrapper, int, int, dolfinx::graph::AdjacencyList<long>], dolfinx::graph::AdjacencyList<int>]) -> dolfinx::mesh::Mesh
Invoked with: <mpi4py.MPI.Intracomm object at 0x7f1e60d27c30>, array([[0, 1],
[1, 2]], dtype=int32), <dolfinx.cpp.fem.CoordinateElement object at 0x7f1e4a8644f0>, array([[1., 0.],
[1., 1.],
[1., 2.]]), <built-in method of PyCapsule object at 0x7f1e4a8823d0>
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "ex.py", line 19, in <module>
Traceback (most recent call last):
File "ex.py", line 19, in <module>
mesh = create_mesh(comm, cells, x, domain)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 181, in create_mesh
mesh = create_mesh(comm, cells, x, domain)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 181, in create_mesh
mesh = _cpp.mesh.create_mesh(comm, _cpp.graph.AdjacencyList_int64(np.cast['int64'](cells)),
RuntimeError: A facet is connected to more than two cells.
mesh = _cpp.mesh.create_mesh(comm, _cpp.graph.AdjacencyList_int64(np.cast['int64'](cells)),
RuntimeError: A facet is connected to more than two cells.
I get the same traceback (RuntimeError: A facet is connected to more than two cells.) using the snippet in this post : MeshEditor equivalent - #6 by dokken
Reminded here for convenience :
from dolfinx.fem import FunctionSpace
from dolfinx.mesh import create_mesh
from mpi4py import MPI
import ufl
import numpy as np
gdim, shape, degree = 1, "interval", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0], [0.2], [0.6], [1], [1.3]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
V = FunctionSpace(mesh, ("P", 1))
print(V.tabulate_dof_coordinates(), mesh.geometry.dim)
print(mesh.topology.index_map(mesh.topology.dim).size_local)
I’m using the docker build of dolfinx, it’s surprising that such a simple mesh fails…
The dtype here should be np.int64 as you are sending in the global topology indices
What version of DOLFINx are you running through docker? I cannot reproduce that error:
xxxx@yyyyy:~/Documents/debug$ docker run -ti --network=host -e DISPLAY=$DISPLAY -v /tmp/.X11-unix:/tmp/.X11-unix -v $(pwd):/root/shared -w /root/shared --rm dolfinx/dolfinx:v0.6.0-r1
root@dokken-XPS-9320:~/shared# ipython
Python 3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.9.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: from dolfinx.fem import FunctionSpace
...: from dolfinx.mesh import create_mesh
...: from mpi4py import MPI
...: import ufl
...: import numpy as np
...:
...: gdim, shape, degree = 1, "interval", 1
...: cell = ufl.Cell(shape, geometric_dimension=gdim)
...: domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
...: x = np.array([[0], [0.2], [0.6], [1], [1.3]])
...: cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
...: mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
...:
...: V = FunctionSpace(mesh, ("P", 1))
...: print(V.tabulate_dof_coordinates(), mesh.geometry.dim)
...: print(mesh.topology.index_map(mesh.topology.dim).size_local)
[[0. 0. 0. ]
[0.2 0. 0. ]
[0.6 0. 0. ]
[1. 0. 0. ]
[1.3 0. 0. ]] 1
4
Version is 0.6.0, problem persists after changing type to int64
Ok I have it ! It fails in parallel only. I’m also not getting an error in serial.
You could have a look at my very experimental new documentation of mesh creation:
https://jsdokken.com/dolfinx_docs/meshes.html
which covers how to input points from numpy arrays in serial and parallel.
If you use python files and not notebooks, you can just ignore the ipyparallel stuff, as it is there so that one can run this with multiple processes in a notebook.
Thanks for the pointer. I preferred the boring workaround of generating the mesh in serial than reading it in parallel again. Not entirely satisfying, but it works.
Question: How to construct meshtags from the pair (facets, values) and not (facet_indices, values)?
I am trying to read mesh data from here and saving it as .xdmf file with meshtag information.
import h5py
import numpy as np
import dolfinx
import ufl
from mpi4py import MPI
filename = "colin27_coarse_boundaries.h5"
file = h5py.File(filename, 'r')
coordinates = np.array(file["mesh/coordinates"])
topology = np.array(file["mesh/topology"])
boundary_coordinates = np.array(file["boundaries/coordinates"])
boundary_facets = np.array(file["boundaries/topology"])
tags = np.array(file["boundaries/values"], dtype=np.int64)
gdim = 3
shape = "tetrahedron"
degree = 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
ufl_tetra = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, topology, coordinates, ufl_tetra)
tdim = mesh.topology.dim
fdim = tdim - 1
facet_tag = dolfinx.mesh.meshtags_from_entities(mesh, fdim, boundary_facets, tags)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "brain.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
Thanks
I should have emphasized that I have boundary_facets which store for each facet it’s vertices. If I had facet indices then that would be easier.
>>> boundary_facets
array([[ 0, 1210, 3570],
[ 0, 1210, 9339],
[ 0, 1210, 15409],
...,
[28959, 28989, 28990],
[28966, 28989, 28990],
[28978, 28995, 28997]])
>>> boundary_facets.shape
(224786, 3)
>>> tags.shape
(224786,)
I had tried meshtags as well. It expects that both boundary_facets and tags be of same shape.
Using locate_exterior_entities gives indices of facets on the exterior.
Then an API for getting the vertices from the facet index will be needed, with which the facets identified in boundary_facets through veritces will get correct facet index.
It might be trivial, but I don’t know about it.
Thanks.
Use: dolfinx.io — DOLFINx 0.6.0 documentation
as done in:
Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken (Long tutorial bit)
and talked about in Finding facet index in imported mesh - #2 by dokken
Hello,
I have query regarding using meshtags once we assign them. I am not finding an MWE where meshtags are incorporated and used in assembly using dx.
import numpy as np
tri_points = np.array([[0,0],[0,1],[1,0],[1,1]], dtype=np.float64)
triangles = np.array([[0,1,3], [0,2,3]], dtype=np.int64)
import ufl
import basix
ufl_tri = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3, )))
import dolfinx
from mpi4py import MPI
mesh = dolfinx.mesh.create_mesh(
MPI.COMM_WORLD, triangles, tri_points, ufl_tri)
# Subdomain data (MeshTags)
cell_map = mesh.topology.index_map(mesh.topology.dim)
num_cells_on_process = cell_map.size_local + cell_map.num_ghosts
cells = np.arange(num_cells_on_process, dtype=np.int64)
subdomains = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, cells, np.array([0,2]))
print(subdomains.values)
V = dolfinx.fem.functionspace(mesh, basix.ufl.element("CG", "triangle", 1, shape=(3, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
dx = ufl.Measure('dx')(domain=mesh, subdomain_data=subdomains)
F2 = sum([(dot(u,v)-6)*dx(i) for i in range(2)])
A= dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(ufl.lhs(F2)))
A.assemble()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[11], line 6
4 dx = ufl.Measure('dx')(domain=mesh, subdomain_data=subdomains)
5 F2 = sum([(dot(u,v)-6)*dx(i) for i in range(2)])
----> 6 A= dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(ufl.lhs(F2)))
7 A.assemble()
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:249, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
246 return list(map(lambda sub_form: _create_form(sub_form), form))
247 return form
--> 249 return _create_form(form)
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:244, in form.<locals>._create_form(form)
241 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
242 return form argument"""
243 if isinstance(form, ufl.Form):
--> 244 return _form(form)
245 elif isinstance(form, collections.abc.Iterable):
246 return list(map(lambda sub_form: _create_form(sub_form), form))
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:222, in form.<locals>._form(form)
219 return [(s[0], np.array(s[1])) for s in subdomain]
221 # Subdomain markers (possibly empty list for some integral types)
--> 222 subdomains = {
223 _ufl_to_dolfinx_domain[key]: get_integration_domains(
224 _ufl_to_dolfinx_domain[key], subdomain_data[0]
225 )
226 for (key, subdomain_data) in sd.get(domain).items()
227 }
229 f = ftype(
230 module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)),
231 V,
(...)
236 mesh,
237 )
238 return Form(f, ufcx_form, code)
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:223, in <dictcomp>(.0)
219 return [(s[0], np.array(s[1])) for s in subdomain]
221 # Subdomain markers (possibly empty list for some integral types)
222 subdomains = {
--> 223 _ufl_to_dolfinx_domain[key]: get_integration_domains(
224 _ufl_to_dolfinx_domain[key], subdomain_data[0]
225 )
226 for (key, subdomain_data) in sd.get(domain).items()
227 }
229 f = ftype(
230 module.ffi.cast("uintptr_t", module.ffi.addressof(ufcx_form)),
231 V,
(...)
236 mesh,
237 )
238 return Form(f, ufcx_form, code)
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:214, in form.<locals>._form.<locals>.get_integration_domains(integral_type, subdomain)
212 subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim)
213 subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1)
--> 214 domains = _cpp.fem.compute_integration_domains(
215 integral_type, subdomain._cpp_object
216 )
217 return [(s[0], np.array(s[1])) for s in domains]
218 except AttributeError:
TypeError: compute_integration_domains(): incompatible function arguments. The following argument types are supported:
1. compute_integration_domains(integral_type: dolfinx.cpp.fem.IntegralType, meshtags: dolfinx.cpp.mesh.MeshTags_int32) -> list[tuple[int, list[int]]]
Invoked with types: dolfinx.cpp.fem.IntegralType, dolfinx.cpp.mesh.MeshTags_int64
Can you share where meshtags used from input data of user and done solving var formulation similar to above what I am getting error in?
Change the dtype here to numpy.int32.
There are multiple demos illustrating the usage, like
and
and in
Thanks for sharing the links.
It worked with below change and other dtype with np.int32. The post helped.
subdomains = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, cells, np.array(matid[ii], dtype=np.int32))