Writes mdpa files, cf.
    <https://github.com/KratosMultiphysics/Kratos/wiki/Input-data>.

def write(filename, mesh, write_binary=False):
    """Writes mdpa files, cf.
    <https://github.com/KratosMultiphysics/Kratos/wiki/Input-data>.
    """
    assert not write_binary
    if mesh.points.shape[1] == 2:
        logging.warning(
            "mdpa requires 3D points, but 2D points given. "
            "Appending 0 third component."
        )
        mesh.points = numpy.column_stack(
            [mesh.points[:, 0], mesh.points[:, 1], numpy.zeros(mesh.points.shape[0])]
        )

    # Kratos cells are mostly ordered like VTK, with a few exceptions:
    cells = mesh.cells.copy()
    if "hexahedron20" in cells:
        cells["hexahedron20"] = cells["hexahedron20"][
            :, [0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 10, 9, 16, 17, 18, 19, 12, 15, 14, 13]
        ]
    if "hexahedron27" in cells:
        cells["hexahedron27"] = cells["hexahedron27"][
            :,
            [
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                11,
                10,
                9,
                16,
                17,
                18,
                19,
                12,
                15,
                14,
                13,
                22,
                24,
                21,
                23,
                20,
                25,
                26,
            ],
        ]

    with open(filename, "wb") as fh:
        # Write some additional info
        fh.write(("Begin ModelPartData\n").encode("utf-8"))
        fh.write(("//  VARIABLE_NAME value\n").encode("utf-8"))
        fh.write(("End ModelPartData\n\n").encode("utf-8"))
        fh.write(("Begin Properties 0\n").encode("utf-8"))
        fh.write(("End Properties\n\n").encode("utf-8"))

        # Split the cell data: gmsh:physical and gmsh:geometrical are tags, the
        # rest is actual cell data.
        tag_data = {}
        other_data = {}
        for cell_type, a in mesh.cell_data.items():
            tag_data[cell_type] = {}
            other_data[cell_type] = {}
            for key, data in a.items():
                if key in ["gmsh:physical", "gmsh:geometrical"]:
                    tag_data[cell_type][key] = data.astype(numpy.int32)
                else:
                    other_data[cell_type][key] = data

        # We identity which dimension are we
        dimension = 2
        for cell_type in cells.keys():
            name_elem = _meshio_to_mdpa_type[cell_type]
            if local_dimension_types[name_elem] == 3:
                dimension = 3
                break

        # We identify the entities
        _write_nodes(fh, mesh.points, write_binary)
        _write_elements_and_conditions(fh, cells, tag_data, write_binary, dimension)
        for name, dat in mesh.point_data.items():
            _write_data(fh, "NodalData", name, dat, write_binary)
        cell_data_raw = raw_from_cell_data(other_data)
        for (
            name,
            dat,
        ) in (
            cell_data_raw.items()
        ):  # NOTE: We will assume always when writing that the components are elements (for now)
            _write_data(fh, "ElementalData", name, dat, write_binary)

    return