11.3.2021, 9:00 - 11:00: Due to updates GitLab may be unavailable for some minutes between 09:00 and 11:00.

Commit f26b2c4a authored by Benedikt Zoennchen's avatar Benedikt Zoennchen

DistMesh on the GPU

parent 7d728459
......@@ -45,6 +45,20 @@
<target>1.8</target>
</configuration>
</plugin>
<plugin>
<groupId>com.googlecode.mavennatives</groupId>
<artifactId>maven-nativedependencies-plugin</artifactId>
<version>0.0.5</version>
<executions>
<execution>
<id>unpacknatives</id>
<phase>generate-resources</phase>
<goals>
<goal>copy</goal>
</goals>
</execution>
</executions>
</plugin>
</plugins>
</build>
......
kernel void computeForces(
__global float2* vertices,
__global int4* edges,
__global float2* lengths,
__global float* scalingFactor,
__global float2* forces)
{
int i = get_global_id(0);
float2 p1 = vertices[edges[i].s0];
float2 p2 = vertices[edges[i].s1];
float2 v = normalize(p1-p2);
// F= ...
float len = lengths[i].s0;
float desiredLen = lengths[i].s1 * (*scalingFactor) * 1.2f;
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
//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);
}
// computation of the scaling factor:
kernel void computeLengths(
__global float2* vertices,
__global int4* edges,
__global float2* lengths,
__global float2* qLengths)
{
int i = get_global_id(0);
float2 p1 = vertices[edges[i].s0];
float2 p2 = vertices[edges[i].s1];
float2 v = p1-p2;
//TODO: desiredLenfunction required
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);
}
// 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);
if(gid < size){
int global_index = gid;
float2 accumulator = (float2)(0, 0);
// Loop sequentially over chunks of input vector
while (global_index < size) {
float2 element = qlengths[global_index];
accumulator += element;
global_index += get_global_size(0);
}
int lid = get_local_id(0);
int group_size = get_local_size(0);
float2 len = qlengths[gid];
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) {
partialSums[lid] += partialSums[lid + i];
}
barrier(CLK_LOCAL_MEM_FENCE);
}
if(lid == 0){
output[get_group_id(0)] = partialSums[0];
}
}
}
// kernel for 1 work-group
kernel void computeCompleteSF(__const int size, __global float2* qlengths, __local float2* partialSums, __global float* 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 = sqrt(partialSums[0].s0 / partialSums[0].s1);
}
}
}
\ No newline at end of file
#ifdef DOUBLE_FP
#ifdef AMD_FP
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#else
#ifndef CL_VERSION_1_2
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#endif
#endif
#define varfloat double
#define _255 255.0
#else
#define varfloat float
#define _255 255.0f
#endif
#ifdef USE_TEXTURE
#define OUTPUT_TYPE __write_only image2d_t
#else
#define OUTPUT_TYPE global uint *
#endif
/**
* For a description of this algorithm please refer to
* http://en.wikipedia.org/wiki/Mandelbrot_set
* @author Michael Bien
*/
kernel void mandelbrot(
const int width, const int height,
const varfloat x0, const varfloat y0,
const varfloat rangeX, const varfloat rangeY,
OUTPUT_TYPE output, global uint *colorMap,
const int colorMapSize, const int maxIterations
) {
unsigned int ix = get_global_id(0);
unsigned int iy = get_global_id(1);
varfloat r = x0 + ix * rangeX / width;
varfloat i = y0 + iy * rangeY / height;
varfloat x = 0;
varfloat y = 0;
varfloat magnitudeSquared = 0;
int iteration = 0;
while ( magnitudeSquared < 4 && iteration < maxIterations ) {
varfloat x2 = x*x;
varfloat y2 = y*y;
y = 2 * x * y + i;
x = x2 - y2 + r;
magnitudeSquared = x2+y2;
iteration++;
}
if ( iteration == maxIterations ) {
#ifdef USE_TEXTURE
write_imagef(output, (int2)(ix, iy), (float4)0);
#else
output[iy * width + ix] = 0;
#endif
} else {
float alpha = (float)iteration / maxIterations;
int colorIndex = (int)(alpha * colorMapSize);
#ifdef USE_TEXTURE
// We could have changed colorMap to a texture + sampler, but the
// unpacking below has minimal overhead and it's kinda interesting.
uint c = colorMap[colorIndex];
write_imageui(output, (int2)(ix, iy), (uint4)(
(c >> 0) & 0xFF,
(c >> 8) & 0xFF,
(c >> 16) & 0xFF,
0xFF
));
#else
output[iy * width + ix] = colorMap[colorIndex];
#endif
// monochrom
//output[iy * width + ix] = 255*iteration/maxIterations;
}
}
\ No newline at end of file
kernel void add(__global float2* a, __global float2* b, __global float2* c) {
int i = get_global_id(0);
c[i] = a[i] * b[i];
}
\ No newline at end of file
package org.vadere.util.geometry.mesh.gen;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.shapes.IPoint;
/**
* Created by bzoennchen on 06.09.17.
*/
public class AFace<P extends IPoint> implements IFace<P> {
/**
* One of the half-edges bordering this face.
*/
private int edge;
private int id;
private boolean border;
private boolean destroyed = false;
/**
* Default constructor. To construct a face where you have already some half-edges
* bordering this face.
*
* @param edge one of the half-edges bordering this face.
*/
AFace(@NotNull final int id, @NotNull final int edge) {
this(id, edge, false);
}
AFace(@NotNull final int id,@NotNull final int edge, boolean border) {
this.border = border;
this.edge = edge;
this.id = id;
}
/**
* This constructor can be used for constructing a new face without having
* constructed the bordering half-edges jet.
*/
AFace(@NotNull final int id, @NotNull boolean border) {
this.border = border;
this.edge = -1;
this.id = id;
}
AFace() {
this.border = false;
}
boolean isBorder() {
return border;
}
void destroy() {
setEdge(-1);
destroyed = true;
}
/**
* Sets one of the half-edges bordering this face.
*
* @param edge half-edge bordering this face
*/
void setEdge(final int edge) {
this.edge = edge;
}
int getEdge() {
return edge;
}
int getId() {
return id;
}
boolean isDestroyed() {
return destroyed;
}
}
package org.vadere.util.geometry.mesh.gen;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.mesh.inter.IHalfEdge;
import org.vadere.util.geometry.shapes.IPoint;
/**
* Created by bzoennchen on 06.09.17.
*/
public class AHalfEdge<P extends IPoint> implements IHalfEdge<P> {
private int id;
/**
* point at the end of the half edge.
*/
private int end;
/**
* next half-edge around the face.
*/
private int next;
/**
* previous half-edge around the face.
*/
private int previous;
private int twin;
/**
* the face the half-edge borders.
*/
private int face;
private boolean destroyed;
protected AHalfEdge(@NotNull final int id, @NotNull final int end, @NotNull final int face) {
this.id = id;
this.end = end;
this.face = face;
this.destroyed = false;
}
protected AHalfEdge(@NotNull final int id, @NotNull final int end) {
this.id = id;
this.end = end;
this.face = -1;
this.destroyed = false;
}
int getFace() {
return face;
}
void setFace(final int face) {
this.face = face;
}
int getEnd() {
return end;
}
boolean hasNext() {
return next != -1;
}
int getNext() {
return next;
}
int getPrevious() {
return previous;
}
int getId() {
return id;
}
int getTwin() { return twin; }
/**
* removes the cyclic pointer structure such that the GC can deleteBoundaryFace these objects.
*/
void destroy() {
setNext(-1);
setPrevious(-1);
setFace(-1);
destroyed = true;
}
boolean isValid() {
return next != -1 && previous != -1 && face != -1;
}
void setPrevious(final int previous) {
this.previous = previous;
}
void setNext(final int next) {
this.next = next;
}
void setTwin(final int twin) {
this.twin = twin;
}
void setEnd(final int end) {
this.end = end;
}
public boolean isDestroyed() {
return destroyed;
}
@Override
public String toString() {
return ""+id;
}
/*
* A half-edge is defined by its end vertex and its face. In a geometry there can not be more than
* one half-edge part of face and ending at end.
*/
/*@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
PHalfEdge<?> halfEdge = (PHalfEdge<?>) o;
if (!end.equals(halfEdge.end)) return false;
return face != null ? face.equals(halfEdge.face) : halfEdge.face == null;
}
@Override
public int hashCode() {
int result = end.hashCode();
result = 31 * result + (face != null ? face.hashCode() : 0);
return result;
}*/
}
package org.vadere.util.geometry.mesh.gen;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.mesh.inter.IMesh;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.triangulation.IPointConstructor;
import java.util.*;
import java.util.function.Predicate;
import java.util.stream.Collectors;
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 List<AFace<P>> faces;
private boolean elementRemoved;
private int numberOfVertices;
private int numberOfEdges;
private int numberOfFaces;
//private List<PFace<P>> borderFaces;
private AFace<P> boundary;
private List<AHalfEdge<P>> edges;
private IPointConstructor<P> pointConstructor;
private List<AVertex<P>> vertices;
public AMesh(final IPointConstructor<P> pointConstructor) {
this.faces = new ArrayList<>();
//this.borderFaces = new ArrayList<>();
this.edges = new ArrayList<>();
this.vertices = new ArrayList<>();
this.boundary = new AFace<>(-1, true);
//this.faces.add(boundary);
this.pointConstructor = pointConstructor;
this.elementRemoved = false;
this.numberOfFaces = 0;
this.numberOfEdges = 0;
this.numberOfVertices = 0;
}
@Override
public AHalfEdge<P> getNext(@NotNull AHalfEdge<P> halfEdge) {
return edges.get(halfEdge.getNext());
}
@Override
public AHalfEdge<P> getPrev(@NotNull AHalfEdge<P> halfEdge) {
return edges.get(halfEdge.getPrevious());
}
@Override
public AHalfEdge<P> getTwin(@NotNull AHalfEdge<P> halfEdge) {
return edges.get(halfEdge.getTwin());
}
@Override
public AFace<P> getFace(@NotNull AHalfEdge<P> halfEdge) {
int edgeId = halfEdge.getFace();
if(edgeId == -1) {
if(halfEdge.isDestroyed()) {
throw new IllegalArgumentException(halfEdge + " is already destroyed.");
}
return boundary;
}
else {
return faces.get(halfEdge.getFace());
}
}
@Override
public AHalfEdge<P> getEdge(@NotNull AVertex<P> vertex) {
return edges.get(vertex.getEdge());
}
@Override
public AHalfEdge<P> getEdge(@NotNull AFace<P> face) {
return edges.get(face.getEdge());
}
@Override
public P getPoint(@NotNull AHalfEdge<P> halfEdge) {
return getVertex(halfEdge).getPoint();
}
@Override
public AVertex<P> getVertex(@NotNull AHalfEdge<P> halfEdge) {
return vertices.get(halfEdge.getEnd());
}
// the vertex should not be contained in vertices, only the up/down
@Override
public AVertex<P> getDown(@NotNull AVertex<P> vertex) {
return vertices.get(vertex.getDown());
}
// the vertex should not be contained in vertices, only the up/down
@Override
public void setDown(@NotNull AVertex<P> up, @NotNull AVertex<P> down) {
up.setDown(down.getId());
}
@Override
public P getPoint(@NotNull AVertex<P> vertex) {
return vertex.getPoint();
}
@Override
public AFace<P> getFace() {
return faces.stream().filter(f -> !isDestroyed(f)).findAny().get();
}
@Override
public boolean isBoundary(@NotNull AFace<P> face) {
return face == boundary;
}
@Override
public boolean isBoundary(@NotNull AHalfEdge<P> halfEdge) {
return halfEdge.getFace() == boundary.getId();
}
@Override
public boolean isDestroyed(@NotNull AFace<P> face) {
return face.isDestroyed();
}
@Override
public boolean isDestroyed(@NotNull AHalfEdge<P> edge) {
return edge.isDestroyed();
}