/************************************************************************************************ * User-defined function for arbitrary polygonal holes * * Called from inside a groundplane definition as: * + hole user1 (x1,y1,x2,y2,x3,y3,...) * * Defines a polygon by a set of three or more tuples giving (x,y) coordinates for the polygon * vertices. Removes all nodes from the groundplane that are inside or on the edge of the polygon. * * The function uses a standard ray-casting algorithm (http://en.wikipedia.org/wiki/Point_in_polygon) * to determine if each groundplane node is inside or outside the polygon: for each node, extend a ray * in any direction (here in the -x direction) and count the number of polygon edges it crosses. If * the count is odd, the point is inside the polygon, if the count is even, the point is outside the * polygon. This algorithm handles self-crossing polygons, like hourglasses, or pentagrams, but is not * recommended for them as it creates groundplane islands, disconnected or linked by a single segment. * * Note: This user function assumes ground plane and hole are both at constant z. Feel free to * generalize it if you want to do the work. * *************************************************************************************************/ hole_user1(holep, gp, relx, rely, relz, units) GROUNDPLANE *gp; HoleList *holep; double relx,rely,relz,units; { int i, j, e; NODES *origin; int gn1, gn2; double gx0, gy0; double gd1, gd2; double ux1, uy1, ux2, uy2; int nedges; double ex0, ey0, ex1, ey1; int inPolygon; if ( ((holep->numvals)%2 != 0) || (holep->numvals < 6) ) // Check for correct number of inputs (expect even number greater than 5) hole_error("Need an even number of inputs for at least three (x,y) pairs.","(x,y)",holep); nedges = holep->numvals / 2; // Number of edges of hole polygon // Extract useful variables from the groundplane structure origin = gp->pnodes[0][0]; gx0 = origin->x / units; // Origin corner of groundplane (x-coord) gy0 = origin->y / units; // Origin corner of groundplane (y-coord) gn1 = gp->num_nodes1; // Number of nodes along first side of groundplane gn2 = gp->num_nodes2; // Number of nodes along second side of groundplane gd1 = gp->d1 / units; // Groundplane segment length along first side gd2 = gp->d2 / units; // Groundplane segment length along second side ux1 = gp->ux1; // Unit vector between segments along first side (x-coord) uy1 = gp->uy1; // Unit vector between segments along first side (y-coord) ux2 = gp->ux2; // Unit vector between segments along second side (x-coord) uy2 = gp->uy2; // Unit vector between segments along second side (y-coord) for(j = 0; j < gn2; j++) { // For each ground plane node for(i = 0; i < gn1; i++) { inPolygon = FALSE; for(e = 0; e < nedges; e++) { // For each edge of the polygon // Get the endpoints of the edge in the groundplane coordinate system (ex0,ey0) and (ex1,ey1) // (Possible confusion: I use a different convention from that of the groundplane struct which has axis 1 and axis 2 // and use x and y in the groundplane coordinate system for later thinking convenience.) ex0 = (holep->vals[0+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[1+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1; ey0 = (holep->vals[0+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[1+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2; if(e+1 == nedges) { ex1 = (holep->vals[0] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[1] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1; ey1 = (holep->vals[0] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[1] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2; } else { ex1 = (holep->vals[2+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux1 + (holep->vals[3+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy1; ey1 = (holep->vals[2+2*e] - (gx0 + j*gd2*ux2 + i*gd1*ux1))*ux2 + (holep->vals[3+2*e] - (gy0 + j*gd2*uy2 + i*gd1*uy1))*uy2; } // Toggle inPolygon if the leftward ray from the groundplane point intersects the polygon edge if(ey0*ey1 < 0) { if((ex0 - ey0*(ex1-ex0)/(ey1-ey0)) <= 0) { // If the polygon edge crosses the horizontal line and does it to the left of the node, toggle inside/outside status inPolygon = !inPolygon; } } else if( ((ey0 == 0) && (ey1 < 0) && (ex0 <= 0)) || ((ey1 == 0) && (ey0 < 0) && (ex1 <= 0))) { // If one endpoint of a polygon edge touches the horizontal to the left of the node and the other node is under the line, toggle inside/outside status inPolygon = !inPolygon; } } if(inPolygon) { // Delete node if inside polygon delete_node(gp->pnodes[i][j]); } } } }