hexahedral mesh

Asked by Lars on 2010-06-23

Hello,
I wonder if I can create hexahedral meshes in Dolfin?
The hexahedron is listed as one of the reference cells in the Dolfin user guide.
But the following python snippet fails with 'RuntimeError: *** Error: Unknown cell type "hexahedron"':

  from dolfin import *
  mesh=Mesh();
  e=MeshEditor()
  e.open(mesh,'hexahedron',3,3)

Diving into dolfin-0.9.7/dolfin/mesh/MeshEditor.cpp I find the open function:

  void MeshEditor::open(Mesh& mesh, std::string type, uint tdim, uint gdim)
  {
  if (type == "point")
    open(mesh, CellType::point, tdim, gdim);
  else if (type == "interval")
    open(mesh, CellType::interval, tdim, gdim);
  else if (type == "triangle")
    open(mesh, CellType::triangle, tdim, gdim);
  else if (type == "tetrahedron")
    open(mesh, CellType::tetrahedron, tdim, gdim);
  else
    error("Unknown cell type \"%s\".", type.c_str());
  }

Apparently hexahedrons and quadrilaterals are not supported.
Can I somehow create a hexahedral mesh?
If not, what is the purpose of the hexahedron cell type?

Thanks
Lars

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Lars
Solved:
2010-06-28
Last query:
2010-06-28
Last reply:
2010-06-28
Anders Logg (logg) said : #1

On Wed, Jun 23, 2010 at 03:11:41PM -0000, Lars wrote:
> New question #115551 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/115551
>
> Hello,
> I wonder if I can create hexahedral meshes in Dolfin?
> The hexahedron is listed as one of the reference cells in the Dolfin user guide.
> But the following python snippet fails with 'RuntimeError: *** Error: Unknown cell type "hexahedron"':
>
> from dolfin import *
> mesh=Mesh();
> e=MeshEditor()
> e.open(mesh,'hexahedron',3,3)
>
> Diving into dolfin-0.9.7/dolfin/mesh/MeshEditor.cpp I find the open function:
>
> void MeshEditor::open(Mesh& mesh, std::string type, uint tdim, uint gdim)
> {
> if (type == "point")
> open(mesh, CellType::point, tdim, gdim);
> else if (type == "interval")
> open(mesh, CellType::interval, tdim, gdim);
> else if (type == "triangle")
> open(mesh, CellType::triangle, tdim, gdim);
> else if (type == "tetrahedron")
> open(mesh, CellType::tetrahedron, tdim, gdim);
> else
> error("Unknown cell type \"%s\".", type.c_str());
> }
>
> Apparently hexahedrons and quadrilaterals are not supported.
> Can I somehow create a hexahedral mesh?
> If not, what is the purpose of the hexahedron cell type?

It's there if someone wants to go ahead and implement it.

It is essentially supported by the mesh data structure, but there is
no MeshEditor for it, no built-in hexahedral meshes, and probably
assumptions on simplex meshes in a few places.

--
Anders

Lars (lars-g-m-x) said : #2

Thank you, Anders. I am interested to work on the implementation.
I read your paper "Efficient representation of computational meshes", which is a good explanation of the data structures.

Can you recommend more documentation on the MeshEditor concepts and implementation, or do I have to study the sources in detail? I think I have read somewhere that without the MeshEditor you can construct meshes of arbitrary cells - can I find some code examples (c++ or python) how to achieve this?

And finally, is there any reason why the MeshEditor handles only meshes of one single CellType? Why not allow a combination of triangles and quadrilaterals, for example, or tets and hexes and wedges in 3D?

Thanks
Lars

Anders Logg (logg) said : #3

On Fri, Jun 25, 2010 at 08:00:58AM -0000, Lars wrote:
> Question #115551 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/115551
>
> Status: Answered => Open
>
> Lars is still having a problem:
> Thank you, Anders. I am interested to work on the implementation.

Good luck! It will be interesting to see whether you can make it work.

Note that even if you get the mesh to work, other parts of FEniCS
(like FIAT and FFC) need some work to get everything running on hex
meshes.

> I read your paper "Efficient representation of computational meshes", which is a good explanation of the data structures.
>
> Can you recommend more documentation on the MeshEditor concepts and
> implementation, or do I have to study the sources in detail? I think I
> have read somewhere that without the MeshEditor you can construct meshes
> of arbitrary cells - can I find some code examples (c++ or python) how
> to achieve this?

You would probably need to read the source in detail.

What MeshEditor does is that it interacts with the low-level data
structures of the Mesh class (MeshTopology and MeshGeometry).

You probably need to rename some of the overloaded functions in
MeshEditor. In particular add_cell which is overloaded on the number
of arguments (2 arguments is an interval, 3 is a triangle, 4 is a
tet) since a tetrahedron has the same number of vertices as a
quadrilateral. Better names would be add_interval, add_triangle etc.

> And finally, is there any reason why the MeshEditor handles only meshes
> of one single CellType? Why not allow a combination of triangles and
> quadrilaterals, for example, or tets and hexes and wedges in 3D?

It's because the Mesh only has one cell type which is specified by the
CellType* _cell_type member in the Mesh class. This is a limitation
that we might also consider changing, but that would also require a
bit of work.

--
Anders

Lars (lars-g-m-x) said : #4

Thank you. After browsing the sources and reading your second reply the task seems much more time consuming than after your first reply. I will check if splitting my hexahedral mesh into tetrahedrons is a better option.

Lars