Commit 54b88621 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen

DistMesh on the GPU

parent 875c7cb7
......@@ -40,11 +40,13 @@ public class DelaunayHierarchy<P extends IPoint, V extends IVertex<P>, E extends
private Supplier<ITriangulation<P, V, E, F>> triangulationSupplier;
// see delaunay-hierarchy paper!
private double alpha = 30;
/*private double alpha = 30;
private int maxLevel = 5;
private int minSize = 20;*/
private int minSize = 20;
private double alpha = 30;
private int maxLevel = 5;
private int minSize = 20;
private LinkedList<F> prevLocationResult;
......@@ -125,6 +127,8 @@ public class DelaunayHierarchy<P extends IPoint, V extends IVertex<P>, E extends
prev = v;
}
}
//System.out.println(this);
}
private ITriangulation<P, V, E, F> getLevel(final int level) {
......@@ -140,7 +144,8 @@ public class DelaunayHierarchy<P extends IPoint, V extends IVertex<P>, E extends
}
private V getDown(V src, int srcLevel) {
return getLevel(srcLevel).getMesh().getDown(src);
// srcLevel-1 since the resulting vertex is contained in the mesh one level below the src vertex!
return getLevel(srcLevel-1).getMesh().getDown(src);
}
private int randomLevel() {
......@@ -230,7 +235,16 @@ public class DelaunayHierarchy<P extends IPoint, V extends IVertex<P>, E extends
}
}
/**
@Override
public String toString() {
StringBuilder builder = new StringBuilder("");
for(int i = 0; i < hierarchySets.size(); i++) {
builder.append("["+i+"]:" + hierarchySets.get(i).getMesh().getNumberOfVertices()+"\n");
}
return builder.toString();
}
/**
* Returns vertex of the triangulation of the face with the smallest distance to point.
* Assumption: The face has to be part of the mesh of the triangulation.
*
......
......@@ -384,7 +384,7 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
*/
@Override
public boolean isIllegal(E edge, V p) {
if(!mesh.isBoundary(mesh.getTwinFace(edge))) {
if(!mesh.isAtBoundary(edge)) {
//assert mesh.getVertex(mesh.getNext(edge)).equals(p);
//V p = mesh.getVertex(mesh.getNext(edge));
E t0 = mesh.getTwin(edge);
......@@ -547,7 +547,7 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
points.add(new VPoint(80,70));*/
Random r = new Random();
for(int i=0; i< 1000000; i++) {
for(int i=0; i< 10000; i++) {
VPoint point = new VPoint(width*r.nextDouble(), height*r.nextDouble());
points.add(point);
}
......@@ -560,11 +560,9 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
IPointLocator.Type.DELAUNAY_HIERARCHY,
points,
pointConstructor);
new IncrementalTriangulation<>(points);
bw.finalize();
Set<VLine> edges = bw.getEdges();
System.out.println(System.currentTimeMillis() - ms);
Set<VLine> edges = bw.getEdges();
JFrame window = new JFrame();
window.setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE);
......@@ -572,6 +570,22 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
window.getContentPane().add(new Lines(edges, points, max));
window.setVisible(true);
ms = System.currentTimeMillis();
ITriangulation<VPoint, AVertex<VPoint>, AHalfEdge<VPoint>, AFace<VPoint>> bw2 = ITriangulation.createATriangulation(
IPointLocator.Type.DELAUNAY_HIERARCHY,
points,
pointConstructor);
bw2.finalize();
System.out.println(System.currentTimeMillis() - ms);
Set<VLine> edges2 = bw2.getEdges();
JFrame window2 = new JFrame();
window2.setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE);
window2.setBounds(0, 0, max, max);
window2.getContentPane().add(new Lines(edges2, points, max));
window2.setVisible(true);
/*ms = System.currentTimeMillis();
BowyerWatsonSlow<VPoint> bw2 = new BowyerWatsonSlow<VPoint>(points, (x, y) -> new VPoint(x, y));
bw2.execute();
......@@ -600,7 +614,7 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
window3.getContentPane().add(new Lines(edges3, edges3.stream().flatMap(edge -> edge.streamPoints()).collect(Collectors.toSet()), max));
window3.setVisible(true);*/
VRectangle bound = new VRectangle(0, 0, width, height);
/*VRectangle bound = new VRectangle(0, 0, width, height);
ITriangulation triangulation = ITriangulation.createVPTriangulation(bound);
VPUniformRefinement uniformRefinement = new VPUniformRefinement(
triangulation,
......@@ -615,7 +629,7 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
window4.setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE);
window4.setBounds(0, 0, max, max);
window4.getContentPane().add(new Lines(edges4, edges4.stream().flatMap(edge -> edge.streamPoints()).collect(Collectors.toSet()), max));
window4.setVisible(true);
window4.setVisible(true);*/
}
private static class Lines extends JComponent{
......
......@@ -180,6 +180,16 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
return faces.size();
}
@Override
public boolean tryLock(@NotNull PVertex<P> vertex) {
return vertex.getLock().tryLock();
}
@Override
public void unlock(@NotNull PVertex<P> vertex) {
vertex.getLock().unlock();
}
@Override
public PHalfEdge<P> createEdge(@NotNull PVertex<P> vertex) {
PHalfEdge<P> edge = new PHalfEdge<>(vertex);
......@@ -273,11 +283,21 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
@Override
public Stream<PHalfEdge<P>> streamEdges() {
return edges.stream();
return edges.stream().filter(e -> !isDestroyed(e));
}
@Override
public Stream<PVertex<P>> streamVertices() { return vertices.stream(); };
public Stream<PHalfEdge<P>> streamEdgesParallel() {
return edges.parallelStream().filter(e -> !isDestroyed(e));
}
@Override
public Stream<PVertex<P>> streamVertices() { return vertices.stream().filter(v -> !isDestroyed(v)); }
@Override
public Stream<PVertex<P>> streamVerticesParallel() {
return vertices.parallelStream().filter(v -> !isDestroyed(v));
}
@Override
public Iterable<PHalfEdge<P>> getEdgeIt() {
......@@ -303,4 +323,8 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
return faces.stream().filter(predicate);
}
@Override
public int getNumberOfEdges() {
return edges.size();
}
}
......@@ -3,12 +3,16 @@ package org.vadere.util.geometry.mesh.gen;
import org.vadere.util.geometry.mesh.inter.IVertex;
import org.vadere.util.geometry.shapes.IPoint;
import java.util.concurrent.locks.Lock;
import java.util.concurrent.locks.ReentrantLock;
/**
* @author Benedikt Zoennchen
* @param <P>
*/
public class PVertex<P extends IPoint> implements IVertex<P> {
private final Lock lock;
private final P point;
private PVertex<P> down;
private PHalfEdge<P> halfEdge;
......@@ -18,6 +22,7 @@ public class PVertex<P extends IPoint> implements IVertex<P> {
this.point = point;
this.destroyed = false;
this.down = null;
this.lock = new ReentrantLock();
}
@Override
......@@ -62,6 +67,10 @@ public class PVertex<P extends IPoint> implements IVertex<P> {
destroyed = true;
}
public Lock getLock() {
return lock;
}
@Override
public int hashCode() {
return point.hashCode();
......
......@@ -5,12 +5,15 @@ import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.Vector2D;
import org.vadere.util.geometry.mesh.iterators.Ring1Iterator;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VPoint;
import java.util.List;
import java.util.Optional;
import java.util.Random;
import java.util.concurrent.locks.Lock;
import java.util.concurrent.locks.ReentrantLock;
import java.util.function.Predicate;
//TODO: check unused methods!
......@@ -38,6 +41,10 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
default void splitEdgeEvent(F original, F f1, F f2) {}
default Iterable<E> getRing1It(V vertex) {
return () -> new Ring1Iterator<>(getMesh(), getMesh().getEdge(vertex));
}
/**
* Inserts a point into the mesh which is contained in the boundary. The point will be
* connected by the start and end point of the edge. Note that the user has to make sure
......@@ -280,6 +287,65 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
return splitEdge(p, halfEdge, true);
}
/*default void flipLock(@NotNull final E edge) {
}*/
default void flipSync(@NotNull final E edge) {
IMesh<P, V, E, F> mesh = getMesh();
E a0 = edge;
E a1 = mesh.getNext(a0);
E a2 = mesh.getNext(a1);
E b0 = mesh.getTwin(edge);
E b1 = mesh.getNext(b0);
V v1 = mesh.getVertex(a0);
V v2 = mesh.getVertex(a1);
V v3 = mesh.getVertex(a2);
V v4 = mesh.getVertex(b1);
// TODO: a very first and simple aquire all locks implementation => improve it
while (true) {
// lock all 4 involved vertices
if (mesh.tryLock(v1)) {
if (mesh.tryLock(v2)) {
if (mesh.tryLock(v3)) {
if (mesh.tryLock(v4)) {
break;
}
else {
mesh.unlock(v3);
mesh.unlock(v2);
mesh.unlock(v1);
}
}
else {
mesh.unlock(v2);
mesh.unlock(v1);
}
} else {
mesh.unlock(v1);
}
}
}
try {
// if everything is locked flip
flip(edge);
}
// unlock all locks
finally {
mesh.unlock(v4);
mesh.unlock(v3);
mesh.unlock(v2);
mesh.unlock(v1);
}
}
/**
* Flips an edge in the triangulation assuming the egdge which will be
* created is not jet there.
......@@ -288,6 +354,7 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
* @param edge the edge which will be flipped.
*/
default void flip(@NotNull final E edge) {
IMesh<P, V, E, F> mesh = getMesh();
// 1. gather all the references required
......@@ -345,6 +412,22 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
flipEdgeEvent(fa, fb);
}
/**
* Tests if the face is counter-clockwise oriented.
* Assumption: The face is a triangle!
*
* @param triangleFace
* @return
*/
default boolean isCCW(F triangleFace) {
E edge = getMesh().getEdge(triangleFace);
P p1 = getMesh().getPoint(edge);
P p2 = getMesh().getPoint(getMesh().getNext(edge));
P p3 = getMesh().getPoint(getMesh().getPrev(edge));
return GeometryUtils.isCCW(p1, p2, p3);
}
default boolean isTriangle(F face) {
IMesh<P, V, E, F> mesh = getMesh();
List<E> edges = mesh.getEdges(face);
......
package org.vadere.util.opencl.examples;
import org.apache.commons.math.util.MathUtils;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.lwjgl.BufferUtils;
import org.lwjgl.PointerBuffer;
import org.lwjgl.opencl.*;
import org.lwjgl.system.MemoryStack;
import org.vadere.util.geometry.mesh.gen.*;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VPoint;
import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.FloatBuffer;
import java.nio.IntBuffer;
import java.util.ArrayList;
import java.util.Collection;
import static org.lwjgl.opencl.CL10.*;
import static org.lwjgl.opencl.CL10.clEnqueueNDRangeKernel;
import static org.lwjgl.system.MemoryUtil.NULL;
import static org.lwjgl.system.MemoryUtil.memUTF8;
/**
* Created by bzoennchen on 08.09.17.
*/
public class CLDistMesh {
// CL ids
private MemoryStack stack;
private long clPlatform;
private long clDevice;
private long clContext;
private long clQueue;
private long clProgram;
// CL kernel ids
private long clKernelForces;
private long clKernelMove;
private long clKernelLengths;
private long clKernelPartialSF;
private long clKernelCompleteSF;
// error code buffer
private IntBuffer errcode_ret;
// CL callbacks
private CLContextCallback contextCB;
private CLProgramCallback programCB;
private static Logger log = LogManager.getLogger(CLDistMesh.class);
// data on the host
private FloatBuffer v;
private FloatBuffer scalingFactor;
private IntBuffer e;
private IntBuffer t;
private float delta = 0.5f;
// addresses to memory on the GPU
private long clVertices;
private long clEdges;
private long clTriangles;
private long clForces;
private long clLengths;
private long clqLengths;
private long clPartialSum;
private long clScalingFactor;
private ArrayList<Long> clSizes = new ArrayList<>();
// size
private int n;
private int numberOfVertices;
private int numberOfEdges;
private int numberOfFaces;
private long maxGroupSize;
private long maxComputeUnits;
private PointerBuffer clGlobalWorkSizeEdges;
private PointerBuffer clGlobalWorkSizeVertices;
public CLDistMesh(@NotNull AMesh<? extends IPoint> mesh) {
this.stack = MemoryStack.stackPush();
this.v = CLGatherer.getVertices(mesh, stack);
this.e = CLGatherer.getEdges(mesh, stack);
this.t = CLGatherer.getFaces(mesh, stack);
this.numberOfVertices = mesh.getNumberOfVertices();
this.numberOfEdges = mesh.getNumberOfEdges();
this.numberOfFaces = mesh.getNumberOfFaces();
}
private void initCallbacks() {
contextCB = CLContextCallback.create((errinfo, private_info, cb, user_data) ->
{
log.warn("[LWJGL] cl_context_callback");
log.warn("\tInfo: " + memUTF8(errinfo));
});
programCB = CLProgramCallback.create((program, user_data) ->
{
log.info("The cl_program [0x"+program+"] was built " + (InfoUtils.getProgramBuildInfoInt(program, clDevice, CL_PROGRAM_BUILD_STATUS) == CL_SUCCESS ? "successfully" : "unsuccessfully"));
});
}
private void initCL() {
// helper for the memory allocation in java
//stack = MemoryStack.stackPush();
errcode_ret = stack.callocInt(1);
IntBuffer numberOfPlatforms = stack.mallocInt(1);
clGetPlatformIDs(null, numberOfPlatforms);
PointerBuffer platformIDs = stack.mallocPointer(numberOfPlatforms.get(0));
clGetPlatformIDs(platformIDs, numberOfPlatforms);
clPlatform = platformIDs.get(0);
IntBuffer numberOfDevices = stack.mallocInt(1);
clGetDeviceIDs(clPlatform, CL_DEVICE_TYPE_GPU, null, numberOfDevices);
PointerBuffer deviceIDs = stack.mallocPointer(numberOfDevices.get(0));
clGetDeviceIDs(clPlatform, CL_DEVICE_TYPE_GPU, deviceIDs, numberOfDevices);
clDevice = deviceIDs.get(0);
printDeviceInfo(clDevice, "CL_DEVICE_NAME", CL_DEVICE_NAME);
PointerBuffer ctxProps = stack.mallocPointer(3);
ctxProps.put(CL_CONTEXT_PLATFORM)
.put(clPlatform)
.put(NULL)
.flip();
clContext = clCreateContext(ctxProps, clDevice, contextCB, NULL, errcode_ret);
InfoUtils.checkCLError(errcode_ret);
clQueue = clCreateCommandQueue(clContext, clDevice, 0, errcode_ret);
PointerBuffer pp = stack.mallocPointer(1);
clGetDeviceInfo(clDevice, CL_DEVICE_MAX_WORK_GROUP_SIZE, pp, null);
maxGroupSize = pp.get(0);
clGetDeviceInfo(clDevice, CL_DEVICE_MAX_COMPUTE_UNITS, pp, null);
maxComputeUnits = pp.get(0);
}
private void buildProgram() {
PointerBuffer strings = BufferUtils.createPointerBuffer(1);
PointerBuffer lengths = BufferUtils.createPointerBuffer(1);
ByteBuffer source;
try {
source = IOUtil.ioResourceToByteBuffer("DistMesh.cl", 4096);
} catch (IOException e) {
throw new RuntimeException(e);
}
strings.put(0, source);
lengths.put(0, source.remaining());
clProgram = clCreateProgramWithSource(clContext, strings, lengths, errcode_ret);
int errcode = clBuildProgram(clProgram, clDevice, "", programCB, NULL);
InfoUtils.checkCLError(errcode);
clKernelLengths = clCreateKernel(clProgram, "computeLengths", errcode_ret);
clKernelPartialSF = clCreateKernel(clProgram, "computePartialSF", errcode_ret);
clKernelCompleteSF = clCreateKernel(clProgram, "computeCompleteSF", errcode_ret);
clKernelForces = clCreateKernel(clProgram, "computeForces", errcode_ret);
clKernelMove = clCreateKernel(clProgram, "moveVertices", errcode_ret);
}
private void createMemory() {
clVertices = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, v, errcode_ret);
clEdges = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, e, errcode_ret);
clTriangles = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, t, errcode_ret);
clForces = clCreateBuffer(clContext, CL_MEM_READ_WRITE, 4 * 2 * numberOfVertices, errcode_ret);
clLengths = clCreateBuffer(clContext, CL_MEM_READ_WRITE, 4 * 2 * numberOfEdges, errcode_ret);
clqLengths = clCreateBuffer(clContext, CL_MEM_READ_WRITE, 4 * 2 * numberOfEdges, errcode_ret);
clScalingFactor = clCreateBuffer(clContext, CL_MEM_READ_WRITE, 4, errcode_ret);
//TODO: smaller
clPartialSum = clCreateBuffer(clContext, CL_MEM_READ_WRITE, 4 * 2 * numberOfEdges, errcode_ret);
}
private void initialKernelArgs() {
clSetKernelArg1p(clKernelLengths, 0, clVertices);
clSetKernelArg1p(clKernelLengths, 1, clEdges);
clSetKernelArg1p(clKernelLengths, 2, clLengths);
clSetKernelArg1p(clKernelLengths, 3, clqLengths);
clSetKernelArg1p(clKernelPartialSF, 1, clqLengths);
clSetKernelArg(clKernelPartialSF, 2, 4 * 2 * numberOfEdges);
clSetKernelArg1p(clKernelPartialSF, 3, clPartialSum);
clSetKernelArg1p(clKernelCompleteSF, 1, clqLengths);
clSetKernelArg(clKernelCompleteSF, 2, 4 * 2 * numberOfEdges);
clSetKernelArg1p(clKernelCompleteSF, 3, clScalingFactor);
clSetKernelArg1p(clKernelForces, 0, clVertices);
clSetKernelArg1p(clKernelForces, 1, clEdges);
clSetKernelArg1p(clKernelForces, 2, clLengths);
clSetKernelArg1p(clKernelForces, 3, clScalingFactor);
clSetKernelArg1p(clKernelForces, 4, clForces);
clSetKernelArg1p(clKernelMove, 0, clVertices);
clSetKernelArg1p(clKernelMove, 1, clForces);
clSetKernelArg1f(clKernelMove, 2, delta);
clGlobalWorkSizeEdges = BufferUtils.createPointerBuffer(1);
clGlobalWorkSizeVertices = BufferUtils.createPointerBuffer(1);
clGlobalWorkSizeEdges.put(0, numberOfEdges);
clGlobalWorkSizeVertices.put(0, numberOfVertices);
}
public void step() {
/*
* DistMesh-Loop
* 1. compute scaling factor
* 2. compute forces;
* 3. update vertices;
*
*/
int potGrpSize = (int)Math.ceil(MathUtils.log(2, maxGroupSize));
int potWorkLoad = (int)Math.ceil(MathUtils.log(2, numberOfEdges));
int rounds = potWorkLoad - potGrpSize;
clEnqueueNDRangeKernel(clQueue, clKernelLengths, 1, null, clGlobalWorkSizeEdges, null, null, null);
int size = numberOfEdges;
int globalWorkSize = (int)(maxGroupSize * maxComputeUnits);
int localWorkSize = (int)maxGroupSize;
PointerBuffer clGlobalWorkSize = BufferUtils.createPointerBuffer(1);
PointerBuffer clLocalWorkSize = BufferUtils.createPointerBuffer(1);
clGlobalWorkSize.put(0, globalWorkSize);
clLocalWorkSize.put(0, localWorkSize);
int index = 0;
while(maxGroupSize < size) {
IntBuffer clSize = BufferUtils.createIntBuffer(1);
clSize.put(0, size);
clSetKernelArg1i(clKernelPartialSF, 0, size);
clEnqueueNDRangeKernel(clQueue, clKernelPartialSF, 1, null, clGlobalWorkSize, clLocalWorkSize, null, null);
if(index == 0) {
size = (int)maxComputeUnits;