#include #include #include #include #include #define maxThreads 85 #define expectedVertices 12 #define maxWidth 400 #define maxHeight 300 #define currentMaxDegree 5 #define sizeOfGraph 73 #define maxFaceSize 10 #define maxNumberOfFaces 20 #define RIGID -1 #define separatorASCIICode 32 int totalCorrectGraphs=0; int totalProcessedGraphs = 0; #define rigidAngleInformation nonRigidVertexDegrees,nonRigidEdges,angleArray,rigidEdges,rigidVertexDegrees #define graphInformation nonRigidVertexDegrees,nonRigidEdges,faceArray,numberOfFaces,angleArray,rigidVertexDegrees,rigidEdges double ACOSONETHIRD = 1.23095942; int legitSemiRigidGraphs = 0; int legitRigidGraphs = 0; double square(double x) { //returns x^2 return (x*x); } double vectorDistance(double vector1[3],double vector2[3]) { //returns the square of the distance between vector 1 and vector 2 double totalDistance = 0; //used to store the total distance between v1 and v2 int i; for(i = 0; i < 3; i++) { totalDistance += square(vector1[i] - vector2[i]); //iterate through each i, adding (v1_i - v2_i)^2 each time } return totalDistance; //return the total distance } void splitGraphToVertices(int graph[sizeOfGraph],int lengthOfGraph,int vertexDegrees[expectedVertices],int vertices[expectedVertices][currentMaxDegree]) { int currentEdge = 0; //used to temporarily store the current edge or a 0 if the end of the vertex is specified int currentVertexPosition = 0; //the current degree of the vertex int currentGraphEntry=0; //the current entry in the graph int vertexNumber = 0; //how many vertices have we currently processed for(currentGraphEntry = 1;currentGraphEntry < lengthOfGraph; currentGraphEntry++) { currentEdge = graph[currentGraphEntry]; if(currentEdge != 0) { vertices[vertexNumber][currentVertexPosition] = currentEdge - 1; //plantri outputs vertices starting from vertex 1, but we would prefer to label them from 0 currentVertexPosition++; } else { //we have finished processing this vertex vertexDegrees[vertexNumber] = currentVertexPosition; currentVertexPosition = 0; vertexNumber++; } } if(vertexNumber != expectedVertices) { fprintf(stderr,"A graph was passed to splitGraphToVertices with the incorrect number of vertices. %i vertices were passed and %i vertices were expected",vertexNumber,expectedVertices); } } int findNextVertex(int vertices[expectedVertices][currentMaxDegree],int previousVertex,int nextVertex,int currentVertexDegree) { /* We consider the list of vertices connected to the vertex we are currently at. Once we have found the edge that we just took to reach this vertex (by finding the previous vertex in the aforementioned list) we take the next edge in the list to get the next vertex (clockwise) */ int i = 0; int returnValue = -1; for(i = 0; i < (currentVertexDegree - 1); i++) { if((vertices[nextVertex][i]) == previousVertex) { //We have found the edge that we just took returnValue = vertices[nextVertex][i+1]; //this is the next edge clockwise break; } } if(i == (currentVertexDegree - 1)) { /*the final edge of the current vertex connects to the previous vertex: thus the next edge to traverse clockwise is the first edge listed*/ returnValue = vertices[nextVertex][0]; } return returnValue; } void retrieveFaces(int vertexDegrees[expectedVertices],int vertices[expectedVertices][currentMaxDegree],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize],int numberOfVertices) { /*Input is vertexDegrees (the degrees of each vertex) and vertices (the degrees of each vertex and its neighbours in clockwise order). This procedure modifies the values of faceArray and numberOfFaces numberOfFaces stores the number of faces in each entry e.g. numberOfFaces[2] is the number of triangles found faceArray stores information about types of faces in the graph e.g. faceArray[3][1][2] is the third vertex in the second face from the collection of quadrilateral faces that the graph contains */ int vertexNumber = 0; //which vertex are we currently processing int nextVertex; //the next vertex in the face int previousVertex; //the previous vertex in the face int faceSides = 0;//stores the number of sides of the face int edgeProcessed; //the current edge that we are processing int currentFace[maxFaceSize]; int i; //used for iteration int alreadyCountedFace=0;// if the face has already been counted then this flag will be triggered // Begin by setting the numberOfFaces to 0 : this way we can increment them for(i = 0; i < maxFaceSize; i++) { numberOfFaces[i] = 0; } /*The strategy for processing faceArray and faceArray lengths is to go through each vertex (except the last two - no face can contain just two vertices) and try to work out which face it is connected to. We must use the assumption of the clockwise ordering of the edges in vertices */ for(vertexNumber = 0; vertexNumber < (numberOfVertices - 2); vertexNumber++) { for(edgeProcessed = 0; edgeProcessed < vertexDegrees[vertexNumber]; edgeProcessed++) { //We traverse the face until we return to the starting vertex previousVertex = vertexNumber; nextVertex = vertices[previousVertex][edgeProcessed]; alreadyCountedFace = 0; /*Avoid faces which have already been counted*/ if(!(nextVertex < previousVertex)) // if nextVertex < previousVertex then all faces attached to this vertex will have already been added { currentFace[0] = vertexNumber; faceSides = 1; while(nextVertex != vertexNumber) { currentFace[faceSides] = nextVertex; //add this vertex to the list of vertices traversed nextVertex = findNextVertex(vertices,previousVertex,nextVertex,vertexDegrees[nextVertex]); previousVertex = currentFace[faceSides]; faceSides++; /*Avoid faces which have already been counted*/ if(nextVertex 2) { // we need to add this face to the list of faces for(i = 0; i < faceSides; i++) { faceArray[(faceSides - 1)][numberOfFaces[(faceSides - 1)]][i] = currentFace[i]; } numberOfFaces[(faceSides - 1)]++; } } } } } int setAngleToRigid(int nonRigidVertexDegrees[expectedVertices],int nonRigidEdges[expectedVertices][currentMaxDegree], double angleArray[expectedVertices][currentMaxDegree],int rigidEdgeArray[expectedVertices][currentMaxDegree], int rigidVertexDegrees[expectedVertices], int centreVertex,int anticlockwiseVertex,double angleValue) { //sets the angle with centre at centreVertex and with anticlockwise vertex to be angleValue int currentEdge = 0; //the current edge (may be rigid) int i; //used to iterate through non rigid edges int edgeLeadingFromCentreVertex; for(i = 0; i < nonRigidVertexDegrees[centreVertex]; currentEdge++) { edgeLeadingFromCentreVertex = nonRigidEdges[centreVertex][currentEdge]; if(edgeLeadingFromCentreVertex == anticlockwiseVertex) { //we have found the correct angle nonRigidVertexDegrees[centreVertex] --; //there is one less edge attached to the centre vertex nonRigidEdges[centreVertex][currentEdge] = RIGID; //mark this edge as rigid in the list of non rigid edges angleArray[centreVertex][currentEdge] = angleValue; //set the angle in angleArray rigidEdgeArray[centreVertex][currentEdge] = anticlockwiseVertex; //set the edge to be rigid in the list of rigid edges rigidVertexDegrees[centreVertex]++; return 1; //we removed the vertex successfully } if(edgeLeadingFromCentreVertex != -1) { i++; //we encountered a non rigid edge so we need to increment } } return 0; //the edge may have already been rigid } int isAngleRigid(int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int centreVertex,int anticlockwiseVertex) { //returns the position in the array of the rigid angle if the angle starting at centre vertex with anticlockwise vertex anticlockwisevertex if that angle is rigid, or -1 if it is not. int currentRigidEdge = 0; //the current edge (may not be rigid - if it is not then rigidEdges[centreVertex][currentRigidEdge] will be -1) int i; //used to iterate through non rigid edges int edgeLeadingFromCentreVertex; for(i = 0; i < rigidVertexDegrees[centreVertex]; currentRigidEdge++) { edgeLeadingFromCentreVertex = rigidEdges[centreVertex][currentRigidEdge]; if(edgeLeadingFromCentreVertex == anticlockwiseVertex) { //we have found the correct angle return currentRigidEdge; //we removed the vertex successfully } if(edgeLeadingFromCentreVertex != -1) { i++; //we encountered a non rigid edge so we need to increment } } return -1; } int fixFinalAnglesInVertices(int nonRigidVertexDegrees[expectedVertices],int nonRigidEdges[expectedVertices][currentMaxDegree],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize], double angleArray[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int rigidEdges[expectedVertices][currentMaxDegree],double totalAngleValue) { int currentVertex; //used to iterate through each vertex to check whether or not it can be made rigid int numberOfNewRigidAngles = 0; //the number of rigid angles we have set through this process int currentRigidAngle; int numberOfRigidAngles; //the number of rigid angles we have already summed double currentRigidAngleSum; //the sum of all angles so far attached to the vertex int remainingNonRigidAngle; //the last angle that remains int currentNonRigidAngle; //used to iterate through possible non rigid angles for(currentVertex = 0; currentVertex < expectedVertices;currentVertex ++) { if(nonRigidVertexDegrees[currentVertex] == 1) { //this vertex can be fixed as rigid - it has one edge remaining. We need to determine what the final angle is currentRigidAngleSum = 0; currentRigidAngle = 0; for(numberOfRigidAngles = 0; numberOfRigidAngles < rigidVertexDegrees[currentVertex]; currentRigidAngle ++) { if(rigidEdges[currentVertex][currentRigidAngle] != -1) //this edge is rigid { currentRigidAngleSum += angleArray[currentVertex][currentRigidAngle]; numberOfRigidAngles ++; } } //now we need to determine which angle remains remainingNonRigidAngle = -1; currentNonRigidAngle = 0; while(remainingNonRigidAngle == -1) { if(nonRigidEdges[currentVertex][currentNonRigidAngle] != -1) { remainingNonRigidAngle = nonRigidEdges[currentVertex][currentNonRigidAngle]; } currentNonRigidAngle++; } numberOfNewRigidAngles += setAngleToRigid(rigidAngleInformation,currentVertex,remainingNonRigidAngle,(totalAngleValue - currentRigidAngleSum)); } } return numberOfNewRigidAngles; } int fixOppositeAnglesEqualInFace(int nonRigidVertexDegrees[expectedVertices],int nonRigidEdges[expectedVertices][currentMaxDegree],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize], double angleArray[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int rigidEdges[expectedVertices][currentMaxDegree],int faceSize) { int currentFace; //used to iterate through faces int currentVertex; int previousVertex; int halfFaceSize = faceSize / 2; int numberOfNewRigidAngles = 0; //used to count how many rigid angles have been fixed by this process int previousOppositeVertex; int currentOppositeVertex; int rigidAnglePosition; //used to store whether or not an angle is rigid, and if it is, the position in the array int currentPositionInFace; //used to iterate through the currentFace if((halfFaceSize * 2) != faceSize) { #pragma omp critical(printerror) { fprintf(stderr,"You cannot refer to opposite angles in faces with an odd number of vertices"); } return 0; } for(currentFace = 0; currentFace < numberOfFaces[(faceSize - 1)]; currentFace++) { previousVertex = faceArray[(faceSize - 1)][currentFace][(faceSize - 1)]; previousOppositeVertex = faceArray[(faceSize - 1)][currentFace][(halfFaceSize - 1)]; //the opposite angle should be half the facesize away from the original angle for(currentPositionInFace = 0; currentPositionInFace < halfFaceSize;currentPositionInFace++) { //work out what the next vertex is currentVertex = faceArray[(faceSize - 1)][currentFace][currentPositionInFace]; //now we need to get the vertex opposite currentOppositeVertex = faceArray[(faceSize - 1)][currentFace][(currentPositionInFace+halfFaceSize)]; rigidAnglePosition = isAngleRigid(rigidEdges,rigidVertexDegrees,currentVertex,previousVertex); if(rigidAnglePosition != - 1) //we have a rigid angle centred at the currentVertex. Set the angle at currentOppositeVertex to be rigid { numberOfNewRigidAngles += setAngleToRigid(rigidAngleInformation,currentOppositeVertex,previousOppositeVertex,angleArray[currentVertex][rigidAnglePosition]); } else { //the first angle was not rigid - see if the opposite angle is rigid rigidAnglePosition = isAngleRigid(rigidEdges,rigidVertexDegrees,currentOppositeVertex,previousOppositeVertex); if(rigidAnglePosition != -1) //a rigid angle centred at the currentOppositeVertex. Set the angle at currentVertex to be rigid { numberOfNewRigidAngles += setAngleToRigid(rigidAngleInformation,currentOppositeVertex,previousOppositeVertex,angleArray[currentVertex][rigidAnglePosition]); } } previousVertex = currentVertex; previousOppositeVertex = currentOppositeVertex; } } return numberOfNewRigidAngles; } int fixAnglesInFace(int nonRigidVertexDegrees[expectedVertices],int nonRigidEdges[expectedVertices][currentMaxDegree],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize], double angleArray[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int rigidEdges[expectedVertices][currentMaxDegree],int faceSize,double value) { //cycle through all triangles int currentFace; //used to iterate through faces int currentPositionInFace; //used to iterate through the vertices in an individual triangle int currentVertex; int previousVertex; int numberOfNewRigidAngles = 0; //used to count how many rigid angles have been fixed by this process for(currentFace = 0;currentFace < numberOfFaces[(faceSize - 1)];currentFace++) { previousVertex = faceArray[(faceSize - 1)][currentFace][(faceSize - 1)]; //set the last vertex in the face to be the first anticlockwise edge for(currentPositionInFace=0; currentPositionInFace < faceSize; currentPositionInFace++) { currentVertex = faceArray[2][currentFace][currentPositionInFace]; numberOfNewRigidAngles += setAngleToRigid(rigidAngleInformation,currentVertex,previousVertex,value); //set the angle between these two vertices to rigid previousVertex = currentVertex; } } return numberOfNewRigidAngles; } void clearRigidAngleVariables(int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices]) { int i; //used to iterate through the arrays int currentVertex; //used to iterate through the rigidEdges array for (i=0;i 0) { //the last vertex was isolated setVertexNonRigid(rigidVertices,positionInRigidVertexArray); //set this vertex to non rigid *numberOfRigidVertices = ((*numberOfRigidVertices) - 1); //reduce the number of rigid vertices verticesSetAsIsolated = 1; returnValue ++; } currentRigidVertex++; } } } return returnValue; } void removeNonRigidEdges(int rigidVertices[expectedVertices],double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int centreVertex, double newAngleArray[expectedVertices][currentMaxDegree],int newRigidEdges[expectedVertices][currentMaxDegree],int newRigidVertexDegrees[expectedVertices]) { //this function removes all edges which are not connected to a rigid vertex centered at centreVertex int currentEdge = -1; //we need to process all rigid edges that this vertex possesses int newRigidVertexDegree = 0; //used to store the degree of the vertex when all the non rigid edges are removed int currentConnectedVertex; int isRigid = 0; int haveGoneThroughAllVertices = 0; int myRigidVertexNumber = rigidVertices[centreVertex]; int rigidConnectedVertexNumber; int returnValue; //returnValue is the number of vertices eliminated double totalAngle; while(!isRigid) //find the first rigid vertex connected to this one (there must be at least one or this vertex would not be considered rigid after removeIsolatedRigidVertices) { currentEdge++; currentConnectedVertex = rigidEdges[centreVertex][currentEdge]; if(rigidVertices[currentConnectedVertex]!=-1 ) { isRigid = 1; //we have found a rigid vertex } } while(!haveGoneThroughAllVertices) { rigidConnectedVertexNumber = rigidVertices[currentConnectedVertex]; newRigidEdges[myRigidVertexNumber][newRigidVertexDegree] = rigidConnectedVertexNumber; isRigid = 0; totalAngle = 0; //we need to work out what the angle is while(!isRigid) //find the first rigid vertex connected to this one (there must be at least one or this vertex would not be considered rigid after removeIsolatedRigidVertices) { totalAngle += angleArray[centreVertex][currentEdge]; //add this angle currentEdge++; if(currentEdge == rigidVertexDegrees[centreVertex]) { //we have gone through all vertices haveGoneThroughAllVertices = 1; currentEdge = 0; } currentConnectedVertex = rigidEdges[centreVertex][currentEdge]; if(rigidVertices[currentConnectedVertex]!=-1 ) { isRigid = 1; //we have found a rigid vertex } } newAngleArray[myRigidVertexNumber][newRigidVertexDegree] = totalAngle; newRigidVertexDegree++; } newRigidVertexDegrees[myRigidVertexNumber] = newRigidVertexDegree; } void eliminateNonRigidVertices(int rigidVertices[expectedVertices],double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int *numberOfRigidVertices,double newAngleArray[expectedVertices][currentMaxDegree],int newRigidEdges[expectedVertices][currentMaxDegree],int newRigidVertexDegrees[expectedVertices]) { int currentVertex; removeIsolatedRigidVertices(rigidVertices,rigidEdges,rigidVertexDegrees,numberOfRigidVertices); //all vertices that are not of degree > 1 are useless to us for(currentVertex = 0; currentVertex < expectedVertices;currentVertex++) { if(rigidVertices[currentVertex]!=-1) { removeNonRigidEdges(rigidVertices,angleArray,rigidEdges,rigidVertexDegrees,currentVertex,newAngleArray,newRigidEdges,newRigidVertexDegrees); //delete all edges that are not attached to a rigid vertex } } } int allowedEdge(int intendedVertex,int disallowedVertices[expectedVertices],int numberOfDisallowedVertices) { int i; //check to see whether an edge is allowed, given a list of vertices that we are not allowed to connect to for(i = 0; i < numberOfDisallowedVertices;i++) { if(intendedVertex == disallowedVertices[i]) return 0; } return 1; } int findPath(int startVertex,int finishVertex,int disallowedVertices[expectedVertices],int numberOfDisallowedVertices,int vertices[expectedVertices][currentMaxDegree],int vertexDegrees[expectedVertices]) { int currentEdge; int nextVertex; for(currentEdge = 0; currentEdge < vertexDegrees[startVertex]; currentEdge++) { nextVertex = vertices[startVertex][currentEdge]; if(nextVertex == finishVertex) { return 1; //we have returned to the intended vertex. This must mean that this path is possible } if(allowedEdge(nextVertex,disallowedVertices,numberOfDisallowedVertices)) { //we are allowed to use this edge - try to get to the finish vertex disallowedVertices[numberOfDisallowedVertices] = startVertex; //ensure that we cannot return to this vertex numberOfDisallowedVertices++; if(findPath(nextVertex,finishVertex,disallowedVertices,numberOfDisallowedVertices,vertices,vertexDegrees)) { return 1; //we managed to get to the intended vertex } else { //we failed to get to the intended vertex by going to nextVertex numberOfDisallowedVertices --; } } } return 0; } int isEdgeAttachedToFace(int firstVertex,int secondVertex,int vertices[expectedVertices][currentMaxDegree],int vertexDegrees[expectedVertices]) { //we see if this edge actually connects a face. This is necessary because retrieveFaces crashes if there is an edge that is not part of a face int disallowedVertices[expectedVertices]; int numberOfDisallowedVertices = 1; int currentEdge; int numberOfEdges = vertexDegrees[secondVertex]; disallowedVertices[0] = secondVertex; for(currentEdge = 0; currentEdge < numberOfEdges;currentEdge++) { if(vertices[secondVertex][currentEdge]!=firstVertex) { if(findPath(vertices[secondVertex][currentEdge],firstVertex,disallowedVertices,numberOfDisallowedVertices,vertices,vertexDegrees)) { //we found a path back to the start vertex return 1; } } } return 0; } void removeEdge(int vertices[expectedVertices][currentMaxDegree],int vertexDegrees[expectedVertices],double angleArray[expectedVertices][currentMaxDegree],int startVertex,int edge) { //removes the edge connected to startVertex int currentVertexDegree = vertexDegrees[startVertex]; int currentEdge; //the current edge that is being processed //we need to fix the angles (they change after removing this edge) if(edge == 0) { angleArray[startVertex][currentVertexDegree - 1] += angleArray[startVertex][0]; } else { angleArray[startVertex][(edge - 1)] += angleArray[startVertex][edge]; } //now that the angles have been corrected, we must delete this edge for(currentEdge = (edge + 1);currentEdge currentVertex) //if this is not the case then the edge has already been processed { //we should check to see whether this edge works if(!isEdgeAttachedToFace(currentVertex,nextVertex,vertices,vertexDegrees)) { //this edge is not attached to a face. We must remove it. removeEdge(vertices,vertexDegrees,angleArray,currentVertex,currentEdge); secondEdge = findConnectingEdge(nextVertex,currentVertex,vertices,vertexDegrees); removeEdge(vertices,vertexDegrees,angleArray,nextVertex,secondEdge); returnValue = 1; } } } } return returnValue; } void findRigidGraph(int graph[sizeOfGraph],int lengthOfGraph,double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],int *numberOfRigidVertices,double newAngleArray[expectedVertices][currentMaxDegree],int newRigidEdges[expectedVertices][currentMaxDegree],int newRigidVertexDegrees[expectedVertices],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize], int numberOfFaces[maxFaceSize]) { int vertexDegrees[expectedVertices]; int vertices[expectedVertices][currentMaxDegree]; int fixedVertices[expectedVertices]; *numberOfRigidVertices =0; //the number of rigid vertices splitGraphToVertices(graph,lengthOfGraph,vertexDegrees,vertices); retrieveFaces(vertexDegrees,vertices,faceArray,numberOfFaces,expectedVertices); findRigidAngles(vertexDegrees,vertices,faceArray,numberOfFaces,angleArray,rigidEdges,rigidVertexDegrees); //we need to produce a list of rigid vertices getListOfRigidVertices(vertexDegrees, fixedVertices, numberOfRigidVertices); if(*numberOfRigidVertices!=expectedVertices) { eliminateNonRigidVertices(fixedVertices,angleArray,rigidEdges,rigidVertexDegrees,numberOfRigidVertices,newAngleArray,newRigidEdges,newRigidVertexDegrees);//remove all vertices from the graph which were not rigid removeEdgesNotAttachedToFace(newRigidEdges,newRigidVertexDegrees,newAngleArray,*numberOfRigidVertices); //remove all edges that were not attached to a face retrieveFaces(newRigidVertexDegrees,newRigidEdges,faceArray,numberOfFaces,*numberOfRigidVertices); //get the faces in the new rigid graph *numberOfRigidVertices = *numberOfRigidVertices; } } void applyMatrixToVector(double matrix[3][3],double previousVector[3],double resultVector[3]) { //sets resultVector = matrix(previousVector) double temporaryValue; //used to store the current entry for each element of the matrix int i; //the row of the matrix int j; //the column of the matrix for(i=0;i<3;i++) { temporaryValue = 0; for(j=0;j<3;j++) { temporaryValue += matrix[i][j] *previousVector[j]; } resultVector[i] = temporaryValue; } } void getCoordinatesGivenAngle(double centreVertexCoordinates[3],double previousVertexCoordinates[3],double newVertexCoordinates[3],double angleOfRotation) { double x = centreVertexCoordinates[0]; double y = centreVertexCoordinates[1]; double z = centreVertexCoordinates[2]; double XZ = sqrt(square(x)+ square(z)); double YZ = sqrt(square(y) + square(z)); double cosTheta = cos(angleOfRotation); double sineTheta = sin(angleOfRotation); double matrix[3][3]; //set the matrix up so that it is a rotation of angle angleOfRotation about the axis centreVertex matrix[0][0] = z*cosTheta/XZ; matrix[0][1] = x*y*cosTheta/(XZ*YZ) + z*sineTheta/YZ; matrix[0][2] = x*z*cosTheta/(XZ*YZ) - y*sineTheta/YZ; matrix[1][0] = -z*sineTheta/(XZ); matrix[1][1] = z*cosTheta/(YZ) - x*y*sineTheta/(XZ*YZ); matrix[1][2] = -y*cosTheta/(YZ) - x*z*sineTheta/(XZ*YZ); matrix[2][0] = -x/(XZ); matrix[2][1] = y*z/(XZ*YZ); matrix[2][2] = square(z)/(XZ*YZ); applyMatrixToVector(matrix,previousVertexCoordinates,newVertexCoordinates); } double getAngle(int anticlockwiseVertex,int centreVertex,double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int centreVertexDegree) { int i; //we go through all vertices to find the correct angle for(i = 0; i < centreVertexDegree; i++) { if(rigidEdges[centreVertex][i] == anticlockwiseVertex) { //we found the edge required return angleArray[centreVertex][i]; } } fprintf(stderr,"There was an error getting the requested angle"); return -1; //error! } void buildArrayOfAnglesInFace(int face[maxFaceSize],int faceSize,double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices],double *faceAngles) { int currentAngleInFace = 0; int currentVertex = face[0]; int previousVertex = face[(faceSize - 1)]; for(currentAngleInFace = 0; currentAngleInFace < faceSize;currentAngleInFace++) //go around the face getting the angles at each vertex { currentVertex = face[currentAngleInFace]; faceAngles[currentAngleInFace] = getAngle(previousVertex,currentVertex,angleArray,rigidEdges,rigidVertexDegrees[currentVertex]); previousVertex = currentVertex; } } int deriveContradictionInFace(int face[maxFaceSize],int sizeOfFace,double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices]) { double clockwiseCoordinates[maxFaceSize][3]; double anticlockwiseCoordinates[maxFaceSize][3]; double startVertexCoordinates[3]; int currentPositionInFace; int previousPositionInFace; double *faceAngles; double delta = 0.001; double sqrt3 = 1.73205080757; faceAngles = (double *)malloc(sizeOfFace * sizeof(double)); buildArrayOfAnglesInFace(face,sizeOfFace,angleArray,rigidEdges,rigidVertexDegrees,faceAngles); //find a list of angles in this face //set the first coordinate to be at the top clockwiseCoordinates[0][0] = 0; clockwiseCoordinates[0][1] = 0; clockwiseCoordinates[0][2] = 1; //we fix the second vertex. This should fix the entire face clockwiseCoordinates[1][0] = sqrt3/(double)2; clockwiseCoordinates[1][1] = 0; clockwiseCoordinates[1][2] = (double)(1)/(double)2; for(currentPositionInFace = 2; currentPositionInFace < sizeOfFace;currentPositionInFace++) { getCoordinatesGivenAngle(clockwiseCoordinates[currentPositionInFace - 1],clockwiseCoordinates[currentPositionInFace - 2],clockwiseCoordinates[currentPositionInFace],faceAngles[(currentPositionInFace - 1)]); } getCoordinatesGivenAngle(clockwiseCoordinates[sizeOfFace - 1],clockwiseCoordinates[sizeOfFace - 2],startVertexCoordinates,faceAngles[sizeOfFace - 1]); if(vectorDistance(startVertexCoordinates,clockwiseCoordinates[0])>delta) //if the position obtained by traversing the face does not agree with the position of the original fixed vertex then there was a contradiction { //contradiction return 1; } anticlockwiseCoordinates[0][0] = 0; anticlockwiseCoordinates[0][1] = 0; anticlockwiseCoordinates[0][2] = 1; anticlockwiseCoordinates[1][0] = sqrt3/(double)2; anticlockwiseCoordinates[1][1] = 0; anticlockwiseCoordinates[1][2] = (double)1/(double)2; getCoordinatesGivenAngle(anticlockwiseCoordinates[0],anticlockwiseCoordinates[1],anticlockwiseCoordinates[sizeOfFace - 1],-faceAngles[0]); if(vectorDistance(anticlockwiseCoordinates[sizeOfFace - 1],clockwiseCoordinates[sizeOfFace - 1])>delta) { return 1; } previousPositionInFace = 0; for(currentPositionInFace = (sizeOfFace - 1); currentPositionInFace > 2;currentPositionInFace--) { getCoordinatesGivenAngle(anticlockwiseCoordinates[currentPositionInFace],anticlockwiseCoordinates[previousPositionInFace],anticlockwiseCoordinates[currentPositionInFace - 1],-faceAngles[(currentPositionInFace)]); previousPositionInFace = currentPositionInFace; if(vectorDistance(anticlockwiseCoordinates[currentPositionInFace - 1],clockwiseCoordinates[currentPositionInFace - 1])>delta) //we traverse the face in two different ways - clockwise or anticlockwise. If any vertices have different coordinates after both traversals then filter this graph { return 1; } } return 0; } int deriveContradiction(double angleArray[expectedVertices][currentMaxDegree],int rigidEdges[expectedVertices][currentMaxDegree],int rigidVertexDegrees[expectedVertices], int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize]) { int currentFaceSize = 3; //look at all faces, including triangles (triangles may not be the same as in the original contact graph, and so contradictions may still arise. int currentFace = 0; //start with the first quadrilateral int contradiction = 0; //there has not been a contradiction so far for(currentFaceSize = 3; currentFaceSize < maxFaceSize; currentFaceSize ++) { for(currentFace = 0; currentFace < numberOfFaces[currentFaceSize - 1]; currentFace++) //test each face with size currentFaceSize { contradiction += deriveContradictionInFace(faceArray[currentFaceSize - 1][currentFace],currentFaceSize,angleArray,rigidEdges,rigidVertexDegrees); } } if(contradiction > 0) { return 1; //there was a contradiction } return 0; } int isGraphLegal(int graph[sizeOfGraph],int lengthOfGraph) { int rigidFaceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize]; //int ***rigidFaceArray; int i,j; int numberOfRigidFaces[maxFaceSize]; double angleArray[expectedVertices][currentMaxDegree]; int rigidEdges[expectedVertices][currentMaxDegree]; int rigidVertexDegrees[expectedVertices]; int numberOfRigidVertices; double newAngleArray[expectedVertices][currentMaxDegree]; int newRigidEdges[expectedVertices][currentMaxDegree]; int newRigidVertexDegrees[expectedVertices]; int legality = 0; /*rigidFaceArray = malloc(maxFaceSize * sizeof(int **)); for(i = 0; i < maxFaceSize; i++) { rigidFaceArray[i] = malloc(maxNumberOfFaces * sizeof(int *)); for(j = 0; j < maxNumberOfFaces;j++) { rigidFaceArray[i][j] = malloc(maxFaceSize * sizeof(int)); } }*/ findRigidGraph(graph,lengthOfGraph,angleArray,rigidEdges,rigidVertexDegrees,&numberOfRigidVertices,newAngleArray,newRigidEdges,newRigidVertexDegrees,rigidFaceArray,numberOfRigidFaces); //get the rigid graph if(numberOfRigidVertices == expectedVertices) { legality = deriveContradiction(angleArray,rigidEdges,rigidVertexDegrees,rigidFaceArray,numberOfRigidFaces); //the rigid graph was the same as the original graph } else if(numberOfRigidVertices == 0) { legality = 0; //the graph had no rigid vertices, and so we should simply allow it } else { legality = deriveContradiction(newAngleArray,newRigidEdges,newRigidVertexDegrees,rigidFaceArray,numberOfRigidFaces); //check for contradictions in the remaining rigid graph } if(legality > 0) { if(numberOfRigidVertices > 3 && numberOfRigidVertices < expectedVertices) { #pragma omp atomic legitSemiRigidGraphs++; //this graph was not completely rigid, but had rigid vertices } else if(numberOfRigidVertices == expectedVertices) { #pragma omp atomic legitRigidGraphs++; //this graph was rigid } return 0; } return 1; } void passGraph(int graph[sizeOfGraph], int length) { int i = 0; int j = 0; if(isGraphLegal(graph,length)) //this graph was permitted { #pragma omp critical(writefile) //avoids OMP race conditions { totalCorrectGraphs++; fprintf(stderr, "Edges of remaining graph %i:\n", totalCorrectGraphs); printf("%c",expectedVertices); // fprintf(stderr, "%i ", expectedVertices); j = 1; for(i = 1; i0) { fprintf(stderr, "(%i,%i) ", j,graph[i]); } else { j++; fprintf(stderr, "\n"); } } fprintf(stderr,"\n"); } } } int main() { int vertexCount; //which vertex are we currently processing int graph[sizeOfGraph]; //the information contained within each graph (per thread) int lengthOfGraph = 0; //the length of the information contained within each graph (per thread) char a; //used to read the file int exitLoop = 0; int threadsToUse = maxThreads; #pragma omp parallel { //this section of code automatically determines the optimal number of threads, constrained by maxThreads if(omp_get_thread_num() == 0) //only the first thread should execute this piece of code { if(threadsToUse>omp_get_num_threads()) { threadsToUse = omp_get_num_threads(); //if the user has asked for too many threads then the number of threads will be reset to the optimal number } fprintf(stderr,"%i threads will be used to execute this program \n",threadsToUse); } } #pragma omp barrier #pragma omp parallel private(graph,lengthOfGraph,a) shared(exitLoop) num_threads(threadsToUse) //each thread has it's own graph, length of graph and input character(a). The flag to exit the loop is shared between all threads { a = 0; while(!exitLoop) { //if exitloop is set then all threads will eventually exit this loop #pragma omp critical(readGraph) //reading a graph needs to be done in serial, not parallel. { //a thread may enter this section of code even though there are no graphs left to process. Consequently we need to check whether or not another thread has set exitloop lengthOfGraph = 0; if(!exitLoop) //we have not necessarily reached the end of file { a = getchar(); //read the first character if(a!=expectedVertices) { if(a!=EOF) { fprintf(stderr,"A graph started with %i vertices when %i vertices were expected",a,expectedVertices); } else { exitLoop = 1; //now we have reached the end of the file. Signal to all other threads to stop processing new graphs } } else { //read in the graph: there are some vertices left graph[0] = a; for(vertexCount=0;vertexCount