# Copyright (C) 2017-2021 Chris N. Richardson, Garth N. Wells, Michal Habera
# and Jørgen S. Dokken
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
"""IO module for input data, post-processing and checkpointing"""
import numpy
import ufl
from dolfinx import cpp, fem
import typing
[docs]class VTKFile(cpp.io.VTKFile):
"""Interface to VTK files
VTK supports arbitrary order Lagrangian finite elements for the
geometry description. XDMF is the preferred format for geometry
order <= 2.
"""
[docs] def write_mesh(self, mesh: cpp.mesh.Mesh, t: float = 0.0) -> None:
"""Write mesh to file for a given time (default 0.0)"""
self.write(mesh, t)
[docs] def write_function(self, u: typing.Union[typing.List[fem.Function], fem.Function], t: float = 0.0) -> None:
"""
Write a single function or a list of functions to file for a given time (default 0.0)
"""
cpp_list = None
if isinstance(u, list):
cpp_list = [getattr(u_, "_cpp_object", u_) for u_ in u]
else:
cpp_list = [getattr(u, "_cpp_object", u)
]
super().write(cpp_list, t)
[docs]class XDMFFile(cpp.io.XDMFFile):
[docs] def write_function(self, u, t=0.0, mesh_xpath="/Xdmf/Domain/Grid[@GridType='Uniform'][1]"):
u_cpp = getattr(u, "_cpp_object", u)
super().write_function(u_cpp, t, mesh_xpath)
[docs] def read_mesh(self, ghost_mode=cpp.mesh.GhostMode.shared_facet, name="mesh", xpath="/Xdmf/Domain"):
# Read mesh data from file
cell_shape, cell_degree = super().read_cell_type(name, xpath)
cells = super().read_topology_data(name, xpath)
x = super().read_geometry_data(name, xpath)
# Construct the geometry map
cell = ufl.Cell(cpp.mesh.to_string(cell_shape), geometric_dimension=x.shape[1])
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, cell_degree))
# Build the mesh
cmap = cpp.fem.CoordinateElement(cell_shape, cell_degree)
mesh = cpp.mesh.create_mesh(self.comm(), cpp.graph.AdjacencyList_int64(cells),
cmap, x, ghost_mode, cpp.mesh.partition_cells_graph)
mesh.name = name
domain._ufl_cargo = mesh
mesh._ufl_domain = domain
return mesh
[docs]def extract_gmsh_topology_and_markers(gmsh_model, model_name=None):
"""Extract all entities tagged with a physical marker
in the gmsh model, and collects the data per cell type.
Returns a nested dictionary where the first key is the gmsh
MSH element type integer. Each element type present
in the model contains the cell topology of the elements
and corresponding markers.
"""
if model_name is not None:
gmsh_model.setCurrent(model_name)
# Get the physical groups from gmsh on the form
# [(dim1, tag1),(dim1, tag2), (dim2, tag3),...]
phys_grps = gmsh_model.getPhysicalGroups()
topologies = {}
for dim, tag in phys_grps:
# Get the entities for a given dimension:
# dim=0->Points, dim=1->Lines, dim=2->Triangles/Quadrilaterals
# etc.
entities = gmsh_model.getEntitiesForPhysicalGroup(dim, tag)
for entity in entities:
# Get data about the elements on a given entity:
# NOTE: Assumes that each entity only have one cell-type
element_data = gmsh_model.mesh.getElements(dim, tag=entity)
element_types, element_tags, node_tags = element_data
assert len(element_types) == 1
# The MSH type of the cells on the element
element_type = element_types[0]
num_el = len(element_tags[0])
# Determine number of local nodes per element to create
# the topology of the elements
properties = gmsh_model.mesh.getElementProperties(element_type)
name, dim, order, num_nodes, local_coords, _ = properties
# 2D array of shape (num_elements,num_nodes_per_element) containing
# the topology of the elements on this entity
# NOTE: GMSH indexing starts with 1 and not zero.
element_topology = node_tags[0].reshape(-1, num_nodes) - 1
# Gather data for each element type and the
# corresponding physical markers
if element_type in topologies.keys():
topologies[element_type]["topology"] = numpy.concatenate(
(topologies[element_type]["topology"], element_topology), axis=0)
topologies[element_type]["cell_data"] = numpy.concatenate(
(topologies[element_type]["cell_data"], numpy.full(num_el, tag)), axis=0)
else:
topologies[element_type] = {"topology": element_topology, "cell_data": numpy.full(num_el, tag)}
return topologies
# Map from Gmsh int to DOLFINx cell type and degree
# http://gmsh.info//doc/texinfo/gmsh.html#MSH-file-format
_gmsh_to_cells = {1: ("interval", 1), 2: ("triangle", 1),
3: ("quadrilateral", 1), 4: ("tetrahedron", 1),
5: ("hexahedron", 1), 8: ("interval", 2),
9: ("triangle", 2), 10: ("quadrilateral", 2),
11: ("tetrahedron", 2), 12: ("hexahedron", 2),
15: ("point", 0), 21: ("triangle", 3),
26: ("interval", 3), 29: ("tetrahedron", 3),
36: ("quadrilateral", 3)}
[docs]def ufl_mesh_from_gmsh(gmsh_cell: int, gdim: int):
"""
Create a UFL mesh from a Gmsh cell identifier and the geometric dimension.
See: # http://gmsh.info//doc/texinfo/gmsh.html#MSH-file-format
"""
shape, degree = _gmsh_to_cells[gmsh_cell]
cell = ufl.Cell(shape, geometric_dimension=gdim)
scalar_element = ufl.FiniteElement("Lagrange", cell, degree, variant="equispaced")
return ufl.Mesh(ufl.VectorElement(scalar_element))