Introduction

About mefikit

mefikit — short for Mesh and Field Kit — is a comprehensive library for generating, manipulating, and analyzing unstructured meshes together with associated scalar, vector, and tensor fields. Its goal is to provide a unified in-memory representation of meshes and fields, along with a set of robust, efficient tools that support numerical simulation workflows in research, engineering, and scientific computing.

mefikit focuses on three core principles:

1. A flexible mesh model capable of representing mixed-element unstructured meshes.
2. A consistent field and group architecture that attaches data to mesh entities across dimensions.
3. Zero-copy, iterator-based access patterns to efficiently navigate and process large meshes.

This design positions mefikit as a modern core library for algorithm development.

mefikit provides

• umesh the unstructured mesh container supporting mixed element types, fields, and groups.
• io modules for reading and writing meshes in various formats (e.g., VTK, serde_json, serde_yaml).
• topology tools for analyzing mesh connectivity, computing descending meshes, neighbours, domain frontier, etc.
• geometry tools for computing element measures, centroids, etc.
• Selector utilities for querying and filtering mesh elements based on geometric or topological criteria.

Definitions

Element

An element is an abstract mesh entity specified by:

• its type (e.g., triangle, hexahedron, polygon, polyhedron)
• its connectivity (indices referencing the global coordinate array)

In mefikit, elements are ephemeral, zero-copy views constructed on the fly when iterating through a mesh block. They do not own memory; instead they refer to underlying index and coordinate buffers. This behavior mirrors VTK’s “cell iterator” pattern while improving cache locality and reducing overhead.

It is only accessible through the rust API. The python bindings do not expose elements directly because python is only intended to be used for high-level scripting. Do not fear of going to the rust side for performance critical code, mefikit was designed to be simple to use for rust newcommers.

Topological Dimension

The topological dimension of an element is the dimension of the mathematical object it represents:

• 0D → vertices
• 1D → edges / segments
• 2D → faces
• 3D → volumes

This categorization is independent of the embedding space and is central for field association and mesh validity checks.

Spatial Dimension

The spatial dimension is the dimension of the coordinate system in which the mesh is embedded, typically 2D or 3D. All element coordinates must be consistent with this spatial dimension. For example, one may have:

• 2D elements in 2D (pure surface mesh)
• 2D elements in 3D (surface embedded in 3D)
• 3D volumetric elements in 3D

mefikit does not assume a fixed coordinate system beyond this requirement.

Mesh

A mesh in mefikit is a container that holds:

• the global node coordinate array,
• a set of blocks, each containing elements of a given type,
• fields (scalar/vector/tensor) attached to elements of a given topological dimension,
• groups, i.e., labeled subsets of elements for boundary conditions, materials, or post-processing.

Properties:

• You may attach any number of fields or groups.
• Groups may overlap arbitrarily.
• A field must cover all elements of its associated topological dimension (similar to Gmsh’s “node data” or “element data” consistency).

If these constraints seem restrictive, mefikit encourages using multiple meshes (e.g., one per partition, or one per physical domain) to better match specialized workflows.

Topologically Valid Mesh

A topologically valid mesh in mefikit satisfies the following conditions:

1. No overlapping higher-dimensional elements Two 3D elements may share faces or edges, but they must not occupy the same region of space. (Duplicate faces in a hex-dominant mesh are allowed and interpreted in context.)

2. Lower-dimensional entities derive consistently from higher-dimensional ones

   • All edges must be edges of some face or volume.
   • All faces must be boundaries of some volume.
   • No “floating” or orphaned lower-dimension elements are allowed. This is consistent with the expectations of many finite-element and finite-volume codes.

3. Non-degenerate geometry Elements must not be geometrically degenerate:

   • triangles must have non-zero area,
   • tetrahedra must have non-zero volume,
   • polygons must be well-defined and non-self-intersecting, etc. These checks are similar to the geometric validation steps performed by MeshGems and various meshing libraries.

Element Conventions

This chapter defines the node ordering and face orientation conventions used throughout mefikit. These conventions follow the VTK convention for element node ordering and face definitions.

Note: TET10 uses MEDFile mid-side node numbering (deferred migration).

Element types

Poly element representation

Poly elements (SPLINE, PGON, PHED) store their connectivity as a flat array of node indices. Faces are delimited by usize::MAX sentinel values.

• PGON: nodes form a single face, no sentinel needed. [n0, n1, n2, n3] — a 4-node polygon.

• PHED: each face is a closed polygon (PGON), separated by sentinels.

  [f0_n0, f0_n1, f0_n2, MAX, f1_n0, f1_n1, f1_n2, f1_n3, MAX, ...]
  

  The last face has no trailing sentinel.

Face orientation convention

All subentities (edges of 2D elements, faces of 3D elements) are defined with consistent counter-clockwise (CCW) winding when viewed from outside the element. This ensures that:

• Each outward-facing normal follows the right-hand rule.
• Shared lower-dimensional entities are traversed in opposite directions by adjacent higher-dimensional elements.

The subentity definitions below are the canonical reference for node orderings.

2D elements

TRI3 / TRI6 / TRI7 — Triangle

Nodes 0, 1, 2 are the three vertices. Edges follow CCW winding:

        2
       / \
      /   \
     /     \
    0-------1

For TRI6/TRI7, mid-side nodes are placed as: node 3 on edge 01, node 4 on edge 12, node 5 on edge 20. TRI7 additionally has a centroid node (node 6).

QUAD4 / QUAD8 / QUAD9 — Quadrilateral

Nodes 0–3 are the four vertices, numbered counter-clockwise:

    3-------2
    |       |
    |       |
    0-------1

3D elements

TET4 — Tetrahedron

Nodes 0–3 are the four vertices. The tetrahedron is defined by its four triangular faces (VTK convention):

          3
         /|\
        / | \
       /  |  \
      /   |   \
     /    |    \
    2-----+-----1
     \    |    /
      \   |   /
       \  |  /
        \ | /
         \|/
          0

Each face lists the three nodes that do not include the opposite vertex, in CCW order when viewed from outside the element.

For TET10, the vertex convention is: nodes 0–3 vertices, mid-side nodes 4 (edge 01), 5 (edge 12), 6 (edge 02), 7 (edge 03), 8 (edge 13), 9 (edge 23).

HEX8 — Hexahedron

Nodes 0–7 are the eight vertices of a unit cube, numbered following the VTK convention:

        7---------6
       /|        /|
      / |       / |
     /  |      /  |
    4---------5   |
    |   3-----|---2
    |  /      |  /
    | /       | /
    |/        |/
    0---------1

Bottom face: nodes 0, 1, 2, 3 (CCW viewed from below). Top face: nodes 4, 5, 6, 7 (CCW viewed from above). Node 4 is directly above node 0, node 5 above node 1, etc.

All faces are wound CCW when viewed from outside the element.

Differences with MEDFile: MEDFile HEX8 has a different node numbering: MED node 0 is top-left-front, while VTK node 0 is bottom-left-front. The MED→VTK node permutation is [4,5,6,7,0,1,2,3] (self-inverse). Do not mix VTK and MEDFile conventions

For HEX21, the vertex convention is: nodes 0–7 vertices, mid-side nodes 8 (edge 01), 9 (edge 12), 10 (edge 23), 11 (edge 30), 12 (edge 45), 13 (edge 56), 14 (edge 67), 15 (edge 74), 16 (edge 04), 17 (edge 15), 18 (edge 26), 19 (edge 37).

PHED — Polyhedron

A polyhedron is stored as a sequence of polygonal faces, each wound CCW when viewed from outside the element. Faces are delimited by usize::MAX sentinels in the flat connectivity array.

The constraints for a valid PHED are:

1. Each face is a closed polygon — its nodes form a simple, non-self-intersecting loop.
2. Consistent orientation — every edge shared by exactly two faces is traversed in opposite directions by those two faces. Equivalently, the outward normal of every face follows the right-hand rule.
3. Closed surface — the collection of faces forms a topologically closed manifold (no dangling edges).

Topological basic operations

Subgraph and neighbours computation

The subgraph computation is done supposing that the mesh is topologically valid. The subgraph of a volume mesh is defined by a set of either faces, edges or vertices elements depending on the subgraph relative dimension.

The subgraph computation and the neighbours computation are tightly linked. The neighbours are defined as the elements that share a common lower dimensional element. For example, two volumes are neighbours if they share a common face, two faces are neighbours if they share a common edge, and so on.

In mefikit, the neighbours computation builds an adjacency graph, having as nodes the elements of the mesh and as edges the shared lower dimensional elements. Then, the neighbours of an element can be found by looking at its adjacent nodes in the graph.

Connectivity equivalence

There is different ways to represent the connectivity of an element. Here are the different class of equivalence:

• exact representation equality: two elements are equivalent if their connectivity is exactly the same, including the order of the nodes.
• rotational equivalence: two elements are equivalent if their connectivity can be made the same by rotating the order of the nodes. The topological shape is strictly equivalent under all operations, preserving the measure orientation.
• chiral equivalence: two elements are equivalent if their connectivity can be made the same by reversing the order of the nodes and rotating the nodes. This equivalence does not preserve the measure orientation (if it is positive or negative).
• node set equivalence: two elements are equivalent if they have the same set of nodes, regardless of the order. This equivalence is theoretical only, as it does not preserve the shape of the element. But it is very useful to detect duplicate elements in a mesh.

Python guide

This page is a compact reference for the Python API. The notebooks under Python Examples show the same features in context; here they are gathered as tables.

Meshes are UMesh objects. Fields and element groups live in two dict-like mappings on the mesh, and selections are lazy views that only evaluate when queried.

The fields mapping

mesh.fields behaves like a dict keyed by field name. Each entry returns a FieldRef, a handle to read, reduce, or write the stored values.

Accepted values for creation and writes:

• a float — broadcast to every row
• an np.ndarray — full column or per-block rows
• a dict {etype: array} — per element type
• an expression: mf.Field("T") * 2 or another field object
• a string naming an existing field (e.g. "T") — copies it

Reductions over all elements carrying the field: min(), max(), sum(), mean(), var(ddof=0), std(ddof=0), integral() (measure-weighted).

Partial reads and writes

ref[selector] gathers the selected rows as {etype: array}; assigning through ref[selector] = value writes them. Selectors are:

• wildcards: ..., : (full slice), or None
• an ids dict: {"QUAD4": [0, 3]}
• any selection expression: mf.sel.rect(...), mf.Field("T") > 1.0, …

The groups mapping

mesh.groups behaves like a dict keyed by group name. Each entry is a GroupRef.

Groups feed back into selections through mf.sel.group("wall") and mf.sel.exclude_group("wall").

Selections

Selection factories live in the mf.sel module:

Selections compose with &, |, ^, -, ~. Field thresholds produce selections too: mf.Field("T") > 1.0.

Two families of spatial selectors are available (showcased in the selection notebook):

• the n* variants (nbbox, nrect, nsphere, ncircle, nids) match node positions and take an all= flag (all vs. any node of the element must match);
• bbox, rect, sphere, circle, ids match element centroids.

Lazy results

mesh.select(expr) does not build a mesh; it returns a lightweight SelectionResult that re-evaluates on every call:

• result.ids() → {etype: array}
• len(result)
• reductions with any field expression: min/max/sum/mean(expr), var/std(expr, ddof=0), integral(expr)
• result.to_mesh(with_fields=True) materializes a sub-mesh when needed

hot = mesh.select(mf.Field("energy") > 1e6)
print(hot.mean("energy"))
submesh = hot.to_mesh()

Mesh modification

Most topological tools return a new UMesh; the *_update variants operate in-place and return a new mesh only when the result displaced elements (otherwise None).

Field expressions (notably mf.M for the on-the-fly measure) can be evaluated without a stored field:

• mesh.eval(expr, dim=None) → {etype: array}, e.g. mf.M or mf.Field("T") * 2
• mesh.eval_update(name, expr, dim=None) stores the result in-place
• mesh.measure() → per-type measures; mesh.measure_update() materializes a "Measure" field (usually unnecessary, prefer mf.M)

Python Examples

UMesh basics

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Building cartesian meshes

volumes = mf.build_cmesh(
    range(2), np.linspace(0.0, 1.0, 3), np.logspace(0.0, 1.0, 4) / 10.0
)

print(volumes)

UMeshBase {
    coords: [[0.0, 0.0, 0.1],
     [1.0, 0.0, 0.1],
     [0.0, 0.5, 0.1],
     [1.0, 0.5, 0.1],
     [0.0, 1.0, 0.1],
     [1.0, 1.0, 0.1],
     [0.0, 0.0, 0.2154434690031884],
     [1.0, 0.0, 0.2154434690031884],
     [0.0, 0.5, 0.2154434690031884],
     [1.0, 0.5, 0.2154434690031884],
     [0.0, 1.0, 0.2154434690031884],
     [1.0, 1.0, 0.2154434690031884],
     [0.0, 0.0, 0.46415888336127786],
     [1.0, 0.0, 0.46415888336127786],
     [0.0, 0.5, 0.46415888336127786],
     [1.0, 0.5, 0.46415888336127786],
     [0.0, 1.0, 0.46415888336127786],
     [1.0, 1.0, 0.46415888336127786],
     [0.0, 0.0, 1.0],
     [1.0, 0.0, 1.0],
     [0.0, 0.5, 1.0],
     [1.0, 0.5, 1.0],
     [0.0, 1.0, 1.0],
     [1.0, 1.0, 1.0]], shape=[24, 3], strides=[3, 1], layout=Cc (0x5), const ndim=2,
    element_blocks: {
        HEX8: ElementBlockBase {
            cell_type: HEX8,
            connectivity: Regular(
                [[0, 1, 3, 2, 6, 7, 9, 8],
                 [2, 3, 5, 4, 8, 9, 11, 10],
                 [6, 7, 9, 8, 12, 13, 15, 14],
                 [8, 9, 11, 10, 14, 15, 17, 16],
                 [12, 13, 15, 14, 18, 19, 21, 20],
                 [14, 15, 17, 16, 20, 21, 23, 22]], shape=[6, 8], strides=[8, 1], layout=Cc (0x5), const ndim=2,
            ),
            fields: {},
            families: [0, 0, 0, 0, 0, 0], shape=[6], strides=[1], layout=CFcf (0xf), const ndim=1,
            groups: ArcGroups(
                {},
            ),
        },
    },
}

The mesh is composed of a coordinates array, and several blocks.

volumes.to_pyvista().plot(show_edges=True)

Building mesh with custom connectivity

x, y = np.meshgrid(np.linspace(0.0, 1.0, 5), np.linspace(0.0, 1.0, 5))
coords = np.c_[x.flatten(), y.flatten()]
conn = np.array(
    [
        [0, 1],
        [1, 6],
        [6, 5],
        [5, 0],
        [6, 7],
        [7, 12],
        [12, 11],
        [11, 6],
        [12, 17],
        [17, 16],
        [16, 11],
    ],
    dtype=np.uint,
)

mesh = mf.UMesh(coords)
mesh.add_regular_block("VERTEX", np.arange(13, 22, dtype=np.uint)[..., np.newaxis])
mesh.add_regular_block("SEG2", conn)
mesh.add_regular_block("QUAD4", np.array([[3, 4, 9, 8]], dtype=np.uint))

mesh.to_pyvista(dim="all").plot(cpos="xy", show_edges=True)

Input / Output

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

volumes = mf.build_cmesh(
    range(2), np.linspace(0.0, 1.0, 5), np.logspace(0.0, 1.0, 5) / 10.0
)

Memory exports

• Through numpy arrays manipulations:
  • medcoupling
  • meshio
  • pyvista
• Through string translation to Python:
  • json

print(volumes.to_mc())

Unstructured mesh with name : "mf_UMesh"
Description of mesh : ""
Time attached to the mesh [unit] : 0 []
Iteration : -1 Order : -1
 Mesh dimension has not been set or is invalid !3
Info attached on space dimension : "" "" ""
Number of nodes : 50
Number of cells : 15
Cell types present : NORM_HEXA8

print(volumes.to_pyvista())

UnstructuredGrid (0x7cf1a2281cc0)
  N Cells:    16
  N Points:   50
  X Bounds:   0.000e+00, 1.000e+00
  Y Bounds:   0.000e+00, 1.000e+00
  Z Bounds:   1.000e-01, 1.000e+00
  N Arrays:   0

volumes.to_pyvista().plot(show_edges=True)

print(volumes.to_meshio())

<meshio mesh object>
  Number of points: 50
  Number of cells:
    hexahedron: 16

print(volumes.to_json())

{"coords":{"v":1,"dim":[50,3],"data":[0.0,0.0,0.1,1.0,0.0,0.1,0.0,0.25,0.1,1.0,0.25,0.1,0.0,0.5,0.1,1.0,0.5,0.1,0.0,0.75,0.1,1.0,0.75,0.1,0.0,1.0,0.1,1.0,1.0,0.1,0.0,0.0,0.17782794100389226,1.0,0.0,0.17782794100389226,0.0,0.25,0.17782794100389226,1.0,0.25,0.17782794100389226,0.0,0.5,0.17782794100389226,1.0,0.5,0.17782794100389226,0.0,0.75,0.17782794100389226,1.0,0.75,0.17782794100389226,0.0,1.0,0.17782794100389226,1.0,1.0,0.17782794100389226,0.0,0.0,0.31622776601683794,1.0,0.0,0.31622776601683794,0.0,0.25,0.31622776601683794,1.0,0.25,0.31622776601683794,0.0,0.5,0.31622776601683794,1.0,0.5,0.31622776601683794,0.0,0.75,0.31622776601683794,1.0,0.75,0.31622776601683794,0.0,1.0,0.31622776601683794,1.0,1.0,0.31622776601683794,0.0,0.0,0.5623413251903491,1.0,0.0,0.5623413251903491,0.0,0.25,0.5623413251903491,1.0,0.25,0.5623413251903491,0.0,0.5,0.5623413251903491,1.0,0.5,0.5623413251903491,0.0,0.75,0.5623413251903491,1.0,0.75,0.5623413251903491,0.0,1.0,0.5623413251903491,1.0,1.0,0.5623413251903491,0.0,0.0,1.0,1.0,0.0,1.0,0.0,0.25,1.0,1.0,0.25,1.0,0.0,0.5,1.0,1.0,0.5,1.0,0.0,0.75,1.0,1.0,0.75,1.0,0.0,1.0,1.0,1.0,1.0,1.0]},"element_blocks":{"HEX8":{"cell_type":"HEX8","connectivity":{"Regular":{"v":1,"dim":[16,8],"data":[0,1,3,2,10,11,13,12,2,3,5,4,12,13,15,14,4,5,7,6,14,15,17,16,6,7,9,8,16,17,19,18,10,11,13,12,20,21,23,22,12,13,15,14,22,23,25,24,14,15,17,16,24,25,27,26,16,17,19,18,26,27,29,28,20,21,23,22,30,31,33,32,22,23,25,24,32,33,35,34,24,25,27,26,34,35,37,36,26,27,29,28,36,37,39,38,30,31,33,32,40,41,43,42,32,33,35,34,42,43,45,44,34,35,37,36,44,45,47,46,36,37,39,38,46,47,49,48]}},"fields":{},"families":{"v":1,"dim":[16],"data":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]},"groups":{}}}}

File read/write

• On rust side, file I/O with the read/write methods, driven by the file extension:
  • vtk (legacy binary vtk 2.0)
  • yaml
  • json
  • vtkhdf / h5 / hdf5 (HDF5-based VTK)
  • medfile

The legacy vtk reader/writer only supports the old binary vtk 2.0 file format (no rust crate is doing better so far). The HDF5-based .vtkhdf reader/writer is the recommended option for a more modern and HPC friendly format. CGNS support is planned.

import pathlib

pathlib.Path("data").mkdir(exist_ok=True)
for ext in ("vtk", "yaml", "json", "vtkhdf", "med"):
    volumes.write(f"data/volumes.{ext}")
    volumes_from_disk = mf.UMesh.read(f"data/volumes.{ext}")
    assert volumes_from_disk
    assert (
        volumes != volumes_from_disk
    )  # this is a new instance, with a different memory adress

Mesh extrusions

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Building mesh with custom connectivity

First let’s build a 2D mesh wich will by used to demonstrate extrusions.

x, y = np.meshgrid(np.linspace(0.0, 1.0, 5), np.linspace(0.0, 1.0, 5))
coords = np.c_[x.flatten(), y.flatten()]
conn = np.array(
    [
        [0, 1],
        [1, 6],
        [6, 5],
        [5, 0],
        [6, 7],
        [7, 12],
        [12, 11],
        [11, 6],
        [12, 17],
        [17, 16],
        [16, 11],
    ],
    dtype=np.uint,
)

mesh = mf.UMesh(coords)
mesh.add_regular_block("VERTEX", np.arange(13, 22, dtype=np.uint)[..., np.newaxis])
mesh.add_regular_block("SEG2", conn)
mesh.add_regular_block("QUAD4", np.array([[3, 4, 9, 8]], dtype=np.uint))

mesh.to_pyvista(dim="all").plot(cpos="xy", show_edges=True)

Extrusion along an existing axis

Build simple extruded mesh

extruded = mesh.extrude(range(3))

extruded.to_pyvista(dim="all").plot(show_edges=True)

Build extruded mesh along a 3d line with parallel z faces

n = 50
x = np.sin(np.linspace(0.0, np.pi, n))
y = np.cos(np.linspace(0.0, np.pi, n))
z = np.linspace(0.0, 4.0, n)
line = np.c_[x, y, z]

extruded_par = mesh.extrude_parallel(line)

extruded_par.to_pyvista(dim="all").plot(show_edges=True)

Build curvilinear extrusion mesh

mesh = mf.build_cmesh(range(2), range(2))
n = 20
x = np.zeros((n,))
y = np.cos(np.linspace(0.0, np.pi, n))
z = np.sin(np.linspace(0.0, np.pi, n))
line = np.c_[x, y, z]

extruded_curv = mesh.extrude_curv(line)

extruded_curv.to_pyvista(dim="all").plot(show_edges=True)

Topological tools

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Submesh functionality

x = range(3)
y = np.linspace(0.0, 3.0, 5, endpoint=True)
z = np.logspace(-0.5, 1, 4, endpoint=True)
volumes = mf.build_cmesh(x, y, z)

Simple descending_mesh

This functionality is able to compute the descending connectivity of the provided mesh. It can act on elements of dimension 1, 2 or 3.

faces = volumes.descend()
edges = faces.descend()
vertex = edges.descend()

plotter = pv.Plotter(shape=(1, 3))
plotter.subplot(0, 0)
plotter.add_mesh(faces.to_pyvista().shrink(0.8), show_edges=True)
plotter.subplot(0, 1)
plotter.add_mesh(edges.to_pyvista().shrink(0.8))
plotter.subplot(0, 2)
plotter.add_mesh(vertex.to_pyvista())
plotter.show()

Submesh in one go

You might want to directly access either the node mesh or the edges mesh. You can ! And going into one step is ever faster than chaining multiple .descend() calls.

edges = volumes.descend(target_dim=1)
vertex = volumes.descend(target_dim=0)

plotter = pv.Plotter(shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_mesh(edges.to_pyvista().shrink(0.8))
plotter.subplot(0, 1)
plotter.add_mesh(vertex.to_pyvista())
plotter.show()

Boundaries computation

As it is very common to compute boundaries on a mesh (for boundary conditions for ex), there is a custom boundaries computation method.

face_bounds = volumes.boundaries()
edge_bounds = volumes.boundaries(target_dim=1)
vertex_bounds = volumes.boundaries(target_dim=0)

plotter = pv.Plotter(shape=(1, 3))
plotter.subplot(0, 0)
plotter.add_mesh(face_bounds.to_pyvista().shrink(0.8), show_edges=True)
plotter.subplot(0, 1)
plotter.add_mesh(edge_bounds.to_pyvista().shrink(0.8))
plotter.subplot(0, 2)
plotter.add_mesh(vertex_bounds.to_pyvista())
plotter.show()

Descend / boundaries update

You can directly update the mesh inplace when computing the descending mesh or the boundaries mesh.

volumes.boundaries_update()
volumes.boundaries_update(target_dim=1)
volumes.to_pyvista(dim="all").shrink(0.8).plot(show_edges=True)

When using the _update version, the elements of the same dimension of the generated mesh are returned as a new mesh.

old_face_mesh = volumes.descend_update()
volumes.to_pyvista(dim="all").shrink(0.8).plot(show_edges=True)

old_face_mesh.to_pyvista().shrink(0.8).plot(show_edges=True)

Connected components

x, y = np.meshgrid(np.linspace(0.0, 1.0, 5), np.linspace(0.0, 1.0, 5))
coords = np.c_[x.flatten(), y.flatten()]
conn = np.array(
    [
        [0, 1, 6, 5],
        [6, 7, 12, 11],
        # [2, 3, 8, 7],
        [11, 12, 17, 16],
    ],
    dtype=np.uint,
)
mesh = mf.UMesh(coords)
mesh.add_regular_block("QUAD4", conn)
mesh.add_regular_block("VERTEX", np.arange(len(coords), dtype=np.uint)[..., np.newaxis])

compos_link_edge = mesh.connected_components(link_dim=1)
compos_link_node = mesh.connected_components(link_dim=0)

print(f"{len(compos_link_edge)=}")
print(f"{len(compos_link_node)=}")

len(compos_link_edge)=2
len(compos_link_node)=1

edges = mesh.descend()

shape = (3, 2)
row_weights = [1.0, 0.5, 0.5]
groups = [
    (0, np.s_[:]),
    (1, 0),
    (2, 0),
    (np.s_[1:], 1),
]

plotter = pv.Plotter(shape=shape, groups=groups, row_weights=row_weights)
plotter.subplot(0, 0)
plotter.add_text("Original mesh")
plotter.add_mesh(mesh.to_pyvista(), show_edges=True)
plotter.camera_position = "xy"

for i, compo in enumerate(compos_link_edge):
    plotter.subplot(i + 1, 0)
    plotter.add_text(f"Compo linked by edge: n°{i}")
    plotter.add_mesh(edges.to_pyvista())
    plotter.add_mesh(compo.to_pyvista(), show_edges=True)
    plotter.camera_position = "xy"

for i, compo in enumerate(compos_link_node):
    plotter.subplot(i + 1, 1)
    plotter.add_text(f"Compo linked by node: n°{i}")
    plotter.add_mesh(edges.to_pyvista())
    plotter.add_mesh(compo.to_pyvista(), show_edges=True)
    plotter.camera_position = "xy"
plotter.show()

Crack

This feature is the contrary of the merge_nodes feature. It duplicates nodes such that the resulting mesh does not connect on the descending_mesh given.

x = range(2)
y = np.linspace(0.0, 3.0, 3, endpoint=True)
z = np.logspace(0.0, 1.0, 3, endpoint=True)
volumes = mf.build_cmesh(x, y, z)
faces = volumes.descend()

cracked = volumes.crack(faces)

edges = faces.descend()
compos_original = volumes.connected_components()
compos_cracked = cracked.connected_components()

assert len(compos_original) == 1

n_compos = len(compos_cracked)

shape = (3, n_compos + 1)
groups = [
    (0, 0),
    (0, np.s_[1:]),
    (np.s_[1:], 0),
    (1, np.s_[1:]),
    *((2, i + 1) for i in range(n_compos)),
]
row_weights = [1.0, 0.1, 1.0]
col_weights = [1.5, *(0.5,) * n_compos]
pv.set_jupyter_backend("static")
plotter = pv.Plotter(
    shape=shape, groups=groups, row_weights=row_weights, col_weights=col_weights
)

plotter.subplot(0, 0)
plotter.add_text("Original mesh")
plotter.add_mesh(volumes.to_pyvista(), show_edges=True)
plotter.subplot(0, 1)
plotter.add_text("Cut mesh used for the crack")
plotter.add_mesh(faces.to_pyvista().shrink(0.8), show_edges=True)

plotter.subplot(1, 0)
plotter.add_text("Compo of original mesh")
plotter.add_mesh(edges.to_pyvista())
plotter.add_mesh(compos_original[0].to_pyvista(), show_edges=True)

plotter.subplot(1, 1)
plotter.add_text("Compos of cracked mesh")

for i, compo in enumerate(compos_cracked):
    plotter.subplot(2, i + 1)
    plotter.add_mesh(edges.to_pyvista())
    plotter.add_mesh(compo.to_pyvista(), show_edges=True)
    plotter.camera.zoom(2)
plotter.show()

Split

This tool is usefull to split cells into smaller cells. It does not change the topology of the domain not the element type. Using it it gives you 2^n times the number of elements where n is the dimension of of the elements.

x = np.linspace(0.0, 3.0, 2, endpoint=True)
y = np.logspace(0.0, 1.0, 2, endpoint=True)
z = range(2)
mesh = mf.build_cmesh(x, y, z)

mesh_splitted = mesh.split()

pt = pv.Plotter()
pt.add_mesh(
    mesh.to_pyvista().shrink(0.95), show_edges=True, edge_color="yellow", line_width=2
)
pt.add_mesh(mesh_splitted.to_pyvista(), style="wireframe", color="red", line_width=2)
pt.show()

Polyze

This functionnality is useful to generate a poly mesh from a regular one. A poly mesh is a mesh of PGON 2d elements and PHED 3d elements.

x = np.linspace(0.0, 3.0, 4, endpoint=True)
y = np.logspace(0.0, 1.0, 4, endpoint=True)
mesh = mf.build_cmesh(x, y)
print(mesh.blocks())

{'QUAD4': array([[ 0,  1,  5,  4],
       [ 1,  2,  6,  5],
       [ 2,  3,  7,  6],
       [ 4,  5,  9,  8],
       [ 5,  6, 10,  9],
       [ 6,  7, 11, 10],
       [ 8,  9, 13, 12],
       [ 9, 10, 14, 13],
       [10, 11, 15, 14]], dtype=uint64)}

mesh_polyzed = mesh.polyze()

print(mesh_polyzed.blocks())

{'PGON': (array([ 0,  1,  5,  4,  1,  2,  6,  5,  2,  3,  7,  6,  4,  5,  9,  8,  5,
        6, 10,  9,  6,  7, 11, 10,  8,  9, 13, 12,  9, 10, 14, 13, 10, 11,
       15, 14], dtype=uint64), array([ 4,  8, 12, 16, 20, 24, 28, 32, 36], dtype=uint64))}

mesh_polyzed.to_pyvista().plot(show_edges=True)

unpolyzed = mesh_polyzed.unpolyze()
print(unpolyzed.blocks())
unpolyzed.to_pyvista().plot(show_edges=True)

{'QUAD4': array([[ 0,  1,  5,  4],
       [ 1,  2,  6,  5],
       [ 2,  3,  7,  6],
       [ 4,  5,  9,  8],
       [ 5,  6, 10,  9],
       [ 6,  7, 11, 10],
       [ 8,  9, 13, 12],
       [ 9, 10, 14, 13],
       [10, 11, 15, 14]], dtype=uint64)}

Geometrical tools

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Snap points

x = np.linspace(0.0, 3.0, 10, endpoint=True)
mesh = mf.build_cmesh(x, x)

eps = 0.1
dec = x[-1] / len(x) + eps
x2 = np.linspace(dec, x[-1] + dec, len(x), endpoint=True)
mesh2 = mf.build_cmesh(x2, x2)

Note: epsilon value used in the following operation is big enough so that there are multiple candidates for some points, but low enough so that there is no degenerated cell created.

snaped = mesh.snap(mesh2, eps=x[-1] / len(x))

pt = pv.Plotter()
pt.add_mesh(mesh.to_pyvista(), show_edges=True)
pt.add_mesh(mesh2.descend(target_dim=0).to_pyvista(), color="red")
pt.show(cpos="xy")

pt = pv.Plotter()
pt.add_mesh(snaped.to_pyvista(), show_edges=True)
pt.add_mesh(mesh2.descend(target_dim=0).to_pyvista(), color="red")
pt.show(cpos="xy")

Merge nodes

x = range(2)
y = np.linspace(0.0, 3.0, 3, endpoint=True)
z = np.logspace(0.0, 1.0, 3, endpoint=True)
volumes = mf.build_cmesh(x, y, z)
faces = volumes.descend()
cracked = volumes.crack(faces)

merged = cracked.merge_nodes()

edges = faces.descend()
compos_merged = merged.connected_components()
compos_cracked = cracked.connected_components()

assert len(compos_merged) == 1

n_compos = len(compos_cracked)

shape = (3, n_compos + 1)
groups = [
    (0, np.s_[:-1]),  # cracked
    (0, n_compos),  # merged
    (1, np.s_[:-1]),  # cracked txt
    (np.s_[1:], n_compos),  # merged compos
    *((2, i) for i in range(n_compos)),  # cracked compos
]
row_weights = [1.0, 0.1, 1.0]
col_weights = [*(0.5,) * n_compos, 1.5]
pv.set_jupyter_backend("static")
plotter = pv.Plotter(
    shape=shape, groups=groups, row_weights=row_weights, col_weights=col_weights
)

plotter.subplot(0, n_compos)
plotter.add_text("Merged mesh")
plotter.add_mesh(merged.to_pyvista(), show_edges=True)
plotter.subplot(0, 0)
plotter.add_text("Cut mesh used for the crack")
plotter.add_mesh(faces.to_pyvista().shrink(0.8), show_edges=True)

plotter.subplot(1, n_compos)
plotter.add_text("Compo of merged mesh")
plotter.add_mesh(edges.to_pyvista())
plotter.add_mesh(compos_merged[0].to_pyvista(), show_edges=True)

plotter.subplot(1, 0)
plotter.add_text("Compos of cracked mesh")

for i, compo in enumerate(compos_cracked):
    plotter.subplot(2, i)
    plotter.add_mesh(edges.to_pyvista())
    plotter.add_mesh(compo.to_pyvista(), show_edges=True)
    plotter.camera.zoom(2)
plotter.show()

Overlay

The intersection is valid in the following conditions :

• mesh1 and mesh2 are valid (no self recovering),
• mesh1 and mesh2 are 2d (xy),
• correctly oriented (CCW element connectivity),
• and fully merged (no unmerged nodes).

Mesh1 and mesh2 may have any number of kind of 2d elements of the first order (TRI3, QUAD4, PGON). A future version of this algorithm will work with quadratic elements (TRI7, QUAD8, QPGON).

x = np.linspace(0.0, 3.0, 4, endpoint=True)
mesh = mf.build_cmesh(x, x)

eps = 0.5
dec = x[-1] / len(x) + eps
x2 = np.linspace(dec, x[-1] + dec, len(x), endpoint=True)
mesh2 = mf.build_cmesh(x2, x2)

imprint2on1 = mesh.overlay(mesh2, mf.OverlayOperation.IMPRINT)
imprint1on2 = mesh2.overlay(mesh, mf.OverlayOperation.IMPRINT)

union = mesh.overlay(mesh2, mf.OverlayOperation.UNION)
diff = mesh.overlay(mesh2, mf.OverlayOperation.DIFFERENCE)
symdiff = mesh.overlay(mesh2, mf.OverlayOperation.SYMMETRIC_DIFFERENCE)
intersection = mesh.overlay(mesh2, mf.OverlayOperation.INTERSECTION)

meshes = [
    [mesh, mesh2],
    [union, intersection],
    [diff, symdiff],
    [imprint2on1, imprint1on2],
]
labels = [
    ["Mesh1", "Mesh2, staggered"],
    ["Union", "Intersection"],
    ["Difference", "Symmetric difference"],
    ["Imprint over 1", "Imprint over 2"],
]

pt = pv.Plotter(shape=(4, 2))
for i in range(4):
    for j in range(2):
        m = meshes[i][j]
        t = labels[i][j]
        pt.subplot(i, j)
        pt.add_text(t)
        pt.add_mesh(m.to_pyvista(), show_edges=True)
        pt.camera_position = "xy"
pt.show(cpos="xy")

The ugly cell in the center in the difference and symmetric difference comes from the plotting of non convex cells in pyvista. It is just a known plotting bug (due to optimisation quirks).

Geometric transforms

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

coords = np.array(
    [
        [0.0, 0.0],
        [1.0, 0.0],
        [0.0, 1.0],
        [1.0, 1.0],
    ]
)
mesh = mf.UMesh(coords)
mesh.add_regular_block("QUAD4", np.array([[0, 1, 3, 2]], dtype=np.uint))

All angles are given in radians. A Transform is an affine transformation represented by a 4x4 homogeneous matrix (last row 0 0 0 1). UMesh objects are transformed out-of-place: each transform returns a new mesh and leaves the source unchanged.

Out-of-place transforms

translated = mesh.translate([10.0, 0.0])
translated.to_pyvista().plot()

rotated = mesh.rotate([0.0, 0.0, 1.0], np.pi / 6)
rotated.to_pyvista().plot()

mirrored = mesh.rotate([0.0, 0.0, 1.0], np.pi / 6).mirror([1.0, 0.0, 0.0])
mirrored.to_pyvista().plot()

scaled = mesh.scale([2.0, 3.0])
scaled.to_pyvista().plot()

uniform = mesh.scale_uniform(4.0)
uniform.to_pyvista().plot()

The input mesh is never modified.

Unidirectional (left-to-right) composition

Transform supports composition in both conventions:

• a @ b (matrix product) applies b first;
• a.then(b) applies a, then b.

tr = mf.Transform.translation([1.0, 0.0]).then(
    mf.Transform.rotation([0.0, 0.0, 1.0], np.pi / 2)
)
out = mesh.transform(tr)
assert np.allclose(out.coords()[1], [0.0, 2.0], atol=1e-9)

matrix = mf.Transform.translation([1.0, 2.0, 3.0]) @ mf.Transform.scaling([2.0])
assert np.allclose(
    matrix.matrix(),
    mf.Transform.translation([1.0, 2.0, 3.0]).matrix()
    @ mf.Transform.scaling([2.0]).matrix(),
)

Rather than a Transform, mesh.transform(...) also accepts a raw 4x4 numpy array, and Transform.from_matrix(...) wraps one.

mat = np.eye(4)
mat[0, 3] = 7.0
assert np.isclose(mesh.transform(mat).coords()[0, 0], 7.0)

Dimensional consistency is enforced: trying to move a 2D (or 1D) mesh out of its plane (or line) raises a ValueError.

try:
    mesh.translate([0.0, 0.0, 0.5])
except ValueError as err:
    print(err)

This transform moves the mesh out of its plane/line: rows beyond the space dimension must leave the extra coordinates unchanged.

Duplicating a mesh

duplicate(step, n) returns n copies, each transformed by the powers step, step @ step, … of step.

ts = mf.Transform.translation([0.0, 3.0, 0.0])
rt = mf.Transform.rotation([0.0, 0.0, 1.0], np.pi / 4)
column = mesh.duplicate(ts @ rt, 3)
assert column.coords().shape[0] == 12

column.to_pyvista().plot(show_edges=True)

Arbitrary arrangements can be built with the module-level aggregate / concat functions, which concatenate meshes while preserving blocks, fields, families and groups (element ids are relabelled so that the resulting mesh stays valid).

line = mf.concat(mesh, mesh.translate([5.0, 0.0]))
three = mf.aggregate([mesh, mesh.translate([5.0, 0.0]), mesh.translate([15.0, 0.0])])
assert line.coords().shape[0] == 8
assert three.coords().shape[0] == 12
line.to_pyvista().plot()
three.to_pyvista().plot()

Transforms preserve the mesh metadata: blocks, fields, families and groups survive untouched.

Reconstructing coordinates

For non-affine coordinate changes, use UMesh.from_mesh. It accepts a same-shaped coordinate array and preserves the selected source connectivity, fields, families and groups. Omitting coords reuses the source coordinates.

coords = mesh.coords()
coords[:, 0] *= 2.0
coords[:, 1] *= coords[:, 0] + 1.0
warped = mf.UMesh.from_mesh(
    mesh,
    coords,
)
# assert np.allclose(warped.coords()[:, 0], mesh.coords()[:, 0] ** 2)
assert np.allclose(warped.coords(), coords)
warped.to_pyvista().plot()

from_mesh can select source element blocks by topological dimension or by element type. The two selectors are mutually exclusive; the source mesh is never changed.

surfaces = mf.UMesh.from_mesh(mesh, dim=2)
quads = mf.UMesh.from_mesh(mesh, element_types=["QUAD4"])

Coordinate arrays are structurally validated when a mesh is reconstructed; geometric quality checks are not performed.

Selection tool

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Element selection expressions

Elements can be selected based on their:

• types
• ids
• dimensions

Elements can be selected based on their centroid position. The selection methods using centroids are :

• bbox
• sphere
• rectangle
• circle

As you can guess, bbox and sphere should be used for 3d meshes and rectangle and circle for 2d meshes.

Elements can be selected based on the nodes position and a boolean : whether to select element with all nodes matching the condition or any node matching the condition.

• nbbox
• nsphere
• nrect
• ncircle
• nids

Elements can be selected based on their group membership :

• group(name)
• exclude_group(name)

Elements can be selected based on their scalar fields values :

• FieldExpr > FieldExpr
• FieldExpr >= FieldExpr
• FieldExpr < FieldExpr
• FieldExpr <= FieldExpr
• FieldExpr == FieldExpr

Wherever a selector is expected, the wildcards None, ... and [:] select every element — the explicit form is mf.sel.all().

sphere = mf.sel.sphere([0.5, 0.5, 0.5], 0.5)
clip_x = mf.sel.bbox([0.5, -np.inf, -np.inf], [np.inf, np.inf, np.inf])  # x > 0.5
compound_sel = ~(sphere & clip_x)

Selections are light objects, there are cheap to create and combine and independent from a support. You should just not mix 2d selectors with 3d selectors (sphere vs circle, bbox vs rectangle).

How does it works ?

Each objects generated by a selection function (a function from the mf.sel module) is of the Selection type. It implements operators so that it knows how to compose in an expression. That way an expression generates a new SelectionExpr which can be interpreted by the .select(expr) method.

print(sphere)

CentroidSelection(
    Sphere {
        center: [
            0.5,
            0.5,
            0.5,
        ],
        r: 0.5,
    },
)

print(compound_sel)

NotExpr(
    NotExpr(
        BinarayExpr(
            BinarayExpr {
                operator: And,
                left: CentroidSelection(
                    Sphere {
                        center: [
                            0.5,
                            0.5,
                            0.5,
                        ],
                        r: 0.5,
                    },
                ),
                right: CentroidSelection(
                    BBox {
                        min: [
                            0.5,
                            -inf,
                            -inf,
                        ],
                        max: [
                            inf,
                            inf,
                            inf,
                        ],
                    },
                ),
            },
        ),
    ),
)

You can see two layers of NotExpr and BinaryExpr. That is perfectly normal, it does not mean that the operation is applied twice, both NotExpr operations and both BinaryExpr op actually comes from different namespaces and it is just a form of encapsulation (first is a variant, second is an enum).

Select to extract part of a mesh

Here the types manipulated are Selection. They are simple objects that know how to compose themselves. When applied on a mesh the mesh selection is computed.

clip = mf.sel.bbox([-np.inf, -np.inf, -np.inf], [np.inf, np.inf, 0.5])
sphere = mf.sel.sphere([0.5, 0.5, 0.5], 0.5)

x = np.linspace(0.0, 1.0, 20, endpoint=True)
volumes = mf.build_cmesh(x, x, x)
volumes.to_pyvista().plot(show_edges=True)

volumes.select(clip).to_mesh().to_pyvista().plot()

volumes.select(sphere).to_mesh().to_pyvista().plot()

Selection composition

One of the great strength of the select method is its composability ! Watch by yourself.

The operators &, |, ^, - and ~ are available.

2D composition examples

x = np.linspace(0.0, 2.0, 100)
y = np.linspace(0.0, 1.0, 50)
faces = mf.build_cmesh(x, y)

circle1 = mf.sel.circle([0.75, 0.5], 0.5)
circle2 = mf.sel.circle([1.25, 0.5], 0.5)

union = faces.select(circle1 | circle2).to_mesh()

pt = pv.Plotter()
pt.add_mesh(faces.descend().to_pyvista())
pt.add_mesh(union.to_pyvista())
pt.camera_position = "xy"
pt.show()

intersection = faces.select(circle1 & circle2).to_mesh()

pt = pv.Plotter()
pt.add_mesh(faces.descend().to_pyvista())
pt.add_mesh(intersection.to_pyvista())
pt.camera_position = "xy"
pt.show()

sym_diff = faces.select(circle1 ^ circle2).to_mesh()

pt = pv.Plotter()
pt.add_mesh(faces.descend().to_pyvista())
pt.add_mesh(sym_diff.to_pyvista())
pt.camera_position = "xy"
pt.show()

diff = faces.select(circle1 - circle2).to_mesh()

pt = pv.Plotter()
pt.add_mesh(faces.descend().to_pyvista())
pt.add_mesh(diff.to_pyvista())
pt.camera_position = "xy"
pt.show()

notsel = faces.select(~circle1).to_mesh()

pt = pv.Plotter()
pt.add_mesh(faces.descend().to_pyvista())
pt.add_mesh(notsel.to_pyvista())
pt.camera_position = "xy"
pt.show()

A 3D complex example

sphere = mf.sel.sphere([0.5, 0.5, 0.5], 0.5)
clip_x = mf.sel.bbox([0.5, -np.inf, -np.inf], [np.inf, np.inf, np.inf])  # x > 0.5
clip_z = mf.sel.bbox([-np.inf, -np.inf, -np.inf], [np.inf, np.inf, 0.5])  # z < 0.5

volumes.select(
    (clip_x & sphere & clip_z) | (sphere & ~clip_x & ~clip_z)
).to_mesh().to_pyvista().plot()

Select API

.select() returns a lazy view: ids(), len(), to_mesh() and reductions evaluate it on demand.

volumes.select(sphere)

SelectionResult(n_elements=3695)

The number of elements referenced here is the number of elements of volumes. The selection was not yet computed.

two_quarters_expr = (clip_x & sphere & clip_z) | (sphere & ~clip_x & ~clip_z)

volumes.select(two_quarters_expr).ids()

{'HEX8': array([ 143,  144,  162, ..., 6715, 6716, 6735],
       shape=(1838,), dtype=uint64)}

len(volumes.select(two_quarters_expr))

1838

Selection to reduction

One might want to compute a reduction (min, max, mean, std, etc) on part of a mesh. This can be done applying a reduction and using a field expression (see fields section).

volumes.select(two_quarters_expr).mean(mf.X)

0.5171525113109214

Groups

Groups API

Selections can be stored on the mesh as named groups. The mesh.groups mapping behaves like a dict: assign a selection expression (or a {etype: ids} dict) to create or replace a group, then manage groups with the usual operations.

volumes.groups["two_quarters"] = two_quarters_expr

volumes.groups

GroupsMapping(["two_quarters"])

volumes.groups["two_quarters"].ids()

{'HEX8': array([ 143,  144,  162, ..., 6715, 6716, 6735],
       shape=(1838,), dtype=uint64)}

tq = volumes.groups["two_quarters"].to_mesh()
tq.to_pyvista().plot()

len(volumes.groups["two_quarters"])

1838

volumes.groups.rename("two_quarters", "quarter_tag")
volumes.groups

GroupsMapping(["quarter_tag"])

del volumes.groups["quarter_tag"]
volumes.groups

GroupsMapping([])

Groups inplace modifications

Groups can be modified inplace, either using SelectionExpr (including existing groups expr):

volumes.groups["two_quarters"] = two_quarters_expr
n = len(volumes.groups["two_quarters"])
added_sel = mf.sel.bbox([-np.inf, -np.inf, -np.inf], [np.inf, np.inf, 0.2])
volumes.groups["two_quarters"].add(added_sel)
print(len(volumes.groups["two_quarters"]), "after adding a slab (was", n, ")")

3089 after adding a slab (was 1838 )

volumes.groups["two_quarters"].to_mesh().to_pyvista().plot()

… or directly with element ids per element type.

print(len(volumes.groups["two_quarters"]))
volumes.groups["two_quarters"].remove({"HEX8": [0, 1]})
print(len(volumes.groups["two_quarters"]))
volumes.groups["two_quarters"].add({"HEX8": [0, 1]})
print(len(volumes.groups["two_quarters"]), "back to the original size")

3089
3087
3089 back to the original size

Fields

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Field expressions

FieldExpr are composition of floats and mf.sel.field(“fieldname”) or custom fields. Available operations on fields include :

• binary expresions + - / *
• unary expr sin(), cos(), abs(), ln(), log10(), exp()
• primitives :
  • M: the measure, ie length/area/volume of an element
  • C: the node centroid of an element, not the volume barycenter
  • X: the X compo of the node centroid
  • Y: the Y compo of the node centroid
  • Z: the Z compo of the node centroid

Scalar binary operations

toto = mf.Field("toto")
tata = mf.Field("tata")
toto + tata
toto * tata
toto - tata
toto / tata
toto + 2.0
toto - 2.0
toto * 2.0

Scalar unary ops

toto.sin() toto.cos() toto.abs() toto.exp() toto.ln() toto.square() toto.sqrt() toto.tan() toto.log10();

Vector ops

toto.dot(tata)
toto @ tata  # same as dot
toto[0]

Primitives

m = mf.M  # measure: length/area/volume
c = mf.C  # node centroids
x = mf.X  # x compo of node centroid
y = mf.Y  # y compo of node centroid
z = mf.Z  # z compo of node centroid
n = mf.Normal  # normal to n-1 dim elements, ie 2d elems in 3d or 1d elems in 2d
nx = mf.Nx  # x compo of element normal
ny = mf.Ny  # y compo of element normal
nz = mf.Nz  # z compo of element normal

How does it work ?

The operations build a binary operation tree structure. Mefikit knows how to interpret this binary tree to compute fields.

print(toto * mf.M + 3.0 * mf.X)

BinaryExpr {
    operator: Add,
    left: BinaryExpr {
        operator: Mul,
        left: Field(
            "toto",
        ),
        right: Measure,
    },
    right: BinaryExpr {
        operator: Mul,
        left: Array(
            3.0, shape=[], strides=[], layout=CFcf (0xf), dynamic ndim=0,
        ),
        right: X,
    },
}

This is quite handfull because it enables two patterns:

• reusability and composition of filters
• evaluation optimizations of selections : some selection filters are evaluated in parallel, some are evaluated first if they are discriminant

Mesh fields mapping

x = np.logspace(-5, 0.0, 50)
z = np.linspace(0.0, 0.1, 3)
mesh2 = mf.build_cmesh(x, x, z)

mesh2.to_pyvista().plot(show_edges=True)

Fields attribute is dictionnary like: fields can be accessed, modified, added, defined through it using field expressions evaluation on the mesh.

Fields expressions are independent from the mesh and light, fields are evaluated field expressions stored alongside the mesh.

mesh2.fields["Measure"] = mf.M

mesh2.to_pyvista().plot()

mesh2.fields["toto"] = mf.X + mf.Y

pvm = mesh2.to_pyvista()
pvm.active_scalars_name = "toto"
pvm.plot()

# List and look up fields by name.
print(mesh2.fields.keys())
print(mesh2.fields.values())

['Measure', 'toto']
[FieldRef("Measure"), FieldRef("toto")]

for n, f in mesh2.fields.items():
    print(n, ":", f)

Measure : FieldRef("Measure")
toto : FieldRef("toto")

del mesh2.fields["toto"]  # remove it
print(mesh2.fields)

FieldsMapping(["Measure"])

Field references

Fields live in a dict-like mapping on the mesh, keyed by name. Each entry is a handle (FieldRef) to read values, reduce them, or write through selectors. Fields reference are always bound to a given mesh.

mes = mesh2.fields["Measure"]
print("shape:", mes.shape, "| elements:", len(mes))

shape: (1,) | elements: 4802

Whole reductions

Field references support reductions evaluated eagerly :

# Whole-domain reductions over every element carrying the field.
print(mes.min(), mes.max(), mes.mean())

3.507414294773286e-13 0.002192327517292864 2.0824239902124112e-05

Numpy input/ouput

# Bulk export as {etype: array} (or a single array via `.numpy()` when the
# mesh has one element type).
vals = mes.values()
print(vals.keys())
shortcut_vals = mes.numpy()
print(shortcut_vals.shape)
assert np.allclose(vals["HEX8"], shortcut_vals)

dict_keys(['HEX8'])
(4802,)

# Bulk import a dict[str, ndarray] as new field
# Field size checks are done so that the field lay on all elements of the same dim.
mesh2.fields["toto"] = {"HEX8": shortcut_vals * 3.0}

Regional filtered fields

Lazy selections (see selection notebook) allow to compute regional reduction with any field expression, including plain existing field names as strings.

rect = mf.sel.bbox([0.25, 0.25, 0.0], [0.7, 0.7, 0.1])
zone = mesh2.select(rect)
print(zone.mean("toto"), zone.max(mf.M * 4))

0.0013598121884658915 0.003426116772966117

Writes accept scalars, arrays, field expressions or existing field names, targeted by wildcards (...) or selectors.

mesh2.fields["Scratch"] = 0.0  # create by broadcast
mesh2.fields["Scratch"][...] = (
    "Measure"  # whole selection, copy an existing field inplace
)

A field can be overwritten on a specific region. The overwrite is done using an expression formulae which can even reference the previous field values.

sel = mf.sel.bbox([0.0, 0.0, 0.0], [0.3, 1.0, 0.1])
mesh2.fields["Scratch"][sel] = mf.Field("Scratch") * 2  # scaled sub-region

sel2 = mf.sel.sphere(center=[0.5, 0.5, 0.05], r=0.3)
m = mesh2.select(sel2).mean("Scratch")  # compute the mean of scratch in a region
mesh2.fields["Scratch"][sel2] = m  # assign this constant value to the whole region

pvm = mesh2.to_pyvista()
pvm.active_scalars_name = "Scratch"
pvm.plot()

Direct field expression evaluation to numpy

It is not really recommended not to use the .fields storing mecanism as it provides complete integration with mefikit, but it is nevertheless possible to evaluate an expression on a field and export it directly as a numpy array. The eval method does exaclty this.

m = mf.Field("Measure")
m2 = mf.Field("4 * M2")
mesh2.fields["4 * M2"] = 4.0 * m * m
mesh2.eval(m2 - 4.0 * mf.M.square())

{'HEX8': array([0., 0., 0., ..., 0., 0., 0.], shape=(4802,))}

Field to Selection

Field expressions can be converted to threshold selections expressions. The available comparisons are <, <=, >, >=, == :

maxM = mesh2.select(mf.sel.all()).max(mf.M)
minM = mesh2.select(mf.sel.all()).min(mf.M)
meanM = mesh2.select(mf.sel.all()).mean(mf.M)
lb = (minM + meanM) / 2.0
hb = (maxM + meanM) / 2.0
m = mf.Field("Measure")

th = (m > lb) & (m <= hb)

m2sel = mesh2.select(th).to_mesh()
pvm2: pv.UnstructuredGrid = m2sel.to_pyvista()
pvm2.active_scalars_name = "Measure"
pvm2.plot()

Those threshold selections can be combined with other selections.

r = mf.sel.bbox([0.5, 0.5, 0.0], [0.8, 0.8, 0.1])
c = mf.sel.sphere([0.12, 0.12, 0.05], 0.05)

mesh2.select(th - r - c).to_mesh().to_pyvista().plot()

Vector/Matrix/Tensor fields

Normals of hyperplane dim elements

Let first add the elements of the boundaries to the current mesh :

mesh2.boundaries_update()

Now normals can be computed on those elements.

mesh2.fields["N"] = mf.Normal
mesh2.fields["N"].values()

{'QUAD4': array([[ 0.,  0., -1.],
        [ 0.,  0., -1.],
        [ 0.,  0.,  1.],
        ...,
        [ 0.,  0.,  1.],
        [ 0.,  0.,  1.],
        [ 0.,  0.,  1.]], shape=(5194, 3))}

Dot product / matrix multiplication is available through a numpy like syntax with the .dot operator or the @ operator :

ev = mesh2.eval(mf.Normal @ np.array([0.0, 1.0, 0.0]))["QUAD4"]
(ev > 0.5).sum()  # number of faces oriented towards y

np.int64(98)

mesh2.groups["top"] = mf.Nz > 0.9
top = mesh2.select(mf.sel.group("top")).to_mesh()

pt = pv.Plotter()
pt.add_mesh(mesh2.descend(target_dim=1).to_pyvista())
pt.add_mesh(top.to_pyvista(), show_edges=True)
pt.show()

Field transfers

This notebook shows how to remap fields between meshes with the mf.transfer operators: interpolation, extrapolation and conservative remapping, all sharing the same prepare / apply split.

import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

# Same mesh and selection context as the Fields notebook.
x = np.logspace(-5, 0.0, 1000)
mesh2 = mf.build_cmesh(x, x)

mesh2.fields["Measure"] = mf.M

m = mf.Field("Measure")
m2 = mf.Field("4 * M2")
mesh2.fields["4 * M2"] = m * m * 4.0

r = mf.sel.rect([0.25, 0.25], [0.7, 0.7])
c = mf.sel.circle([0.875, 0.875], 0.05)

Transferring fields

The transfer function can be

• interpolation
• extrapolation
• conservative
• non-conservative
• using cells
• using cell centers and point clouds methods
• etc

They are many.

ConstantPiecewise Transfer

The transfer is very simple. It is based on the cells of src_mesh and the cells center of the target mesh. It assigns to a cell from the target the value of the cell in which the center is located in. This is a point location based value assignment. By default the centroid (mean of cell nodes) is used because it is fast to comupute and is accurate with regular cells.

m_src = mesh2.select((m2 > 4e-9) - r - c).to_mesh()
m_tgt = mf.build_cmesh(np.linspace(0.0, 1.5, 20), np.linspace(0.0, 1.5, 20))

# The transfer is computed between source and target geometry.
# This step is computationnaly heavy, but done once.
cpt = mf.transfer.ConstantPiecewise(m_src, m_tgt)

# The transfer is applied. This step is much faster.
cpt.apply_update(m_src, "Measure", m_tgt, tgt_field_name="Projection", def_val=np.nan)

m_tgt.to_pyvista().plot(show_edges=True)

As you can see the "Measure" field from m_src was used to compute the "Projection" field on m_tgt. Both mesh are not completly overlapping but that is not an issue. Cells from m_tgt whose center is not in a cell from m_src take a default value def_val. Default is 0.0 but any floating point value, such as np.nan is accepted.

This interpolation is good when coarseing a mesh and you do not need conservation. It might be useful in other circumstances I do not know of. It is quite fast but not that much because of the is_in_cell exact geometrical query.

MovingMean Transfer

This Transfer is based on m_src cell center positions and m_tgt cell centers positions. It is a “meshless” operation as it does not care about connectivity. There are several options :

• normal mean
• weighted mean

Pros :

• it is extremly fast to compute
• it does not overshoot / undershoot

Cons :

• It lacks precision

MovingLeastSquare Transfer

This Transfer is based on m_src cell center positions and m_tgt cell centers positions. It is a “meshless” operation as it does not care about connectivity. There are several options :

• linear least square : the projection is the least square linear approx of the solution (can be an extrapolation)
• weighted least square : the projection is the weighted least square linear approx of the solution, there are several possibilities for the weighting function but it depends on the relative distance to the target interpolation point.

m_src = mesh2.select((m2 > 4e-9) - r - c).to_mesh()
m_tgt = mf.build_cmesh(np.linspace(0.0, 1.5, 20), np.linspace(0.0, 1.5, 20))

# The transfer is computed between source and target geometry.
# This step is computationnaly heavy, but done once.
mlsqt = mf.transfer.MovingLeastSquares(m_src, m_tgt, k=10)

# The transfer is applied. This step is much faster.
mlsqt.apply_update(m_src, "Measure", m_tgt, tgt_field_name="Projection", def_val=np.nan)

m_tgt.to_pyvista().plot(show_edges=True)

As you can see the projection gets a value everywhere, even outside the initial domain. This can lead to bad extrapolations, so it is to your responsability. Here as an example the extrapolated Measure can have negative values.

Inside the domain the interpolation works like a charm.

Transfer methods comparison

def compare_src_tgt(m_src, m_tgt):
    pt = pv.Plotter(shape=(1, 2))
    pt.subplot(0, 0)
    pt.add_text("Source")
    pt.add_mesh(m_src.to_pyvista(), show_edges=True, clim=[0.0, 0.06])
    pt.camera_position = "xy"
    pt.subplot(0, 1)
    pt.add_text("Target")
    pt.add_mesh(
        m_tgt.to_pyvista(), clim=[0.0, 0.06], below_color="pink", above_color="red"
    )
    pt.add_mesh(m_src.descend().to_pyvista(), show_edges=True, line_width=1)
    pt.camera_position = "xy"
    pt.show()

import time

transfers = (
    mf.transfer.ConstantPiecewise,
    lambda src, tgt: mf.transfer.MovingLeastSquares(src, tgt, k=5),
    mf.transfer.MovingLeastSquares,
    lambda src, tgt: mf.transfer.MovingLeastSquares(src, tgt, k=20),
    mf.transfer.MovingLeastSquares,
    lambda src, tgt: mf.transfer.MovingLeastSquares(
        src, tgt, weighting=mf.transfer.DistanceWeighting.Gaussian()
    ),
    lambda src, tgt: mf.transfer.MovingLeastSquares(
        src, tgt, weighting=mf.transfer.DistanceWeighting.InverseDistance(1.0)
    ),
    lambda src, tgt: mf.transfer.InverseDistance(src, tgt, k=3),
    lambda src, tgt: mf.transfer.InverseDistance(src, tgt, k=5),
    lambda src, tgt: mf.transfer.InverseDistance(src, tgt, k=10),
    mf.transfer.ConservativeP0,
)
trasfers_labels = (
    "CPW",
    "MLS k5",
    "MLS k10",
    "MLS k20",
    "MLS",
    "MLS gaussian",
    "MLS inv_dist",
    "ID k3",
    "ID k5",
    "ID k10",
    "ConservativeP0",
)
prepare_times = []
apply_times = []

for T, label in zip(transfers, trasfers_labels):
    m_src = mf.build_cmesh(np.logspace(-2.0, 0.0, 20), np.logspace(-2.0, 0.0, 20))
    m_tgt = mf.build_cmesh(np.linspace(-0.05, 1.1, 40), np.linspace(-0.05, 1.1, 40))
    m_src.fields["Measure"] = mf.M
    t0 = time.time()
    tr = T(m_src, m_tgt)
    t1 = time.time()
    tr.apply_update(m_src, "Measure", m_tgt, label + " Transfered Measure")
    t2 = time.time()

    prepare_times.append((t1 - t0) * 1000.0)
    apply_times.append((t2 - t1) * 1000.0)

    compare_src_tgt(m_src, m_tgt)

import matplotlib.pyplot as plt

chart_data = {
    "Prepare": prepare_times,
    "Apply": apply_times,
}

fig, ax = plt.subplots(figsize=(10, 5))

res = ax.grouped_bar(chart_data, tick_labels=trasfers_labels, group_spacing=1)
for container in res.bar_containers:
    ax.bar_label(container, padding=3)

# Add some text for labels, title, etc.
ax.set_ylabel("Time (ms)")
ax.set_title("Time per step")
ax.legend(loc="upper left", ncols=3)
fig.tight_layout()
plt.show()

Bubbles example

import numpy as np
import pyvista as pv

import mefikit as mf

rng = np.random.default_rng(seed=123)
pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

Setup

xmax = 5.0
ymax = 1.0
r = 0.17
nb = 15
nr = 12.5

nx = int(xmax / r * nr)
ny = int(ymax / r * nr)
print(f"Number of elements : {nx * ny * ny:,}")

Number of elements : 1,955,743

xc = rng.uniform(r, xmax - r, nb)
yc = rng.uniform(r, ymax - r, nb)
zc = rng.uniform(r, ymax - r, nb)
spheres = [mf.sel.sphere([x, y, z], r) for x, y, z in zip(xc, yc, zc)]
sphere_union = spheres[0]
for s in spheres[1:]:
    sphere_union = sphere_union | s

x = np.linspace(0.0, xmax, nx)
y = np.linspace(0.0, ymax, ny)
volumes = mf.build_cmesh(x, y, y)

volumes.boundaries().to_pyvista().plot(opacity=0.4)

Selecting bubbles

# `select` returns a lazy view: materialize it into a sub-mesh with `.to_mesh()`.
inner_bubbles = volumes.select(sphere_union).to_mesh()
interface = inner_bubbles.boundaries()
cracked = volumes.crack(interface)

# Named groups live in a dict-like mapping on the mesh:
# assign any selection expression (or {etype: ids} dict) to tag elements.
volumes.groups["bubbles"] = sphere_union
print(len(volumes.groups["bubbles"]), "elements tagged in group 'bubbles'")

105990 elements tagged in group 'bubbles'

Cracking and connected components

cracked.boundaries().to_pyvista().plot(opacity=0.4)

bubble_groups = inner_bubbles.connected_components()

pv.global_theme.color_cycler = "default"
pl = pv.Plotter()
for c in bubble_groups:
    compo = c.to_pyvista()
    pl.add_mesh(compo)
pl.add_mesh(volumes.boundaries(target_dim=1).to_pyvista())
pl.show()
pv.global_theme.color_cycler = None

clip1 = mf.sel.bbox([-np.inf] * 3, [np.inf, ymax / 3.0, np.inf])
pl = pv.Plotter()
pl.add_mesh(volumes.select(clip1 & ~sphere_union).to_mesh().to_pyvista())
pl.add_mesh(interface.to_pyvista(), opacity=0.4)
pl.show()

Computing statistics

bubble_volumes = volumes.select("bubbles").sum(mf.M)
print(nb * 4.0 / 3.0 * np.pi * r**3.0)
bubble_volumes

0.3086928941417331





0.2793115007084435

bubbles_mean_pos = volumes.select("bubbles").mean(mf.C)
print(bubbles_mean_pos)

[2.49923348 0.4928932  0.5314482 ]

mefikit vs. medcoupling

A friendly, side-by-side look at two libraries that solve the same class of problems, in slightly different ways.

medcoupling is the reference implementation of the MED world: mature, broad and battle-tested in industrial workflows. It set the bar for mesh interchange — mefikit itself happily reads and writes the very same .med files.

mefikit is a younger library, written in Rust and exposed to Python. It pursues the same goal — manipulate unstructured meshes and their fields — with a shorter, more expressive API and a performance-oriented core.

This notebook is neither a verdict nor a “rip and replace”. It is an invitation: we build the same meshes, run the same operations, and compare the interface and the timings step by step. If you know medcoupling, you should find your bearings immediately — and, hopefully, be tempted.

Every mesh below is built from the exact same geometry on both sides, and every value is cross-checked against the other library before we say one word about speed.

import os
import time
from pathlib import Path

import matplotlib.pyplot as plt
import medcoupling as mc
import numpy as np
import pyvista as pv

import mefikit as mf

pv.set_plot_theme("dark")
pv.set_jupyter_backend("static")

print("mefikit     :", mf.__file__)
print("medcoupling :", mc.__file__)
print("medcoupling version:", mc.__version__)
print("numpy       :", np.__version__)

mefikit     : /home/asonolet/Codes/mefikit/src/mefikit/__init__.py
medcoupling : /home/asonolet/Codes/mefikit/.venv/lib/python3.13/site-packages/medcoupling.py
medcoupling version: V9_15_0
numpy       : 2.5.2

A couple of plumbing helpers

Nothing to learn here — two small utilities used to cross-check node counts and surface areas between the libraries below.

def used_nodes(mesh):
    ids = np.concatenate([np.asarray(b) for b in mesh.blocks().values()])
    return int(np.unique(ids).size)


def area_2d(mesh):
    return float(sum(np.asarray(v).sum() for v in mesh.measure().values()))

Building the same mesh

Every mesh workflow starts with creating a grid. In mefikit, build_cmesh builds a structured (cartesian) mesh from the coordinate axes in one call.

In medcoupling the standard route goes through a structured MEDCouplingCMesh, materialised as an unstructured mesh with buildUnstructured(). Three lines, still very readable.

The good news: the result is the same mesh — same nodes, same cells — because a cartesian grid is just a very regular unstructured mesh.

def medcoupling_cmesh(*axes):
    # medcoupling counterpart of mf.build_cmesh(*axes): build a structured
    # MEDCouplingCMesh, then materialise it as an unstructured mesh with the
    # standard buildUnstructured() call.
    cmesh = mc.MEDCouplingCMesh()
    cmesh.setCoords(
        *[mc.DataArrayDouble(np.ascontiguousarray(a, np.float64)) for a in axes]
    )
    umesh = cmesh.buildUnstructured()
    umesh.setMeshDimension(len(axes))
    return umesh


x = np.linspace(0.0, 1.0, 6)
mf_mesh = mf.build_cmesh(x, x)  # mefikit    : one call
mc_mesh = medcoupling_cmesh(x, x)  # medcoupling: structured -> unstructured

print(
    "mefikit     :",
    mf_mesh.num_elements(),
    "QUAD4 cells,",
    mf_mesh.coords().shape[0],
    "nodes",
)
print(
    "medcoupling :",
    mc_mesh.getNumberOfCells(),
    "QUAD4 cells,",
    mc_mesh.getNumberOfNodes(),
    "nodes",
)

assert mc_mesh.getNumberOfCells() == mf_mesh.num_elements()
assert np.allclose(mc_mesh.getCoords().toNumPyArray(), mf_mesh.coords())
print("identical geometry (nodes and cells): OK")

mefikit     : 25 QUAD4 cells, 36 nodes
medcoupling : 25 QUAD4 cells, 36 nodes
identical geometry (nodes and cells): OK

# The bridge in the other direction is a one-liner: mefikit meshes
# export to medcoupling with to_mc(), geometry untouched.
twin = mf_mesh.to_mc()
twin.setMeshDimension(2)
print(
    "to_mc() bridge:",
    twin.getNumberOfCells(),
    "cells,",
    twin.getNumberOfNodes(),
    "nodes",
)

to_mc() bridge: 25 cells, 36 nodes

def mc_to_pyvista(mmesh):
    # Render a medcoupling mesh with pyvista, in-memory (no temp files).
    vtk_type = {
        mc.NORM_SEG2: pv.CellType.LINE,
        mc.NORM_TRI3: pv.CellType.TRIANGLE,
        mc.NORM_QUAD4: pv.CellType.QUAD,
        mc.NORM_POLYGON: pv.CellType.POLYGON,
        mc.NORM_HEXA8: pv.CellType.HEXAHEDRON,
    }
    n = mmesh.getNumberOfCells()
    coords = mmesh.getCoords().toNumPyArray()
    if coords.shape[1] == 2:
        coords = np.c_[coords, np.zeros(len(coords))]
    conn = mmesh.getNodalConnectivity().toNumPyArray()
    off = mmesh.getNodalConnectivityIndex().toNumPyArray()
    cells_arr = np.concatenate(
        [
            np.r_[off[i + 1] - off[i] - 1, conn[off[i] + 1 : off[i + 1]]]
            for i in range(n)
        ]
    ).astype(np.int64)
    cell_types = np.array(
        [vtk_type[mmesh.getTypeOfCell(i)] for i in range(n)], np.uint8
    )
    return pv.UnstructuredGrid(cells_arr, cell_types, coords)


pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("mefikit - mf.build_cmesh(x, x)")
pt.add_mesh(mf_mesh.to_pyvista(), show_edges=True)
pt.camera_position = "xy"
pt.subplot(0, 1)
pt.add_text("medcoupling - MEDCouplingCMesh.buildUnstructured()")
pt.add_mesh(mc_to_pyvista(mc_mesh), show_edges=True)
pt.camera_position = "xy"
pt.show()

Fields, measures, and one-line expressions

Everyone needs per-cell measures in remapping workflows. mefikit exposes a symbolic field mf.M — the measure — evaluated on demand, together with a small expression DSL on mf.Field: you write the computation, mefikit evaluates it on the mesh, and combinations of select / mean / sum become one-liners.

medcoupling takes the more explicit route: you build a DataArrayDouble, wrap it into a MEDCouplingFieldDouble, and attach it to the mesh. Both compute exactly the same numbers.

# --- mefikit ----------------------------------------------------------------
mf_mesh.fields["Measure"] = mf.M  # symbolic measure
mf_mesh.fields["T"] = 1.0 + mf.X**2 + 0.5 * mf.Y

print("measure sum  :", mf_mesh.fields["Measure"].sum())
hot = mf_mesh.select(mf.Field("T") > 1.5)
print("cells T > 1.5:", len(hot))
print("mean T above :", hot.mean("T"))

measure sum  : 1.0
cells T > 1.5: 13
mean T above : 1.8338461538461535

# --- medcoupling -------------------------------------------------------------
mc_measure = mc_mesh.getMeasureField(True)
mc_measure.setNature(mc.IntensiveConservation)

centers = np.asarray(mc_mesh.computeCellCenterOfMass().toNumPyArray())
mc_T_vals = 1.0 + centers[:, 0] ** 2 + 0.5 * centers[:, 1]

mc_T = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
_arr = mc.DataArrayDouble(np.ascontiguousarray(mc_T_vals, np.float64))
_arr.setName("T")
mc_T.setArray(_arr)
mc_T.setMesh(mc_mesh)
mc_T.setNature(mc.IntensiveConservation)

print("measure sum  :", float(mc_measure.getArray().toNumPyArray().sum()))

measure sum  : 1.0

# identical numbers, whichever library computed them
mf_T_vals = np.asarray(mf_mesh.fields["T"].numpy()).ravel()
mc_T_vals = mc_T.getArray().toNumPyArray()
print("max |mf_T - mc_T|:", float(np.abs(mf_T_vals - mc_T_vals).max()))
assert np.allclose(mf_T_vals, mc_T_vals, atol=1e-12)

max |mf_T - mc_T|: 1.3322676295501878e-15

pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("mefikit - mf.Field / select / eval")
g1 = mf_mesh.to_pyvista()
g1["T"] = mf_T_vals
pt.add_mesh(g1, scalars="T", show_edges=True)
pt.camera_position = "xy"
pt.subplot(0, 1)
pt.add_text("medcoupling - MEDCouplingFieldDouble")
g2 = mc_to_pyvista(mc_mesh)
g2["T"] = mc_T_vals
pt.add_mesh(g2, scalars="T", show_edges=True)
pt.camera_position = "xy"
pt.show()

Common operations, side by side

Descending connectivity

The faces of a volume mesh — the bread and butter of boundary conditions and surface integrals. mefikit: one call, descend(). medcoupling: buildDescendingConnectivity(). Both report the same face count.

axes3 = [np.linspace(0.0, 1.0, 5)] * 3
m3 = mf.build_cmesh(*axes3)
mc3 = medcoupling_cmesh(*axes3)

mf_faces = m3.descend()
mc_desc = mc3.buildDescendingConnectivity()
mc_faces = mc_desc[0] if isinstance(mc_desc, tuple) else mc_desc

print("faces (mefikit)    :", mf_faces.num_elements())
print("faces (medcoupling):", mc_faces.getNumberOfCells())
assert mf_faces.num_elements() == mc_faces.getNumberOfCells()
print("identical face count: OK")

faces (mefikit)    : 240
faces (medcoupling): 240
identical face count: OK

pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("descend() - faces in black")
pt.add_mesh(m3.to_pyvista(), show_edges=True, opacity=0.35)
pt.add_mesh(mf_faces.to_pyvista(), show_edges=True, color="black", line_width=2)
pt.subplot(0, 1)
pt.add_text("buildDescendingConnectivity()")
pt.add_mesh(mc_to_pyvista(mc3), show_edges=True, opacity=0.35)
pt.add_mesh(mc_to_pyvista(mc_faces), show_edges=True, color="black", line_width=2)
pt.show()

Merging duplicated nodes

Meshes often carry duplicated nodes (split or intersected interfaces). The reference cure is a node merge. mefikit’s merge_nodes() and medcoupling’s mergeNodes() collapse the same duplicates.

A small nuance worth knowing: mefikit keeps the coordinate array untouched and rewires the connectivity — cheap and zero-copy friendly — while medcoupling physically compacts nodes. The number of used nodes is identical.

# crack a small hex stack so its shared interface is duplicated
volumes = mf.build_cmesh([0.0, 1.0], np.linspace(0.0, 1.0, 5), np.linspace(0.0, 1.0, 5))
faces = volumes.descend()
cracked = volumes.crack(faces)
cracked_mc = cracked.to_mc()
cracked_mc.setMeshDimension(3)

merged = cracked.merge_nodes()
merged_mc = cracked_mc.deepCopy()
merged_mc.mergeNodes(1e-12)

print("components before merge:", len(cracked.connected_components()))
print("components after merge :", len(merged.connected_components()))
print(
    "used nodes  before/after (mefikit)    :",
    used_nodes(cracked),
    "->",
    used_nodes(merged),
)
print(
    "nodes       before/after (medcoupling):",
    cracked_mc.getNumberOfNodes(),
    "->",
    merged_mc.getNumberOfNodes(),
)
assert len(merged.connected_components()) == 1
assert used_nodes(merged) == merged_mc.getNumberOfNodes()
print("identical result: OK")

components before merge: 16
components after merge : 1
used nodes  before/after (mefikit)    : 128 -> 50
nodes       before/after (medcoupling): 128 -> 50
identical result: OK

pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text(f"cracked - {len(cracked.connected_components())} components")
for compo in cracked.connected_components():
    pt.add_mesh(compo.to_pyvista(), show_edges=True, opacity=0.6)
pt.subplot(0, 1)
pt.add_text("merge_nodes() - 1 component")
pt.add_mesh(merged.to_pyvista(), show_edges=True)
pt.show()

2D overlay / imprint

Boolean overlay of two 2D meshes (insert an embedded mesh into a background one). mefikit: overlay(). medcoupling: Intersect2DMeshes(). Both imprint the embedded grid into the background and preserve the total area.

g1 = mf.build_cmesh(np.linspace(0.0, 1.0, 7), np.linspace(0.0, 1.0, 7))
g2 = mf.build_cmesh(np.linspace(0.25, 0.75, 5), np.linspace(0.25, 0.75, 5))

imprint = g1.overlay(g2, mf.OverlayOperation.IMPRINT)

g1m, g2m = g1.to_mc(), g2.to_mc()
g1m.setMeshDimension(2)
g2m.setMeshDimension(2)
mc_imprint = mc.MEDCouplingUMesh.Intersect2DMeshes(g1m, g2m, 1e-12)[0]

a_mf = area_2d(imprint)
a_mc = float(np.asarray(mc_imprint.getMeasureField(True).getArray().getValues()).sum())
print("imprint area (mefikit)    :", a_mf)
print("imprint area (medcoupling):", a_mc)
assert abs(a_mf - 1.0) < 1e-9 and abs(a_mc - 1.0) < 1e-9
print("identical (unit) area: OK")

imprint area (mefikit)    : 1.0000000000000002
imprint area (medcoupling): 1.0
identical (unit) area: OK

pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("overlay(IMPRINT)")
pt.add_mesh(g1.to_pyvista(), show_edges=True, opacity=0.25, color="grey")
pt.add_mesh(imprint.to_pyvista().shrink(0.8), show_edges=True, line_width=2)
pt.camera_position = "xy"
pt.subplot(0, 1)
pt.add_text("Intersect2DMeshes()")
pt.add_mesh(mc_to_pyvista(g1m), show_edges=True, opacity=0.25, color="grey")
pt.add_mesh(mc_to_pyvista(mc_imprint).shrink(0.8), show_edges=True, line_width=2)
pt.camera_position = "xy"
pt.show()

Conservative P0/P0 field transfer (QUAD4)

The flagship operation: transfer a cell field from a source mesh to a target mesh, conserving mass. Same prepare / apply split on both sides:

• mefikit: build the operator once (mf.transfer.ConservativeP0), then call apply_update.
• medcoupling: MEDCouplingRemapper, prepare("P0P0"), then transferField.

The transferred fields match to machine precision.

# --- mefikit ----------------------------------------------------------------
src2 = mf.build_cmesh(np.linspace(0.0, 1.0, 20), np.linspace(0.0, 1.0, 20))
tgt2 = mf.build_cmesh(np.linspace(0.0, 1.0, 24), np.linspace(0.0, 1.0, 24))
src2.fields["T"] = 1.0 + (mf.X - 0.5) ** 2 + 0.5 * mf.Y

op = mf.transfer.ConservativeP0(src2, tgt2)  # prepare once
op.apply_update(src2, "T", tgt2, "T", def_val=0.0)
mf_tgt_vals = tgt2.fields["T"].numpy()

# --- medcoupling -------------------------------------------------------------
sc = src2.to_mc()
sc.setMeshDimension(2)
tc = tgt2.to_mc()
tc.setMeshDimension(2)

f_src = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
_arr = mc.DataArrayDouble(
    np.ascontiguousarray(np.asarray(src2.fields["T"].numpy()).ravel())
)
_arr.setName("T")
f_src.setArray(_arr)
f_src.setMesh(sc)
f_src.setNature(mc.IntensiveConservation)

remap = mc.MEDCouplingRemapper()
remap.prepare(sc, tc, "P0P0")  # prepare once
f_tgt = remap.transferField(f_src, 0.0)
mc_tgt_vals = f_tgt.getArray().toNumPyArray()

print(
    "max |mefikit - medcoupling| after P0P0:",
    float(np.abs(mf_tgt_vals - mc_tgt_vals).max()),
)
assert np.allclose(mf_tgt_vals, mc_tgt_vals, atol=1e-9)
print("transferred fields match: OK")

max |mefikit - medcoupling| after P0P0: 3.552713678800501e-15
transferred fields match: OK

vmin, vmax = 0.5, 1.75
pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("Source")
s_pv = src2.to_pyvista()
s_pv["T"] = np.asarray(src2.fields["T"].numpy()).ravel()
pt.add_mesh(s_pv, scalars="T", clim=[vmin, vmax], show_edges=True)
pt.camera_position = "xy"
pt.subplot(0, 1)
pt.add_text("Target (transferred)")
t_pv = tgt2.to_pyvista()
t_pv["T"] = mf_tgt_vals
pt.add_mesh(t_pv, scalars="T", clim=[vmin, vmax], show_edges=True)
pt.camera_position = "xy"
pt.show()

Polyhedra, where mefikit shines the most

Real-world meshes are very often polyhedral: Voronoi meshes, cell-centred finite volumes, unrolled CAD… Both libraries store them; the difference shows when it is time to use them.

We take two different 2000-cell polyhedral meshes shipped with mefikit’s test suite — mesh_36.med and mesh_27.med — and feed the same .med files to both libraries.

# locate the test meshes whatever the current working directory is
root = Path(os.getcwd())
while root != root.parent and not (root / "tests" / "data" / "mesh_36.med").exists():
    root = root.parent
data = root / "tests" / "data"

mf_src = mf.UMesh.read(str(data / "mesh_36.med"))  # mefikit : 1 line
mf_tgt = mf.UMesh.read(str(data / "mesh_27.med"))
mc_src = mc.ReadMeshFromFile(str(data / "mesh_36.med"), 0)  # medcoupling
mc_tgt = mc.ReadMeshFromFile(str(data / "mesh_27.med"), 0)

print("mefikit  mesh_36 :", mf_src.block_types(), mf_src.num_elements(), "cells")
print("mefikit  mesh_27 :", mf_tgt.block_types(), mf_tgt.num_elements(), "cells")
print(
    "medcoupl mesh_36 :",
    mc_src.getNumberOfCells(),
    "cells,",
    mc_src.getNumberOfNodes(),
    "nodes",
)
print("medcoupl mesh_27 :", mc_tgt.getNumberOfCells(), "cells")

# measure + a field + a selection, poly edition, all in mefikit
mf_src.fields["Measure"] = mf.M
mf_src.fields["T"] = 1.0 + 2.0 * mf.M
print("volume(mesh_36):", mf_src.fields["Measure"].sum())
big = mf_src.select(mf.M > mf_src.fields["Measure"].mean())
print("cells above mean volume:", len(big))

mefikit  mesh_36 : ['PHED'] 2000 cells
mefikit  mesh_27 : ['PHED'] 2000 cells
medcoupl mesh_36 : 2000 cells, 12404 nodes
medcoupl mesh_27 : 2000 cells


volume(mesh_36): 0.99999999869571
cells above mean volume: 898

# boundaries and descending connectivity on polyhedra
mf_bnd = mf_src.boundaries()
mf_desc = mf_src.descend()

mc_bnd = mc_src.buildBoundaryMesh(True)
mc_d = mc_src.buildDescendingConnectivity()
mc_desc = mc_d[0] if isinstance(mc_d, tuple) else mc_d

print(
    "boundary faces (mefikit)     :",
    mf_bnd.num_elements(),
    "| medcoupling:",
    mc_bnd.getNumberOfCells(),
)
print(
    "descending    (mefikit)     :",
    mf_desc.num_elements(),
    "| medcoupling:",
    mc_desc.getNumberOfCells(),
)
assert mf_bnd.num_elements() == mc_bnd.getNumberOfCells()
assert mf_desc.num_elements() == mc_desc.getNumberOfCells()
print("identical counts: OK")

boundary faces (mefikit)     : 888 | medcoupling: 888
descending    (mefikit)     : 14401 | medcoupling: 14401
identical counts: OK

pt = pv.Plotter(shape=(1, 2))
pt.subplot(0, 0)
pt.add_text("mefikit - boundaries() of mesh_36 (PHED)")
pt.add_mesh(mf_bnd.to_pyvista(), color="grey", show_edges=True)
pt.camera_position = "xy"
pt.subplot(0, 1)
pt.add_text("medcoupling - buildBoundaryMesh() (PHED)")
pt.add_mesh(mc_to_pyvista(mc_bnd), color="grey", show_edges=True)
pt.camera_position = "xy"
pt.show()

Polyhedral-to-polyhedral remap

Same P0/P0 transfer, but between the two real polyhedral meshes. Prepare once, apply many times is the name of the game for unsteady runs, so both sides are timed separately. We scan a few mesh sizes (always the same geometry on both sides) and check that the actual transferred fields coincide.

def poly_subset(path, n, mc_side):
    if mc_side:
        m = mc.ReadMeshFromFile(str(path), 0)
        m = m.buildPartOfMySelf(np.arange(n, dtype=np.int64).tolist())
        m.mergeNodes(1e-12)
        return m
    m = mf.UMesh.read(str(path))
    return m.select(mf.sel.ids({"PHED": np.arange(n)})).to_mesh()


def mc_p0_field(mesh, vals):
    f = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
    a = mc.DataArrayDouble(np.ascontiguousarray(vals, np.float64))
    a.setName("T")
    f.setArray(a)
    f.setMesh(mesh)
    f.setNature(mc.IntensiveConservation)
    return f


def median_ms(fn, n=3):
    times = []
    for _ in range(n):
        t0 = time.perf_counter()
        fn()
        times.append(time.perf_counter() - t0)
    return float(np.median(times) * 1e3)


poly_res = []
for n in (100, 200, 400):
    p_src = poly_subset(data / "mesh_36.med", n, False)
    p_tgt = poly_subset(data / "mesh_27.med", n, False)
    p_src.fields["T"] = 1.0 + 2.0 * mf.M

    mf_prepare = median_ms(
        lambda p_src=p_src, p_tgt=p_tgt: mf.transfer.ConservativeP0(p_src, p_tgt)
    )
    op = mf.transfer.ConservativeP0(p_src, p_tgt)
    mf_apply = median_ms(
        lambda op=op, p_src=p_src, p_tgt=p_tgt: op.apply_update(
            p_src, "T", p_tgt, "T", def_val=0.0
        )
    )

    c_src = poly_subset(data / "mesh_36.med", n, True)
    c_tgt = poly_subset(data / "mesh_27.med", n, True)
    f_src = mc_p0_field(c_src, np.asarray(p_src.fields["T"].numpy()).ravel())

    remap = mc.MEDCouplingRemapper()
    mc_prepare = median_ms(
        lambda remap=remap, c_src=c_src, c_tgt=c_tgt: remap.prepare(
            c_src, c_tgt, "P0P0"
        )
    )
    mc_apply = median_ms(
        lambda remap=remap, f_src=f_src: remap.transferField(f_src, 0.0)
    )

    op.apply_update(p_src, "T", p_tgt, "T", def_val=0.0)
    out_mf = np.asarray(p_tgt.fields["T"].numpy()).ravel()
    remap.prepare(c_src, c_tgt, "P0P0")
    out_mc = remap.transferField(f_src, 0.0).getArray().toNumPyArray()

    diff = float(np.abs(out_mf - out_mc).max())
    poly_res.append((n, mf_prepare, mf_apply, mc_prepare, mc_apply, diff))
    print(
        f"poly n={n:4d} | mefikit {mf_prepare:7.2f} ms (prepare) / {mf_apply:5.2f} ms (apply)"
        f" | medcoupling {mc_prepare:8.1f} ms / {mc_apply:5.2f} ms | {mc_prepare / mf_prepare:4.0f}x faster"
    )
    print(
        f"            | max |mefikit - medcoupling| on the transferred field: {diff:.2e}"
    )

assert all(r[5] < 1e-9 for r in poly_res)
print("transferred fields match at every size: OK")

poly n= 100 | mefikit    2.59 ms (prepare) /  0.01 ms (apply) | medcoupling    153.5 ms /  0.20 ms |   59x faster
            | max |mefikit - medcoupling| on the transferred field: 5.88e-15


poly n= 200 | mefikit    7.14 ms (prepare) /  0.01 ms (apply) | medcoupling    568.9 ms /  0.38 ms |   80x faster
            | max |mefikit - medcoupling| on the transferred field: 6.38e-15


poly n= 400 | mefikit   23.26 ms (prepare) /  0.02 ms (apply) | medcoupling   2347.4 ms /  0.74 ms |  101x faster
            | max |mefikit - medcoupling| on the transferred field: 7.22e-15
transferred fields match at every size: OK

ns = [r[0] for r in poly_res]
fig, ax = plt.subplots(figsize=(8, 5))
ax.loglog(ns, [r[1] for r in poly_res], "o-", label="mefikit prepare")
ax.loglog(ns, [r[3] for r in poly_res], "s-", label="medcoupling prepare")
ax.set_xlabel("polyhedral cells")
ax.set_ylabel("prepare time (ms)")
ax.set_title("P0/P0 remap prepare on real polyhedral meshes (mesh_36 -> mesh_27)")
ax.grid(True, which="both", ls="--", alpha=0.4)
ax.legend()
fig.tight_layout()
plt.show()


r = poly_res[-1]
print(
    f"At {r[0]} cells mefikit prepares the polyhedral remap {r[3] / r[1]:.0f} times faster "
    "than medcoupling, and the gap grows with the mesh size."
)

At 400 cells mefikit prepares the polyhedral remap 101 times faster than medcoupling, and the gap grows with the mesh size.

Performance on common operations

Now the same benchmark spirit on the structured meshes (quad / hexa) of daily life. Timings are medians of several runs on this machine; tiny absolute values should be read with perspective. Both libraries always work on the exact same geometry: this is what fairness looks like.

Workloads are the same ones used in mefikit’s tests/bench_vs_medcoupling.py.

N_ITER = 10
MC_INTENSIVE = 37
MC_EXTENSIVE = 35

N2D = 96  # 96x96  = 9216 QUAD4 cells
N3D = 16  # 16^3   = 4096 HEX8 cells
NPOLY = 16  # poly remap target, 16^3 source
MERGE_N = 24  # 2 stacked 24x24 HEX8 layers, duplicated interface
DESCEND_N = 24  # 24^3 hexa grid -> faces
OVERLAY_N = 32  # 32x32 grid overlayed by an embedded 8x8 block
CRACK_N = 20  # 20^3 hexa grid, cracked along all its faces

RTOL = 1e-9
ATOL = 1e-9


def median_time(fn, n=N_ITER):
    # median wall time of fn over n runs, in milliseconds
    times = []
    for _ in range(n):
        t0 = time.perf_counter()
        fn()
        times.append(time.perf_counter() - t0)
    return float(np.median(times) * 1e3)


def mc_mesh(mesh, dim):
    m = mesh.to_mc()
    m.setMeshDimension(dim)
    return m


def mc_field(mmesh, vals, nature, name="T"):
    f = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
    a = mc.DataArrayDouble(np.ascontiguousarray(vals.ravel()))
    a.setName(name)
    f.setArray(a)
    f.setMesh(mmesh)
    f.setNature(nature)
    return f


def field_2d(nx):
    i, j = np.meshgrid(np.arange(nx), np.arange(nx), indexing="ij")
    xc, yc = (i + 0.5) / nx, (j + 0.5) / nx
    return (1.0 + 0.5 * np.sin(2 * np.pi * xc) * np.cos(np.pi * yc)).reshape(-1, 1)


def field_3d(nx):
    i, j, k = np.meshgrid(np.arange(nx), np.arange(nx), np.arange(nx), indexing="ij")
    xc, yc, zc = (i + 0.5) / nx, (j + 0.5) / nx, (k + 0.5) / nx
    return (
        1.0 + 0.5 * np.sin(2 * np.pi * xc) * np.cos(np.pi * yc) * np.cos(np.pi * zc)
    ).reshape(-1, 1)


def dump_merged_mesh(nx):
    # two stacked HEX8 layers; the shared interface is duplicated (2 node sets)
    gx = np.linspace(0.0, 1.0, nx + 1)
    px, py = np.meshgrid(gx, gx, indexing="ij")
    z0 = np.c_[px.ravel(), py.ravel(), np.zeros((nx + 1) ** 2)]
    z1 = np.c_[px.ravel(), py.ravel(), np.ones((nx + 1) ** 2)]
    coords = np.ascontiguousarray(np.vstack([z0, z1, z1, z0 + 2.0]), np.float64)

    def nid(i, j, layer):
        return layer * (nx + 1) ** 2 + i * (nx + 1) + j

    conn = []
    for i in range(nx):
        for j in range(nx):
            conn += [
                [
                    nid(i, j, 0),
                    nid(i + 1, j, 0),
                    nid(i + 1, j + 1, 0),
                    nid(i, j + 1, 0),
                    nid(i, j, 1),
                    nid(i + 1, j, 1),
                    nid(i + 1, j + 1, 1),
                    nid(i, j + 1, 1),
                ],
                [
                    nid(i, j, 2),
                    nid(i + 1, j, 2),
                    nid(i + 1, j + 1, 2),
                    nid(i, j + 1, 2),
                    nid(i, j, 3),
                    nid(i + 1, j, 3),
                    nid(i + 1, j + 1, 3),
                    nid(i, j + 1, 3),
                ],
            ]
    mesh = mf.UMesh(coords)
    mesh.add_regular_block("HEX8", np.ascontiguousarray(np.array(conn), np.uintp))
    return mesh

def bench_remap(dim, n, build_iter, poly=False):
    x = np.linspace(0.0, 1.0, n + 1)
    axes = [x] * dim
    shift = [0.5 / n] + [0.0] * (dim - 1)
    et = "QUAD4" if dim == 2 else "HEX8"
    vals = field_2d(n) if dim == 2 else field_3d(n)

    src = mf.build_cmesh(*axes)
    tgt = mf.build_cmesh(*[a + s for a, s in zip(axes, shift)])
    if poly:
        tgt = tgt.polyze()

    sm = mc_mesh(src, dim)
    tm = mc_mesh(mf.build_cmesh(*[a + s for a, s in zip(axes, shift)]), dim)
    if poly:
        tm.convertAllToPoly()

    mf_build = median_time(lambda: mf.ConservativeP0(src, tgt), build_iter)

    vt = mc.MEDCouplingRemapper()
    mc_prepare = median_time(lambda: vt.prepare(sm, tm, "P0P0"), build_iter)

    src.set_field("T", {et: np.ascontiguousarray(vals)})
    op = mf.ConservativeP0(src, tgt)
    mf_apply = median_time(lambda: op.apply_update(src, "T", tgt, "T", def_val=0.0))

    field = mc_field(sm, vals, MC_INTENSIVE)
    mc_transfer = median_time(lambda: vt.transferField(field, 0.0))

    checks = {}

    def read_mf():
        parts = [np.asarray(v).ravel() for v in tgt.fields["T"].values().values()]
        return np.concatenate(parts)

    op.apply_update(src, "T", tgt, "T", def_val=0.0)
    out_mf = read_mf()
    out_mc = np.asarray(vt.transferField(field, 0.0).getArray().getValues())

    if not poly:
        diff = np.max(np.abs(out_mf - out_mc))
        checks["intensive match (mf == mc)"] = (
            np.allclose(out_mf, out_mc, rtol=RTOL, atol=ATOL),
            diff,
        )

    vol = 1.0 / n**dim
    mass_analytic = float(vals.sum() * vol)
    cs = mf.build_cmesh(*axes)
    ct = cs.polyze() if poly else mf.build_cmesh(*axes)
    cs.set_field("T", {et: np.ascontiguousarray(vals)})
    opc = mf.ConservativeP0(cs, ct)
    opc.apply_update(cs, "T", ct, "T", def_val=0.0, extensive=True)
    parts = [np.asarray(v).ravel() for v in ct.fields["T"].values().values()]
    mass_mf = float(np.concatenate(parts).sum())
    csm = mc_mesh(cs, dim)
    ctm = mc_mesh(mf.build_cmesh(*axes), dim)
    if poly:
        ctm.convertAllToPoly()
    mass_mc = float(
        np.asarray(
            vt.transferField(mc_field(csm, vals * vol, MC_EXTENSIVE), 0.0)
            .getArray()
            .getValues()
        ).sum()
    )
    tol = max(1e-9, mass_analytic * 1e-9)
    checks["mass (mf == analytic)"] = (
        abs(mass_mf - mass_analytic) <= tol,
        (round(mass_analytic, 10), round(mass_mf, 10)),
    )
    checks["mass (mc == mf)"] = (
        abs(mass_mc - mass_mf) <= tol,
        (round(mass_mf, 10), round(mass_mc, 10)),
    )

    return {
        "mf_build": mf_build,
        "mc_prepare": mc_prepare,
        "mf_apply": mf_apply,
        "mc_transfer": mc_transfer,
        "checks": checks,
    }


def bench_merge():
    mesh = dump_merged_mesh(MERGE_N)
    mm = mesh.to_mc()
    mf_t = median_time(lambda: mesh.merge_nodes(1e-12))
    mc_t = median_time(lambda: mm.mergeNodes(1e-12))
    used_after = used_nodes(mesh.merge_nodes(1e-12))
    mc_ref = mc_mesh(mesh, 3)
    mc_ref.mergeNodes(1e-12)
    nodes_after = mc_ref.getNumberOfNodes()
    return {
        "mf": mf_t,
        "mc": mc_t,
        "checks": {
            "used nodes == 1875": (used_after == 1875, used_after),
            "mc nodes == mf": (nodes_after == used_after, (used_after, nodes_after)),
        },
    }


def bench_descend():
    n = DESCEND_N
    axes = [np.linspace(0.0, 1.0, n + 1)] * 3
    mesh = mf.build_cmesh(*axes)
    mm = mc_mesh(mesh, 3)
    mf_t = median_time(mesh.descend)
    mc_t = median_time(mm.buildDescendingConnectivity)
    f_mf = int(mesh.descend().blocks()["QUAD4"].shape[0])
    f_mc = int(mm.buildDescendingConnectivity()[0].getNumberOfCells())
    expected = 3 * n * n * (n + 1)
    return {
        "mf": mf_t,
        "mc": mc_t,
        "checks": {
            "faces == 3 n^2 (n+1)": (
                (f_mf == expected) and (f_mc == expected),
                (expected, f_mf, f_mc),
            )
        },
    }


def bench_overlay():
    n = OVERLAY_N
    m1 = mf.build_cmesh(np.linspace(0.0, 1.0, n + 1), np.linspace(0.0, 1.0, n + 1))
    m2 = mf.build_cmesh(np.linspace(0.2, 0.7, 9), np.linspace(0.2, 0.7, 9))
    m1m = mc_mesh(m1, 2)
    m2m = mc_mesh(m2, 2)
    mf_t = median_time(lambda: m1.overlay(m2), 5)
    mc_t = median_time(
        lambda: mc.MEDCouplingUMesh.Intersect2DMeshes(m1m, m2m, 1e-12), 5
    )
    a_mf = area_2d(m1.overlay(m2))
    a_mc = float(
        np.asarray(
            mc.MEDCouplingUMesh.Intersect2DMeshes(m1m, m2m, 1e-12)[0]
            .getMeasureField(True)
            .getArray()
            .getValues()
        ).sum()
    )
    return {
        "mf": mf_t,
        "mc": mc_t,
        "checks": {
            "area == 1 (both)": (
                abs(a_mf - 1.0) < ATOL and abs(a_mc - 1.0) < ATOL,
                (a_mf, a_mc),
            )
        },
    }


def bench_crack():
    n = CRACK_N
    axes = [np.linspace(0.0, 1.0, n + 1)] * 3
    mesh = mf.build_cmesh(*axes)
    faces = mesh.descend()

    # medcoupling cracks an MEDFileUMesh along a group of M1 faces
    vm = mc_mesh(mesh, 3)
    fm = vm.buildDescendingConnectivity()[0]
    fm.setName(vm.getName())
    grp = mc.DataArrayInt(np.arange(fm.getNumberOfCells(), dtype=np.int64))
    grp.setName("crack-line")
    fmu = mc.MEDFileUMesh.New()
    fmu.setMeshAtLevel(0, vm)
    fmu.setMeshAtLevel(-1, fm)
    fmu.setGroupsAtLevel(-1, [grp])

    def run_mc():
        box = fmu.deepCopy()
        box.crackAlong("crack-line")
        return box

    mf_t = median_time(lambda: mesh.crack(faces))
    mc_t = median_time(run_mc, n=3)

    n_mf = used_nodes(mesh.crack(faces))
    n_mc = run_mc().getNumberOfNodes()
    return {
        "mf": mf_t,
        "mc": mc_t,
        "checks": {"node count (mf == mc)": (n_mf == n_mc, (n_mf, n_mc))},
    }

bench = {}
bench["remap-2d"] = bench_remap(2, N2D, build_iter=10, poly=False)
print("remap-2d done")
bench["remap-3d"] = bench_remap(3, N3D, build_iter=10, poly=False)
print("remap-3d done")
bench["remap-3d-poly"] = bench_remap(3, NPOLY, build_iter=5, poly=True)
print("remap-3d-poly done")
bench["merge-nodes"] = bench_merge()
print("merge-nodes done")
bench["descend"] = bench_descend()
print("descend done")
bench["overlay"] = bench_overlay()
print("overlay done")
bench["crack"] = bench_crack()
print("crack done")

remap-2d done


remap-3d done


remap-3d-poly done
merge-nodes done


descend done


overlay done


crack done

rows = [
    (
        "remap-2d",
        "build/prepare",
        bench["remap-2d"]["mf_build"],
        bench["remap-2d"]["mc_prepare"],
    ),
    (
        "remap-2d",
        "transfer",
        bench["remap-2d"]["mf_apply"],
        bench["remap-2d"]["mc_transfer"],
    ),
    (
        "remap-3d",
        "build/prepare",
        bench["remap-3d"]["mf_build"],
        bench["remap-3d"]["mc_prepare"],
    ),
    (
        "remap-3d",
        "transfer",
        bench["remap-3d"]["mf_apply"],
        bench["remap-3d"]["mc_transfer"],
    ),
    (
        "remap-3d-poly",
        "build/prepare",
        bench["remap-3d-poly"]["mf_build"],
        bench["remap-3d-poly"]["mc_prepare"],
    ),
    (
        "remap-3d-poly",
        "transfer",
        bench["remap-3d-poly"]["mf_apply"],
        bench["remap-3d-poly"]["mc_transfer"],
    ),
    ("merge-nodes", "merge", bench["merge-nodes"]["mf"], bench["merge-nodes"]["mc"]),
    ("descend", "run", bench["descend"]["mf"], bench["descend"]["mc"]),
    ("overlay", "run", bench["overlay"]["mf"], bench["overlay"]["mc"]),
    ("crack", "crack", bench["crack"]["mf"], bench["crack"]["mc"]),
]

print(
    f"{'case':<14s} {'step':<14s} {'mefikit ms':>12s} {'medcoup ms':>12s} {'mc/mf':>9s}"
)
print("-" * 62)
for case, step, mf_t, mc_t in rows:
    print(f"{case:<14s} {step:<14s} {mf_t:>12.3f} {mc_t:>12.3f} {mc_t / mf_t:>8.1f}x")

print()
print("correctness cross-checks:")
all_ok = True
for tag in bench:
    for label, (ok, info) in bench[tag]["checks"].items():
        all_ok &= ok
        print(f"  [{'OK' if ok else 'FAIL'}] {tag}: {label}  {info}")
assert all_ok, "a cross-check failed"

case           step             mefikit ms   medcoup ms     mc/mf
--------------------------------------------------------------
remap-2d       build/prepare        27.045       24.817      0.9x
remap-2d       transfer              0.181        3.905     21.5x
remap-3d       build/prepare       162.564      692.317      4.3x
remap-3d       transfer              0.094        1.235     13.2x
remap-3d-poly  build/prepare       356.550     8746.085     24.5x
remap-3d-poly  transfer              0.099        2.171     22.0x
merge-nodes    merge                 0.538        4.019      7.5x
descend        run                  56.470      103.696      1.8x
overlay        run                   2.071       68.392     33.0x
crack          crack               235.632     1624.416      6.9x

correctness cross-checks:
  [OK] remap-2d: intensive match (mf == mc)  7.752687380957468e-13
  [OK] remap-2d: mass (mf == analytic)  (1.0, 1.0)
  [OK] remap-2d: mass (mc == mf)  (1.0, 1.0)
  [OK] remap-3d: intensive match (mf == mc)  0.0
  [OK] remap-3d: mass (mf == analytic)  (1.0, 1.0)
  [OK] remap-3d: mass (mc == mf)  (1.0, 1.0)
  [OK] remap-3d-poly: mass (mf == analytic)  (1.0, 1.0)
  [OK] remap-3d-poly: mass (mc == mf)  (1.0, 1.0)
  [OK] merge-nodes: used nodes == 1875  1875
  [OK] merge-nodes: mc nodes == mf  (1875, 1875)
  [OK] descend: faces == 3 n^2 (n+1)  (43200, 43200, 43200)
  [OK] overlay: area == 1 (both)  (1.0, 1.0)
  [OK] crack: node count (mf == mc)  (64000, 64000)

def twin_bars(labels, mf_vals, mc_vals, ylabel, title, rot=0):
    x = np.arange(len(labels))
    w = 0.36
    fig, ax = plt.subplots(figsize=(10, 5))
    b1 = ax.bar(x - w / 2, mf_vals, w, label="mefikit")
    b2 = ax.bar(x + w / 2, mc_vals, w, label="medcoupling")
    ax.set_xticks(x)
    ax.set_xticklabels(labels, rotation=rot, ha="right")
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    ax.legend()
    ax.bar_label(b1, fmt="%.1f", padding=1, fontsize=8)
    ax.bar_label(b2, fmt="%.1f", padding=1, fontsize=8)
    fig.tight_layout()
    plt.show()


# --- prepare / build: every operation has exactly one ---
codes = ["remap-2d", "remap-3d", "remap-3d-poly"]
singles = ["merge-nodes", "descend", "overlay", "crack"]
labels = codes + singles
twin_bars(
    labels,
    [bench[c]["mf_build"] for c in codes] + [bench[c]["mf"] for c in singles],
    [bench[c]["mc_prepare"] for c in codes] + [bench[c]["mc"] for c in singles],
    "time (ms)",
    "Prepare / build time — all operations",
    rot=25,
)


# --- transfer / apply: only the P0/P0 remaps have a separate apply step ---
twin_bars(
    codes,
    [bench[c]["mf_apply"] for c in codes],
    [bench[c]["mc_transfer"] for c in codes],
    "time (ms)",
    "Transfer / apply time — P0/P0 remaps",
)

labels = [f"{c}\n{s}" for c, s, _, _ in rows]
ratios = [mc_t / mf_t for _, _, mf_t, mc_t in rows]

fig, ax = plt.subplots(figsize=(10, 6))
colors = ["tab:blue" if r >= 1 else "tab:red" for r in ratios]
yloc = np.arange(len(ratios))[::-1]
ax.barh(yloc, ratios, color=colors)
ax.axvline(1.0, color="white", ls="--", lw=1)
ax.set_yticks(yloc)
ax.set_yticklabels(labels)
ax.set_xscale("log")
ax.set_xlabel("medcoupling time / mefikit time  (>1 means mefikit is faster)")
ax.set_title("How many times longer medcoupling takes than mefikit, per operation")
for y, r in zip(yloc, ratios):
    ax.text(r * 1.25, y, f" {r:.1f}x", va="center", ha="left", fontsize=8)
fig.tight_layout()
plt.show()

Feature comparison at a glance

Both libraries cover the common operations shown above. Beyond that, mefikit adds ergonomics of its own — most of it already exercised in this notebook.

Let us be perfectly honest: medcoupling remains far broader than mefikit — decades of API surface, advanced field machinery, spline remapping and a huge ecosystem around the MED format. mefikit does not try to replace that. It tries to be direct for the operations above, and fast where it counts: two propositions you have just measured.

Closing thoughts

A fair question: should a medcoupling user switch to mefikit? As always, it depends.

• You already have a solid, mature medcoupling pipeline: nothing here is broken, and mefikit reads the very same .med files — you can adopt it as a complement. mesh.to_mc() hands you the medcoupling twin of a mefikit mesh in one line.
• You are starting a new project, work a lot with polyhedral meshes, value a short readable API and care about remap performance: give mefikit a try, you should feel at home very quickly.

Both libraries agree on the most important thing: the meshes are the same, the physics is the same, and the results match to machine precision. Thank you medcoupling for setting such a high bar — mefikit simply tries to reach it with fewer keystrokes and a bit more speed.

If you want to go further: browse the other notebooks of this book, have a look at the roadmap and try the operations above on your own meshes.