I have produced a mesh by converting the .msh and .med files to a .h5 file which also spits out the .xml files too. The mesh is of a simple cylinder but I have grouped the surfaces and volumes seperately as I wish to use this to model a surface energy source. The .msh and .xml files each recognise that different groups are present. I want to just do a sanity check to see if the .h5 file is importing the files correctly by asking it to tell me the area of one of the surfaces.
However, this consistently gives me an answer of zero, even though it should be non-zero. My minimal working example is here:
from __future__ import print_function
import numpy as np
from dolfin import *
import ufl
import math
from mpi4py import MPI
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
set_log_level(50)
#--Import the required mesh files produced in Salome and gmsh-----
mesh_file = 'Mesh_convergence/'
filename = 'mesh_3mm_9'
data = mesh_file+filename
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells_f = MeshFunction('size_t', mesh, 3)
hdf.read(cells_f, '/cells')
facets_f = MeshFunction('size_t', mesh, 2)
hdf.read(facets_f, '/facets')
#-----------------------------------------------
# Mark the tags for 'top' and 'surface' groups
#-----------------------------------------------
tag_s = 2 #Tags for 'surface' group in .msh file
tag_t = 3 # Tag for 'top' group in .msh file
i,j,k = ufl.indices(3)
x = SpatialCoordinate(mesh)
#------------------------------------------------
da = Measure('ds', domain=mesh, subdomain_data = facets_f, metadata = {'quadrature_degree': 2}) #area element
#-------------------------------------------
print(assemble(1*da(tag_t)))
I am wondering if there is an issue with my installation of dolfin that I am running through a conda environment as the function HDF5File comes from dolfin. Does anyone have any recommendations?
It seems like you have written your own medtoh5 converter (or have you gotten it from someone?). If you have written it yourself. are you sure that the data is added to the h5 file?
If you haven’t written it yourself, what is the original source of the file?
I would also prefer that such a script is embedded on the forum with
```python
# add code here
```
rather than in a google drive, that might disappear in the future.
Thank you for your reply. What package do I need to install/use in order to test the h5dump command that you recommended? For the converter, this was provided by a colleague. This converter has been very successful in the past and has previously worked well with prior meshes produced by Salome. The script itself is very long (around 800 lines) so embedding this all into the forum may not be that practical, hence the Google Drive.
Are these, which I do not believe should cause any issues as all variables are continuous across internal facets.
So I guess you have another calculation somewhere, where you would like to compute dot(grad(T),n) across one of these interior facets.
Then you need to restrict the integral to be from a given side, as grad(T) is discontinuous across the facet.
Whereas the current mesh is connected to two cells, like you mentioned. Which would cause the calculation to not compute correctly. So perhaps there is an issue with either how the original .med/.msh file is being output or how the file is being converted to .h5
Hmm I see. I made the mesh via Salome (output as the .med file in the Google Drive) and then converted it to a .msh file via gmsh (again, this is in the Google Drive). You may attempt to convert this yourself into a h5 file to see if the issue lies with the converter or how the mesh was produced in Salome.
That is unfortunately not something I have the bandwidth to do.
As I’ve shown above, you seem to have meshed elements on both sides of the interface in Salome.
That is not something that I can easily circumvent in a meshing conversion script, as the elements are already present in the mesh.
So you are aware, I have mostly solved the issue now. The problem was arising form my Fenics installation in docker, in that whenever I restarted the container, the packages did not fully carry over. Now when the code you provided runs, I obtain the following:
How the long export script you got from your collaborator does the conversion.
Without some effort put into identifying which of these are the issue (and what section of your collaborators code causes the issue), I do not have the bandwidth to give further guidance.
That being said, this is an open source forum, and maybe someone else has the time an expertise to dig into this further.
gmsh writes out msh files, which you can convert to xml or xdmf for DOLFINx with meshio.
There are many tutorials on this on the forum (search for “DOLFIN meshio”)