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

Commit 6540e923 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

improve Eikmesh to handle extreme slim obstacles.

parent 60aeb959
Pipeline #138079 passed with stages
in 103 minutes and 58 seconds
# nVertices dimension nAttributes boundaryMarker
20 2 0 0
# vertexId x y
1 0.500000 0.500000
2 61.903197 8.064521
3 38.841269 14.606642
4 40.029782 19.726267
5 0.500000 99.500000
6 44.390559 12.250599
7 43.789522 11.865840
8 45.212895 18.026666
9 40.404106 19.127630
10 31.586493 15.912108
11 31.682452 90.622759
12 39.213217 13.997772
13 38.261243 89.048447
14 62.005890 8.553862
15 31.179013 91.257350
16 99.500000 0.500000
17 45.813087 18.402190
18 32.084388 16.302644
19 99.500000 99.500000
20 38.377608 89.534717
#
# nSegments boundaryMarker
20 0
# lineId vertexId1 vertexId2
1 16 19
2 19 5
3 5 1
4 1 16
5 3 4
6 4 17
7 17 6
8 6 14
9 14 2
10 2 7
11 7 8
12 8 9
13 9 12
14 12 10
15 10 15
16 15 20
17 20 13
18 13 11
19 11 18
20 18 3
#
# nHoles
1
# vertexId x y (of a vertex which lies inside the hole)
1 31.184013 90.334521
\ No newline at end of file
......@@ -4,6 +4,7 @@ import org.jetbrains.annotations.NotNull;
import org.vadere.meshing.mesh.gen.MeshPanel;
import org.vadere.meshing.mesh.impl.PMeshPanel;
import org.vadere.meshing.mesh.impl.PSLG;
import org.vadere.meshing.mesh.triangulation.DistanceFunctionApproxBF;
import org.vadere.meshing.mesh.triangulation.EdgeLengthFunctionApprox;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.mesh.triangulation.triangulator.gen.GenRuppertsTriangulator;
......@@ -11,6 +12,7 @@ import org.vadere.meshing.mesh.triangulation.triangulator.impl.PContrainedDelaun
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PDelaunayTriangulator;
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PRuppertsTriangulator;
import org.vadere.meshing.utils.io.poly.PSLGGenerator;
import org.vadere.util.math.IDistanceFunction;
import java.io.IOException;
import java.io.InputStream;
......@@ -25,12 +27,23 @@ public class BackgroundMeshExamples {
//localFeatureSize("/poly/narrowCorridor.poly");
//localFeatureSize("/poly/bridge.poly");
localFeatureSize("/poly/mf_small_very_simple.poly");
//distance("/poly/mf_small_very_simple.poly");
//distance("/poly/mf_small_very_simple.poly");
}
public static void localFeatureSize(@NotNull final String fileName) throws IOException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream(fileName);
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg);
//edgeLengthFunctionApprox.smooth(0.2);
edgeLengthFunctionApprox.printPython();
}
public static void distance(@NotNull final String fileName) throws IOException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream(fileName);
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
DistanceFunctionApproxBF distFunctionApprox = new DistanceFunctionApproxBF(pslg, IDistanceFunction.create(pslg.getSegmentBound(), pslg.getHoles()));
distFunctionApprox.printPython();
}
}
......@@ -57,9 +57,11 @@ public class EikMeshPlots {
kaiserslautern();
ruppertsAndEikMeshKaiserslautern();*/
bridge();
//bridge();
//roomLFS();
//cornerLFS();
uniformRing(0.3);
}
public static void randomDelaunay() throws IOException {
......@@ -432,9 +434,9 @@ public class EikMeshPlots {
private static void write(final String string, final String filename) throws IOException {
File outputFile = new File("./eikmesh/"+filename+".tex");
try(FileWriter fileWriter = new FileWriter(outputFile)) {
fileWriter.write(string);
}
// try(FileWriter fileWriter = new FileWriter(outputFile)) {
//fileWriter.write(string);
// }
}
private static String toTexDocument(final String tikz) {
......
......@@ -5,6 +5,7 @@ import org.vadere.meshing.mesh.impl.PMeshPanel;
import org.vadere.meshing.mesh.impl.PSLG;
import org.vadere.meshing.mesh.triangulation.EdgeLengthFunctionApprox;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PRuppertsTriangulator;
import org.vadere.meshing.utils.io.poly.PSLGGenerator;
import org.vadere.meshing.utils.io.tex.TexGraphGenerator;
import org.vadere.util.geometry.shapes.VPolygon;
......@@ -27,13 +28,14 @@ public class EikMeshPoly {
//meshPoly("/poly/bridge.poly");
//meshPoly("/poly/room.poly");
//meshPoly("/poly/corner.poly");
//meshPoly("/poly/railing.poly");
}
public static void meshPoly(@NotNull final String fileName) throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream(fileName);
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg);
edgeLengthFunctionApprox.smooth(0.3);
edgeLengthFunctionApprox.smooth(0.4);
edgeLengthFunctionApprox.printPython();
......@@ -41,26 +43,37 @@ public class EikMeshPoly {
VPolygon segmentBound = pslg.getSegmentBound();
IDistanceFunction distanceFunction = IDistanceFunction.create(segmentBound, holes);
var ruppert = new PRuppertsTriangulator(
pslg,
p -> Double.POSITIVE_INFINITY,
0,
true
);
ruppert.generate();
// (3) use EikMesh to improve the mesh
double h0 = 1.0;
var meshImprover = new PEikMesh(
distanceFunction,
edgeLengthFunctionApprox,
h0,
new VRectangle(segmentBound.getBounds2D()),
pslg.getHoles()
pslg.getBoundingBox(),
pslg.getAllPolygons()
);
var meshPanel = new PMeshPanel(meshImprover.getMesh(), 1000, 800);
var meshPanel = new PMeshPanel(meshImprover.getMesh()/*, f -> distanceFunction.apply(meshImprover.getMesh().toTriangle(f).midPoint()) > 0*/, 1000, 800);
meshPanel.display("Combined distance functions " + h0);
meshImprover.improve();
while (!meshImprover.isFinished()) {
meshImprover.improve();
Thread.sleep(20);
synchronized (meshImprover.getMesh()) {
meshImprover.improve();
}
//Thread.sleep(200);
meshPanel.repaint();
}
//meshImprover.generate();
write(toTexDocument(TexGraphGenerator.toTikz(meshImprover.getMesh(), f-> lightBlue, 1.0f)), "mesh.tex");
write(toTexDocument(TexGraphGenerator.toTikz(meshImprover.getMesh(), f-> lightBlue, 1.0f)), "mesh");
System.out.println(meshImprover.getMesh().getNumberOfVertices());
// display the mesh
......
package org.vadere.meshing.examples;
import org.jetbrains.annotations.NotNull;
import org.vadere.meshing.mesh.impl.PMeshPanel;
import org.vadere.meshing.mesh.impl.PSLG;
import org.vadere.meshing.mesh.triangulation.IEdgeLengthFunction;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PRuppertsTriangulator;
import org.vadere.meshing.utils.io.poly.PSLGGenerator;
import org.vadere.meshing.utils.io.tex.TexGraphGenerator;
import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPolygon;
import org.vadere.util.math.IDistanceFunction;
import java.io.IOException;
import java.io.InputStream;
import java.util.Collection;
public class RuppertExamples {
public static void main(String... args) throws InterruptedException, IOException {
rupperts("/poly/mf_small_very_simple.poly");
}
public static void rupperts(@NotNull final String fileName) throws IOException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream(fileName);
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
Collection<VLine> lines = pslg.getAllSegments();
Collection<VPolygon> holes = pslg.getHoles();
VPolygon segmentBound = pslg.getSegmentBound();
//IEdgeLengthFunction h = p -> 0.01 /*+ 0.2*Math.abs(distanceFunction.apply(p))*/;
var ruppert = new PRuppertsTriangulator(
pslg,
p -> Double.POSITIVE_INFINITY,
0,
false
);
//(mesh, f -> false, width, height, colorFunction
PMeshPanel panel = new PMeshPanel(ruppert.getMesh(), f -> ruppert.getMesh().getBooleanData(f, "boundary"),1000, 1000);
panel.display("Ruppert's Algorithm");
while (!ruppert.isFinished()) {
try {
Thread.sleep(200);
} catch (InterruptedException e) {
e.printStackTrace();
}
synchronized (ruppert.getMesh()) {
ruppert.step();
}
panel.repaint();
}
System.out.println(TexGraphGenerator.toTikz(ruppert.getMesh()));
System.out.println("finished");
}
}
......@@ -23,7 +23,8 @@ import javax.swing.*;
*/
public class MeshPanel<V extends IVertex, E extends IHalfEdge, F extends IFace> extends Canvas {
private MeshRenderer<V, E, F> meshRenderer;
private MeshRenderer<V, E, F> meshRenderer;
public static Function defaultFaceColors = f -> new Color(0.8584083044982699f, 0.9134486735870818f, 0.9645674740484429f);
/**
* The width of the canvas.
......@@ -121,7 +122,7 @@ public class MeshPanel<V extends IVertex, E extends IHalfEdge, F extends IFace>
@NotNull final IMesh<V, E, F> mesh,
final double width,
final double height) {
this(mesh, f -> false, width, height, f -> new Color(0.8584083044982699f, 0.9134486735870818f, 0.9645674740484429f), e -> Color.GRAY);
this(mesh, f -> false, width, height, defaultFaceColors, e -> Color.GRAY);
}
@Override
......
......@@ -172,10 +172,11 @@ public class MeshRenderer<V extends IVertex, E extends IHalfEdge, F extends IFac
for(F face : faces) {
VPolygon polygon = mesh.toPolygon(face);
if(faceColorFunction != null) {
if(alertPred.test(face)) {
graphics.setColor(new Color(200, 0, 0));
}
else if(faceColorFunction != null) {
graphics.setColor(faceColorFunction.apply(face));
} else if(alertPred.test(face)) {
graphics.setColor(new Color(100, 0, 0));
} else {
graphics.setColor(Color.GRAY);
}
......
......@@ -368,11 +368,11 @@ public interface IMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
<CF> Optional<CF> getData(@NotNull F face, @NotNull final String name, @NotNull final Class<CF> clazz);
default boolean getBooleanData(@NotNull F face, @NotNull final String name) {
return getData(face, name, Boolean.class).or(null).get();
return getData(face, name, Boolean.class).orElse(false);
}
default double getDoubleData(@NotNull F face, @NotNull final String name) {
return getData(face, name, Double.class).or(null).get();
return getData(face, name, Double.class).orElse(0.0);
}
/**
......@@ -1162,7 +1162,7 @@ public interface IMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
* @return a parallel stream {@link Stream} of (all alive) interior faces
*/
default Stream<F> streamFaces() {
return streamFaces(f -> true);
return streamFaces(f -> !isBoundary(f));
}
/**
......@@ -1316,7 +1316,7 @@ public interface IMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
*/
default VTriangle toTriangle(@NotNull final F face) {
List<V> vertices = getVertices(face);
assert vertices.size() == 3;
assert vertices.size() == 3 : "number of vertices of " + face + " is " + vertices.size();
return new VTriangle(new VPoint(vertices.get(0)), new VPoint(vertices.get(1)), new VPoint(vertices.get(2)));
}
......@@ -1939,6 +1939,8 @@ public interface IMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
return Optional.empty();
}
//default Optional<E> getBoundaryEdge()
/**
* Returns a deep clone of this mesh.
*
......
......@@ -9,6 +9,7 @@ import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.logging.Logger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;
import java.util.Optional;
......@@ -39,6 +40,7 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
*/
Logger log = Logger.getLogger(IPolyConnectivity.class);
/**
* <p>Returns the mesh of this poly-connectivity {@link IPolyConnectivity}.</p>
*
......@@ -553,6 +555,66 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
return Optional.of(currentFace);
}
// TODO: improve performance by remembering faces
/**
* <p>A virus like working algorithm which searches for neighbouring faces by starting at the face until
* the <tt>markCondition</tt> does no longer hold or the maximal dept is reached.
* This requires in the worst case O(n), where n is the number of edges of all involved faces
* (i.e. the face and the merged faces).</p>
*
* <p>Does not changes the connectivity.</p>
*
* @param face the face
* @param markCondition the mark condition.
* @param maxDept the maximum dept / neighbouring distance at which faces can be marked / found
*
* @return the merge result i.e. the resulting face.
*/
default List<F> findFaces(
@NotNull final F face,
@NotNull final Predicate<F> markCondition,
final int maxDept) {
int dept = 0;
if(!markCondition.test(face)) {
return Collections.EMPTY_LIST;
}
Set<F> markedFaces = new HashSet<>();
List<F> result = new ArrayList<>();
result.add(face);
markedFaces.add(face);
List<E> toProcess = getMesh().getEdges(face);
while (!toProcess.isEmpty()) {
dept++;
List<E> newToProcess = new ArrayList<>();
//assert toProcess.isEmpty() || toProcess.stream().anyMatch(e -> getMesh().getFace(e))) : "each edge of both faces is a link to the other face";
for(E edge : toProcess) {
// the face might be destroyed by an operation before
F candidate = getMesh().getTwinFace(edge);
if(!markedFaces.contains(candidate) && markCondition.test(candidate)) {
result.add(candidate);
markedFaces.add(candidate);
for(E e : getMesh().getEdgeIt(candidate)) {
newToProcess.add(e);
}
}
}
if(maxDept > 0 && dept >= maxDept) {
return result;
}
toProcess = newToProcess;
}
return result;
}
/**
* Creates a new hole or extends an existing hole by removing neighbouring faces by a virus algorithm
* which consumes faces as long as the merge condition holds.
......
......@@ -2813,7 +2813,7 @@ public interface ITriConnectivity<V extends IVertex, E extends IHalfEdge, F exte
};
//log.debug(getMesh().streamFaces().filter(f -> !getMesh().isDestroyed(f)).filter(f -> !getMesh().isBoundary(f)).filter(e -> !orientationPredicate.test(e)).count() + " invalid triangles");
return getMesh().streamFaces().filter(f -> !getMesh().isDestroyed(f)).filter(f -> !getMesh().isBoundary(f)).allMatch(orientationPredicate);
return getMesh().streamFaces().filter(f -> !getMesh().isDestroyed(f)).allMatch(orientationPredicate);
}
/**
......
package org.vadere.meshing.mesh.triangulation;
import org.jetbrains.annotations.NotNull;
import org.vadere.meshing.mesh.gen.PFace;
import org.vadere.meshing.mesh.gen.PHalfEdge;
import org.vadere.meshing.mesh.gen.PVertex;
import org.vadere.meshing.mesh.impl.PSLG;
import org.vadere.meshing.mesh.inter.IIncrementalTriangulation;
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PRuppertsTriangulator;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.math.IDistanceFunction;
import org.vadere.util.math.InterpolationUtil;
import java.util.function.Function;
public class DistanceFunctionApproxBF implements IDistanceFunction {
private IIncrementalTriangulation<PVertex, PHalfEdge, PFace> triangulation;
private static final String propName = "distObs";
public DistanceFunctionApproxBF(
@NotNull final PSLG pslg,
@NotNull final Function<IPoint, Double> circumRadiusFunc,
@NotNull final IDistanceFunction exactDistanceFunc) {
//IPointConstructor<DataPoint<Double>> pointConstructor = (x, y) -> new DataPoint<>(x, y);
/**
* Add a bound around so the edge function is also defined outside.
*/
VRectangle bound = GeometryUtils.boundRelative(pslg.getSegmentBound().getPoints(), 0.3);
PSLG boundedPSLG = pslg.conclose(bound);
var ruppertsTriangulator = new PRuppertsTriangulator(boundedPSLG, circumRadiusFunc, 10, false);
triangulation = ruppertsTriangulator.generate();
//TODO: maybe transform into an immutable triangulation / mesh!
triangulation.setCanIllegalPredicate(e -> true);
// compute and set the local feature size
var vertices = triangulation.getMesh().getVertices();
for(var v : vertices) {
double distance = exactDistanceFunc.apply(v);
triangulation.getMesh().setDoubleData(v, propName, distance);
}
}
public DistanceFunctionApproxBF(@NotNull final PSLG pslg, @NotNull final IDistanceFunction exactDistanceFunc) {
this(pslg, p -> Double.POSITIVE_INFINITY, exactDistanceFunc);
//IPointConstructor<DataPoint<Double>> pointConstructor = (x, y) -> new DataPoint<>(x, y);
}
@Override
public Double apply(@NotNull final IPoint p) {
var face = triangulation.locateFace(p.getX(), p.getY()).get();
var mesh = triangulation.getMesh();
if(mesh.isBoundary(face)) {
return Double.POSITIVE_INFINITY;
}
else {
double x[] = new double[3];
double y[] = new double[3];
double z[] = new double[3];
triangulation.getTriPoints(face, x, y, z, propName);
double totalArea = GeometryUtils.areaOfPolygon(x, y);
return InterpolationUtil.barycentricInterpolation(x, y, z, totalArea, p.getX(), p.getY());
}
}
public void printPython() {
System.out.println(triangulation.getMesh().toPythonTriangulation(v -> triangulation.getMesh().getDoubleData(v, propName)));
/*var points = triangulation.getMesh().getPoints();
System.out.print("[");
for(var dataPoint : points) {
System.out.print("["+dataPoint.getX()+","+dataPoint.getY()+"],");
}
System.out.println("]\n\n");
System.out.print("[");
for(var dataPoint : points) {
System.out.print(dataPoint.getData()+",");
}
System.out.println("]");*/
}
}
......@@ -43,17 +43,21 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
// compute and set the local feature size
var vertices = triangulation.getMesh().getVertices();
var mesh = triangulation.getMesh();
for(var v : vertices) {
double minEdgeLen = Double.MAX_VALUE;
for(var u : triangulation.getMesh().getAdjacentVertexIt(v)) {
double len = v.distance(u) * 0.9;
if(len < minEdgeLen) {
minEdgeLen = len;
for(var e : triangulation.getMesh().getEdges(v)) {
if(!mesh.getBooleanData(mesh.getFace(e), "boundary")
|| !mesh.getBooleanData(mesh.getTwinFace(e), "boundary")) {
var u = triangulation.getMesh().getTwinVertex(e);
double len = v.distance(u) * 0.9;
if(len < minEdgeLen) {
minEdgeLen = len;
}
}
}
triangulation.getMesh().setData(v, propName, minEdgeLen);
triangulation.getMesh().setDoubleData(v, propName, minEdgeLen);
}
}
......@@ -67,21 +71,21 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
// smooth the function based such that it is g-Lipschitz
var mesh = triangulation.getMesh();
PriorityQueue<PVertex> heap = new PriorityQueue<>(
Comparator.comparingDouble(v1 -> mesh.getData(v1, propName, Double.class).get())
Comparator.comparingDouble(v1 -> mesh.getDoubleData(v1, propName))
);
heap.addAll(mesh.getVertices());
while (!heap.isEmpty()) {
var v = heap.poll();
double hv = mesh.getData(v, propName, Double.class).get();
double hv = mesh.getDoubleData(v, propName);
for (var u : mesh.getAdjacentVertexIt(v)) {
double hu = mesh.getData(u, propName, Double.class).get();
double hu = mesh.getDoubleData(u, propName);
double min = Math.min(hu, hv + g * v.distance(u));
// update heap
if (min < hu) {
heap.remove(u);
mesh.setData(u, propName, min);
mesh.setDoubleData(u, propName, min);
heap.add(u);
}
}
......@@ -110,7 +114,7 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
}
public void printPython() {
System.out.println(triangulation.getMesh().toPythonTriangulation(v -> triangulation.getMesh().getData(v, propName, Double.class).get()));
System.out.println(triangulation.getMesh().toPythonTriangulation(v -> triangulation.getMesh().getDoubleData(v, propName)));
/*var points = triangulation.getMesh().getPoints();
System.out.print("[");
for(var dataPoint : points) {
......
......@@ -59,7 +59,7 @@ public class PEikMesh extends GenEikMesh<PVertex, PHalfEdge, PFace> {
double initialEdgeLen,
@NotNull VRectangle bound
) {
super(distanceFunc, edgeLengthFunc, fixPoints, initialEdgeLen, bound, null, Collections.EMPTY_LIST,() -> new PMesh());