#!/usr/bin/env python3
"""
Comparison: Conforming vs Non-Conforming 1D-in-3D Meshes

This script demonstrates:
1. netgen.occ.Glue() produces NON-CONFORMING meshes (pipe floats inside)
2. gmsh.fragment() produces CONFORMING meshes (pipe nodes shared with tets)

The key difference is topological: with Glue(), pipe vertices are NOT
vertices of any tetrahedron (except at endpoints). With fragment(),
ALL pipe vertices are shared with surrounding tetrahedra.

Conclusion: Use gmsh for geometry fragmentation, then transfer mesh data
to NGSolve. The mesh format conversion is non-trivial but the conformity
is verified.
"""

import numpy as np


def analyze_glue_mesh():
    """Analyze mesh created with netgen.occ.Glue()."""
    from netgen.occ import Box, Pnt, Segment, Glue, OCCGeometry
    from ngsolve import Mesh, VOL, BBND

    print("="*70)
    print("METHOD 1: netgen.occ.Glue() - NON-CONFORMING")
    print("="*70)

    box = Box(Pnt(0, 0, 0), Pnt(1, 1, 1))
    pipe = Segment(Pnt(0, 0.5, 0.5), Pnt(1, 0.5, 0.5))
    pipe.name = "gamma"

    geo = Glue([box, pipe])
    mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.3))

    print(f"\nMesh: nv={mesh.nv}, ne={mesh.ne}")

    # Get VOL vertices
    vol_vertices = set()
    for el in mesh.Elements(VOL):
        for v in el.vertices:
            vol_vertices.add(v.nr)

    # Get pipe (gamma) vertices
    pipe_vertices = set()
    gamma_idx = list(mesh.GetBBoundaries()).index("gamma")
    for el in mesh.Elements(BBND):
        if el.index == gamma_idx:
            for v in el.vertices:
                pipe_vertices.add(v.nr)

    # Check conformity
    shared = pipe_vertices & vol_vertices
    isolated = pipe_vertices - vol_vertices

    print(f"\nConformity analysis:")
    print(f"  Total pipe vertices: {len(pipe_vertices)}")
    print(f"  Shared with VOL elements: {len(shared)}")
    print(f"  ISOLATED (not in any tet): {len(isolated)}")

    if len(isolated) > 0:
        print(f"\n  *** NON-CONFORMING: {len(isolated)} pipe vertices float inside volume ***")
    else:
        print(f"\n  *** CONFORMING: All pipe vertices shared with tets ***")

    return len(isolated)


def analyze_gmsh_fragment():
    """Analyze mesh created with gmsh.fragment()."""
    import gmsh

    print("\n" + "="*70)
    print("METHOD 2: gmsh.model.occ.fragment() - CONFORMING")
    print("="*70)

    gmsh.initialize()
    gmsh.option.setNumber("General.Verbosity", 0)
    gmsh.model.add("pipe_in_box")

    # Create geometry
    box = gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
    p0 = gmsh.model.occ.addPoint(0, 0.5, 0.5)
    p1 = gmsh.model.occ.addPoint(1, 0.5, 0.5)
    pipe = gmsh.model.occ.addLine(p0, p1)

    # FRAGMENT: This is the key operation
    gmsh.model.occ.fragment([(3, box)], [(1, pipe)])
    gmsh.model.occ.synchronize()

    # Find interior edge (pipe)
    edges = gmsh.model.getEntities(dim=1)
    pipe_tags = []
    for dim, tag in edges:
        com = gmsh.model.occ.getCenterOfMass(dim, tag)
        if abs(com[0] - 0.5) < 0.01 and abs(com[1] - 0.5) < 0.01 and abs(com[2] - 0.5) < 0.01:
            pipe_tags.append(tag)

    if pipe_tags:
        gmsh.model.addPhysicalGroup(1, pipe_tags, tag=1, name="gamma")

    gmsh.model.addPhysicalGroup(3, [t for d, t in gmsh.model.getEntities(dim=3)], name="soil")

    # Mesh
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.3)
    gmsh.model.mesh.generate(3)

    # Get mesh data
    node_tags, node_coords, _ = gmsh.model.mesh.getNodes()

    # Get tet nodes
    tet_types, _, tet_nodes = gmsh.model.mesh.getElements(dim=3)
    tet_node_set = set()
    n_tets = 0
    for i, etype in enumerate(tet_types):
        if etype == 4:
            tet_node_set.update(tet_nodes[i].astype(int))
            n_tets = len(tet_nodes[i]) // 4

    # Get pipe nodes
    pipe_node_set = set()
    for tag in pipe_tags:
        seg_types, _, seg_nodes = gmsh.model.mesh.getElements(dim=1, tag=tag)
        for nodes in seg_nodes:
            pipe_node_set.update(nodes.astype(int))

    gmsh.finalize()

    # Check conformity
    shared = pipe_node_set & tet_node_set
    isolated = pipe_node_set - tet_node_set

    print(f"\nMesh: nodes={len(node_tags)}, tets={n_tets}")
    print(f"\nConformity analysis:")
    print(f"  Total pipe nodes: {len(pipe_node_set)}")
    print(f"  Shared with tetrahedra: {len(shared)}")
    print(f"  ISOLATED: {len(isolated)}")

    if len(isolated) > 0:
        print(f"\n  *** NON-CONFORMING ***")
    else:
        print(f"\n  *** CONFORMING: All pipe nodes are tet vertices ***")

    return len(isolated)


def main():
    print("\n" + "#"*70)
    print("# CONFORMING vs NON-CONFORMING 1D-in-3D MESH COMPARISON")
    print("#"*70)

    isolated_glue = analyze_glue_mesh()
    isolated_gmsh = analyze_gmsh_fragment()

    print("\n" + "="*70)
    print("SUMMARY")
    print("="*70)

    print(f"""
Method              | Isolated Nodes | Conforming?
--------------------|----------------|------------
netgen.occ.Glue()   | {isolated_glue:>14} | {'NO' if isolated_glue > 0 else 'YES'}
gmsh.occ.fragment() | {isolated_gmsh:>14} | {'NO' if isolated_gmsh > 0 else 'YES'}

KEY INSIGHT:
- Glue() creates a COMPOUND shape where pipe and box are geometrically
  coincident but topologically SEPARATE. The mesher treats them independently.

- fragment() computes all intersections and creates a conforming topology
  where the pipe is actually embedded in the solid's geometric structure.

IMPLICATION FOR NGSolve:
- With Glue() mesh, Trace().Trace() fails on isolated pipe vertices because
  they have no connection to the 3D field's DOFs.

- With fragment() mesh, pipe nodes ARE tet vertices, so the 3D field's
  values can be correctly traced to the 1D pipe elements.

NEXT STEPS:
1. Use gmsh for geometry construction and meshing
2. Export mesh data (nodes, tets, pipe segments)
3. Build NGSolve mesh manually from this data
   (requires proper CD2Name setup for 1D elements)
""")


if __name__ == "__main__":
    main()