tt

VadereUtils/resources/DistMesh.cl

+
+//IDistanceFunction distanceFunc = p -> Math.abs(6 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 4;
+// abs(6 - length(v))-4
+
+inline void atomicAdd_g_f(volatile __global float2 *addr, float2 val) {
+    union{
+        unsigned int u32;
+        float2 f32;
+    } next, expected, current;
+   	current.f32 = *addr;
+    do{
+   	    expected.f32 = current.f32;
+        next.f32 = expected.f32 + val;
+   		current.u32 = atomic_cmpxchg( (volatile __global unsigned int *) addr, expected.u32, next.u32);
+    } while(current.u32 != expected.u32);
+}
+
 kernel void computeForces(
     __global float2* vertices,
     __global int4* edges,
@@ -6,6 +23,7 @@ kernel void computeForces(
     __global float2* forces)
 {
     int i = get_global_id(0);
+    int p1Index = edges[i].s0;
     float2 p1 = vertices[edges[i].s0];
     float2 p2 = vertices[edges[i].s1];
     float2 v = normalize(p1-p2);
@@ -13,23 +31,47 @@ kernel void computeForces(
 
     // F= ...
     float len = lengths[i].s0;
-    float desiredLen = lengths[i].s1 * (*scalingFactor) * 1.2f;
+    float desiredLen = lengths[i].s1 * (*scalingFactor) * 1.1f;
     float lenDiff = (desiredLen - len);
     lenDiff = lenDiff > 0 ? lenDiff : 0;
 
 
     float2 partialForce = v * lenDiff; // TODO;
 
-    forces[edges[i].s0] = forces[edges[i].s0] + partialForce; // TODO sync
+    atomicAdd_g_f(&forces[p1Index], partialForce.s0);
+    //atomicAdd_g_f(&forces[p1Index].s1, partialForce.s1);
+
+    //forces[edges[i].s0] = forces[edges[i].s0] + partialForce; // TODO sync
     //forces[edges[i].s0] = (float2) (1.0f, 1.0f);
     //forces[edges[i].s1] = forces[edges[i].s1] - partialForce; // TODO sync
 }
 
 kernel void moveVertices(__global float2* vertices, __global float2* forces, const float delta) {
     int i = get_global_id(0);
+
     //float2 force = (float2)(1.0, 1.0);
     float2 force = forces[i];
-    vertices[i] = vertices[i] + (force * delta);
+    float2 v = vertices[i];
+    v = v + (force * 0.01f);
+
+    // project back if necessary
+    float distance = fabs(6.0f - length(v))-4.0f;
+    if(distance <= 0) {
+
+    }
+    else {
+        float deps = 1.4901e-8f * 1.0f;
+        float2 dX = (deps, 0.0f);
+        float2 dY = (0.0f, deps);
+        float dGradPX = ((fabs(6 - length(v + dX))-4)-distance)/deps;
+        float dGradPY = ((fabs(6 - length(v + dY))-4)-distance)/deps;
+        float2 projection = (dGradPX * distance, dGradPY * distance);
+        v = v - projection;
+    }
+
+    // set force to 0.
+    forces[i] = (0.0f, 0.0f);
+    vertices[i] = v;
 }
 
 // computation of the scaling factor:
@@ -48,13 +90,14 @@ kernel void computeLengths(
     float desiredLen = 1.0f;
     float2 len = (float2) (length(v), desiredLen);
     lengths[i] = len;
-    qLengths[i] = (float2) (v.s0*v.s0 + v.s1*v.s1, desiredLen*desiredLen);
+    qLengths[i] = (float2) (length(v)*length(v), desiredLen*desiredLen);
 }
 
 
 // kernel for multiple work-groups
 kernel void computePartialSF(__const int size, __global float2* qlengths, __local float2* partialSums, __global float2* output) {
     int gid = get_global_id(0);
+    int lid = get_local_id(0);
 
     if(gid < size){
         int global_index = gid;
@@ -66,15 +109,16 @@ kernel void computePartialSF(__const int size, __global float2* qlengths, __loca
             global_index += get_global_size(0);
         }
 
-        int lid = get_local_id(0);
         int group_size = get_local_size(0);
-        float2 len = qlengths[gid];
+        float2 len = accumulator;
+
+        //float2 len = (float2)(1.f, 1.0f);
         partialSums[lid] = len;
 
         barrier(CLK_LOCAL_MEM_FENCE);
         // group_size has to be a power of 2!
         for(int i = group_size/2; i > 0; i>>=1){
-            if(lid < i && lid + i < size) {
+            if(lid < i && gid + i < size) {
                 partialSums[lid] += partialSums[lid + i];
             }
             barrier(CLK_LOCAL_MEM_FENCE);
@@ -84,6 +128,11 @@ kernel void computePartialSF(__const int size, __global float2* qlengths, __loca
             output[get_group_id(0)] = partialSums[0];
         }
     }
+    else {
+        if(lid == 0){
+            output[get_group_id(0)] = (float2)(0.0f, 0.0f);
+        }
+    }
 }
 
 // kernel for 1 work-group
@@ -106,6 +155,8 @@ kernel void computeCompleteSF(__const int size, __global float2* qlengths, __loc
 
 
         if(lid == 0) {
+          //*scaleFactor =  qlengths[10].s0;
+          //*scaleFactor =  partialSums[0].s0;
           *scaleFactor = sqrt(partialSums[0].s0 / partialSums[0].s1);
         }
     }

VadereUtils/resources/DistMeshDouble.cl

+
+//IDistanceFunction distanceFunc = p -> Math.abs(6 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 4;
+// abs(6 - length(v))-4
+
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
+#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
+#define LOCK(a) atomic_cmpxchg(a, 0, 1)
+#define UNLOCK(a) atomic_xchg(a, 0)
+
+inline void atomicAdd_g_f(volatile __global double2 *addr, double2 val) {
+    union{
+        unsigned int u32;
+        double2 f32;
+    } next, expected, current;
+   	current.f32 = *addr;
+    do{
+   	    expected.f32 = current.f32;
+        next.f32 = expected.f32 + val;
+   		current.u32 = atomic_cmpxchg( (volatile __global unsigned int *) addr, expected.u32, next.u32);
+    } while(current.u32 != expected.u32);
+}
+
+kernel void computeForces(
+    __global double2* vertices,
+    __global int4* edges,
+    __global double2* lengths,
+    __global double* scalingFactor,
+    __global double2* forces,
+    __global int* mutexes)
+{
+    int i = get_global_id(0);
+    int p1Index = edges[i].s0;
+    double2 p1 = vertices[edges[i].s0];
+    double2 p2 = vertices[edges[i].s1];
+    double2 v = normalize(p1-p2);
+
+    // F= ...
+    double len = lengths[i].s0;
+    double desiredLen = lengths[i].s1 * (*scalingFactor) * 1.2;
+    double lenDiff = (desiredLen - len);
+    lenDiff = lenDiff > 0 ? lenDiff : 0;
+    double2 partialForce = v * lenDiff;
+
+    volatile __global int* addr = &mutexes[p1Index];
+
+    // TODO does this sync work properly? This syncs too much?
+    int waiting = 1;
+    while (waiting) {
+        while (LOCK(addr)) {}
+        forces[p1Index] = forces[p1Index] + partialForce;
+        UNLOCK(addr);
+        waiting = 0;
+    }
+}
+
+inline double dabs(double d) {return d < 0 ? -d : d;}
+
+kernel void moveVertices(__global double2* vertices, __global double2* forces, const double delta) {
+    int i = get_global_id(0);
+
+    //double2 force = (double2)(1.0, 1.0);
+    double2 force = forces[i];
+    double2 v = vertices[i];
+
+
+    // project back if necessary
+    double distance = dabs(6.0 - length(v))-4.0;
+    if(distance <= 0) {
+        v = v + (force * 0.08);
+    }
+    else {
+        double deps = 0.0001;
+        double2 dX = (deps, 0.0);
+        double2 dY = (0.0, deps);
+        double dGradPX = ((dabs(6.0 - length(v + dX))-4.0)-distance) / deps;
+        double dGradPY = ((dabs(6.0 - length(v + dY))-4.0)-distance) / deps;
+        double2 projection = (dGradPX * distance, dGradPY * distance);
+        v = v - projection;
+    }
+
+    // set force to 0.
+    forces[i] = (0.0, 0.0);
+    vertices[i] = v;
+}
+
+// computation of the scaling factor:
+kernel void computeLengths(
+    __global double2* vertices,
+    __global int4* edges,
+    __global double2* lengths,
+    __global double2* qLengths)
+{
+    int i = get_global_id(0);
+    double2 p1 = vertices[edges[i].s0];
+    double2 p2 = vertices[edges[i].s1];
+    double2 v = p1-p2;
+
+    //TODO: desiredLenfunction required
+    double desiredLen = 1.0f;
+    double2 len = (double2) (length(v), desiredLen);
+    lengths[i] = len;
+    qLengths[i] = (double2) (length(v)*length(v), desiredLen*desiredLen);
+}
+
+
+// kernel for multiple work-groups
+kernel void computePartialSF(__const int size, __global double2* qlengths, __local double2* partialSums, __global double2* output) {
+    int gid = get_global_id(0);
+    int lid = get_local_id(0);
+
+    if(gid < size){
+        int global_index = gid;
+        double2 accumulator = (double2)(0, 0);
+        // Loop sequentially over chunks of input vector
+        while (global_index < size) {
+            double2 element = qlengths[global_index];
+            accumulator += element;
+            global_index += get_global_size(0);
+        }
+
+        int group_size = get_local_size(0);
+        double2 len = accumulator;
+
+        //double2 len = (double2)(1.0, 1.0);
+        partialSums[lid] = len;
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+        // group_size has to be a power of 2!
+        for(int i = group_size/2; i > 0; i>>=1){
+            if(lid < i && gid + i < size) {
+                partialSums[lid] += partialSums[lid + i];
+            }
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+
+        if(lid == 0){
+            output[get_group_id(0)] = partialSums[0];
+        }
+    }
+    else {
+        if(lid == 0){
+            output[get_group_id(0)] = (double2)(0.0, 0.0);
+        }
+    }
+}
+
+// kernel for 1 work-group
+kernel void computeCompleteSF(__const int size, __global double2* qlengths, __local double2* partialSums, __global double* scaleFactor) {
+    int gid = get_global_id(0);
+    int lid = get_local_id(0);
+    if(lid < size){
+        int group_size = get_local_size(0);
+        partialSums[lid] = qlengths[lid];
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        // group_size has to be a power of 2!
+        for(int i = group_size/2; i > 0; i>>=1){
+            if(lid < i && lid + i < size) {
+                partialSums[lid] += partialSums[lid + i];
+            }
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+
+
+        if(lid == 0) {
+          //*scaleFactor =  qlengths[10].s0;
+          //*scaleFactor =  partialSums[0].s0;
+          *scaleFactor = sqrt(partialSums[0].s0 / partialSums[0].s1);
+        }
+    }
+}
\ No newline at end of file

VadereUtils/src/org/vadere/util/geometry/mesh/gen/AFace.java

@@ -5,7 +5,7 @@ import org.vadere.util.geometry.mesh.inter.IFace;
 import org.vadere.util.geometry.shapes.IPoint;
 
 /**
- * Created by bzoennchen on 06.09.17.
+ * @author Benedikt Zoennchen
  */
 public class AFace<P extends IPoint> implements IFace<P> {
 
@@ -46,10 +46,6 @@ public class AFace<P extends IPoint> implements IFace<P> {
 		this.id = id;
 	}
 
-	AFace() {
-		this.border = false;
-	}
-
 	boolean isBorder() {
 		return border;
 	}
@@ -64,10 +60,18 @@ public class AFace<P extends IPoint> implements IFace<P> {
 	 *
 	 * @param edge half-edge bordering this face
 	 */
-	void setEdge(final int edge) {
+	void setEdge(@NotNull final int edge) {
 		this.edge = edge;
 	}
 
+	/**
+	 * This method should only be called by the garbage collector in AMesh.
+	 * @param id
+	 */
+	void setId(@NotNull final int id) {
+		this.id = id;
+	}
+
 	int getEdge() {
 		return edge;
 	}

VadereUtils/src/org/vadere/util/geometry/mesh/gen/AHalfEdge.java

@@ -114,6 +114,14 @@ public class AHalfEdge<P extends IPoint> implements IHalfEdge<P> {
 		return destroyed;
 	}
 
+	/**
+	 * This method should only be called by the garbage collector in AMesh.
+	 * @param id
+	 */
+	void setId(@NotNull final int id) {
+		this.id = id;
+	}
+
 	@Override
 	public String toString() {
 		return ""+id;

VadereUtils/src/org/vadere/util/geometry/mesh/gen/AMesh.java

@@ -357,4 +357,89 @@ public class AMesh<P extends IPoint> implements IMesh<P, AVertex<P>, AHalfEdge<P
 	public void unlock(@NotNull AVertex<P> vertex) {
         vertex.getLock().unlock();
 	}
+
+	/**
+	 * removes all destroyed object from this mesh and re-arranges all indices.
+	 */
+	public void garbageCollection() {
+		Map<Integer, Integer> faceIdMap = new HashMap<>();
+		Map<Integer, Integer> edgeIdMap = new HashMap<>();
+		Map<Integer, Integer> vertexIdMap = new HashMap<>();
+
+		int i = 0;
+		int j = 0;
+		for(AFace<P> face : faces) {
+			if(face.isDestroyed()) {
+				j--;
+			}
+			else {
+				faceIdMap.put(i, j);
+			}
+			i++;
+			j++;
+		}
+
+		i = 0;
+		j = 0;
+		for(AHalfEdge<P> edge : edges) {
+			if(edge.isDestroyed()) {
+				j--;
+			}
+			else {
+				edgeIdMap.put(i, j);
+			}
+			i++;
+			j++;
+		}
+
+		i = 0;
+		j = 0;
+		for(AVertex<P> vertex : vertices) {
+			if(vertex.isDestroyed()) {
+				j--;
+			}
+			else {
+				vertexIdMap.put(i, j);
+			}
+			i++;
+			j++;
+		}
+
+		faces = faces.stream().filter(f -> !f.isDestroyed()).collect(Collectors.toList());
+		edges = edges.stream().filter(e -> !e.isDestroyed()).collect(Collectors.toList());
+		vertices = vertices.stream().filter(v -> !v.isDestroyed()).collect(Collectors.toList());
+
+		i = 0;
+		for(AFace<P> face : faces) {
+			face.setId(faceIdMap.get(face.getId()));
+			face.setEdge(edgeIdMap.get(face.getEdge()));
+			assert face.getId() == i;
+			i++;
+		}
+
+		i = 0;
+		for(AVertex<P> vertex : vertices) {
+			vertex.setId(vertexIdMap.get(vertex.getId()));
+			vertex.setEdge(edgeIdMap.get(vertex.getEdge()));
+			assert vertex.getId() == i;
+			i++;
+		}
+
+		i = 0;
+		for(AHalfEdge<P> edge : edges) {
+			edge.setId(edgeIdMap.get(edge.getId()));
+			edge.setEnd(vertexIdMap.get(edge.getEnd()));
+			edge.setNext(edgeIdMap.get(edge.getNext()));
+			edge.setPrevious(edgeIdMap.get(edge.getPrevious()));
+			edge.setTwin(edgeIdMap.get(edge.getTwin()));
+			if(edge.getFace() != boundary.getId()) {
+				edge.setFace(faceIdMap.get(edge.getFace()));
+			}
+
+			assert edge.getId() == i;
+			i++;
+		}
+
+		assert (getNumberOfVertices() == vertices.size()) && (getNumberOfEdges() == edges.size()) && (getNumberOfFaces() == faces.size());
+	}
 }

VadereUtils/src/org/vadere/util/geometry/mesh/gen/AVertex.java

@@ -8,7 +8,7 @@ import java.util.concurrent.locks.Lock;
 import java.util.concurrent.locks.ReentrantLock;
 
 /**
- * Created by bzoennchen on 06.09.17.
+ * @author Benedikt Zoennchen
  */
 public class AVertex<P extends IPoint> implements IVertex<P> {
 
@@ -16,7 +16,7 @@ public class AVertex<P extends IPoint> implements IVertex<P> {
 	private final P point;
 	private int down;
 	private int halfEdge;
-	private final int id;
+	private int id;
 	private boolean destroyed;
 
 
@@ -64,6 +64,14 @@ public class AVertex<P extends IPoint> implements IVertex<P> {
 		return point.equals(((AVertex<P>)obj).getPoint());
 	}
 
+	/**
+	 * This method should only be called by the garbage collector in AMesh.
+	 * @param id
+	 */
+	void setId(@NotNull final int id) {
+		this.id = id;
+	}
+
 	public boolean isDestroyed() {
 		return destroyed;
 	}

VadereUtils/src/org/vadere/util/geometry/mesh/gen/CLGatherer.java

@@ -2,8 +2,10 @@ package org.vadere.util.geometry.mesh.gen;
 
 import org.jetbrains.annotations.NotNull;
 import org.lwjgl.system.MemoryStack;
+import org.lwjgl.system.MemoryUtil;
 import org.vadere.util.geometry.shapes.IPoint;
 
+import java.nio.DoubleBuffer;
 import java.nio.FloatBuffer;
 import java.nio.IntBuffer;
 import java.util.Collection;
@@ -13,11 +15,26 @@ import java.util.Collection;
  */
 public class CLGatherer {
 
-    public static <P extends IPoint> FloatBuffer getVertices(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
+    public static <P extends IPoint> DoubleBuffer getVerticesD(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
         Collection<AVertex<P>> vertices = mesh.getVertices();
-        FloatBuffer vertexBuffer = stack.mallocFloat(vertices.size()*2);
+        DoubleBuffer vertexBuffer = MemoryUtil.memAllocDouble(vertices.size()*2);
         int index = 0;
         for(AVertex<P> vertex : vertices) {
+            assert index/2.0 == vertex.getId();
+            vertexBuffer.put(index, vertex.getX());
+            index++;
+            vertexBuffer.put(index, vertex.getY());
+            index++;
+        }
+        return vertexBuffer;
+    }
+
+    public static <P extends IPoint> FloatBuffer getVerticesF(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
+        Collection<AVertex<P>> vertices = mesh.getVertices();
+        FloatBuffer vertexBuffer = MemoryUtil.memAllocFloat(vertices.size()*2);
+        int index = 0;
+        for(AVertex<P> vertex : vertices) {
+            assert index/2.0 == vertex.getId();
             vertexBuffer.put(index, (float)vertex.getX());
             index++;
             vertexBuffer.put(index, (float)vertex.getY());
@@ -28,7 +45,7 @@ public class CLGatherer {
 
     public static <P extends IPoint> IntBuffer getEdges(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
         Collection<AHalfEdge<P>> edges = mesh.getEdges();
-        IntBuffer edgeBuffer = stack.mallocInt(edges.size()*4);
+        IntBuffer edgeBuffer =  MemoryUtil.memAllocInt(edges.size()*4);
         int index = 0;
         for(AHalfEdge<P> edge : edges) {
             edgeBuffer.put(index, mesh.getPrev(edge).getEnd());
@@ -45,7 +62,7 @@ public class CLGatherer {
 
     public static <P extends IPoint> IntBuffer getFaces(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
         Collection<AFace<P>> faces = mesh.getFaces();
-        IntBuffer faceBuffer = stack.mallocInt(faces.size()*4);
+        IntBuffer faceBuffer = MemoryUtil.memAllocInt(faces.size()*4);
         int index = 0;
         for(AFace<P> face : faces) {
             AHalfEdge<P> edge = mesh.getEdge(face);

VadereUtils/src/org/vadere/util/geometry/mesh/inter/IFace.java

@@ -2,6 +2,7 @@ package org.vadere.util.geometry.mesh.inter;
 
 import org.vadere.util.geometry.mesh.gen.*;
 import org.vadere.util.geometry.shapes.IPoint;
+import org.vadere.util.geometry.shapes.MPoint;
 import org.vadere.util.geometry.shapes.VPoint;
 
 /**
@@ -10,22 +11,22 @@ import org.vadere.util.geometry.shapes.VPoint;
  */
 public interface IFace<P extends IPoint> {
 
-    static AMesh<VPoint> createSimpleTriMesh() {
-        AMesh<VPoint> mesh;
-        AFace<VPoint> face1;
-        AFace<VPoint> face2;
-        AFace<VPoint> border;
-        AVertex<VPoint> x, y, z, w;
-        AHalfEdge<VPoint> zx ;
-        AHalfEdge<VPoint> xy;
-        AHalfEdge<VPoint> yz;
-
-        AHalfEdge<VPoint> wx;
-        AHalfEdge<VPoint> xz;
-        AHalfEdge<VPoint> yw;
-        AHalfEdge<VPoint> zy;
-
-        mesh = new AMesh<>((x1, y1) -> new VPoint(x1, y1));
+    static AMesh<MPoint> createSimpleTriMesh() {
+        AMesh<MPoint> mesh;
+        AFace<MPoint> face1;
+        AFace<MPoint> face2;
+        AFace<MPoint> border;
+        AVertex<MPoint> x, y, z, w;
+        AHalfEdge<MPoint> zx ;
+        AHalfEdge<MPoint> xy;
+        AHalfEdge<MPoint> yz;
+
+        AHalfEdge<MPoint> wx;
+        AHalfEdge<MPoint> xz;
+        AHalfEdge<MPoint> yw;
+        AHalfEdge<MPoint> zy;
+
+        mesh = new AMesh<>((x1, y1) -> new MPoint(x1, y1));
         border = mesh.createFace(true);