Isoparametric facet function not loading correctly

mesh
1

I am reposting a previous question, with a slightly modified code and a full working example.

I have one issue and one question…

Issue: I am trying to generate a mesh with triangle6 isoparametric elements, and the second-last line of this code yields an error, something which didn’t happen in the regular order 1 mesh:

import pygmsh
import gmsh
import meshio
import matplotlib.pyplot as plt

from dolfin import *
from dolfin_adjoint import *

geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()

resolution = .1
path = "/home/leonardo/PycharmProjects/masters_thesis/meshes/annulus/"

# A circle centered at the origin and radius 1
circle = model.add_circle([0.0, 0.0, 0.0], 1, mesh_size=2.5 * resolution)

# A hole
hole = model.add_circle([0.0, 0.0, 0.0], .5, mesh_size=1 * resolution)

# My surface
plane_surface = model.add_plane_surface(circle.curve_loop, [hole.curve_loop])

model.synchronize()

model.add_physical([plane_surface], "volume")
model.add_physical(circle.curve_loop.curves, "outer_ring")
model.add_physical(hole.curve_loop.curves, "inner_ring")

geometry.generate_mesh(dim=2, order=2)

gmsh.write(path + "mesh_" + str(resolution) + ".msh")
gmsh.clear()
geometry.__exit__()

mesh_from_file = meshio.read(path + "mesh_" + str(resolution) + ".msh")


def mesh_to_meshio(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)  # get the cells of some type: it will change!
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh


line_mesh = mesh_to_meshio(mesh_from_file, "line3", prune_z=True)
meshio.write(path + "facet_mesh_" + str(resolution) + ".xdmf", line_mesh)
triangle_mesh = mesh_to_meshio(mesh_from_file, "triangle6", prune_z=True)
meshio.write(path + "mesh_" + str(resolution) + ".xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile(path + "mesh_" + str(resolution) + ".xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(path + "facet_mesh_" + str(resolution) + ".xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = MeshFunction("size_t", mesh, mvc)

The error is:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  13fd4d8cac9cac928b724b3540e1aa8bfb50b34d
*** -------------------------------------------------------------------------

Question: are the boundary edges curved in this formulation or are they just straight-edge polygons again? How can I visualize this in Fenics?

2

There is very limited support for higher order meshes (geometries) in legacy DOLFIN. This support has been improved in DOLFINx, see for instance: dolfinx/demo_gmsh.py at main · FEniCS/dolfinx · GitHub

2 Likes
3

Thanks for the answer. As a short follow-up question to your comment, is dolfin-x compatible with dolfin-adjoint? This would be my main reason to/not to go to dolfin-x…

4

DOLFINx is not compatible with dolfin_adjoint, as none of the developers (including me) has not had time to do this.

1 Like