Commit 87921029 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen Committed by Benedikt Zoennchen

working on the flipping alg on the gpu

parent ffbb4ce0
......@@ -32,23 +32,22 @@ inline double2 getCircumcenter(double2 p1, double2 p2, double2 p3) {
return (double2) (x, y);
}
inline lock void(int3 ta, int3 tb) {}
// TODO: lots of if/else which is bad for the GPU
kernel void flip(__global double2* vertices,
__global int4* edges,
__global int3* triangles,
__global int* mutexes,
__global int* triMutexes,
__global int* R)
{
int edgeId = get_global_id(0);
// check if edge is illegal
eps = 0.00001;
double eps = 0.00001;
int v0 = edges[edgeId][0];
int v1 = edges[edgeId][1];
int ta = edges[edgeId][2];
int tb = edges[edgeId][3];
int v0 = edges[edgeId].s0;
int v1 = edges[edgeId].s1;
int ta = edges[edgeId].s2;
int tb = edges[edgeId].s3;
// edge is a non-boundary edge
if(tb != -1) {
......@@ -59,38 +58,20 @@ kernel void flip(__global double2* vertices,
int tbU2 = -1;
// search for v2 in ta
if(triangles[ta][0] != v0 && triangles[ta][0] != v1) {
u1 = triangles[ta][0];
} else if(triangles[ta][1] != v0 && triangles[ta][1] != v1) {
u1 = triangles[ta][1];
if(triangles[ta].s0 != v0 && triangles[ta].s0 != v1) {
u1 = triangles[ta].s0;
} else if(triangles[ta].s1 != v0 && triangles[ta].s1 != v1) {
u1 = triangles[ta].s1;
} else {
u1 = triangles[ta][2];
u1 = triangles[ta].s2;
}
if(triangles[tb][0] != v0 && triangles[tb][0] != v1) {
u2 = triangles[tb][0];
} else if(triangles[ta][1] != v0 && triangles[ta][1] != v1) {
u2 = triangles[tb][1];
if(triangles[tb].s0!= v0 && triangles[tb].s0 != v1) {
u2 = triangles[tb].s0;
} else if(triangles[ta].s1 != v0 && triangles[ta].s1 != v1) {
u2 = triangles[tb].s1;
} else {
u2 = triangles[tb][2];
}
if(triangles[tb][0] == v0) {
tbV0 = 0;
} else if(triangles[tb][1] == v0) {
tbV0 = 1;
}
else {
tbV0 = 2;
}
if(triangles[ta][0] == v1) {
taV1 = 0;
} else if(triangles[ta][1] == v1) {
taV1 = 1;
}
else {
taV1 = 2;
u2 = triangles[tb].s2;
}
......@@ -101,55 +82,93 @@ kernel void flip(__global double2* vertices,
// flip
if(rad + eps < length(c-vertices[u2])) {
// 1. lock ta and tb
lock(ta, tb);
triangles[ta][taV1] = u2;
triangles[tb][tbV0] = u1;
edges[0] = u1;
edges[1] = u2;
// save the relation of the changes
R[ta] = tb;
R[tb] = ta;
volatile __global int* addr1 = &triMutexes[ta];
volatile __global int* addr2 = &triMutexes[tb];
// we dont have both locks
if(!LOCK(addr1)) {
if(!LOCK(addr2)) {
if(triangles[ta].s0 == v1) {
triangles[ta].s0 = u2;
} else if(triangles[ta].s1 == v1) {
triangles[ta].s1 = u2;
}
else {
triangles[ta].s2 = u2;
}
if(triangles[tb].s0 == v0) {
triangles[tb].s0 = u1;
} else if(triangles[tb].s1 == v0) {
triangles[tb].s1 = u1;
}
else {
triangles[tb].s2 = u1;
}
edges[0] = u1;
edges[1] = u2;
// save the relation of the changes
R[ta] = tb;
R[tb] = ta;
UNLOCK(addr1);
UNLOCK(addr2);
} else {
UNLOCK(addr1);
}
}
}
}
}
kernel void repair(
__global int4* edges,
kernel void repair(__global int4* edges,
__global int3* triangles,
__global int* R)
{
int edgeId = get_global_id(0);
int v0 = edges[edgeId][0];
int v1 = edges[edgeId][1];
int ta = edges[edgeId][2];
int tb = edges[edgeId][3];
int v0 = edges[edgeId].s0;
int v1 = edges[edgeId].s1;
int ta = edges[edgeId].s2;
int tb = edges[edgeId].s3;
int c = 0;
for(int i = 0; i <=3; i++) {
if(triangles[ta][i] == v0 || triangles[ta][i] == v1){
c++;
}
if(triangles[ta].s0 == v0 || triangles[ta].s0 == v1) {
c++;
}
if(triangles[ta].s1 == v0 || triangles[ta].s1 == v1) {
c++;
}
if(triangles[ta].s2 == v0 || triangles[ta].s2 == v1) {
c++;
}
// invalid
if(c != 2) {
edges[edgeId][2] = R[ta];
edges[edgeId].s2 = R[ta];
}
c = 0;
for(int i = 0; i <=3; i++) {
if(triangles[tb][i] == v0 || triangles[tb][i] == v1){
c++;
}
if(triangles[tb].s0 == v0 || triangles[tb].s0 == v1) {
c++;
}
if(triangles[tb].s1 == v0 || triangles[tb].s1 == v1) {
c++;
}
if(triangles[tb].s2 == v0 || triangles[tb].s2 == v1) {
c++;
}
// invalid
if(c != 2) {
edges[edgeId][2] = R[tb];
edges[edgeId].s3 = R[tb];
}
}
}
kernel void computeForces(
......
package org.vadere.util.geometry.mesh.gen;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.mesh.inter.IMesh;
import org.vadere.util.geometry.shapes.IPoint;
......@@ -14,7 +16,7 @@ import java.util.stream.Stream;
* Created by bzoennchen on 06.09.17.
*/
public class AMesh<P extends IPoint> implements IMesh<P, AVertex<P>, AHalfEdge<P>, AFace<P>> {
private final static Logger log = LogManager.getLogger(AMesh.class);
private List<AFace<P>> faces;
private boolean elementRemoved;
private int numberOfVertices;
......@@ -305,6 +307,7 @@ public class AMesh<P extends IPoint> implements IMesh<P, AVertex<P>, AHalfEdge<P
@Override
public Stream<AHalfEdge<P>> streamEdges() {
log.info(edges.stream().filter(e -> !e.isDestroyed()).count());
return edges.stream().filter(e -> !e.isDestroyed());
}
......
......@@ -24,6 +24,7 @@ public class AVertex<P extends IPoint> implements IVertex<P> {
this.point = point;
this.id = id;
this.lock = new ReentrantLock();
this.destroyed = false;
}
@Override
......
......@@ -15,7 +15,7 @@ import java.util.Collection;
*/
public class CLGatherer {
public static <P extends IPoint> DoubleBuffer getVerticesD(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
public static <P extends IPoint> DoubleBuffer getVerticesD(@NotNull final AMesh<P> mesh) {
Collection<AVertex<P>> vertices = mesh.getVertices();
DoubleBuffer vertexBuffer = MemoryUtil.memAllocDouble(vertices.size()*2);
int index = 0;
......@@ -29,7 +29,7 @@ public class CLGatherer {
return vertexBuffer;
}
public static <P extends IPoint> FloatBuffer getVerticesF(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
public static <P extends IPoint> FloatBuffer getVerticesF(@NotNull final AMesh<P> mesh) {
Collection<AVertex<P>> vertices = mesh.getVertices();
FloatBuffer vertexBuffer = MemoryUtil.memAllocFloat(vertices.size()*2);
int index = 0;
......@@ -43,7 +43,7 @@ public class CLGatherer {
return vertexBuffer;
}
public static <P extends IPoint> IntBuffer getEdges(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
public static <P extends IPoint> IntBuffer getEdges(@NotNull final AMesh<P> mesh) {
Collection<AHalfEdge<P>> edges = mesh.getEdges();
IntBuffer edgeBuffer = MemoryUtil.memAllocInt(edges.size()*4);
int index = 0;
......@@ -60,7 +60,7 @@ public class CLGatherer {
return edgeBuffer;
}
public static <P extends IPoint> IntBuffer getFaces(@NotNull final AMesh<P> mesh, @NotNull MemoryStack stack) {
public static <P extends IPoint> IntBuffer getFaces(@NotNull final AMesh<P> mesh) {
Collection<AFace<P>> faces = mesh.getFaces();
IntBuffer faceBuffer = MemoryUtil.memAllocInt(faces.size()*4);
int index = 0;
......
......@@ -42,7 +42,7 @@ public class CLPSMeshing implements IPSMeshing {
private Object gobalAcessSynchronizer = new Object();
private CLDistMesh clDistMesh;
private CLDistMesh<MeshPoint> clDistMesh;
public CLPSMeshing(
final IDistanceFunction distanceFunc,
......@@ -103,7 +103,8 @@ public class CLPSMeshing implements IPSMeshing {
//log.info(scalingFactor);
clDistMesh.step();
retriangulate();
flipEdges();
//retriangulate();
// get new cooridnates
......@@ -134,6 +135,43 @@ public class CLPSMeshing implements IPSMeshing {
log.info("#points: " + getMesh().getVertices().size());*/
}
public boolean flipEdges() {
refresh();
boolean anyFlip = false;
// Careful, iterate over all half-edges means iterate over each "real" edge twice!
/*for(AHalfEdge<MeshPoint> edge : getMesh().getEdgeIt()) {
if(triangulation.isIllegal(edge)) {
//triangulation.flip(edge);
anyFlip = true;
}
}*/
if(clDistMesh != null) {
clDistMesh.finish();
}
clDistMesh = new CLDistMesh<>((AMesh<MeshPoint>) triangulation.getMesh());
clDistMesh.init();
return anyFlip;
}
public void finish() {
if(clDistMesh != null) {
clDistMesh.finish();
}
triangulation = ITriangulation.createATriangulation(IPointLocator.Type.DELAUNAY_HIERARCHY, clDistMesh.getResult(), (x, y) -> new MeshPoint(x, y, false));
removeTrianglesInsideObstacles();
triangulation.finalize();
}
public void refresh() {
if(clDistMesh != null) {
clDistMesh.refresh();
}
triangulation = ITriangulation.createATriangulation(IPointLocator.Type.DELAUNAY_HIERARCHY, clDistMesh.getResult(), (x, y) -> new MeshPoint(x, y, false));
removeTrianglesInsideObstacles();
triangulation.finalize();
}
public void retriangulate() {
//Set<MeshPoint> points = getMesh().getVertices().stream().map(vertex -> getMesh().getPoint(vertex)).collect(Collectors.toSet());
//removeLowQualityTriangles();
......@@ -144,7 +182,7 @@ public class CLPSMeshing implements IPSMeshing {
if(clDistMesh != null) {
clDistMesh.finish();
}
clDistMesh = new CLDistMesh((AMesh<MeshPoint>) triangulation.getMesh());
clDistMesh = new CLDistMesh<>((AMesh<MeshPoint>) triangulation.getMesh());
clDistMesh.init();
}
......
......@@ -44,7 +44,7 @@ public class TestEnhancedVersion4 extends JFrame {
while (counter <= 1000) {
while (counter <= 400) {
//obscuteTriangles = meshGenerator.getTriangles().stream().filter(tri -> tri.isNonAcute()).count();
//PriorityQueue<PFace<MeshPoint>> priorityQueue = meshGenerator.getQuailties();
//avgQuality = priorityQueue.stream().reduce(0.0, (aDouble, meshPointPFace) -> aDouble + meshGenerator.faceToQuality(meshPointPFace), (d1, d2) -> d1 + d2) / priorityQueue.size();
......@@ -60,11 +60,14 @@ public class TestEnhancedVersion4 extends JFrame {
ms = System.currentTimeMillis() - ms;
time += ms;
System.out.println("Step-Time: " + ms);
meshGenerator.refresh();
distmeshPanel.update();
distmeshPanel.repaint();
counter++;
}
meshGenerator.finish();
distmeshPanel.update();
distmeshPanel.repaint();
System.out.print("overall time: " + time);
//System.out.print("finished:" + meshGenerator.getMesh().getVertices().stream().filter(v -> !meshGenerator.getMesh().isDestroyed(v)).count());
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment