Commit 9b5f4ae8 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

starting to try to improve the distmesh algorithm

parent 9b5c6caa
......@@ -229,17 +229,24 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
if(!finalized) {
// we have to use other halfedges than he1 and he2 since they might be deleted
// if we deleteBoundaryFace he0!
List<F> faces1 = mesh.getFaces(p0);
faces1.removeIf(f -> mesh.isBoundary(f));
faces1.forEach(f -> removeFace(f, true));
List<F> faces2 = mesh.getFaces(p1);
faces2.removeIf(f -> mesh.isDestroyed(f) || mesh.isBoundary(f));
faces2.forEach(f -> removeFace(f, true));
if(!mesh.isDestroyed(p0)) {
List<F> faces1 = mesh.getFaces(p0);
faces1.removeIf(f -> mesh.isBoundary(f));
faces1.forEach(f -> removeFace(f, true));
}
List<F> faces3 = mesh.getFaces(p2);
faces3.removeIf(f -> mesh.isDestroyed(f) || mesh.isBoundary(f));
faces3.forEach(f -> removeFace(f, true));
if(!mesh.isDestroyed(p1)) {
List<F> faces2 = mesh.getFaces(p1);
faces2.removeIf(f -> mesh.isDestroyed(f) || mesh.isBoundary(f));
faces2.forEach(f -> removeFace(f, true));
}
if(!mesh.isDestroyed(p2)) {
List<F> faces3 = mesh.getFaces(p2);
faces3.removeIf(f -> mesh.isDestroyed(f) || mesh.isBoundary(f));
faces3.forEach(f -> removeFace(f, true));
}
finalized = true;
}
......
......@@ -143,7 +143,7 @@ public class PHalfEdge<P extends IPoint> implements IHalfEdge<P> {
* one half-edge part of face and ending at end.
*/
@Override
/*@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
......@@ -159,5 +159,5 @@ public class PHalfEdge<P extends IPoint> implements IHalfEdge<P> {
int result = end.hashCode();
result = 31 * result + (face != null ? face.hashCode() : 0);
return result;
}
}*/
}
......@@ -135,11 +135,23 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
@Override
public void setEdge(@NotNull PVertex<P> vertex, @NotNull PHalfEdge<P> edge) {
assert edge.getEnd().equals(vertex);
if(!edge.getEnd().equals(vertex)) {
System.out.println("end of the edge is not equals to the vertex.");
}
vertex.setEdge(edge);
}
@Override
public void setVertex(@NotNull PHalfEdge<P> halfEdge, @NotNull PVertex<P> vertex) {
/*if(halfEdge.getEnd().getEdge() == halfEdge) {
System.out.println("error44");
}
if(!vertex.getEdge().getEnd().equals(vertex)) {
System.out.println("error2");
}*/
halfEdge.setEnd(vertex);
}
......@@ -266,6 +278,9 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
return edges.stream();
}
@Override
public Stream<PVertex<P>> streamVertices() { return vertices.stream(); };
@Override
public Iterable<PHalfEdge<P>> getEdgeIt() {
return () -> edges.iterator();
......@@ -273,6 +288,6 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
@Override
public List<PFace<P>> getFaces() {
return streamFaces().filter(face -> !face.isDestroyed()).collect(Collectors.toList());
return streamFaces().filter(face -> !face.isBorder()).filter(face -> !face.isDestroyed()).collect(Collectors.toList());
}
}
......@@ -13,12 +13,14 @@ import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.geometry.shapes.VShape;
import org.vadere.util.geometry.shapes.VTriangle;
import org.vadere.util.triangulation.adaptive.IDistanceFunction;
import org.vadere.util.triangulation.adaptive.IEdgeLengthFunction;
import java.util.Collection;
import java.util.Collections;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.Set;
......@@ -37,12 +39,15 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
private Set<P> points;
private IMesh<P, V, E, F> mesh;
private static final Logger logger = LogManager.getLogger(UniformRefinementTriangulation.class);
private final IDistanceFunction distFunc;
public UniformRefinementTriangulation(
final ITriangulation<P, V, E, F> triangulation,
final VRectangle bound,
final Collection<VShape> boundary,
final IEdgeLengthFunction lenFunc) {
final IEdgeLengthFunction lenFunc,
final IDistanceFunction distFunc) {
this.distFunc = distFunc;
this.triangulation = triangulation;
this.mesh = triangulation.getMesh();
this.boundary = boundary;
......@@ -79,10 +84,31 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
}
private void removeTrianglesOutsideBBox() {
//TODO:
boolean removedSome = true;
while (removedSome) {
removedSome = false;
List<F> candidates = mesh.getFaces(mesh.getBoundary());
for(F face : candidates) {
if(!mesh.isDestroyed(face) && mesh.streamVertices(face).anyMatch(v -> !bbox.contains(v))) {
triangulation.removeFace(face, true);
removedSome = true;
}
}
}
}
private void removeTrianglesInsideObstacles() {
List<F> faces = triangulation.getMesh().getFaces();
for(F face : faces) {
if(!triangulation.getMesh().isDestroyed(face) && distFunc.apply(triangulation.getMesh().toTriangle(face).midPoint()) > 0) {
triangulation.removeFace(face, true);
}
}
}
/*private void removeTrianglesInsideObstacles() {
for(VShape shape : boundary) {
// 1. find a triangle inside the boundary
......@@ -110,7 +136,7 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
logger.warn("no face found");
}
}
}
}*/
private boolean intersectShape(final VLine line, final VShape shape) {
return shape.intersects(line) || shape.contains(line.getP1()) || shape.contains(line.getP2());
......
......@@ -21,6 +21,6 @@ public class VPUniformRefinement extends UniformRefinementTriangulation<VPoint,
final VRectangle bound,
final Collection<VShape> boundary,
final IEdgeLengthFunction lenFunc) {
super(triangulation, bound, boundary, lenFunc);
super(triangulation, bound, boundary, lenFunc, p -> -1.0);
}
}
package org.vadere.util.geometry.mesh.inter;
import org.apache.commons.collections.IteratorUtils;
import org.apache.commons.collections.ListUtils;
import org.apache.commons.lang3.tuple.Triple;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.GeometryUtils;
......@@ -184,6 +185,8 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
Stream<E> streamEdges();
Stream<V> streamVertices();
default VPolygon toPolygon(F face) {
Path2D path2D = new Path2D.Double();
E edge = getEdge(face);
......@@ -308,6 +311,10 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
*/
default Iterable<F> getFaceIt(F face) { return () -> new SurroundingFaceIterator<>(this, face);}
default Iterable<F> getFaceIt(V vertex) { return () -> new AdjacentFaceIterator(this, getEdge(vertex));}
default List<F> getFaces(F face) { return IteratorUtils.toList(new SurroundingFaceIterator<>(this, face)); }
/**
* Returns a Stream consisting of all surrounding faces of the face.
*
......
......@@ -291,6 +291,10 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
if(getMesh().isBoundary(getMesh().getTwin(edge))) {
delEdges.add(edge);
}
else {
// update the edge of the boundary since it might be deleted!
getMesh().setEdge(boundary, edge);
}
vertices.add(getMesh().getVertex(edge));
}
......@@ -310,6 +314,9 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
next1 = getMesh().getNext(h1);
prev1 = getMesh().getPrev(h1);
boolean isolated0 = isSimpleConnected(v0);
boolean isolated1 = isSimpleConnected(v1);
// adjust next and prev half-edges
getMesh().setNext(prev0, next1);
getMesh().setNext(prev1, next0);
......@@ -317,11 +324,11 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
//boolean isolated0 = getMesh().getNext(prev1).equals(getMesh().getTwin(prev1));
//boolean isolated1 = getMesh().getNext(prev0).equals(getMesh().getTwin(prev0));
boolean isolated0 = getMesh().getTwin(h0).equals(getMesh().getNext(h0)) || getMesh().getTwin(h0).equals(getMesh().getPrev(h0));
boolean isolated1 = getMesh().getTwin(h1).equals(getMesh().getNext(h1)) || getMesh().getTwin(h1).equals(getMesh().getPrev(h1));
//boolean isolated0 = getMesh().getTwin(h0) == getMesh().getNext(h0) || getMesh().getTwin(h0) == getMesh().getPrev(h0);
//boolean isolated1 = getMesh().getTwin(h1) == getMesh().getNext(h1) || getMesh().getTwin(h1) == getMesh().getPrev(h1);
// adjust vertices
if(getMesh().getEdge(v0).equals(h0) && !isolated0) {
if(getMesh().getEdge(v0) == h0 && !isolated0) {
getMesh().setEdge(v0, prev1);
}
......@@ -329,7 +336,7 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
getMesh().destroyVertex(v0);
}
if(getMesh().getEdge(v1).equals(h1) && !isolated1) {
if(getMesh().getEdge(v1) == h1 && !isolated1) {
getMesh().setEdge(v1, prev0);
}
......@@ -342,12 +349,31 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
getMesh().destroyEdge(h1);
// TODO: do we need this?
vertices.stream().filter(getMesh()::isAlive).forEach(this::adjustVertex);
//vertices.stream().filter(getMesh()::isAlive).forEach(this::adjustVertex);
}
}
getMesh().destroyFace(face);
}
/**
* Tests whether the vertex has degree smaller or equals 2.
* If an edge gets deleted and the vertex is simple connected
* the vertex becomes isolated.
*
* @param vertex the vertex
* @return true if the vertex has degree smaller or equals 2, false otherwise.
*/
default boolean isSimpleConnected(@NotNull final V vertex) {
if(getMesh().isDestroyed(vertex)) {
return true;
}
// test if degree of the vertex is <= 2
E edge0 = getMesh().getEdge(vertex);
E edge1 = getMesh().getTwin(getMesh().getNext(edge0));
E edge2 = getMesh().getTwin(getMesh().getNext(edge1));
return edge0 == edge1 || edge0 == edge2;
}
/**
* Returns a half-edge such that it is part of face1 and the twin of this half-edge
* is part of face2.
......
......@@ -60,6 +60,10 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
*/
boolean isIllegal(E edge, V p);
default boolean isIllegal(E edge) {
return isIllegal(edge, getMesh().getVertex(getMesh().getNext(edge)));
}
/**
* Splits the half-edge at point p, preserving a valid triangulation.
* Assumption: p is located on the edge!
......@@ -111,7 +115,17 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
mesh.setEdge(v, t1);
mesh.setTwin(e1, t1);
/*
* These two operations are strongly connected.
* Before these operations the vertex of o0 is v2.
* If the edge of v2 is equal to o0, the edge becomes
* invalid after calling mesh.setVertex(o0, v);
*/
mesh.setVertex(o0, v);
if(mesh.getEdge(v2).equals(o0)) {
mesh.setEdge(v2, e1);
}
if(!mesh.isBoundary(h0)) {
F f1 = mesh.createFace();
......@@ -249,6 +263,11 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
}
// TODO: maybe without if, just do it? its faster?
assert mesh.getVertex(b2) == va0;
assert mesh.getVertex(a2) == vb0;
if(mesh.getEdge(va0).equals(a0)) {
mesh.setEdge(va0, b2);
}
......
package org.vadere.util.geometry.mesh.iterators;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.mesh.inter.IHalfEdge;
import org.vadere.util.geometry.mesh.inter.IMesh;
......@@ -20,10 +22,12 @@ import java.util.Iterator;
*/
public class IncidentEdgeIterator<P extends IPoint, V extends IVertex<P>, E extends IHalfEdge<P>, F extends IFace<P>> implements Iterator<E> {
private static Logger log = LogManager.getLogger(IncidentEdgeIterator.class);
private IMesh<P, V, E, F> mesh;
private E current;
private E edge;
private boolean first;
int count = 0;
public IncidentEdgeIterator(final IMesh<P, V, E, F> mesh, final E edge) {
this.mesh = mesh;
......@@ -42,6 +46,9 @@ public class IncidentEdgeIterator<P extends IPoint, V extends IVertex<P>, E exte
E result = current;
current = mesh.getNext(mesh.getTwin(result));
first = false;
count++;
//log.info(count);
return result;
}
}
......@@ -15,11 +15,11 @@ public class VPolygon extends Path2D.Double implements VShape {
public VPolygon(Path2D.Double path) {
this.reset();
this.append(path, false);
this.closePath();
/*if (!path.getBounds().isEmpty()) {
if (!path.getBounds().isEmpty()) {
this.append(path, false);
this.closePath();
}
}*/
}
public VPolygon() {
......
package org.vadere.util.geometry.shapes;
import java.awt.geom.PathIterator;
import java.util.Arrays;
import java.util.stream.Stream;
......@@ -137,6 +138,22 @@ public class VTriangle extends VPolygon {
return orthocenter;
}
public VPoint closestPoint(final IPoint point) {
VPoint currentClosest = null;
double currentMinDistance = java.lang.Double.MAX_VALUE;
for(VLine line : lines) {
VPoint p = GeometryUtils.closestToSegment(line, point);
if(p.distance(point) < currentMinDistance) {
currentMinDistance = p.distance(point);
currentClosest = p;
}
}
return currentClosest;
}
public VPoint getCircumcenter(){
if(center == null) {
center = GeometryUtils.getCircumcenter(p1, p2, p3);
......
......@@ -159,7 +159,10 @@ public class PSDistmesh {
double lenDiff = Math.max(desiredLen - len, 0);
//double lenDiff = desiredLen - len;
//double force = desiredLen - len;
VPoint directedForce = new VPoint((lenDiff / len) * (line.getX1() - line.getX2()), (lenDiff / len) * (line.getY1() - line.getY2()));
IPoint forceDirection = new VPoint(line.getP1()).subtract(new VPoint(line.getP2())).norm();
VPoint directedForce = new VPoint(forceDirection.scalarMultiply((lenDiff / len)));
//VPoint directedForce = new VPoint((lenDiff / len) * (line.getX1() - line.getX2()), (lenDiff / len) * (line.getY1() - line.getY2()));
//VPoint directedForce = new VPoint(3, 3);
line.setVelocity(directedForce);
//line.setVelocity(line.getVelocity().add(new VPoint(Math.random(), Math.random())));
......
package org.vadere.util.triangulation.adaptive;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.vadere.util.geometry.mesh.gen.PFace;
import org.vadere.util.geometry.mesh.gen.PHalfEdge;
import org.vadere.util.geometry.mesh.gen.PVertex;
import org.vadere.util.geometry.mesh.inter.IMesh;
import org.vadere.util.geometry.mesh.inter.IPointLocator;
import org.vadere.util.geometry.mesh.inter.ITriangulation;
import org.vadere.util.geometry.mesh.gen.UniformRefinementTriangulation;
......@@ -15,16 +17,28 @@ import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.geometry.shapes.VShape;
import org.vadere.util.geometry.shapes.VTriangle;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Comparator;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.PriorityQueue;
import java.util.Set;
import java.util.function.BiFunction;
import java.util.function.Function;
import java.util.stream.Collectors;
import java.util.stream.DoubleStream;
/**
* @author Benedikt Zoennchen
*/
public class PSMeshing {
private BiFunction<Double, Double, Double> f = (l, l_0) -> Math.max(l_0 - l, 0.0);
private static final Logger log = LogManager.getLogger(PSMeshing.class);
private IDistanceFunction distanceFunc;
......@@ -33,9 +47,14 @@ public class PSMeshing {
private Collection<? extends VShape> obstacleShapes;
private VRectangle bound;
private double scalingFactor;
private PriorityQueue<PHalfEdge<MeshPoint>> heap;
private PriorityQueue<Pair<PFace<MeshPoint>, Double>> heap;
private double initialEdgeLen = 10.0;
private double deps;
private double sumOfqDesiredEdgeLength;
private double sumOfqLengths;
private final int retrangulationRate = 10;
private int retrangulationCounter = 0;
public PSMeshing(
final IDistanceFunction distanceFunc,
......@@ -44,11 +63,11 @@ public class PSMeshing {
final Collection<? extends VShape> obstacleShapes) {
this.bound = bound;
this.heap = new PriorityQueue<>(new EdgeComperator());
this.heap = new PriorityQueue<>(new FaceComparator());
this.distanceFunc = distanceFunc;
this.edgeLengthFunc = edgeLengthFunc;
this.obstacleShapes = obstacleShapes;
this.deps = 1.4901e-8 * 1.0;
this.deps = 1.4901e-8 * initialEdgeLen;
this.triangulation = ITriangulation.createPTriangulation(IPointLocator.Type.DELAUNAY_HIERARCHY, bound, (x, y) -> new MeshPoint(x, y, false));
}
......@@ -57,132 +76,296 @@ public class PSMeshing {
*/
public void initialize() {
log.info("##### (start) compute a uniform refined triangulation #####");
UniformRefinementTriangulation uniformRefinementTriangulation = new UniformRefinementTriangulation(triangulation, bound, obstacleShapes, edgeLengthFunc);
UniformRefinementTriangulation uniformRefinementTriangulation = new UniformRefinementTriangulation(triangulation, bound, obstacleShapes, p -> 0.7, distanceFunc);
uniformRefinementTriangulation.compute();
for(PHalfEdge<MeshPoint> edge : triangulation.getMesh().getEdgeIt()) {
MeshPoint point = triangulation.getMesh().getPoint(edge);
point.setVelocity(computeDirectedForces(edge));
heap.add(edge);
//computeScalingFactor();
for(PFace<MeshPoint> face : getMesh().getFaces()) {
heap.add(Pair.of(face, faceToQuality(face)));
}
triangulation.getMesh().streamEdges().forEach(halfEdge -> heap.add(halfEdge));
scalingFactor = computeScalingFactor();
log.info("##### (end) compute a uniform refined triangulation #####");
}
public void improve() {
PHalfEdge<MeshPoint> halfEdge = heap.poll();
while (movePoint(halfEdge) < Parameters.TOL);
computeScalingFactor();
int steps = heap.size();
List<PFace<MeshPoint>> updated = new ArrayList<>(steps);
retrangulationCounter++;
for(int i = 0; i < steps / 100; i++) {
Pair<PFace<MeshPoint>, Double> pair = heap.poll();
PFace<MeshPoint> face = pair.getKey();
for(PHalfEdge<MeshPoint> edge : getMesh().getEdgeIt(face)) {
MeshPoint p1 = getMesh().getPoint(edge);
MeshPoint p2 = getMesh().getPoint(getMesh().getPrev(edge));
VLine line = new VLine(p2.toVPoint(), p1.toVPoint());
double len = line.length();
double desiredLen = edgeLengthFunc.apply(line.midPoint()) * Parameters.FSCALE * scalingFactor;
double lenDiff = Math.max(desiredLen - len, 0);
VPoint forceDirection = p1.toVPoint().subtract(p2.toVPoint()).norm();
VPoint partlyForce = forceDirection.scalarMultiply((lenDiff / len));
updatePoint(p1, partlyForce);
updatePoint(p2, partlyForce.scalarMultiply(-1.0));
}
updated.add(face);
//heap.add(Pair.of(face, faceToQuality(face)));
}
if(retrangulationCounter >= retrangulationRate) {
retriangulate();
for(PFace<MeshPoint> face : getMesh().getFaces()) {
heap.add(Pair.of(face, faceToQuality(face)));
}
retrangulationCounter = 0;
}
else {
for(PFace<MeshPoint> face : updated){
heap.add(Pair.of(face, faceToQuality(face)));
}
}
//retriangulate();
}
public double step() {
PHalfEdge<MeshPoint> halfEdge = heap.poll();
return movePoint(halfEdge);
computeScalingFactor();
computeForces();
//updateVertices();
retriangulate();
return 10;
}
public Collection<VTriangle> getTriangles() {
return triangulation.streamTriangles().collect(Collectors.toList());
}
public void computeForces() {
for(PHalfEdge<MeshPoint> edge : getMesh().getEdgeIt()) {
MeshPoint p1 = getMesh().getPoint(edge);
private double movePoint(final PHalfEdge<MeshPoint> halfEdge) {
//recompute!
VPoint directedForce = computeDirectedForces(halfEdge);
MeshPoint point = triangulation.getMesh().getPoint(halfEdge);
IPoint movement = directedForce.scalarMultiply(Parameters.DELTAT);
for(PHalfEdge<MeshPoint> neighbour : getMesh().getIncidentEdges(edge)) {