Adding sidewall angle parameter to polygon object

[smartalecH]

Several planar lithography applications expect a non-normal (i.e. 90 degree) sidewall angle. This can significantly affect the performance of integrated waveguide structures. Fortunately, many foundries can predict the sidewall angle within their process to a relatively high order of accuracy (compared to other metrics). It would be nice to add this optional sidewall angle parameter to the existing polygon object within libctl.

The implementation wouldn't be too difficult. Right now the polygon class projects a 2D polygon into 3D space. We could generalize this step to account for a narrowing (or widening) in the third dimension.

We should first try to resolve the point in polygon methods, however, since the current bugs would also drastically affect this new feature.

[stevengj]

As long as the top polygon is just an affine transformation of the bottom polygon, I think you can just modify the m_c2p matrix.

This should be sufficient to support a single sidewall angle per prism object.

[danielwboyce]

The transformation between the bottom polygon and the top one isn't the one we're interested in here; m_c2p is the matrix of interest but that transforms a three-dimensional point in Cartesian coordinates to the coordinates of the prism. This was previously an affine transformation in all cases because the transformation was only a rotation of the coordinate system. However in this case, we're attempting the transformation of a prism that looks like this in its internal coordinate system (with normal sidewalls)

to a prism that looks like this in the coordinate system of the simulation:

with non-normal sidewalls (note that in the figure the sidewall angle is much greater than what it is likely to be in practice). This transformation, of course, is not affine but rather is a projective (or perpsective) transformation (that is, the transformation maps lines to lines but doesn't necessary preserve parallelism). What this means for us is that we'll have to figure out some way for the scaling factors in m_c2p on the x- and y-coordinates to be dependent on z-coordinate. We may wish to add another helper function that finds a matrix that transforms m_c2p/m_p2c' at a given point and then matrix multiply m_c2p/m_p2c' and the new matrix in prism_coordinate_p2c, prism_vector_p2c, prism_coordinate_c2p, and prism_vector_c2p.

Most of our other methods relating to prisms, like point_in_prism, shouldn't need modification with this scheme. The main exception that comes to mind is geom_object_volume, which calculates the volume of the prism by finding area of the base polygon and multiplying by the prism's height, and so won't work. Perspective transformations are a common problem in animation/computer graphics, though, so I'll also take a look to see if there are any simple solutions to this problem. If not a derivation of volume wouldn't be too difficult even if just the height, base area, and sidewall volume are known. We'll also need to test normal_to_prism, intersect_line_with_prism, and the like to see how they respond to non-affine transformations between the prism coordinate system and the main coordinate system.

[stevengj]

It’s worse than that. , I thought about it last night and it’s not a simple scaling that we want. We’ll need a more general shape with different polygons for the top and bottom, linearly interpolating in between

[danielwboyce]

I attempted to add in a matrix that could transform the m_c2p and m_p2c matrices for a given coordinate. The math behind this was that the matrix operating m_c2p would be defined

with the scaling coordinates defined

and

and the matrix operating on m_p2c being the inverse of the matrix operating on m_c2p. This was implemented

matrix3x3 prism_projective_transformation_for_p2c(prism *prsm, vector3 pp) {
    matrix3x3 c2p;
    double cx;
    double cy;
    double theta = (K_PI / 2) - prsm->sidewall_angle;
    if (prsm->sidewall_angle == 0) {
        cx = 1;
        cy = 1;

    }
    if (pp.x == 0) {
        cx = 0;
    }
    else {
        cx = 1 - pp.z / (pp.x * tan(theta));
    }
    vector3 c0vector = {cx, 0, 0};
    c2p.c0 = c0vector;
    vector3 c1vector = {0, cy, 0};
    c2p.c1 = c1vector;
    vector3 c2vector = {0, 0, 1};
    c2p.c2 = c2vector;
    return c2p;
}

matrix3x3 prism_projective_transformation_for_c2p(prism *prsm, vector3 pc) {
    matrix3x3 p2c = matrix3x3_inverse(prism_projective_transformation_for_c2p(prsm, pc));
    return p2c;
}

I'm dubious that any of this is correct, however, because

1. I might have the matrix operating on m_c2p and the matrix operating on m_p2c backwards, and
2. this only has a chance of working as a transformation on a single point, but couldn't possibly work on a vector transformation because basic vectors don't carry a start or end point with them (and they shouldn't need to).

I've checked out some books on projective geometry and transformations and am still trying to dig through them, but what I'm coming up with so far seems to require implementing homogeneous coordinates, which isn't something I'd love to do.

Any thoughts, hints, or other on how to implement this transformation would be greatly appreciated.

[danielwboyce]

Another wrinkle that came to my attention is the question of how to deal with hollow polygons as they are extruded (a ring resonator, for example). Most scaling transformations that I've come across, and I think projective transformations, as well, would cause the inner radius of the top polygon to be smaller than the inner radius of the bottom polygon.

What we want, though, is for the inner radius of the top polygon to be slightly larger than the inner radius of the bottom polygon, with the outer radius of the top polygon remaining slightly smaller than the outer radius of the bottom polygon (which is the result we would naturally receive with basic scaling transformations and projective transformations).

I think this needed result might throw us into a new kind of transformation other than a basic scaling or projective transformation. I'm going to do some literature searching to see what the term of this would be.

[stevengj]

Basically I think we need one polygon for the bottom and one polygon for the top and then just interpolate. (Or instead of storing the top polygon we can just store the difference vectors between each bottom–top pair of vertices… that would only save us one subtraction at best during point_in_prism checks, though.)

See also https://scholar.google.com/scholar?hl=en&as_sdt=0%2C22&q=piecewise+linear+interpolation+polygonal+slices&btnG= … I suspect that our problem will be easier because we have the same number of vertices and the same topology between top and bottom

[danielwboyce]

I've got the linear interpolation working for a basic case (extruded square with sidewall angle of 1 degree). We can see here the prism's vertices plotted in gnu plot with a normal sidewall and with the 1 degree sidewall:

[danielwboyce]

I've also already updated point_in_prism to work with any sidewall. Now I'll work updating some prism helper functions in geom.c (normal_to_prism, intersect_line_with_prism, get_prism_bounding_box, geom_object_volume, intersect_line_segment_with_prism, mind_distance_to_prism_roof_or_ceiling) and then I'll work on testing more prisms.

[danielwboyce]

For reference the gnuplot vertices were created with a new unit test in utils/test-prism.c

int test_sidewall_prisms_to_gnuplot() {
  void *m = NULL;

  int num_nodes = 4;
  vector3 nodes[num_nodes];
  nodes[0] = make_vector3(-1.0, -1.0, 0.0);
  nodes[1] = make_vector3(-1.0, 1.0, 0.0);
  nodes[2] = make_vector3(1.0, 1.0, 0.0);
  nodes[3] = make_vector3(1.0, -1.0, 0.0);

  double height = 10;
  vector3 zhat = make_vector3(0, 0, 1);

  double normal_sidewall = 0;
  geometric_object normal_sidewall_geom_object = make_prism(m, nodes, num_nodes, height, zhat, normal_sidewall);
  prism *normal_sidewall_prism = normal_sidewall_geom_object.subclass.prism_data;

  double one_degree_sidewall = 1.0 * 2 * K_PI / 360.0;
  geometric_object one_degree_sidewall_geom_object = make_prism(m, nodes, num_nodes, height, zhat, one_degree_sidewall);
  prism *one_degree_sidewall_prism = one_degree_sidewall_geom_object.subclass.prism_data;

  prism2gnuplot(normal_sidewall_prism, "normal_sidewall_gnu_plot.dat");
  prism2gnuplot(one_degree_sidewall_prism, "one_degree_sidewall_gnu_plot.dat");

  return 0;
}

with a new version of prism2gnuplot

void prism2gnuplot(prism *prsm, char *filename) {
  int num_vertices = prsm->vertices_p.num_items;
  double height = prsm->height;
  vector3_list vertices_bottom;
  vertices_bottom.num_items = num_vertices;
  vertices_bottom.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
  memcpy(vertices_bottom.items, prsm->vertices_p.items, num_vertices * sizeof(vector3));
  vector3_list vertices_top;
  vertices_top.num_items = num_vertices;
  vertices_top.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
  int nv;
  for (nv = 0; nv < num_vertices; nv++) {
    vertices_top.items[nv] = vector3_plus(prsm->vertices_p.items[nv], prsm->top_polygon_diff_vectors_p.items[nv]);
  }

  FILE *f = fopen(filename, "w");
  for (nv = 0; nv < num_vertices; nv++) {
    vector3 vap = vertices_bottom.items[nv];
    vap.z = 0.0;
    vector3 vbp = vertices_top.items[nv];
    vbp.z = height;
    vector3 vcp = vertices_top.items[(nv + 1) % num_vertices];
    vcp.z = height;
    vector3 vdp = vertices_bottom.items[(nv + 1) % num_vertices];
    vdp.z = 0.0;
    vector3 vac = prism_coordinate_p2c(prsm, vap);
    vector3 vbc = prism_coordinate_p2c(prsm, vbp);
    vector3 vcc = prism_coordinate_p2c(prsm, vcp);
    vector3 vdc = prism_coordinate_p2c(prsm, vdp);
    
    fprintf(f, "%e %e %e \n", vac.x, vac.y, vac.z);
    fprintf(f, "%e %e %e \n", vbc.x, vbc.y, vbc.z);
    fprintf(f, "%e %e %e \n", vcc.x, vcc.y, vcc.z);
    fprintf(f, "%e %e %e \n", vdc.x, vdc.y, vdc.z);
    fprintf(f, "%e %e %e \n", vac.x, vac.y, vac.z);
    fprintf(f, "\n\n");
  }
  fclose(f);
}

and the following commands fed to gnuplot:

reset
set size square
set view 66,55,1,1
set offset 1,1,1,1
set style textbox opaque border lc "blue"
unset key

set title "normal sidewall angle"
splot 'normal_sidewall_gnu_plot.dat' using 1:2:3 title column with lines lw 3,  \
     ''          using 1:2:3:(sprintf("(%0.2f, %0.2f, %0.2f)", $1, $2, $3)) with labels center boxed

set title "one degree sidewall angle"
splot 'one_degree_sidewall_gnu_plot.dat' using 1:2:3 title column with lines lw 3,  \
     ''          using 1:2:3:(sprintf("(%0.2f, %0.2f, %0.2f)", $1, $2, $3)) with labels center boxed

[danielwboyce]

As of now all of the math in geom.c has been updated so that it should work with the sidewall feature (in a basic case). I'm a little unsure on how we should find the correct average between the top polygon area and bottom polygon area when finding the prism volume. I think the volume should be the RMS area multiplied by the height, so if that's correct my volume calculation would be correct.

[danielwboyce]

Right now I believe this only works for a basic case, where the base polygon is solid. I'm unsure how to make this work correctly for a case like a ring resonator or a big C, where the outer points taper in and the inner points taper out.

I feel like I've read somewhere that it's considered standard practice for outer vertices when listed sequentially to be listed in a clockwise direction and then for inner vertices to be listed in a counterclockwise direction. If this is true I could utilize this to decide which direction to taper the sidewall in. Do we follow this practice?

[danielwboyce]

After some more thinking about it, I realized that a simple test of clockwise or counter clockwise to determine whether the point will taper in or out won't suffice, because if you get a shape like an tapered, extruded C

not all transforms will point directly towards or directly away from the origin ("in" or "out"), shown here in a top view of the extruded C:

so the issue is more complicated.

[oskooi]

from @stevengj

[stevengj]

In the diagram above, (a,b,c) are vertices in the top polygon and (a',b',c') are vertices in the bottom polygon. The sidewall angle is presumably defined as the slope of the surface perpendicular to each edge, e.g. along the line p–p' shown above.

To get the volume, we need to add up the volumes of all the little polyhedral wedges that make up the difference between the top and bottom polygons, i.e. between the green arrows in the diagram above. In general, the green arrows (the intersections of adjacent faces) will not be parallel, which means that there won't be a simple function like height×base/2 for the volume of the wedge.

In principle, it's just a Small Matter of Algebra™ to work out the the volumes of these wedges (as well as the coordinates of the green arrows etcetera), with no advanced math — just intersections of planes, so just linear algebra basically. But it is a bit tedious and will require care.

[stevengj]

Note that the prism objects only support simply-connected polygons at the moment (no holes). Holes are added by adding another prism object.

[danielwboyce]

To get the volume, we need to add up the volumes of all the little polyhedral wedges that make up the difference between the top and bottom polygons, i.e. between the green arrows in the diagram above. In general, the green arrows (the intersections of adjacent faces) will not be parallel, which means that there won't be a simple function like height×base/2 for the volume of the wedge.

Even if the triangular faces on each end of the wedge are not parallel, I believe the three lines connecting the vertices of the triangles (a->b, a'->b', and the projected line of a'->b' to the z=0 plane, using the points in the above diagram) to form the wedge will always be parallel. This is because the point of the tapering transform is that it yields lines a->b and a'->b' that are parallel to each other, and the lines a'->b' and a'->b' projected to the z=0 plane will be parallel by definition. Therefore we should be able to use the method proved in this blog post to find the volume of the wedge, so long as we can find the cross-sectional area A of the wedge (which will be one of the end triangles projected into some plane perpendicular to a->b, a'->b' and a'->b' projected to the z=0).

Being able to properly implement this is now more a matter of figuring out how to implement the taper transformation properly.

[danielwboyce]

I made the connection that because the tapering transform yields edges that are parallel to one another with a constant normal distance, it's not right to try and calculate the new vertices of the top polygon directly but instead the new edges of the top polygon should be calculated and the intersections of these edges can be found to store as the vertices of the top polygon. I've implemented a rough draft of this algorithm in geom.c, which I'm in the process of debugging.

The rough draft code is here (as a part of init_prism):

  // Calculate difference vertices of the top polygon and vectors between bottom
  // polygon and the top polygon, where:
  //  * the bottom polygon is the one passed in to the the make_prism() function,
  //    stored in vertices_bottom and vertices_bottom_p,
  //  * the top polygon is the top surface (parallel to the bottom polygon) resulting
  //    from the extrusion of the bottom polygon. Whether or not the extrusion tapers
  //    is dependent on the value of sidewall_angle.
  //
  // The top polygon is calculated by first copying the values of vertices_bottom_p into
  // vertices_top_p, except z=prsm->height for all top vertices. If prsm->sidewall_angle
  // is equal to zero, then no further calculations are performed on the top vertices.
  // If not, we know that all EDGES of the the top polygon will be offset so that in the
  // xy plane they are parallel to the edges of the bottom polygon. The offset amount is
  // determined by the sidewall angle and the height of the prism. To perform the
  // calculation, each of the edges of the top polygon (without an offset) are stored in
  // an array of edges (edge is a struct defined if prsm->sidewall_angle!=0 containing
  // the endpoints a1 a2, with a third vector v defined a2-a1). Then the vector normal to
  // v is calculated, and the offset vector. A test is performed to determine in which
  // direction (the direction of +offset or -offset) from the edge we can find points
  // inside the polygon by performing a node_in_or_on_polygon test at
  //          edge.a1 + 0.5*edge.v + 1e-3*offset.
  // This information is used to determine in which direction the offset of the edge is
  // applied, in conjunction with whether prsm->sidewall_angle is positive or negative
  // (if positive, the offset will be applied in towards the points where 
  // node_in_or_on_polygon is true, else the offset will be applied out away from those 
  // points). After the offsets are applied to the edges, the intersections between the 
  // new edges are calculated, which are the new values of vertices_top_p.
  //
  // Some side notes on the difference vectors:
  //   * The value of each of the top polygon vertices can be found
  //             vertices_bottom_p + top_polygon_diff_vectors_p
  //             vertices_bottom   + top_polygon_diff_vectors
  //   * A linearly interpolated value of the polygon vertices between the bottom
  //     polygon and the top can be found
  //             vertices_bottom_p + top_polygon_diff_vectors_scaled_p * z
  if (isnan(prsm->sidewall_angle)) {
      prsm->sidewall_angle = 0.0;
  }
  number theta = (K_PI/2) - abs(prsm->sidewall_angle);

  prsm->vertices_top_p.num_items = num_vertices;
  prsm->vertices_top_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
  CHECK(prsm->vertices_top_p.items, "out of memory");
  memcpy(prsm->vertices_top_p.items, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3));
  for (nv = 0; nv < num_vertices; nv++) {
    prsm->vertices_top_p.items[nv].z = prsm->height;
  }

  if (prsm->sidewall_angle != 0.0) {
    struct edge {
      vector3 a1, a2, v;  // v will be defined as a2 - a1
    } edge;

    edge *top_polygon_edges;
    top_polygon_edges = (edge *)malloc(num_vertices * sizeof(edge));
    number w = prsm->height / tan(theta);

    for (nv = 0; nv < num_vertices; nv++) {
      top_polygon_edges[nv].a1 = prsm->vertices_top_p.items[(nv - 1 == -1 ? num_vertices - 1 : nv - 1)];
      top_polygon_edges[nv].a2 = prsm->vertices_top_p.items[nv];
      top_polygon_edges[nv].v = vector3_minus(top_polygon_edges[nv].a2, top_polygon_edges[nv].a1);

      vector3 normal_vector;
      normal_vector.x = top_polygon_edges[nv].v.y;
      normal_vector.y = -1 * top_polygon_edges[nv].v.x;
      normal_vector.z = 0;
      normal_vector = unit_vector3(normal_vector);
      vector3 offset = vector3_scale(w, normal_vector);

      vector3 midpoint = vector3_plus(top_polygon_edges[nv].a1, vector3_scale(0.5, top_polygon_edges[nv].v));
      boolean midpoint_on_side_of_postive_offset_is_in_polygon = node_in_or_on_polygon(vector3_plus(midpoint, vector3_scale(1e-3, offset)), prsm->vertices_top_p.items, prsm->vertices_top_p.num_items, 1);

      if (midpoint_on_side_of_postive_offset_is_in_polygon) {
        // positive sidewall angles means the prism tapers in towards to the rest of the prism body
        if (prsm->sidewall_angle > 0) {
          top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset);
        }
        // negative sidewall angles means the prism tapers out away from the rest of the prism body
        else {
          top_polygon_edges[nv].a1 = vector3_minus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_minus(top_polygon_edges[nv].a2, offset);
        }
      }
      else {
        // positive sidewall angles means the prism tapers in towards to the rest of the prism body
        if (prsm->sidewall_angle > 0) {
          top_polygon_edges[nv].a1 = vector3_minus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_minus(top_polygon_edges[nv].a2, offset);
        }
        // negative sidewall angles means the prism tapers out away from the rest of the prism body
        else {
          top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset);
        }
      }
    }

    for (nv = 0; nv < num_vertices; nv++) {
      number x1 = top_polygon_edges[nv].a1.x;
      number y1 = top_polygon_edges[nv].a1.y;
      number x2 = top_polygon_edges[nv].a2.x;
      number y2 = top_polygon_edges[nv].a2.y;
      number x3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.x;
      number y3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.y;
      number x4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.x;
      number y4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.y;

      // Intersection point calculated with https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
      number px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
      number py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
      prsm->vertices_top_p.items[nv].x = px;
      prsm->vertices_top_p.items[nv].y = py;
    }
  }

[danielwboyce]

I updated init_prism so that it now properly initializes simply-connected polygons (both concave and convex) with and without nonzero sidewall angles. The relevant section now reads:

  // Calculate difference vertices of the top polygon and vectors between bottom
  // polygon and the top polygon, where:
  //  * the bottom polygon is the one passed in to the the make_prism() function,
  //    stored in vertices_bottom and vertices_bottom_p,
  //  * the top polygon is the top surface (parallel to the bottom polygon) resulting
  //    from the extrusion of the bottom polygon. Whether or not the extrusion tapers
  //    is dependent on the value of sidewall_angle.
  //
  // The top polygon is calculated by first copying the values of vertices_bottom_p into
  // vertices_top_p, except z=prsm->height for all top vertices. If prsm->sidewall_angle
  // is equal to zero, then no further calculations are performed on the top vertices.
  // If not, we know that all EDGES of the the top polygon will be offset so that in the
  // xy plane they are parallel to the edges of the bottom polygon. The offset amount is
  // determined by the sidewall angle and the height of the prism. To perform the
  // calculation, each of the edges of the top polygon (without an offset) are stored in
  // an array of edges (edge is a struct defined if prsm->sidewall_angle!=0 containing
  // the endpoints a1 a2, with a third vector v defined a2-a1). Then the vector normal to
  // v is calculated, and the offset vector. A test is performed to determine in which
  // direction (the direction of +offset or -offset) from the edge we can find points
  // inside the polygon by performing a node_in_or_on_polygon test at a finite distance
  // away from the midpoint of the edge:
  //          edge.a1 + 0.5*edge.v + 1e-3*offset.
  // This information is used to determine in which direction the offset of the edge is
  // applied, in conjunction with whether prsm->sidewall_angle is positive or negative
  // (if positive, the offset will be applied in towards the points where
  // node_in_or_on_polygon is true, else the offset will be applied out away from those
  // points). After the offsets are applied to the edges, the intersections between the
  // new edges are calculated, which are the new values of vertices_top_p.
  //
  // Some side notes on the difference vectors:
  //   * The value of each of the top polygon vertices can be found
  //             vertices_bottom_p + top_polygon_diff_vectors_p
  //             vertices_bottom   + top_polygon_diff_vectors
  //   * A linearly interpolated value of the polygon vertices between the bottom
  //     polygon and the top can be found
  //             vertices_bottom_p + top_polygon_diff_vectors_scaled_p * z
  if (isnan(prsm->sidewall_angle)) {
      prsm->sidewall_angle = 0.0;
  }
  number theta = (K_PI/2) - fabs(prsm->sidewall_angle);
  prsm->vertices_top_p.num_items = num_vertices;
  prsm->vertices_top_p.items = (vector3 *)malloc(num_vertices * sizeof(vector3));
  CHECK(prsm->vertices_top_p.items, "out of memory");
  memcpy(prsm->vertices_top_p.items, prsm->vertices_bottom_p.items, num_vertices * sizeof(vector3));
  for (nv = 0; nv < num_vertices; nv++) {
    prsm->vertices_top_p.items[nv].z = prsm->height;
  }

  if (prsm->sidewall_angle != 0.0) {
    typedef struct {
      vector3 a1, a2, v;  // v will be defined as a2 - a1
    } edge;

    edge *top_polygon_edges;
    top_polygon_edges = (edge *)malloc(num_vertices * sizeof(edge));
    number w = prsm->height / tan(theta);

    for (nv = 0; nv < num_vertices; nv++) {
      top_polygon_edges[nv].a1 = prsm->vertices_top_p.items[(nv - 1 == -1 ? num_vertices - 1 : nv - 1)];
      top_polygon_edges[nv].a2 = prsm->vertices_top_p.items[nv];
      top_polygon_edges[nv].v = vector3_minus(top_polygon_edges[nv].a2, top_polygon_edges[nv].a1);

      vector3 normal_vector;
      normal_vector.x = top_polygon_edges[nv].v.y;
      normal_vector.y = -1 * top_polygon_edges[nv].v.x;
      normal_vector.z = 0;
      normal_vector = unit_vector3(normal_vector);
      vector3 offset = vector3_scale(w, normal_vector);

      vector3 midpoint = vector3_plus(top_polygon_edges[nv].a1, vector3_scale(0.5, top_polygon_edges[nv].v));
      boolean midpoint_on_side_of_postive_offset_is_in_polygon = node_in_or_on_polygon(vector3_plus(midpoint, vector3_scale(1e-3, offset)), prsm->vertices_top_p.items, prsm->vertices_top_p.num_items, 1);

      if (midpoint_on_side_of_postive_offset_is_in_polygon) {
        // positive sidewall angles means the prism tapers in towards to the rest of the prism body
        if (prsm->sidewall_angle > 0) {
          top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset);
        }
        // negative sidewall angles means the prism tapers out away from the rest of the prism body
        else {
          top_polygon_edges[nv].a1 = vector3_minus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_minus(top_polygon_edges[nv].a2, offset);
        }
      }
      else {
        // positive sidewall angles means the prism tapers in towards to the rest of the prism body
        if (prsm->sidewall_angle > 0) {
          top_polygon_edges[nv].a1 = vector3_minus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_minus(top_polygon_edges[nv].a2, offset);
        }
        // negative sidewall angles means the prism tapers out away from the rest of the prism body
        else {
          top_polygon_edges[nv].a1 = vector3_plus(top_polygon_edges[nv].a1, offset);
          top_polygon_edges[nv].a2 = vector3_plus(top_polygon_edges[nv].a2, offset);
        }
      }
    }

    for (nv = 0; nv < num_vertices; nv++) {
      number x1 = top_polygon_edges[nv].a1.x;
      number y1 = top_polygon_edges[nv].a1.y;
      number x2 = top_polygon_edges[nv].a2.x;
      number y2 = top_polygon_edges[nv].a2.y;
      number x3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.x;
      number y3 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a1.y;
      number x4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.x;
      number y4 = top_polygon_edges[(nv + 1 == num_vertices ? 0 : nv + 1)].a2.y;

      // Intersection point calculated with https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
      number px = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
      number py = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
      prsm->vertices_top_p.items[nv].x = px;
      prsm->vertices_top_p.items[nv].y = py;
    }
  }

The rectangle prism and tapered rectangular prism render the same as before, but I also tested with a concave octagonal c shape, shown here:

The prism with the 2.5-degree sidewall angle has exactly the same vertices as the vertices in this FreeCAD file, which was created by drawing the same base vertices and then extruding it with a taper angle of 2.5 degrees. (FreeCAD is an open-source CAD system, with pre-compiled binaries available for Linux, Mac, and Windows.)

[danielwboyce]

The most recent version of danielwboyce/libctl branch:meep_iss_1067 has updated routines that appear to be accurately measuring the volumes of simply-connected convex and concave polygons both with and without nonzero sidewall angles (with the volumes tested against results taken from CAD softwares). That should mean the update to utils/geom.c is completed, pending additional testing (see NanoComp/libctl#52, where the implementation of a much more prism testing suite is discussed).