PVPlant/Civil/PVPlantEarthWorks.py

# /**********************************************************************
# *                                                                     *
# * Copyright (c) 2026 Javier Braña <javier.branagutierrez@gmail.com>  *
# *                                                                     *
# * EarthWorks - Cálculo de movimiento de tierras                       *
# *                                                                     *
# * Calcula volúmenes de desmonte (cut) y terraplén (fill) entre una   *
# * superficie diseñada y el terreno natural, usando operaciones        *
# * booleanas OCC + mallas para visualización.                          *
# *                                                                     *
# ***********************************************************************

import FreeCAD
import Part
import Mesh
import math
import numpy as np
import ArchComponent

if FreeCAD.GuiUp:
    import FreeCADGui, os
    from PySide import QtCore, QtGui
    from PySide.QtCore import QT_TRANSLATE_NOOP
else:
    def translate(ctxt, txt): return txt
    def QT_TRANSLATE_NOOP(ctxt, txt): return txt

import PVPlantResources
from PVPlantResources import DirIcons as DirIcons

VOLUME_TYPES = ["Fill", "Cut"]


def makeEarthWorksVolume(vtype=0):
    """Crea un objeto de volumen (Fill=0, Cut=1) en el documento activo."""
    obj = FreeCAD.ActiveDocument.addObject("Part::FeaturePython", VOLUME_TYPES[vtype])
    EarthWorksVolume(obj)
    ViewProviderEarthWorksVolume(obj.ViewObject)
    return obj


# =========================================================================
# Funciones de cálculo principales
# =========================================================================

def get_tracker_rows(frames):
    """
    Agrupa trackers en filas y columnas usando la lógica de Placement.

    Returns:
        (rows, columns): listas de listas agrupadas
    """
    try:
        import PVPlantPlacement
        return PVPlantPlacement.getRows(frames)
    except Exception:
        return None, None


def compute_earthworks(frames, terrain_mesh, slope_tolerance=10.0):
    """
    Calcula el movimiento de tierras para una lista de frames/trackers.

    Args:
        frames: lista de objetos tracker
        terrain_mesh: Mesh del terreno
        slope_tolerance: pendiente máxima E-W (grados)

    Returns:
        (mesh_cut, mesh_fill, volume_cut, volume_fill)
    """
    import MeshPart
    rows, columns = get_tracker_rows(frames)
    if rows is None or not rows:
        return None, None, 0, 0

    # Generar superficie diseñada (loft de bordes de filas)
    design_shape = _build_design_surface(rows, columns, slope_tolerance)
    if design_shape is None:
        return None, None, 0, 0

    # Convertir a mesh para booleanos
    design_mesh = MeshPart.meshFromShape(Shape=design_shape, LinearDeflection=500, AngularDeflection=0.5)

    # Corte: diseño por encima del terreno → material a remover
    cut_mesh = _mesh_above(design_mesh, terrain_mesh)

    # Relleno: diseño por debajo del terreno → material a aportar
    fill_mesh = _mesh_below(design_mesh, terrain_mesh)

    volume_cut = _mesh_volume(cut_mesh)
    volume_fill = _mesh_volume(fill_mesh)

    return cut_mesh, fill_mesh, volume_cut, volume_fill


def _build_design_surface(rows, columns, slope_tolerance):
    """
    Construye la superficie diseñada entre filas de trackers.

    Para cada fila, genera líneas de borde (izquierda/derecha) teniendo
    en cuenta la pendiente transversal (E-W). Luego lofting para crear
    la superficie continua entre filas.
    """
    from DraftGeomUtils import isPlanar

    all_lines = []
    tools = []  # (frame, line_left, line_right)

    for group in rows:
        lines = []
        for i, frame in enumerate(group):
            aw = _angle_to_prev(group, i)
            ae = _angle_to_next(group, i)
            anf = (aw + ae) / 2
            if anf > slope_tolerance:
                anf = slope_tolerance

            wdt = int(frame.Setup.Width / 2) if hasattr(frame.Setup, 'Width') else 0
            zz = wdt * math.sin(math.radians(anf))

            # Línea base a lo largo del tracker
            line = _get_tracker_line(frame)

            # Borde izquierdo (sur)
            li = line.copy()
            li.Placement = frame.Placement
            li.Placement.Rotation = frame.Placement.Rotation
            li.Placement.Base.x -= wdt
            li.Placement.Base.z -= zz
            lines.append(li)

            # Borde derecho (norte)
            ld = line.copy()
            ld.Placement = frame.Placement
            ld.Placement.Rotation = frame.Placement.Rotation
            ld.Placement.Base.x += wdt
            ld.Placement.Base.z += zz
            lines.append(ld)

            tools.append([frame, li, ld])

        if len(lines) >= 2:
            try:
                loft = Part.makeLoft(lines, False, True, False)
                if loft and not loft.isNull():
                    all_lines.append(loft)
            except Exception:
                pass

    # Rellenar huecos entre columnas
    for group in rows:
        for frame in group:
            col, idx = _find_in_columns(frame, columns)
            tool = _find_tool(frame, tools)
            if tool is None:
                continue

            if idx < len(col) - 1:
                next_frame = col[idx + 1]
                next_tool = _find_tool(next_frame, tools)
                if next_tool:
                    try:
                        l1 = Part.LineSegment(tool[1].Vertexes[1].Point,
                                              next_tool[1].Vertexes[0].Point).toShape()
                        l2 = Part.LineSegment(tool[2].Vertexes[1].Point,
                                              next_tool[2].Vertexes[0].Point).toShape()
                        if len([l1, l2]) >= 2:
                            loft = Part.makeLoft([l1, l2], False, True, False)
                            if loft and not loft.isNull():
                                all_lines.append(loft)
                    except Exception:
                        pass

    if not all_lines:
        return None

    # Unir todas las caras en un solo sólido
    try:
        faces = []
        for item in all_lines:
            faces.extend(item.Faces)
        if faces:
            shell = Part.makeShell(faces)
            solid = Part.makeSolid(shell)
            return solid
    except Exception:
        return None


def _get_tracker_line(frame):
    """Obtiene la línea base longitudinal de un tracker."""
    try:
        lng = int(frame.Setup.Length / 2)
        return Part.LineSegment(FreeCAD.Vector(-lng, 0, 0),
                                FreeCAD.Vector(lng, 0, 0)).toShape()
    except Exception:
        # Fallback: usar el shape del setup
        try:
            setup_shape = frame.Setup.Shape
            bb = setup_shape.BoundBox
            return Part.LineSegment(FreeCAD.Vector(bb.XMin, 0, 0),
                                    FreeCAD.Vector(bb.XMax, 0, 0)).toShape()
        except Exception:
            return Part.LineSegment(FreeCAD.Vector(-2000, 0, 0),
                                    FreeCAD.Vector(2000, 0, 0)).toShape()


def _angle_to_prev(group, i):
    """Ángulo con el tracker anterior en la fila."""
    if i <= 0:
        return 0
    p0 = FreeCAD.Vector(group[i - 1].Placement.Base)
    p1 = FreeCAD.Vector(group[i].Placement.Base)
    return _angle_between(p0, p1)


def _angle_to_next(group, i):
    """Ángulo con el tracker siguiente en la fila."""
    if i >= len(group) - 1:
        return 0
    p1 = FreeCAD.Vector(group[i].Placement.Base)
    p2 = FreeCAD.Vector(group[i + 1].Placement.Base)
    return _angle_between(p1, p2)


def _angle_between(v1, v2):
    """Ángulo en el plano XZ entre dos vectores (grados)."""
    dx = v2.x - v1.x
    dz = v2.z - v1.z
    return math.degrees(math.atan2(dz, dx))


def _find_in_columns(frame, columns):
    """Busca un frame en la estructura de columnas."""
    for col in columns:
        for group in col:
            if frame in group:
                return group, group.index(frame)
    return [], -1


def _find_tool(frame, tools):
    """Busca el tool [frame, line_left, line_right] asociado a un frame."""
    for t in tools:
        if t[0] == frame:
            return t
    return None


def _mesh_above(design_mesh, terrain_mesh):
    """
    Devuelve el mesh de diseño que está POR ENCIMA del terreno (Cut).
    """
    try:
        common = design_mesh.common(terrain_mesh)
        if common and common.countPoints() > 0:
            return common
    except Exception:
        pass
    return None


def _mesh_below(design_mesh, terrain_mesh):
    """
    Devuelve el mesh de diseño que está POR DEBAJO del terreno (Fill).
    """
    try:
        # Mesh booleano: diseño - terreno = parte del diseño fuera del terreno
        diff = design_mesh.cut(terrain_mesh)
        if diff and diff.countPoints() > 0:
            return diff
    except Exception:
        pass
    return None


def _mesh_volume(mesh):
    """Calcula el volumen aproximado de un mesh (mm³)."""
    if mesh is None or mesh.countPoints() == 0:
        return 0
    try:
        return mesh.Volume
    except Exception:
        return 0


# =========================================================================
# Objeto EarthWorksVolume (FeaturePython)
# =========================================================================

class EarthWorksVolume(ArchComponent.Component):
    def __init__(self, obj):
        ArchComponent.Component.__init__(self, obj)
        self.obj = obj
        self.setProperties(obj)
        obj.Proxy = self

    def setProperties(self, obj):
        pl = obj.PropertiesList

        if "VolumeType" not in pl:
            obj.addProperty("App::PropertyEnumeration",
                            "VolumeType", "Volume",
                            "Fill o Cut").VolumeType = VOLUME_TYPES

        if "VolumeMesh" not in pl:
            obj.addProperty("Mesh::PropertyMeshKernel",
                            "VolumeMesh", "Volume", "Mesh del volumen")
        obj.setEditorMode("VolumeMesh", 2)

        if "Volume" not in pl:
            obj.addProperty("App::PropertyVolume",
                            "Volume", "Volume",
                            "Volumen calculado (mm³)")
        obj.setEditorMode("Volume", 1)

        obj.IfcType = "Civil Element"
        obj.setEditorMode("IfcType", 1)

    def onDocumentRestored(self, obj):
        ArchComponent.Component.onDocumentRestored(self, obj)
        self.setProperties(obj)

    def onChange(self, obj, prop):
        if prop == "VolumeMesh":
            if obj.VolumeMesh:
                obj.Volume = obj.VolumeMesh.Volume

    def execute(self, obj):
        pass

    def __getstate__(self):
        return None

    def __setstate__(self, state):
        return None


# =========================================================================
# ViewProvider (Coin3D visualization)
# =========================================================================

class ViewProviderEarthWorksVolume:
    def __init__(self, vobj):
        pl = vobj.PropertiesList
        is_cut = vobj.Object.VolumeType == "Cut"
        r, g, b = (1.0, 0.0, 0.0) if is_cut else (0.0, 0.0, 1.0)

        if "Transparency" not in pl:
            vobj.addProperty("App::PropertyIntegerConstraint",
                             "Transparency", "Surface Style",
                             "Transparencia de la superficie")
        vobj.Transparency = (50, 0, 100, 1)

        if "ShapeColor" not in pl:
            vobj.addProperty("App::PropertyColor",
                             "ShapeColor", "Surface Style",
                             "Color de superficie")
        vobj.ShapeColor = (r, g, b, vobj.Transparency / 100)

        if "ShapeMaterial" not in pl:
            vobj.addProperty("App::PropertyMaterial",
                             "ShapeMaterial", "Surface Style",
                             "Material de superficie")
        vobj.ShapeMaterial = FreeCAD.Material()

        vobj.Proxy = self
        vobj.ShapeMaterial.DiffuseColor = vobj.ShapeColor

    def onChanged(self, vobj, prop):
        if prop == "ShapeColor" or prop == "Transparency":
            if hasattr(vobj, "ShapeColor") and hasattr(vobj, "Transparency"):
                color = vobj.getPropertyByName("ShapeColor")
                t = vobj.getPropertyByName("Transparency")
                vobj.ShapeMaterial.DiffuseColor = (color[0], color[1], color[2], t / 100)

        if prop == "ShapeMaterial":
            if hasattr(vobj, "ShapeMaterial"):
                mat = vobj.getPropertyByName("ShapeMaterial")
                if hasattr(self, 'face_material'):
                    self.face_material.diffuseColor.setValue(mat.DiffuseColor[:3])
                    self.face_material.transparency = mat.DiffuseColor[3]

    def attach(self, vobj):
        from pivy import coin

        self.geo_coords = coin.SoGeoCoordinate()
        self.triangles = coin.SoIndexedFaceSet()
        self.face_material = coin.SoMaterial()
        self.edge_material = coin.SoMaterial()
        self.edge_style = coin.SoDrawStyle()
        self.edge_style.style = coin.SoDrawStyle.LINES

        shape_hints = coin.SoShapeHints()
        shape_hints.vertex_ordering = coin.SoShapeHints.COUNTERCLOCKWISE
        mat_binding = coin.SoMaterialBinding()
        mat_binding.value = coin.SoMaterialBinding.PER_FACE
        offset = coin.SoPolygonOffset()
        offset.styles = coin.SoPolygonOffset.LINES
        offset.factor = -2.0

        highlight = coin.SoType.fromName('SoFCSelection').createInstance()
        highlight.style = 'EMISSIVE_DIFFUSE'
        highlight.addChild(shape_hints)
        highlight.addChild(mat_binding)
        highlight.addChild(self.geo_coords)
        highlight.addChild(self.triangles)

        face = coin.SoSeparator()
        face.addChild(self.face_material)
        face.addChild(highlight)

        edge = coin.SoSeparator()
        edge.addChild(self.edge_material)
        edge.addChild(self.edge_style)
        edge.addChild(highlight)

        surface_root = coin.SoSeparator()
        surface_root.addChild(face)
        surface_root.addChild(offset)
        surface_root.addChild(edge)
        vobj.addDisplayMode(surface_root, "Surface")

        wireframe_root = coin.SoSeparator()
        wireframe_root.addChild(edge)
        vobj.addDisplayMode(wireframe_root, "Wireframe")

        self.onChanged(vobj, "ShapeColor")

    def updateData(self, obj, prop):
        if prop == "VolumeMesh":
            mesh = obj.VolumeMesh
            if mesh is None or mesh.countPoints() == 0:
                return
            try:
                geo_system = ["UTM", FreeCAD.ActiveDocument.Site.UtmZone, "FLAT"]
            except Exception:
                geo_system = ["UTM", "30N", "FLAT"]
            self.geo_coords.geoSystem.setValues(geo_system)

            copy_mesh = mesh.copy()
            triangles = []
            for i in copy_mesh.Topology[1]:
                triangles.extend(list(i))
                triangles.append(-1)

            self.geo_coords.point.setValues(copy_mesh.Topology[0])
            self.triangles.coordIndex.setValues(triangles)

    def getIcon(self):
        return str(os.path.join(DirIcons, "googleearth.svg"))

    def getDisplayModes(self, vobj):
        return ["Surface", "Wireframe"]

    def getDefaultDisplayMode(self):
        return "Surface"

    def setDisplayMode(self, mode):
        return mode

    def __getstate__(self):
        return None

    def __setstate__(self, state):
        return None


# =========================================================================
# TaskPanel
# =========================================================================

class EarthWorksTaskPanel:
    def __init__(self):
        self.form = FreeCADGui.PySideUic.loadUi(
            os.path.join(PVPlantResources.__dir__, "PVPlantEarthworks.ui"))
        self.form.setWindowIcon(
            QtGui.QIcon(os.path.join(PVPlantResources.DirIcons, "convert.svg")))

    def accept(self):
        import MeshPart
        land = FreeCAD.ActiveDocument.Terrain.Mesh.copy()

        frames = []
        for obj in FreeCADGui.Selection.getSelection():
            if hasattr(obj, "Proxy"):
                proxy_type = getattr(obj.Proxy, "Type", None)
                if proxy_type == "Tracker":
                    if obj not in frames:
                        frames.append(obj)
                elif proxy_type == "FrameArea":
                    for fr in obj.Frames:
                        if fr not in frames:
                            frames.append(fr)

        if not frames:
            FreeCAD.Console.PrintWarning("Selecciona trackers o un FrameArea\n")
            return False

        FreeCAD.ActiveDocument.openTransaction("Movimiento de tierras")

        try:
            cut_mesh, fill_mesh, vol_cut, vol_fill = compute_earthworks(
                frames, land,
                slope_tolerance=getattr(FreeCAD.ActiveDocument,
                                        'MaximumWestEastSlope', 10.0)
            )

            if cut_mesh and cut_mesh.countPoints() > 0:
                vol_obj = makeEarthWorksVolume(1)  # Cut
                vol_obj.VolumeMesh = cut_mesh
                FreeCAD.Console.PrintMessage(f"Volumen corte: {vol_cut:.0f} mm³\n")

            if fill_mesh and fill_mesh.countPoints() > 0:
                vol_obj = makeEarthWorksVolume(0)  # Fill
                vol_obj.VolumeMesh = fill_mesh
                FreeCAD.Console.PrintMessage(f"Volumen relleno: {vol_fill:.0f} mm³\n")

        except Exception as e:
            FreeCAD.Console.PrintError(f"Error en movimiento de tierras: {e}\n")
        finally:
            FreeCAD.ActiveDocument.commitTransaction()

        FreeCADGui.Control.closeDialog()
        return True

    def reject(self):
        FreeCADGui.Control.closeDialog()
        return True


# Comando registrado desde Init.py o InitGui.py