There is currently no vtk output for 1D elements implemented. I asked chatgpt and it gave some surprisingly good code (with 2 small fixes from me this runs):
"""
Export boundary data from a 2D NGSolve mesh to a legacy VTK PolyData (.vtk) file.
- Writes boundary edges as 1D line elements (POLYDATA LINES)
- Attaches a scalar defined on the boundary as POINT_DATA (values at boundary vertices)
- Works from either a volume GridFunction (its trace is projected) or a CoefficientFunction
- Uses first-order H1 on the boundary to get vertex values (piecewise linear along edges)
Limitations / notes:
- Exports scalar data; extend `collect_point_values` to handle vector fields if needed.
- For high-order fields, this linearizes to vertex values. Increase resolution by post-sampling
more points per edge and writing them as polylines (left as an optional extension below).
Usage example:
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
# Example scalar field on the 2D domain
fes = H1(mesh, order=2, dirichlet=[]) # any order
u = GridFunction(fes)
x, y = x, y # ngsolve's coordinate CFs
u.Set( sin(pi*x)*cos(pi*y) )
# Export trace of u on the boundary to VTK
export_boundary_to_vtk(mesh, u, filename="boundary_trace.vtk", name="u_boundary")
Open the resulting `boundary_trace.vtk` in ParaView or VisIt.
"""
from typing import Dict, Iterable, Optional, Tuple
from ngsolve import (
Mesh,
GridFunction,
CoefficientFunction,
H1,
BND,
)
def _collect_boundary_topology(mesh: Mesh, definedon: Optional[Iterable[int]] = None):
"""Return (vertex_coords, edges) for boundary elements.
vertex_coords: dict[vnr] = (x, y, z)
edges: list[(v0, v1)] using the mesh's *global* vertex numbers
If `definedon` is provided, restrict to boundary regions whose material numbers are in it.
"""
vertex_coords: Dict[int, Tuple[float, float, float]] = {}
edges = []
for el in mesh.Elements(BND):
# Optionally filter by boundary region(s)
if definedon is not None and el.mat not in definedon:
continue
vs = el.vertices
assert len(vs) == 2, "Expected 1D boundary elements with 2 vertices"
v0, v1 = vs[0].nr, vs[1].nr
# store coords
for v in (v0, v1):
p = mesh.ngmesh[v+1]
vertex_coords[v] = (float(p.p[0]), float(p.p[1]), 0.0)
edges.append((v0, v1))
return vertex_coords, edges
def _project_to_boundary_vertex_values(
mesh: Mesh,
field: "GridFunction|CoefficientFunction",
definedon: Optional[Iterable[int]] = None,
) -> Dict[int, float]:
"""Project a scalar field's trace to first-order H1 on the boundary and
return a mapping vertex_id -> value at that vertex.
This assumes that for order=1 H1 on a 1D boundary, the DOFs coincide with vertices.
"""
# Build an H1 space *on the boundary only*, order 1 so DOFs == vertices
fes_bnd = H1(mesh, order=1, definedon=mesh.Boundaries(".*"))
# Interpolate / project the trace onto the boundary space
gfb = GridFunction(fes_bnd)
gfb.Set(field, definedon=mesh.Boundaries(".*")) # NGSolve takes the boundary trace automatically for H1
# Collect values per vertex from boundary elements' local dofs
v2val: Dict[int, float] = {}
for el in mesh.Elements(BND):
if definedon is not None and el.mat not in definedon:
continue
dofs = fes_bnd.GetDofNrs(el)
vs = el.vertices
# Sanity: one dof per vertex at order=1
if len(dofs) != 2 or len(vs) != 2:
raise RuntimeError("Unexpected number of boundary dofs/vertices; ensure order=1 H1.")
for ldof, v in zip(dofs, vs):
v2val[v.nr] = float(gfb.vec[ldof])
return v2val
def export_boundary_to_vtk(
mesh: Mesh,
field: "GridFunction|CoefficientFunction",
filename: str,
name: str = "u",
boundary_regions: Optional[Iterable[int]] = None,
):
"""Export the trace of a scalar field on the 2D mesh boundary to a legacy VTK .vtk file.
Parameters
----------
mesh : ngsolve.Mesh
2D mesh to export the boundary from.
field : GridFunction | CoefficientFunction
Scalar field defined in the volume; its boundary trace is exported.
filename : str
Output path ending with .vtk
name : str
Name for the scalar in VTK (POINT_DATA array)
boundary_regions : optional iterable of material numbers
If given, restrict export to these boundary region material numbers.
"""
# 1) gather boundary topology and coordinates
vtx_coords, edges = _collect_boundary_topology(mesh, boundary_regions)
if not edges:
raise ValueError("Mesh has no boundary edges to export (or none matched boundary_regions).")
# 2) compute vertex values via H1 order=1 projection on boundary
vtx_values = _project_to_boundary_vertex_values(mesh, field, boundary_regions)
print(vtx_values)
# VTK needs a compact 0..N-1 indexing for points; build a reindex map
vtx_list = sorted(vtx_coords.keys())
vtx_to_compact = {v: i for i, v in enumerate(vtx_list)}
# Ensure we have a value for every exported vertex
missing = [v for v in vtx_list if v not in vtx_values]
if missing:
raise RuntimeError(
f"Missing boundary values for {len(missing)} vertex/vertices.\n"
"Ensure the field is scalar and the trace was set correctly."
)
# 3) write legacy VTK POLYDATA
with open(filename, "w", encoding="utf-8") as f:
f.write("# vtk DataFile Version 3.0\n")
f.write(f"NGSolve 2D boundary export: {name}\n")
f.write("ASCII\n")
f.write("DATASET POLYDATA\n")
# POINTS
f.write(f"POINTS {len(vtx_list)} float\n")
for v in vtx_list:
x, y, z = vtx_coords[v]
f.write(f"{x} {y} {z}\n")
# LINES: each edge is a 2-vertex polyline. VTK LINES section needs total connectivity size = M + sum(len(line))
m = len(edges)
connectivity_size = m + 2 * m
f.write(f"LINES {m} {connectivity_size}\n")
for (v0, v1) in edges:
f.write(f"2 {vtx_to_compact[v0]} {vtx_to_compact[v1]}\n")
# POINT_DATA (scalar at vertices)
f.write(f"POINT_DATA {len(vtx_list)}\n")
f.write(f"SCALARS {name} float 1\n")
f.write("LOOKUP_TABLE default\n")
for v in vtx_list:
f.write(f"{vtx_values[v]}\n")
# ---------------- Optional extension: denser sampling per edge ---------------- #
# If you want smoother visualization for high-order fields, you can replace the
# vertex-based export with sub-sampling along each boundary edge. That requires
# evaluating the field at physical points; a robust path is to build a 1D H1
# space of higher order on the boundary, interpolate there, then sample that
# GridFunction at reference points via element mappings. This is left as an
# advanced exercise to keep the main implementation compact and dependable.
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = H1(mesh, order=2)
u = GridFunction(fes)
u.Set(sin(pi*x)*cos(pi*y))
export_boundary_to_vtk(mesh, u, "boundary_trace.vtk", name="u_boundary")
it is not the most efficient implementation and no subdivision, but I think a good starting point 
if you want subdivision I’d rather evaluate with something like MapToAllElements instead of the suggested high order interpolation of GF.
best,
Christopher
