Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Triangular Prism Geometry for LSC #56

Open
tommyflynn13 opened this issue Jul 16, 2022 · 20 comments
Open

Triangular Prism Geometry for LSC #56

tommyflynn13 opened this issue Jul 16, 2022 · 20 comments

Comments

@tommyflynn13
Copy link

Hi Daniel,

Firstly just wanted to say thanks a million for the work you've done on this code, the software is very impressive and easy to follow for a beginner like myself! I was shown your code by a fellow university student, and am hoping to manipulate it to simulate a triangular LSC panel being fabricated for use in a geodesic dome installation as part of my Masters dissertation.

Not an issue as such, but I was hoping to pick your brain as to how you would edit the LSC.py code to work with a triangular prism rather than a box? As I mentioned I am a complete beginner with Python, I have figured out how to extrude a triangular polygon using shapely and trimesh, and have passed light through it using one of your tutorials. I am struggling however to fully adapt the LSC.py code to work with this mesh geometry, specifically as the "size" and "self.size" parameters don't seem to apply to the (polygon, height) input using trimesh.creation.extrude_polygon.

Any help would be greatly appreciated.

Regards,
Thomas.

@danieljfarrell
Copy link
Owner

danieljfarrell commented Jul 17, 2022

Hi Tommy,

That's a very interesting project, with lots of angles to investigate. I remember Goetzberger studied triangular geometries in one of his early papers, might be worth a look.

I don't know how amenable the LSC class is to extension to other geometries, but definitely worth looking in to. It was not something I considered when writing it.

The usual pvtrace trace approach is to, define a geometry, material and a light source and to throw rays into the screen and count where they come out.

I think you are on the right track starting with Mesh, later down the line you might want to write your own Triangle geometry class to handle the intersections yourself.

(Edit)

OK, so you have stuff more or less working. But probably I need to see the code to provide useful comments. If you have a fork, post a link to the line in question and the error is gives you.

@tommyflynn13
Copy link
Author

tommyflynn13 commented Jul 17, 2022

No bother at all, I've been working through Jupyter and I'll attach the code below. Apologies if it's a bit messy, I'm picking it up on the fly!

Firstly, to create the polygon:

import trimesh
import shapely
from shapely.geometry import Point, Polygon
polygon = Polygon([[0, 0], [1, 2], [2, 0]]) # polygon

To extrude and visualise the polygon:

import logging
logging.getLogger('pvtrace').setLevel(logging.CRITICAL)
logging.getLogger('trimesh').setLevel(logging.CRITICAL)
from pvtrace import *
import trimesh
import shapely


# always need a world node
world = Node(
    name="world",
    geometry=Sphere(
        radius=100.0
    ),
)

# Triangular Prism
mesh = Node(
    name="mesh (triangular prism)",
    geometry=Mesh(
        trimesh=trimesh.creation.extrude_polygon(polygon, 0.5)
    ),
    parent=world
)
mesh.translate((0.0, 0.0, 0.5))
scene = Scene(world)
vis = MeshcatRenderer(wireframe=True)
vis.render(scene)
vis.vis.jupyter_cell()

Finally, the adapted LSC.py code:

from pvtrace.material.component import Absorber, Luminophore
from pvtrace.light.light import Light, rectangular_mask
from pvtrace.light.event import Event
from pvtrace.scene.node import Node
from pvtrace.material.material import Material
from pvtrace.material.utils import isotropic, cone, lambertian
from pvtrace.scene.scene import Scene
from pvtrace.geometry.box import Box
from pvtrace.geometry.utils import EPS_ZERO
from pvtrace.data import lumogen_f_red_305, fluro_red
from pvtrace.scene.renderer import MeshcatRenderer
from pvtrace.material.surface import Surface, FresnelSurfaceDelegate
from pvtrace.material.distribution import Distribution
from pvtrace.algorithm import photon_tracer
from dataclasses import asdict
import numpy as np
import pandas as pd
import functools
import time


class OptionalMirrorAndSolarCell(FresnelSurfaceDelegate):
    """ A delegate adds an ideal specular mirror to the bottom surface and 
        perfectly indexed matched and perfectly absorbing solar cells to the edges.
    """

    def __init__(self, lsc):
        super(OptionalMirrorAndSolarCell, self).__init__()
        self.lsc = lsc

    def reflectivity(self, surface, ray, geometry, container, adjacent):
        normal = geometry.normal(ray.position)
        cell_locations = self.lsc._solar_cell_surfaces
        back_surface_mirror = self.lsc._back_surface_mirror_info
        want_back_surface_mirror = back_surface_mirror["want_back_surface_mirror"]
        if np.allclose((0, 0, -1), normal) and want_back_surface_mirror:
            return 1.0  # perfect mirror
        elif np.allclose((-1, 0, 0), normal) and "left" in cell_locations:  # left
            return 0.0  # perfect absorption
        elif np.allclose((1, 0, 0), normal) and "right" in cell_locations:  # right
            return 0.0
        elif np.allclose((0, -1, 0), normal) and "near" in cell_locations:  # near
            return 0.0
        elif np.allclose((0, 1, 0), normal) and "far" in cell_locations:  # far
            return 0.0
        return super(OptionalMirrorAndSolarCell, self).reflectivity(
            surface, ray, geometry, container, adjacent
        )  # opt-out of handling custom reflection

    def transmitted_direction(self, surface, ray, geometry, container, adjacent):
        cell_locations = self.lsc._solar_cell_surfaces
        normal = geometry.normal(ray.position)
        if (
            (np.allclose((-1, 0, 0), normal) and "left" in cell_locations)
            or (np.allclose((1, 0, 0), normal) and "right" in cell_locations)
            or (np.allclose((0, -1, 0), normal) and "near" in cell_locations)
            or (np.allclose((0, 1, 0), normal) and "far" in cell_locations)
        ):
            return ray.direction  #  solar cell is perfectly index matched
        return super(OptionalMirrorAndSolarCell, self).transmitted_direction(
            surface, ray, geometry, container, adjacent
        )  # opt-out of handling custom reflection


class AirGapMirror(FresnelSurfaceDelegate):
    def __init__(self, lsc):
        super(AirGapMirror, self).__init__()
        self.lsc = lsc

    def reflectivity(self, surface, ray, geometry, container, adjacent):
        return 1.0  # perfect reflector. NB don't reduce.

    def reflected_direction(self, surface, ray, geometry, container, adjacent):
        if not self.lsc._air_gap_mirror_info["lambertian"]:
            # Specular reflection
            return super(AirGapMirror, self).transmitted_direction(
                surface, ray, geometry, container, adjacent
            )
        else:
            normal = geometry.normal(ray.position)
            if not np.allclose((0.0, 0.0, 1.0), normal):  # top surface
                raise NotImplementedError("Not yet generalised to other surfaces.")
            # Currently this return lambertian direction along +z axis and is not
            # generalised to other orientations. This is simple to do using a transform
            # which first moves into to the z+ frame and then back out.
            return tuple(lambertian().tolist())


class LSC(object):
    """Abstraction of a luminescent solar concentrator.
    
       This is intended to be a high-level API to easy use.
    """

    def __init__(self, wavelength_range=None, n0=1.0, n1=1.5):
        super(LSC, self).__init__()
        if wavelength_range is None:
            self.wavelength_range = np.arange(400, 800)

        self.n0 = n0
        self.n1 = n1

        self._solar_cell_surfaces = set()
        self._back_surface_mirror_info = {"want_back_surface_mirror": False}
        self._air_gap_mirror_info = {"want_air_gap_mirror": False, "lambertian": False}
        self._scene = None
        self._renderer = None
        self._store = None
        self._df = None
        self._counts = None
        self._user_lights = []
        self._user_components = []

    def _make_default_components(self):
        """ Default LSC contains Lumogen F Red 305. With concentration such that
            the absorption coefficient at peak is 10 cm-1.
        """
        x = self.wavelength_range
        coefficient = lumogen_f_red_305.absorption(x) * 10.0  # cm-1
        emission = lumogen_f_red_305.emission(x)
        coefficient = np.column_stack((x, coefficient))
        emission = np.column_stack((x, emission))
        lumogen = {
            "cls": Luminophore,
            "name": "Lumogen F Red 305",
            "coefficient": coefficient,
            "emission": emission,
            "quantum_yield": 1.0,
            "phase_function": None,  # will select isotropic
        }
        background = {"cls": Absorber, "coefficient": 0.1, "name": "Background"}  # cm-1
        return [lumogen, background]

    def _make_default_lights(self):
        """ Default light is a spotlight (cone of 20-deg) of single wavelength 555nm.
        """
        light = {
            "name": "Light",
            "location": (0.0, 0.0, self.size[-1] * 5),
            "rotation": (np.radians(180), (1, 0, 0)),
            "direction": functools.partial(cone, np.radians(20)),
            "wavelength": None,
            "position": None,
        }
        return [light]

    def _make_scene(self):
        """ Creates the scene based on configuration values.
        """
        # Make world node
        world = Node(
            name="World",
            geometry=Box(
                (100, 100, 100), material=Material(refractive_index=self.n0)
            ),
        )

        # Create components (Absorbers, Luminophores and Scatteres)
        if len(self._user_components) == 0:
            self._user_components = self._make_default_components()
        components = []
        for component_data in self._user_components:
            cls = component_data.pop("cls")
            coefficient = component_data.pop("coefficient")
            component = cls(coefficient, **component_data)
            components.append(component)

        # Create LSC node
        lsc = Node(
            name="LSC",
            geometry=Mesh(
                trimesh=trimesh.creation.extrude_polygon(polygon, 0.5),
                material=Material(
                    refractive_index=self.n1,
                    components=components,
                    surface=Surface(delegate=OptionalMirrorAndSolarCell(self)),
                ),
            ),
            parent=world,
        )

        if self._air_gap_mirror_info["want_air_gap_mirror"]:
            sheet_thickness = 0.1  # make it appear thinner than the LSC
            air_gap_mirror = Node(
                name="Air Gap Mirror",
                geometry=Mesh(
                    trimesh=trimesh.creation.extrude_polygon(polygon, 0.1),  # same surface air but very thin
                    material=Material(
                        refractive_index=self.n0,
                        components=[],
                        surface=Surface(delegate=AirGapMirror(self)),
                    ),
                ),
                parent=world,
            )
            # Move adjacent to bottom surface with a small air gap
            air_gap_mirror.translate((0.0, 0.0, -(0.5 * 0.5 + 0.1)))

        # Use user light if any have been given, otherwise use default values.
        if len(self._user_lights) == 0:
            self._user_lights = self._make_default_lights()

        # Create light nodes
        for light_data in self._user_lights:
            name = light_data["name"]
            light = Light(
                name=name,
                direction=light_data["direction"],
                wavelength=light_data["wavelength"],
                position=light_data["position"],
            )
            light_node = Node(name=name, light=light, parent=world)
            light_node.location = light_data["location"]
            if light_data["rotation"]:
                light_node.rotate(*light_data["rotation"])

        self._scene = Scene(world)

    def component_names(self):
        if self._scene is None:
            raise ValueError("Run a simulation before calling this method.")
        return {c["name"] for c in self._user_components}

    def light_names(self):
        if self._scene is None:
            raise ValueError("Run a simulation before calling this method.")
        return {l["name"] for l in self._user_lights}

    def add_luminophore(
        self, name, coefficient, emission, quantum_yield, phase_function=None
    ):
        self._user_components.append(
            {
                "cls": Luminophore,
                "name": name,
                "coefficient": coefficient,
                "emission": emission,
                "quantum_yield": quantum_yield,
                "phase_function": phase_function,
            }
        )

    def add_absorber(self, name, coefficient):
        self._user_components.append(
            {"cls": Absorber, "name": name, "coefficient": coefficient}
        )

    def add_scatterer(self, name, coefficient, phase_function=None):
        self._user_components.append(
            {
                "cls": Scatterer,
                "name": name,
                "coefficient": coefficient,
                "phase_function": phase_function,
            }
        )

    def add_light(
        self,
        name,
        location,  # node location in parent
        rotation=None,  # node rotation in parent frame
        direction=None,  # direction delegate callable
        wavelength=None,  # wavelength delegate callable
        position=None,  # position delegate callable
    ):
        self._user_lights.append(
            {
                "name": name,
                "location": location,
                "rotation": rotation,
                "direction": direction,
                "wavelength": wavelength,
                "position": position,
            }
        )

    def add_solar_cell(self, facets):
        if not isinstance(facets, (list, tuple, set)):
            raise ValueError("Facets should be a set. e.g. `{'left', 'right'}`")
        facets = set(facets)
        allowed = {"left", "near", "far", "right"}
        if not facets.issubset(allowed):
            raise ValueError("Solar cell have allowed surfaces", allowed)

        self._solar_cell_surfaces = facets.union(self._solar_cell_surfaces)

    def add_back_surface_mirror(self):
        self._back_surface_mirror_info = {"want_back_surface_mirror": True}

    def add_air_gap_mirror(self, lambertian=False):
        self._air_gap_mirror_info = {
            "want_air_gap_mirror": True,
            "lambertian": lambertian,
        }

    # Simulate

    def show(
        self,
        wireframe=True,
        baubles=True,
        bauble_radius=None,
        world_segment="short",
        short_length=None,
        open_browser=False,
    ):

        if bauble_radius is None:
            bauble_radius = np.min(self.size) * 0.05

        if short_length is None:
            short_length = np.min(self.size) * 0.1

        self._add_history_kwargs = {
            "bauble_radius": bauble_radius,
            "baubles": baubles,
            "world_segment": world_segment,
            "short_length": short_length,
        }

        if self._scene is None:
            self._make_scene()

        self._renderer = MeshcatRenderer(
            open_browser=open_browser,
            transparency=False,
            opacity=0.5,
            wireframe=wireframe,
            max_histories=50,
        )
        self._renderer.render(self._scene)
        time.sleep(1.0)
        return self._renderer

    def simulate(self, n, progress=None, emit_method="kT"):
        if self._scene is None:
            self._make_scene()
        scene = self._scene

        # Simulate can be called multiple time to append rays to the store
        if self._store is None:
            store = {"entrance_rays": [], "exit_rays": []}

        vis = self._renderer
        count = 0
        for ray in scene.emit(n):
            history = photon_tracer.follow(scene, ray, emit_method=emit_method)
            rays, events = zip(*history)
            store["entrance_rays"].append((rays[1], events[1]))
            if events[-1] in (Event.ABSORB, Event.KILL):
                # final event is a lost store path information at final event
                store["exit_rays"].append((rays[-1], events[-1]))
            elif events[-1] == Event.EXIT:
                # final event hits the world node. Store path information at
                # penultimate location
                store["exit_rays"].append((rays[-2], events[-2]))

            # Update visualiser
            if vis:
                vis.add_history(history, **self._add_history_kwargs)

            # Progress callback
            if progress:
                count += 1
                progress(count)

        self._store = store
        print("Tracing finished.")
        print("Preparing results.")
        df = self._make_dataframe()
        df = self.expand_coords(df, "direction")
        df = self.expand_coords(df, "position")
        df = self.label_facets(df, *self.size)
        self._df = df

    def _make_dataframe(self):

        # to-do: Only need to process additional rays not whole dataframe! Optimise!
        df = pd.DataFrame()

        # Rays entering the scene
        for ray, event in self._store["entrance_rays"]:
            rep = asdict(ray)
            rep["kind"] = "entrance"
            rep["event"] = event.name.lower()
            df = df.append(rep, ignore_index=True)

        # Rays exiting the scene
        for ray, event in self._store["exit_rays"]:
            rep = asdict(ray)
            rep["kind"] = "exit"
            rep["event"] = event.name.lower()
            df = df.append(rep, ignore_index=True)

        self._df = df
        return df

    def expand_coords(self, df, column):
        """ Returns a dataframe with coordinate column expanded into components.
    
            Parameters
            ----------
            df : pandas.DataFrame
                The dataframe
            column : str
                The column label
        
            Returns
            -------
            df : pandas.DataFrame
                The dataframe with the column expanded.
        
            Example
            -------
            Given the dataframe::
        
                df = pd.DataFrame({'position': [(1,2,3)]})
        
            the function will return a new dataframe::
        
                edf = expand_coords(df, 'position')
                edf == pd.DataFrame({'position_x': [1], 'position_y': [2], 'position_z': [3]})
        
        """
        coords = np.stack(df[column].values)
        df["{}_x".format(column)] = coords[:, 0]
        df["{}_y".format(column)] = coords[:, 1]
        df["{}_z".format(column)] = coords[:, 2]
        df.drop(columns=column, inplace=True)
        return df

    def label_facets(self, df, length, width, height):
        """ Label rows with facet names for a box LSC.
    
            Notes
            -----
            This function only works if the coordinates in the dataframe
            are in the local frame of the box. If the coordinates are in the
            world frame then this will still work provided the box is axis
            aligned with the world node and centred at the origin.
        """
        xmin, xmax = -0.5 * length, 0.5 * length
        ymin, ymax = -0.5 * width, 0.5 * width
        zmin, zmax = -0.5 * height, 0.5 * height
        df.loc[(np.isclose(df["position_x"], xmin, atol=EPS_ZERO)), "facet"] = "left"
        df.loc[(np.isclose(df["position_x"], xmax, atol=EPS_ZERO)), "facet"] = "right"
        df.loc[(np.isclose(df["position_y"], ymin, atol=EPS_ZERO)), "facet"] = "far"
        df.loc[(np.isclose(df["position_y"], ymax, atol=EPS_ZERO)), "facet"] = "near"
        df.loc[(np.isclose(df["position_z"], zmin, atol=EPS_ZERO)), "facet"] = "bottom"
        df.loc[(np.isclose(df["position_z"], zmax, atol=EPS_ZERO)), "facet"] = "top"
        return df

    def _make_counts(self, df):

        if self._counts is not None:
            return self._counts

        components = self._scene.component_nodes
        lights = self._scene.light_nodes
        all_components = {component.name for component in components}
        all_lights = {light.name for light in lights}

        # Count solar photons exiting
        solar_out = dict()
        for facet in {"left", "right", "near", "far", "top", "bottom"}:
            solar_out[facet] = self.spectrum(
                facets={facet}, source=all_lights, kind="last"
            ).shape[0]

        # Count solar photons entering
        solar_in = dict()
        for facet in {"left", "right", "near", "far", "top", "bottom"}:
            solar_in[facet] = self.spectrum(
                facets={facet}, source=all_lights, kind="first"
            ).shape[0]

        # Count luminescent photons exiting
        lum_out = dict()
        for facet in {"left", "right", "near", "far", "top", "bottom"}:
            lum_out[facet] = self.spectrum(
                facets={facet}, source=all_components, kind="last"
            ).shape[0]

        # Count luminescent photons entering
        lum_in = dict()
        for facet in {"left", "right", "near", "far", "top", "bottom"}:
            lum_in[facet] = self.spectrum(
                facets={facet}, source=all_components, kind="first"
            ).shape[0]

        self._counts = counts = pd.DataFrame(
            {
                "Solar In": pd.Series(solar_in),
                "Solar Out": pd.Series(solar_out),
                "Luminescent Out": pd.Series(lum_out),
                "Luminescent In": pd.Series(lum_in),
            },
            index=["left", "right", "near", "far", "top", "bottom"],
        )
        return counts

    def spectrum(self, facets=set(), kind="last", source="all", events=None):
        if self._df is None:
            raise ValueError("Run a simulation before calling this method.")

        df = self._df

        if kind is not None:
            if not kind in {"first", "last"}:
                raise ValueError("Direction must be either `'first'` or `'last'.`")

        if kind is None:
            want_kind = True  # Opt-out
        else:
            if kind == "first":
                want_kind = df["kind"] == "entrance"
            else:
                want_kind = df["kind"] == "exit"

        all_sources = self.component_names() | self.light_names()
        if source == "all":
            want_sources = all_sources
        else:
            if isinstance(source, str):
                source = {source}
            if not set(source).issubset(all_sources):
                unknown_source_set = set(source).difference(all_sources)
                raise ValueError("Unknown source requested.", unknown_source_set)

        if source == "all":
            want_source = df["source"].isin(all_sources)
        else:
            want_source = df["source"].isin(set(source))

        if isinstance(facets, (list, tuple, set)):
            if len(facets) > 0:
                want_facets = df["facet"].isin(set(facets))
            else:
                want_facets = True  # don't filter by facets
        else:
            raise ValueError(
                "`'facets'` should be a set `{'left', 'right'}`", {"got": facets}
            )

        if events is None:
            want_events = True  # Don't filter by events
        else:
            all_events = {e.name.lower() for e in Event}
            if isinstance(events, (list, tuple, set)):
                events = set(events)
                if events.issubset(all_events):
                    want_events = df["event"].isin(events)
                else:
                    raise ValueError(
                        "Contained some unknown events",
                        {"got": events, "expected": all_events},
                    )
            else:
                raise ValueError(
                    "Events must be set of event strings", {"allowed": all_events}
                )

        return df.loc[want_kind & want_source & want_facets & want_events]["wavelength"]

    def counts(self):
        df = self._df
        if df is None:
            df = self._make_dataframe()
            df = self.expand_coords(df, "direction")
            df = self.expand_coords(df, "position")
            df = self.label_facets(df, *self.size)

        counts = self._make_counts(df)
        return counts

    def summary(self):
        counts = self._make_counts(self._df)
        all_facets = {"left", "right", "near", "far", "top", "bottom"}

        lum_collected = 0
        for facet in self._solar_cell_surfaces:
            lum_collected += counts["Luminescent Out"][facet]

        lum_escaped = 0
        for facet in all_facets - self._solar_cell_surfaces:
            lum_escaped += counts["Luminescent Out"][facet]

        incident = 0
        for facet in all_facets:
            incident += counts["Solar In"][facet]

        lost = self.spectrum(source="all", events={"absorb"}, kind="last").shape[0]

        optical_efficiency = lum_collected / incident
        waveguide_efficiency = lum_collected / (lum_collected + lum_escaped)

        (l, w, d) = self.size
        a1 = w * l
        a2 = 2 * l * d + 2 * w * d
        Cg = a1 / a2
        n = self.n1
        s = pd.Series(
            {
                "Optical Efficiency": optical_efficiency,
                "Waveguide Efficiency": waveguide_efficiency,
                "Waveguide Efficiency (Thermodynamic Prediction)": (
                    n ** 2 / (Cg + n ** 2)
                ),
                "Non-radiative Loss (fraction):": lost / incident,
                "Incident": incident,
                "Geometric Concentration": Cg,
                "Refractive Index": n,
                "Cell Surfaces": self._solar_cell_surfaces,
                "Components": self.component_names(),
                "Lights": self.light_names(),
            }
        )
        return s

    def report(self):
        print()
        print("Simulation Report")
        print("-----------------")
        print()
        print("Surface Counts:")
        print(self.counts())
        print()
        print("Summary:")
        print(self.summary())

The main changes made to the LSC.py script were in the LSC() class. I removed the size parameter in def_init_, updated the world node as it involved a multiple of the (l, w, d) parameters which no longer exist, and updated the geometry of the LSC node to trimesh.creation.extrude_polygon instead of Box().

lsc = LSC((polygon, 1.0))  # size in cm
lsc.show()  # open visualiser
lsc._renderer.vis.jupyter_cell()

The errors appear when I run the above script to run the LSC report. The errors mention things like the self.size parameter in the bauble_radius code, which I assume means it is attached to the size of the LSC box which no longer exists. I've tried tinkering with this but my python knowledge is very limited so I'm not sure where to begin.

@danieljfarrell
Copy link
Owner

Hi Tommy,

I think this an error from the MeshcatRenderer. Baubles are the little spheres that appear on the ray paths when it hits something, they look like Christmas tree decorations hence the (bad) name.

Can you put the actual traceback here?

But I think the problem is that you are, perhaps, giving pvtrace the trimesh mesh object. You will need to wrap the actual trimesh polygon object with using pvtrace's Mesh class (https://github.com/danieljfarrell/pvtrace/blob/master/pvtrace/geometry/mesh.py). This implement the Geometry protocol so that pvtrace and query the geometry object.

So rather than

LSC((polygon, 1.0))

try,

LSC((Mesh(polygon), 1.0))

@tommyflynn13
Copy link
Author

No problem at all, I'll past the traceback here below:

AttributeError                            Traceback (most recent call last)
Input In [7], in <cell line: 2>()
      1 lsc = LSC((polygon, 1.0))  # size in cm
----> 2 lsc.show()  # open visualiser
      3 lsc._renderer.vis.jupyter_cell()

Input In [6], in LSC.show(self, wireframe, baubles, bauble_radius, world_segment, short_length, open_browser)
    299 def show(
    300     self,
    301     wireframe=True,
   (...)
    306     open_browser=False,
    307 ):
    309     if bauble_radius is None:
--> 310         bauble_radius = np.min(self.size) * 0.05
    312     if short_length is None:
    313         short_length = np.min(self.size) * 0.1

AttributeError: 'LSC' object has no attribute 'size'



@danieljfarrell
Copy link
Owner

And if you change line 1 to be

LSC((Mesh(polygon), 1.0))

@tommyflynn13
Copy link
Author

tommyflynn13 commented Jul 19, 2022

Then it gives this:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 lsc = LSC((Mesh(polygon), 1.0))  # size in cm
      2 lsc.show()  # open visualiser
      3 lsc._renderer.vis.jupyter_cell()

File ~\miniconda3\envs\pvtrace-env\lib\site-packages\pvtrace\geometry\mesh.py:17, in Mesh.__init__(self, trimesh, material)
     15 def __init__(self, trimesh, material=None):
     16     super(Mesh, self).__init__()
---> 17     trimesh.vertices -= trimesh.center_mass
     18     self.trimesh = trimesh
     19     self._material = material

AttributeError: 'Polygon' object has no attribute 'vertices'

Edit: I also tried this input:

lsc = LSC(Mesh(trimesh=trimesh.creation.extrude_polygon(polygon, 1.0)))  # size in cm
lsc.show()  # open visualiser
lsc._renderer.vis.jupyter_cell()

Which gave the traceback:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [11], in <cell line: 2>()
      1 lsc = LSC(Mesh(trimesh=trimesh.creation.extrude_polygon(polygon, 1.0)))  # size in cm
----> 2 lsc.show()  # open visualiser
      3 lsc._renderer.vis.jupyter_cell()

Input In [6], in LSC.show(self, wireframe, baubles, bauble_radius, world_segment, short_length, open_browser)
    299 def show(
    300     self,
    301     wireframe=True,
   (...)
    306     open_browser=False,
    307 ):
    309     if bauble_radius is None:
--> 310         bauble_radius = np.min(self.size) * 0.05
    312     if short_length is None:
    313         short_length = np.min(self.size) * 0.1

AttributeError: 'LSC' object has no attribute 'size'

@danieljfarrell
Copy link
Owner

What's the type of the polygon?

type(polygon)

@danieljfarrell
Copy link
Owner

It should be a trimesh.Trimesh or subclass.

@tommyflynn13
Copy link
Author

It says shapely.geometry.polygon.Polygon

@danieljfarrell
Copy link
Owner

🤔 try and convert that to a trimesh.Trimesh object and it "should just work".

@tommyflynn13
Copy link
Author

No bother, I'll give that a go and let you know how I get on! Should I keep using the LSC.py I updated and removed the 'size' parameters, or would I be best to revert back to your original LSC.py script?

@tommyflynn13
Copy link
Author

Hi Daniel,

Just to give you an update I tried a new approach and it seems to be coming together! Updated the LSC node script to be

 # Create LSC node
        polygon = Polygon([[0, 0], [(w/2), l], [w, 0]])
        lsc = Node(
            name="LSC",
            geometry=Mesh(
                trimesh=trimesh.creation.extrude_polygon(polygon, d),
                material=Material(
                    refractive_index=self.n1,
                    components=components,
                    surface=Surface(delegate=OptionalMirrorAndSolarCell(self)),
                ),
            ),
            parent=world,
        )

I've pasted what appeared in the visualiser after running the new LSC.py script

Screenshot 2022-07-19 203303
.

@danieljfarrell
Copy link
Owner

Oh looks like it's doing something. Think you solved the initial problem. Nice!

@tommyflynn13
Copy link
Author

After updating the air gap to be triangular as well, the visualiser worked however the report gave a new traceback error:

GeometryError                             Traceback (most recent call last)
Input In [20], in <cell line: 754>()
    751     lsc.report()
    754 if __name__ == "__main__":
--> 755     test2()

Input In [20], in test2()
    748 lsc.show()
    749 throw = 300
--> 750 lsc.simulate(throw, emit_method="redshift")
    751 lsc.report()

Input In [20], in LSC.simulate(self, n, progress, emit_method)
    350 count = 0
    351 for ray in scene.emit(n):
--> 352     history = photon_tracer.follow(scene, ray, emit_method=emit_method)
    353     rays, events = zip(*history)
    354     store["entrance_rays"].append((rays[1], events[1]))

File ~\miniconda3\envs\pvtrace-env\lib\site-packages\pvtrace\algorithm\photon_tracer.py:194, in follow(scene, ray, maxsteps, maxpathlength, emit_method)
    192 ray = ray.representation(scene.root, hit)
    193 if surface.is_reflected(ray, hit.geometry, container, adjacent):
--> 194     ray = surface.reflect(ray, hit.geometry, container, adjacent)
    195     ray = ray.representation(hit, scene.root)
    196     history.append((ray, Event.REFLECT))

File ~\miniconda3\envs\pvtrace-env\lib\site-packages\pvtrace\material\surface.py:245, in Surface.reflect(self, ray, geometry, container, adjacent)
    242 def reflect(self, ray, geometry, container, adjacent):
    243     """ Returns ray which is reflected from the interface.
    244     """
--> 245     direction = self.delegate.reflected_direction(
    246         self, ray, geometry, container, adjacent
    247     )
    248     if not isinstance(direction, tuple):
    249         raise ValueError(
    250             "Delegate method `reflected_direction` should return a tuple."
    251         )

Input In [20], in AirGapMirror.reflected_direction(self, surface, ray, geometry, container, adjacent)
     74 def reflected_direction(self, surface, ray, geometry, container, adjacent):
     75     if not self.lsc._air_gap_mirror_info["lambertian"]:
     76         # Specular reflection
---> 77         return super(AirGapMirror, self).transmitted_direction(
     78             surface, ray, geometry, container, adjacent
     79         )
     80     else:
     81         normal = geometry.normal(ray.position)

File ~\miniconda3\envs\pvtrace-env\lib\site-packages\pvtrace\material\surface.py:173, in FresnelSurfaceDelegate.transmitted_direction(self, surface, ray, geometry, container, adjacent)
    171 n2 = adjacent.geometry.material.refractive_index
    172 # Be tolerance with definition of surface normal
--> 173 normal = geometry.normal(ray.position)
    174 if np.dot(normal, ray.direction) < 0.0:
    175     normal = flip(normal)

File ~\miniconda3\envs\pvtrace-env\lib\site-packages\pvtrace\geometry\mesh.py:75, in Mesh.normal(self, surface_point)
     71     raise GeometryError(
     72         "Mesh must have a single closest point to calculate normal."
     73     )
     74 if not np.any(np.absolute(distances) < EPS_ZERO):
---> 75     raise GeometryError(
     76         "Point is not on surface.",
     77         {
     78             "point": surface_point,
     79             "geometry": self,
     80             "distances": distances,
     81             "threshold": EPS_ZERO,
     82         },
     83     )
     84 normal = tuple(mesh.face_normals[triangle_id[0]])
     85 return normal

GeometryError: ('Point is not on surface.', {'point': (-0.48078641173750625, 0.635805804066633, -0.006357771365071763), 'geometry': <pvtrace.geometry.mesh.Mesh object at 0x000001581DF88BE0>, 'distances': array([6.07280425e-13]), 'threshold': 2.220446049250313e-13})

@danieljfarrell
Copy link
Owner

This seems to be a common issue with Mesh geometry. It's related to the precision of the intersection points.

Cab you define your scene in cm or mm rather than meters? This will make everything bigger so thresholding needs less precision.

@tommyflynn13
Copy link
Author

Everything seems to be in working order now, I've attached an image sjowing a quick simulation with 100 rays. My (hopefully!) last question, I was just wondering would you have any idea how I would alter the surfaces in the simulation report? Currently they seem to be (Left, right, top, bottom) as it would be for a box geometry, I imagine I should change the coordinates to have just (left, right, bottom) for mine, with the left and right both being sloped surfaces rather than straight as they would have been in the box.

Simulation Report
-----------------

Surface Counts:
        Solar In  Solar Out  Luminescent Out  Luminescent In
left           0          0                0               0
right          0          0                0               0
near           0          0                0               0
far            0          0                0               0
top          100          3               16               0
bottom         0          1               14               0

Summary:
Optical Efficiency                                                             0.0
Waveguide Efficiency                                                           0.0
Waveguide Efficiency (Thermodynamic Prediction)                           0.473684
Non-radiative Loss (fraction):                                                0.28
Incident                                                                       100
Geometric Concentration                                                        2.5
Refractive Index                                                               1.5
Cell Surfaces                                                                   {}
Components                                         {Lumogen F Red 305, Background}
Lights                                                                     {Light}
dtype: object

Thanks a million for all your help this far also!

@danieljfarrell
Copy link
Owner

Take a look at label_facets function in LSC.py.

This assigns labels to the points sitting on the surface. For example, if x-coordinate of the of the point is on the box plane with x=xmin we definite that as being the left surface.

Remember you are comparing floats so need to have a tolerance.

You will have to come up with a scheme for doing this labelling with your geometry.

It will be the same for the edge surfaces and top and bottom, but they hypotenuse will need more of a complicated comparison to evaluate.

@tommyflynn13
Copy link
Author

I've taken a look at that piece of code, and updated it to the following. Not sure if it makes sense as I don't recognise the formula or attributes. just been playing around with it until it ran successfully.

 def label_facets(self, df, length, width, height):
        """ Label rows with facet names for a box LSC.
    
            Notes
            -----
            This function only works if the coordinates in the dataframe
            are in the local frame of the box. If the coordinates are in the
            world frame then this will still work provided the box is axis
            aligned with the world node and centred at the origin.
        """
        xmin, xmax = -0.5 * width, 0.5 * width
        ymin, ymax = -0.3334 * length, 0.6666 * length
        zmin, zmax = -0.5 * height, 0.5 * height
        df.loc[((df["position_x"] < (width/2))), "facet"] = "left"
        df.loc[((df["position_x"] > (width/2))), "facet"] = "right"
        df.loc[(np.isclose(df["position_y"], ymin, atol=EPS_ZERO)), "facet"] = "near"
        df.loc[(np.isclose(df["position_z"], zmin, atol=EPS_ZERO)), "facet"] = "bottom"
        df.loc[(np.isclose(df["position_z"], zmax, atol=EPS_ZERO)), "facet"] = "top"
        return df

The output was as follows:

Simulation Report
-----------------

Surface Counts:
        Solar In  Solar Out  Luminescent Out  Luminescent In
left         127        163               34               0
right         89         89                0               0
near           0          0                0               0
top           84          2               12               0
bottom         0          0                0               0

Seems to be giving "Solar In" values for the left and right surfaces, so something must be incorrect in the way I currently have it written for the "left" and "right" facet labels.

@danieljfarrell
Copy link
Owner

danieljfarrell commented Jul 23, 2022

I'd suggest continuing with your notebook approach, as this it probably a bit too complicated for now.

Look at the the Quick Start notebook and others in the examples folder. They show how to throw a ray into the scene and record where it escapes.

It too early to make these neat for your geometry. Just get it to work first and then worry about making it simple to use.

The LSC.py module was actually the last thing I coded!

Assigning positions to facets it actually quite tricky if just relying on the coordinates. I think pvtrace should probably make this easier by retuning an ID.

There is a refactor in the works to improve this, but I've been a bit busy the last few years 🥴

#38

@tommyflynn13
Copy link
Author

tommyflynn13 commented Jul 23, 2022

You might be right! Just giving it a final shot as I feel I could be close enough. At the minute, I'm attmepting to define the equation of the line (y1 and y2) on the left and right hand sides, and defining "left" and "right" as being close to that line.

def label_facets(self, df, length, width, height):
     """ Label rows with facet names for a box LSC.
 
         Notes
         -----
         This function only works if the coordinates in the dataframe
         are in the local frame of the box. If the coordinates are in the
         world frame then this will still work provided the box is axis
         aligned with the world node and centred at the origin.
     """
     y1 = ((width/length) * df["position_x"]) + ((2/3) * length)
     y2 = -((width/length) * df["position_x"]) + ((2/3) * length)
     xmin, xmax = -0.5 * width, 0.5 * width
     ymin, ymax = -(1/3) * length, (2/3) * length
     zmin, zmax = -0.5 * height, 0.5 * height
     df.loc[(np.isclose(df["position_y"], y1, atol=EPS_ZERO)), "facet"] = "left"
     df.loc[(np.isclose(df["position_y"], y2, atol=EPS_ZERO)), "facet"] = "right"
     df.loc[(np.isclose(df["position_y"], ymin, atol=EPS_ZERO)), "facet"] = "near"
     df.loc[(np.isclose(df["position_y"], ymax, atol=EPS_ZERO)), "facet"] = "far"
     df.loc[(np.isclose(df["position_z"], zmin, atol=EPS_ZERO)), "facet"] = "bottom"
     df.loc[(np.isclose(df["position_z"], zmax, atol=EPS_ZERO)), "facet"] = "top"
     return df

It doesn't seem to be working as the results aren't adding up, but do you think I could be on the right track with this? Would I be better off trying to define an array of points along the line (the equation of the line is y = 1.4x + 4.666 for the dimensions of my LSC) and using that, and if so how would I do that?

I've attached an image from Desmos of the equations of the three lines which make up my triangular geometry.

image

Thanks again,
Tommy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants