Refactor custom surface properties using facets

[danieljfarrell]

Surface reflection is currently handled with a delegate object approach. A node owns a geometry, the geometry owns material, the material owns a surface object, the surface object owns the delegate.

lsc = Node(
    name="LSC",
    geometry=Box(
      (l, w, d),
      material=Material(
          refractive_index=1.5, components=[],
          surface=Surface(delegate=OptionalMirrorAndSolarCell(self))
      ),
  ),
  parent=world,
)

Here OptionalMirrorAndSolarCell is an object that handles all surface interactions for the full surface area of the object. As you can see from the naming of the delegate it seems to have too much responsibility! The code for the delegate is similarly complicated. The chief problem is that the level of abstraction is too high.

Instead of providing one object to handle custom surface interactions we could supply multiple Facet objects which opt-in an arbitrarily defined sub-surface of the object. For example, for a the box geometry this could be full top surface or a part of the top surface.

lsc = Node(
    name="LSC",
    geometry=Box(
        (l, w, d),
        material=Material(
            refractive_index=1.5, components=...,
        ),
        facets=[
            Facet(
                label="left",
                where=lambda point, geometry: np.isclose(point[0], -geometry.size[0]/2, atol=EPS_ZERO),
                kind=SolarCell()
            ),
        ]
    ),
    parent=world,
)

Define a node with Fresnel reflections based on the material's refractive,

lsc = Node(
    name="LSC",
    geometry=Box(
        (l, w, d),
        material=Material(
            refractive_index=1.5, components=...,
        ),
    ),
    parent=world,
)

Modify so that the left surface has a solar cell,

lsc = Node(
    name="LSC",
    geometry=Box(
        (l, w, d),
        material=Material(
            refractive_index=1.5, components=...,
        ),
        facets=[
            Facet(
                label="left",
                where=lambda point, geometry: np.isclose(point[0], -geometry.size[0]/2, atol=EPS_ZERO),
                kind=SolarCell()
            ),
        ]
    ),
    parent=world,
)

Modify so that the bottom surface has a mirror,

lsc = Node(
    name="LSC",
    geometry=Box(
        (l, w, d),
        material=Material(
            refractive_index=1.5, components=...,
        ),
        facets=[
            Facet(
                label="left",
                where=lambda point, geometry: np.isclose(point[0], -geometry.size[0]/2, atol=EPS_ZERO),
                kind=SolarCell()
            ),
            Facet(
                label="left",
                where=lambda point, geometry: np.isclose(point[0], -geometry.size[0]/2, atol=EPS_ZERO),
                kind=Mirror()
            ),
        ]
    ),
    parent=world,
)

Modify so that the top surface has a coating with spectral dependence,

lsc = Node(
    name="LSC",
    geometry=Box(
        (l, w, d),
        material=Material(
            refractive_index=1.5, components=...,
        ),
        facets=[
            Facet(
                label="left",
                where=lambda point, geometry: np.isclose(point[0], -geometry.size[0]/2, atol=EPS_ZERO),
                kind=SolarCell()
            ),
            Facet(
                label="bottom",
                where=lambda point, geometry: np.isclose(point[2], -geometry.size[2]/2, atol=EPS_ZERO),
                kind=Mirror()
            ),
            Facet(
                label="top",
                where=lambda point, geometry: np.isclose(point[2],  geometry.size[2]/2, atol=EPS_ZERO),
                kind=Coating()
            ),
        ]
    ),
    parent=world,
)

The Facet object has a user defined label. The where parameter is a function which the facet will call to determine if a point is on the facet and therefore should be handled or not. The kind parameter is an object which contains custom reflection, transmission, scattering, diffraction or absorption information.

Here the "kind" object (need a better name) must be provided with simply an angle (with respect to the surface normal at the hit point) and a wavelength. This way it can be very simple -- almost just a look-up function of angular and spectral reflectivity. The facet object is responsible for handling any coordinate system transformations and providing the data to the "kind".

[danieljfarrell]

Feature roadmap

@shomikverma has made a good suggestion on his pvtrace fork that it would be easier to describe facets by surface normal.

In addition to the where keyword we can also include a normal keyword

Facet(
    label="top",
    normal=[0,0,1]
    kind=Coating()
)

which accepts a normal vector describing the facet.

Box

This works well for facets that have a single normal such as the box geometry. A user-friendly way by adding a class variables facets to the concrete Geometry subclasses,

class BoxNormals(object):
    left =   [-1, 0,  0]
    right =  [1,  0,  0]
    near =   [0, -1,  0]
    far =    [0,  1,  0]
    bottom = [0,  0, -1]
    top =    [0,  0,  1]

class Box(Geometry):
    
    normals = BoxNormals()
    ...

When creating facet object for a specific geometry they can now be easily specified,

Facet(
    label="top",
    normal=Box.normals.top
    kind=Coating()
)

This can be used in combination with the where keyword to define a sub-surface,

Facet(
    label="top",
    normal=Box.normals.top,
    where=Box.coordinate.x > L/2 and Box.coordinate.y > W/2
    kind=Coating()
)

Mesh

This approach should also work for mesh geometries. All faces (i.e. a single triangle) of meshes have a unique identifier. This would be a fairly general way of allowing users to define arbitrary facets.

# Let's say a facet of interest on the mesh is described by 4 faces of with IDs
face_ids = [99, 100, 101, 102]

Facet(
    label="edge",
    faces=face_ids,
    kind=Coating()
)

Here the Facet object has an additional keyword faces which is only to be used with Mesh geometry. This approach is nice because it is very general, the downside is the user will have to be familiar with CAD software and understand how to inspect meshes to extract this metadata.

A better approach would be able to "tag" certain faces with a single identifier in the CAD software and have pvtrace read those values.

It is also possible to calculate the surface normal of any mesh face. Therefore the user should also also be able to use where and normal as with the prior case. However, the usefulness of this very much depends on details of the mesh itself. Moreover, this approach is not useful for a mesh which his describes a curved facet because the surface normal will vary continuously across the facet.

Sphere and cylinder

We also need this approach to work for geometry with curved surfaces, such as spheres and cylinders, and with meshes which describe curves. This is more difficult because the surface normals need to be described by ranges in a coordinate system.

To describe a curved facet on a sphere we can define angle ranges,

Facet(
    label="top quadrant",
    normal=(
        Sphere.coordinate.phi > 0 and 
        Sphere.coordinate.phi < pi/2 and 
        Sphere.coordinate.theta > 0 and 
        Sphere.coordinate.theta < pi/2
    ),
    kind=Coating()
)

For a cylinder,

Facet(
    label="top",
    normal=(
        Cylinder.coordinate.z > -D/2 and 
        Cylinder.coordinate.z < D/2 and 
        Cylinder.coordinate.theta > -pi/2 and 
        Cylinder.coordinate.theta < pi/2
    ),
    kind=Coating()
)

Delayed evaluation of coordinates

Some Python magic is needed to implement lines like these,

where=Box.coordinate.x > L/2 and Box.coordinate.y > W/2

This describes a delayed evaluation of a ray hit location with box in the box's local coordinate system.

At runtime the value placeholder Box.coordinate.x will be replaced with the ray's x coordinate and the expression ray_x > L/2 will be evaluated.

This would be nice to have. Essentially it is syntaxic sugar for creating the lambda function,

where=lambda box, point: point[0] > box.size[0]/2 and point[1] > box[1]/2

Counter facet

We should also introduce a facet kind of type Counter which simply counts the rays crossing the area is describes. This will be very useful for labelling where rays end-up and generating useful statistics at the end of the simulation.

Facet(
    label="top",
    normal=Box.normals.top,
    where=Box.coordinate.x > L/2 and Box.coordinate.y > W/2
    kind=Counter()
)