Fix #4848 Added STL Importing (#4967)

sirepo/package_data/static/js/radia.js

@@ -3699,7 +3699,7 @@ SIREPO.app.directive('shapeButton', function(appState, geometry, panelState, plo
 
 SIREPO.app.directive('shapeSelector', function(appState, panelState, plotting, radiaService, utilities) {
 
-    const availableShapes = ['cuboid', 'cylinder', 'ell', 'cee', 'jay', 'extrudedPoints',];
+    const availableShapes = ['cuboid', 'cylinder', 'ell', 'cee', 'jay', 'extrudedPoints', 'stl'];
     const sel = new SIREPO.DOM.UISelect('', [
         new SIREPO.DOM.UIAttribute('data-ng-model', 'field'),
     ]);

sirepo/package_data/static/json/radia-schema.json

@@ -107,6 +107,14 @@
                     [0, 0], [24, 0], [24, 12], [18, 12], [18, 24], [6, 24], [6, 16], [6, 12], [0, 12]
                 ],
                 "stroke": "black"
+            },
+            "stl": {
+                "doClose": true,
+                "fill": "#97B6D5",
+                "points": [
+                    [0, 0], [24, 0], [24, 24], [0, 24]
+                ],
+                "stroke": "black"
             }
         },
         "geomTypeLines": "lines",
@@ -180,6 +188,13 @@
                         "model": "geomGroup",
                         "title": "Group",
                         "type": "geomGroup"
+                    },
+                    {
+                        "isButton": false,
+                        "layoutShape": "polygon",
+                        "model": "stl",
+                        "title": "STL",
+                        "type": "stl"
                     }
                 ],
                 "name": "Objects"
@@ -346,7 +361,8 @@
             ["extrudedPoints", "Extruded Points"],
             ["jay", "Extruded J-Bend"],
             ["racetrack", "RaceTrack"],
-            ["tee", "Extruded Tee"]
+            ["tee", "Extruded Tee"],
+            ["stl", "STL"]
         ],
         "PathType": [
             ["axis", "Axis"],
@@ -934,6 +950,11 @@
             "sbatchCores": ["Cores", "Integer", 1],
             "sbatchProject": ["Project", "OptionalString", ""]
         },
+        "stl": {
+            "_super": ["_", "model", "geomObject"],
+            "type": ["_", "ObjectType", "stl"],
+            "file": ["STL File", "InputFile", ""]
+        },
         "fieldLineoutAnimation": {
             "fieldPaths": ["Path", "FieldPaths", []],
             "fieldType": ["Field", "FieldType", "B"],
@@ -1375,7 +1396,8 @@
                     "stemmed.stemWidth",
                     "stemmed.stemPosition",
                     "jay.hookHeight",
-                    "jay.hookWidth"
+                    "jay.hookWidth",
+                    "stl.file"
                 ]],
                 ["Point Editor", [
                     "extrudedPoints.pointsFile",

sirepo/package_data/template/radia/geometryReport.py.jinja

@@ -20,9 +20,10 @@ _MODEL_GROUP = 'geomGroup'
 _MODEL_OBJECT = 'geomObject'
 _MODEL_EXTRUDED_POLY = 'extrudedPoly'
 _MODEL_RACETRACK = 'racetrack'
+_MODEL_STL = 'stl'
 _EXTRUDED_MODELS = [_MODEL_CEE, _MODEL_ELL, _MODEL_JAY, _MODEL_EXTRUDED_POLY]
 _GROUP_MODELS = [_MODEL_GROUP]
-_OBJ_MODELS = [_MODEL_CUBOID, _MODEL_CYLINDER, _MODEL_RACETRACK, _MODEL_OBJECT] + _EXTRUDED_MODELS
+_OBJ_MODELS = [_MODEL_CUBOID, _MODEL_CYLINDER, _MODEL_RACETRACK, _MODEL_OBJECT, _MODEL_STL] + _EXTRUDED_MODELS
 _MODELS = _OBJ_MODELS + _GROUP_MODELS
 
 sirepo_objs = {{ objects }}
@@ -85,6 +86,22 @@ def _add_object(o, radia_objs, id_map):
             num_segs=o.numSegments,
             curr_density=o.currentDensity,
         )
+    if t == _MODEL_STL:
+        g_id = radia_util.build_stl(
+            filename = o.get('file'),
+            #TODO(BG) Not currently reflective
+            center=ctr,
+            #TODO(BG) Not currently reflective
+            size=sz,
+            vertices=o.stlVertices,
+            faces=o.stlFaces,
+            magnetization=m,
+            material=o.material,
+            rem_mag=o.remanentMag,
+            segments=segs, #TODO: Not implemented
+            h_m_curve=o.h_m_curve,
+            #slices=o.stlSlices TODO: Not implemented
+        )
     if t in _OBJ_MODELS:
         for b in o.get('bevels', []):
             g_id = radia_util.apply_bevel(

sirepo/template/radia.py

@@ -7,8 +7,8 @@
 :copyright: Copyright (c) 2017-2018 RadiaSoft LLC.  All Rights Reserved.
 :license: http://www.apache.org/licenses/LICENSE-2.0.html
 """
-import inspect
 
+import inspect
 from pykern import pkcompat
 from pykern import pkinspect
 from pykern import pkio
@@ -30,6 +30,7 @@
 import sirepo.csv
 import sirepo.sim_data
 import sirepo.util
+import trimesh
 import uuid
 
 _AXES_UNIT = [1, 1, 1]
@@ -346,6 +347,18 @@ def python_source_for_model(data, model):
     return _generate_parameters_file(data, False, for_export=True)
 
 
+def validate_file(file_path, path):
+    if path.ext == ".stl":
+        mesh = _create_stl_trimesh(path)
+        if trimesh.convex.is_convex(mesh) == False:
+            return f"not convex model: {path.basename}"
+        elif len(mesh.faces) > 600:
+            return f"too many faces({len(mesh.faces)}): {path.basename}"
+    else:
+        return f"invalid file type: {path.ext}"
+    return None
+
+
 def write_parameters(data, run_dir, is_parallel):
     pkio.write_text(
         run_dir.join(template_common.PARAMETERS_PYTHON_FILE),
@@ -577,6 +590,29 @@ def _build_undulator_objects(geom_objs, model, **kwargs):
     return _update_geom_from_undulator(geom_objs, model, **kwargs)
 
 
+def _is_binary(file_path):
+    textchars = bytearray({7, 8, 9, 10, 12, 13, 27} | set(range(0x20, 0x100)) - {0x7F})
+    is_binary_string = lambda bytes: bool(bytes.translate(None, textchars))
+    return is_binary_string(open(file_path, "rb").read(1024))
+
+
+def _create_stl_trimesh(file_path):
+    readParam = "r"
+    keyType = "ascii"
+    if _is_binary(file_path):
+        readParam = "rb"
+        keyType = "binary"
+    with open(file_path, readParam) as f:
+        m = trimesh.exchange.stl.load_stl(file_obj=f)
+        if "geometry" in m:
+            return trimesh.Trimesh(
+                vertices=m["geometry"][keyType]["vertices"],
+                faces=m["geometry"][keyType]["faces"],
+                process=True,
+            )
+        return trimesh.Trimesh(vertices=m["vertices"], faces=m["faces"], process=True)
+
+
 # deep copy of an object, but with a new id
 def _copy_geom_obj(o):
     import copy
@@ -1620,6 +1656,10 @@ def _poly_area(pts):
         magnetization=[0.0, 0.0, 0.0],
         segments=[1, 1, 1],
         size=[1.0, 1.0, 1.0],
+        stlVertices=[],
+        stlFaces=[],
+        # TODO(BG) Not implemented
+        # stlSlices = [],
     )
     for k in d:
         v = kwargs.get(k)
@@ -1644,6 +1684,42 @@ def _poly_area(pts):
         )
     if "points" in o:
         o.area = _poly_area(o.points)
+    if o.type == "stl":
+        path = str(
+            _SIM_DATA.lib_file_abspath(
+                _SIM_DATA.lib_file_name_with_type(o.file, "stl-file")
+            )
+        )
+        mesh = _create_stl_trimesh(path)
+        for v in list(mesh.vertices):
+            d.stlVertices.append(list(v))
+        for f in list(mesh.faces):
+            d.stlFaces.append(list(f))
+        o.stlVertices = d.stlVertices
+        o.stlFaces = d.stlFaces
+
+        # TODO(BG) Mesh slicing implementation, option for meshes with 400+ faces although will be approximation
+        """
+        z_extents = mesh.bounds[:,2]
+        z_levels  = numpy.arange(*z_extents, step=1)
+        meshSlices = trimesh.intersections.mesh_multiplane(mesh=mesh, plane_origin=mesh.bounds[0], plane_normal=[0,0,1], heights=z_levels)[0]
+        formattedSlices = []
+        index = 0
+        for s in meshSlices:
+            slicePoints = []
+            for l in s:
+                for p in l:
+                    #Remove redundant points by rounding
+                    p[0] = round(p[0],5)
+                    p[1] = round(p[1],5)
+                    if list(p) not in slicePoints:
+                        slicePoints.append(list(p))
+            formattedSlices.append([list(slicePoints), z_levels[index]])
+            index += 1
+        for s in formattedSlices:
+            s[0] = sort_points_clockwise(s[0])
+        o.stlSlices = formattedSlices
+        """
     return o
 
 
@@ -1695,6 +1771,32 @@ def _update_kickmap(km, und, beam_axis):
     km.periodLength = und.periodLength
 
 
+# TODO(BG) Necessary helper function to implement object slicing with radia.radObjMltExtPgn()
+# Edge Case: Need to remove linear points along same vecter before returning
+"""
+def _sort_points_clockwise(points):
+    angles = []
+    center = (sum([p[0] for p in points]) / len(points), sum([p[1] for p in points]) / len(points))
+    for p in points:
+        vector = [p[0] - center[0], p[1] - center[1]]
+        vlength = math.sqrt(pow(vector[0], 2) + pow(vector[1], 2))
+        if vlength == 0:
+            angles.append(-numpy.pi)
+        else:
+            normalized = [vector[0] / vlength, vector[1] / vlength]
+            angle = math.atan2(normalized[0], normalized[1])
+            # function checks against x-positive, if negative angle add to 2pi for mirror
+            if angle < 0:
+                angle = 2 * numpy.pi + angle
+            angles.append(angle)
+    for i in range(len(angles)):
+        angles[i] = [angles[i], points[i]]
+    angles.sort()  
+    # might have to check by lengths as well if angles are the same
+    return [x[1] for x in angels]
+"""
+
+
 def _validate_objects(objects):
     import numpy.linalg
 

sirepo/template/radia_util.py

@@ -49,7 +49,6 @@
     }
 )
 
-
 _MU_0 = 4 * numpy.pi / 1e7
 _ZERO = [0, 0, 0]
 
@@ -300,6 +299,15 @@ def build_cuboid(**kwargs):
     return g_id
 
 
+def build_stl(**kwargs):
+    d = PKDict(kwargs)
+    g_id = radia.ObjPolyhdr(
+        d.vertices, (numpy.array(d.faces) + 1).tolist(), d.magnetization
+    )
+    radia.MatApl(g_id, _radia_material(d.material, d.rem_mag, d.h_m_curve))
+    return g_id
+
+
 def build_cylinder(**kwargs):
     d = PKDict(kwargs)
     axis = d.extrusion_axis