# Python Script, API Version = V261

import math

#method that evaluates the coordinates of the centroid of a face
def face_coord(face):
    
    # Measure    x
    selection = Selection.Create([face,
        GetRootPart().CoordinateSystems[0]])
    gap = MeasureHelper.MinDistanceBetweenObjects(selection)
    deltaX = gap.DeltaX
    # EndBlock

    # Measure    y
    selection = Selection.Create([face,
        GetRootPart().CoordinateSystems[0]])
    gap = MeasureHelper.MinDistanceBetweenObjects(selection)
    deltaY = gap.DeltaY
    # EndBlock

    # Measure    z
    selection = Selection.Create([face,
        GetRootPart().CoordinateSystems[0]])
    gap = MeasureHelper.MinDistanceBetweenObjects(selection)
    deltaZ = gap.DeltaZ
    # EndBlock

    return deltaX,deltaY,deltaZ 
#method that evaluates the coordiante limitations (x_min,x_max_y_min etc...) of the region from which we want to extract the solution data
def box_limits(body):
    
    faces=body.Faces
    
    x_max=face_coord(faces[0])[0]
    x_min=face_coord(faces[0])[0]

    y_max=face_coord(faces[0])[1]
    y_min=face_coord(faces[0])[1]

    z_max=face_coord(faces[0])[2]
    z_min=face_coord(faces[0])[2]

    

    for face in faces:
        
        if face_coord(face)[0] > x_max:
            x_max=face_coord(face)[0]
            
        if (face_coord(face)[1]) > (y_max):
            y_max=face_coord(face)[1]
            
        if (face_coord(face)[2]) > (z_max):
            z_max=face_coord(face)[2]
            
        if (face_coord(face)[0]) < (x_min):
            x_min=face_coord(face)[0]
            
        if (face_coord(face)[1]) < (y_min):
            y_min=face_coord(face)[1]
            
        if (face_coord(face)[2]) < (z_min):
            z_min=face_coord(face)[2]
    return(x_max,x_min,y_max,y_min,z_max,z_min)

def generate_tetrahedral_points_box(xmin, xmax, ymin, ymax, zmin, zmax, h):
    """
    Generate a tetrahedral-style point distribution inside a box.

    Parameters
    ----------
    xmin, xmax, ymin, ymax, zmin, zmax : float
        Box limits
    h : float
        Characteristic spacing (controls density)

    Returns
    -------
    List of [x, y, z] points
    """

    points = []

    # Vertical spacing for regular tetrahedral packing
    hz = math.sqrt(2.0 / 3.0) * h
    hy = math.sqrt(3.0) / 2.0 * h

    k = 0
    z = zmin
    while z <= zmax:

        j = 0
        y = ymin
        while y <= ymax:

            i = 0
            x = xmin
            while x <= xmax:

                # Alternating offsets (tetrahedral arrangement)
                dx = 0.5 * h if (j + k) % 2 == 1 else 0.0
                dy = hy / 3.0 if k % 2 == 1 else 0.0

                px = x + dx
                py = y + dy
                pz = z

                # Keep strictly inside box
                if (xmin <= px <= xmax and
                    ymin <= py <= ymax and
                    zmin <= pz <= zmax):
                    pt=Point.Create(px,py,pz)
                    points.Add(pt)

                i += 1
                x = xmin + i * h

            j += 1
            y = ymin + j * hy

        k += 1
        z = zmin + k * hz

    return points

#Entire Fluid Domain
fluid=GetRootPart().GetAllBodies()[0]

#QuantityofInterest
resultVariable = Results.ResultVariable.Density

#Discretisation of this region (in meters) 
h=0.0005

for box in GetRootPart().Components[1].Content.Bodies:
    #Create a list of points where results are extracted
    points_list=List[Point]()

    x_max=box_limits(box)[0]
    x_min=box_limits(box)[1]

    y_max=box_limits(box)[2]
    y_min=box_limits(box)[3]

    z_max=box_limits(box)[4]
    z_min=box_limits(box)[5]

    points_list = generate_tetrahedral_points_box(x_min, x_max,y_min, y_max,z_min, z_max,h)

    #Extract results at given points list

    result=ResultsExport.GetResultsForPoints(fluid,points_list,resultVariable) 
    value=result.Values #Values

    print(value)

    #Write in  a csv in the same file directory as the discovery file
    import csv
    import os

    file=os.path.join(os.path.dirname(GetActiveWindow().Document.Path), str(resultVariable)+"_"+box.GetName()+".csv")
    k=0
    with open(file, "w") as f:
        f.write("X,Y,Z,Variable,Value,Unit\n")
        
        for i in range(len(points_list)):
            f.write("{},{},{},Density,{},kg/m3\n".format(points_list[i].X,points_list[i].Y,points_list[i].Z,value[i]))
        print("done")
