Different behaviour of a mesh coming from generate_mesh() resp. UnitSquareMesh()?

Gmsh creates triangles and tetrahedrons by default. Using recombine, see: Files · master · gmsh / gmsh · GitLab

The create mesh function is part of dolfinx not dolfin.

One is to save the mesh, while the other one is saving the facet data. To read it into dolfin, these needs to be separate.

You can load it in the same way as you load in mesh data. See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #3 by dokken

Still, how can I tell which one I have in a mesh so I can read it appropriately?

This is a user defined function. In the thread mentioned earlier, the components of this function are used separate. So I found it easier to make it a function. Why/What do I need to differently?

I do not see how in this comment is presented how to read the named boundaries from gmsh to fenics. I expected to have something reading the name boundaries and then using something like bcu_noslip = DirichletBC(V, Constant((0, 0)), mf, inflow) in order to apply the BC. Also, what do cpp, MeshFunctionSizet do?

Not sure what you mean by this. As that tutorial explains, if you do not use recombine, your mesh will consist of triangles when you call

gmsh.model.mesh.generate(2)

and tetrahedrons if you use

gmsh.model.mesh.generate(3)

As far as I can tell, you are referring to:

from dolfinx.mesh import create_mesh
from mpi4py import MPI
mesh = create_mesh(MPI.COMM_WORLD, cells, x, ufl_domain)

As I said earlier, if you are using fenics/dolfin, you should stop following the tutorial after

gmsh.write("mesh3D.msh")

cpp.MeshFunctionSizet creates the MeshFunction you can use in your snippet. Dolfin does not read in strings, so inflow has to be the integer you tagged your boundary with.

You wrote the same command for either triangles or tetrahedrons. So how can I tell if my mesh is with triangles or tetrahedrons?

No, I am using

def create_mesh(mesh, cell_type): cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type]) data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key] for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type]) mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[data]}) return mesh

So in order to get the inlet from the gmsh file should I write something like:

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
          infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
bcu_noslip = DiriclhetBC(V, Constant((0,0)),mf,3)

where 3 is chosen as this boundary has defined 3rd on gmsh?

There is still no documentation to explain what these commands and their arguments are so I still don’t know what the MeshValueCollection does and what are the input both for that and the MeshFunctionSize correspond to.

I also get this error:

mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
AttributeError: module 'mshr.cpp' has no attribute 'mesh'

Sorry, typo on my side, corrected in the precious post to be 2 for triangles,3 for tetrahedrons.

Then you simply write the string for the cell type you are using. You can expect mesh.cells_dict to see the different typos.

Inlet should match the inlet_marker you defined when generating your mesh.

Are you expecting me to write out the documentation of every function? FEniCS has a Python API
and a C++ API.
For instance, you can find what a MeshValueCollection is in the C++ API:
https://fenicsproject.org/olddocs/dolfin/latest/cpp/d0/db6/classdolfin_1_1MeshValueCollection.html
Under detailed description.

Dolfin in general has plenty of documented demos Which cover the base usage of dolfin: https://fenicsproject.org/olddocs/dolfin/latest/python/demos.html
and there are many posts on this forum on how to read in meshes.

I am sorry if I came across as rude, I meant that I cannot find any documentation and therefore I am asking all these questions here.

quote=“EvaAnton, post:25, topic:1901”]
I also get this error:

mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

AttributeError: module ‘mshr.cpp’ has no attribute ‘mesh’

[/quote]

Any ideas about that?

Thank you these! I will read them carefully.

What do you mean expect?

Remove from mshr import * from your program. Mshr is no longer maintained and should not be imported.

I meant inspect.
For instance by using Ipython’s function embed
(IPython reference — IPython 8.17.2 documentation) you can inspect the data read in by meshio before writing them to XDMF

Hi,

I have imported the mesh in fenics now correctly, thank you for your help.

In the weak formulation, I need dx and ds. When I was generating the mesh using fenics functions, these two were computed “automatically” (I think so at least, don’t know how to say it). Do I need to compute them separately here?

Thank you!

P.S.: I saw the updated tutorials and the new tutorial for dolfin-x on your webpage. Should I move to dolfin-x? I thought that wasn’t supporting fenics.

Any measure, dx or ds In dolfin tried to extract the domain (the mesh) from the arguments in the form. If you want to make a measure with a given mesh, you can do the following:

V = FunctionSpace(mesh, "CG", 1)
u = interpolate(Expression("x[0]", degree=1), V)
dx_C = Measure("dx", domain=mesh)
print(assemble(u*dx_C), assemble(u*dx))

Note that the standard measure dx and dx_C gives the same result. A custom measure (dx_C) should be used whenever you are integrating over a subset of a volume/surface, by supplying the keyword argument: subdomain_id=meshfunction_for_entity, subdomain_id=value_in_mesh_function.

I would not suggest moving to dolfin-X just yet. Dolfin-x is the current development version of dolfin (no active development or bug fixing in the traditional dolfin/fenics).
However, as it is under development, the user interface might change rapidly.
The tutorial I’ve added on my webpage is a work in progress, which I update every time the interface of dolfin-x changes.

Thank you for that. So for now, that I just have a mesh for the whole domain (no sudbomains), I don’t need to compute any new dx or ds, right?

Makes sense! Thanks!

That is correct (you can verify this by using the snippet from my previous post)

1 Like