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

Commit 1d09c2cb authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

fix a bug in Rupperts algorithm and improve the edge length approximation function.

parent d3774a1b
# nVertices dimension nAttributes boundaryMarker
10 2 0 0
# vertexId x y
1 0.500000 0.500000
2 9.500000 0.500000
3 7.000000 6.000000
4 5.000000 8.000000
5 3.000000 8.000000
6 3.000000 4.000000
7 7.000000 4.000000
8 5.000000 6.000000
9 0.500000 9.500000
10 9.500000 9.500000
#
# nSegments boundaryMarker
10 0
# lineId vertexId1 vertexId2
1 2 10
2 10 9
3 9 1
4 1 2
5 5 6
6 6 7
7 7 3
8 3 8
9 8 4
10 4 5
#
# nHoles
1
# vertexId x y (of a vertex which lies inside the hole)
1 4.666667 5.666667
\ No newline at end of file
# nVertices dimension nAttributes boundaryMarker
12 2 0 0
# vertexId x y
1 0.500000 0.500000
2 7.000000 5.000000
3 7.000000 5.300000
4 3.000000 4.300000
5 9.500000 0.500000
6 3.000000 5.300000
7 3.000000 5.000000
8 7.000000 4.300000
9 3.000000 4.000000
10 7.000000 4.000000
11 0.500000 9.500000
12 9.500000 9.500000
#
# nSegments boundaryMarker
12 0
# lineId vertexId1 vertexId2
1 5 12
2 12 11
3 11 1
4 1 5
5 2 3
6 3 6
7 6 7
8 7 2
9 10 8
10 8 4
11 4 9
12 9 10
#
# nHoles
2
# vertexId x y (of a vertex which lies inside the hole)
1 5.000000 5.150000
2 5.000000 4.150000
\ No newline at end of file
# nVertices dimension nAttributes boundaryMarker
33 2 0 0
# vertexId x y
1 0.500000 0.500000
2 2.300000 7.000000
3 9.500000 0.500000
4 2.800000 7.300000
5 2.800000 6.500000
6 4.800000 6.500000
7 0.900000 6.600000
8 1.200000 8.100000
9 2.300000 7.200000
10 2.800000 7.200000
11 4.100000 7.000000
12 5.000000 2.000000
13 1.000000 6.500000
14 1.300000 8.000000
15 0.800000 6.600000
16 7.000000 2.000000
17 0.500000 9.500000
18 1.100000 8.200000
19 9.500000 9.500000
20 4.000000 6.900000
21 4.900000 6.600000
22 1.400000 7.900000
23 2.500000 6.700000
24 2.700000 7.400000
25 0.700000 6.700000
26 4.900000 7.000000
27 7.000000 5.000000
28 2.200000 6.700000
29 2.500000 7.400000
30 5.000000 5.000000
31 2.700000 6.500000
32 2.100000 6.900000
33 4.000000 6.100000
#
# nSegments boundaryMarker
33 0
# lineId vertexId1 vertexId2
1 3 19
2 19 17
3 17 1
4 1 3
5 25 18
6 18 8
7 8 15
8 15 25
9 7 14
10 14 22
11 22 13
12 13 7
13 2 9
14 9 29
15 29 24
16 24 4
17 4 10
18 10 5
19 5 31
20 31 23
21 23 28
22 28 32
23 32 2
24 16 27
25 27 30
26 30 12
27 12 16
28 33 20
29 20 11
30 11 26
31 26 21
32 21 6
33 6 33
#
# nHoles
5
# vertexId x y (of a vertex which lies inside the hole)
1 0.950000 7.400000
2 1.150000 7.250000
3 2.538272 6.968724
4 6.000000 3.500000
5 4.401111 6.647778
\ No newline at end of file
# nVertices dimension nAttributes boundaryMarker
4 2 0 0
# vertexId x y
1 0.500000 0.500000
2 9.500000 0.500000
3 0.500000 9.500000
4 9.500000 9.500000
#
# nSegments boundaryMarker
4 0
# lineId vertexId1 vertexId2
1 2 4
2 4 3
3 3 1
4 1 2
#
# nHoles
0
# vertexId x y (of a vertex which lies inside the hole)
\ No newline at end of file
package org.vadere.meshing.examples;
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.EdgeLengthFunctionApprox;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.mesh.triangulation.triangulator.gen.GenRuppertsTriangulator;
import org.vadere.meshing.mesh.triangulation.triangulator.impl.PContrainedDelaunayTriangulator;
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 java.io.IOException;
import java.io.InputStream;
import java.util.stream.Collectors;
public class BackgroundMeshExamples {
public static void main(String ... args) throws IOException, InterruptedException {
greenlandLFS();
corridorLFS();
}
public static void greenlandLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/kaiserslautern.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
double theta = 10;
PRuppertsTriangulator ruppert = new PRuppertsTriangulator(
pslg,
p -> Double.POSITIVE_INFINITY,
theta
);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg);
edgeLengthFunctionApprox.printPython();
}
public static void roomLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/room.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg, p -> 2.0);
edgeLengthFunctionApprox.printPython();
}
public static void cornerLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/corner.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg, p -> 2.0);
edgeLengthFunctionApprox.printPython();
}
public static void corridorLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/narrowCorridor.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg, p -> 2.0);
edgeLengthFunctionApprox.printPython();
}
}
......@@ -4,6 +4,7 @@ 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.inter.IPointConstructor;
import org.vadere.meshing.mesh.triangulation.EdgeLengthFunctionApprox;
import org.vadere.meshing.mesh.triangulation.IEdgeLengthFunction;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.EikMeshPoint;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
......@@ -56,7 +57,9 @@ public class EikMeshPlots {
kaiserslautern();
ruppertsAndEikMeshKaiserslautern();*/
bridge();
//bridge();
//roomLFS();
cornerLFS();
}
public static void randomDelaunay() throws IOException {
......@@ -120,16 +123,84 @@ public class EikMeshPlots {
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
Collection<VPolygon> holes = pslg.getHoles();
VPolygon segmentBound = pslg.getSegmentBound();
IDistanceFunction distanceFunction = IDistanceFunction.create(segmentBound, holes);
IPointConstructor<EikMeshPoint> pointConstructor = (x, y) -> new EikMeshPoint(x, y);
// (3) use EikMesh to improve the mesh
double h0 = 0.5;
var meshImprover = new PEikMesh(
distanceFunction,
p -> h0 + 0.5 * Math.abs(distanceFunction.apply(p)),
h0,
new VRectangle(segmentBound.getBounds2D()),
pslg.getHoles()
);
var meshPanel = new PMeshPanel(meshImprover.getMesh(), 1000, 800);
meshPanel.display("Combined distance functions " + h0);
while (!meshImprover.isFinished()) {
meshImprover.improve();
Thread.sleep(20);
meshPanel.repaint();
}
//meshImprover.generate();
write(toTexDocument(TexGraphGenerator.toTikz(meshImprover.getMesh(), f-> lightBlue, 1.0f)), "eikmesh_kaiserslautern_"+ Double.toString(h0).replace('.', '_'));
// display the mesh
meshPanel.display("Combined distance functions " + h0);
}
public static void roomLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/room.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg, p -> 2.0);
edgeLengthFunctionApprox.printPython();
Collection<VPolygon> holes = pslg.getHoles();
VPolygon segmentBound = pslg.getSegmentBound();
IDistanceFunction distanceFunction = IDistanceFunction.create(segmentBound, holes);
// (3) use EikMesh to improve the mesh
double h0 = 0.5;
var meshImprover = new PEikMesh(
distanceFunction,
edgeLengthFunctionApprox,
h0,
new VRectangle(segmentBound.getBounds2D()),
pslg.getHoles()
);
var meshPanel = new PMeshPanel(meshImprover.getMesh(), 1000, 800);
meshPanel.display("Combined distance functions " + h0);
while (!meshImprover.isFinished()) {
meshImprover.improve();
Thread.sleep(20);
meshPanel.repaint();
}
//meshImprover.generate();
write(toTexDocument(TexGraphGenerator.toTikz(meshImprover.getMesh(), f-> lightBlue, 1.0f)), "eikmesh_kaiserslautern_"+ Double.toString(h0).replace('.', '_'));
// display the mesh
meshPanel.display("Combined distance functions " + h0);
}
public static void cornerLFS() throws IOException, InterruptedException {
final InputStream inputStream = MeshExamples.class.getResourceAsStream("/poly/corner.poly");
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
EdgeLengthFunctionApprox edgeLengthFunctionApprox = new EdgeLengthFunctionApprox(pslg, p -> 1.0);
edgeLengthFunctionApprox.printPython();
Collection<VPolygon> holes = pslg.getHoles();
VPolygon segmentBound = pslg.getSegmentBound();
IDistanceFunction distanceFunction = IDistanceFunction.create(segmentBound, holes);
// (3) use EikMesh to improve the mesh
double h0 = 1.0;
var meshImprover = new PEikMesh(
distanceFunction,
p -> h0 + 0.5 * Math.abs(distanceFunction.apply(p)),
edgeLengthFunctionApprox,
h0,
new VRectangle(segmentBound.getBounds2D()),
pslg.getHoles()
......
......@@ -101,6 +101,18 @@ public class PSLG {
this(segmentBound, holes, segments, points, new HashSet<>());
}
public PSLG conclose(@NotNull final VRectangle rectangle) {
if(!rectangle.contains(segmentBound.getBounds2D())) {
throw new IllegalArgumentException();
}
ArrayList<VLine> newLines = new ArrayList<>(segments.size() + segmentBound.getLinePath().size());
newLines.addAll(segments);
newLines.addAll(segmentBound.getLinePath());
PSLG pslg = new PSLG(new VPolygon(rectangle), holes, newLines, points);
pslg.protectionDiscs.addAll(this.protectionDiscs);
return pslg;
}
public PSLG addLines(@NotNull final Collection<VLine> lines) {
ArrayList<VLine> newLines = new ArrayList<>(segments.size() + lines.size());
newLines.addAll(segments);
......
......@@ -11,10 +11,12 @@ import org.vadere.meshing.mesh.inter.IPointConstructor;
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.InterpolationUtil;
import java.util.Comparator;
import java.util.PriorityQueue;
import java.util.function.Function;
public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
......@@ -22,10 +24,18 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
private static final String propName = "edgeLength";
public EdgeLengthFunctionApprox(@NotNull final PSLG pslg) {
public EdgeLengthFunctionApprox(
@NotNull final PSLG pslg,
@NotNull final Function<IPoint, Double> circumRadiusFunc) {
//IPointConstructor<DataPoint<Double>> pointConstructor = (x, y) -> new DataPoint<>(x, y);
var ruppertsTriangulator = new PRuppertsTriangulator(pslg, 10);
/**
* 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!
......@@ -37,7 +47,7 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
for(var v : vertices) {
double minEdgeLen = Double.MAX_VALUE;
for(var u : triangulation.getMesh().getAdjacentVertexIt(v)) {
double len = v.distance(u);
double len = v.distance(u) * 0.9;
if(len < minEdgeLen) {
minEdgeLen = len;
}
......@@ -47,6 +57,11 @@ public class EdgeLengthFunctionApprox implements IEdgeLengthFunction {
}
}
public EdgeLengthFunctionApprox(@NotNull final PSLG pslg) {
this(pslg, p -> Double.POSITIVE_INFINITY);
//IPointConstructor<DataPoint<Double>> pointConstructor = (x, y) -> new DataPoint<>(x, y);
}
public void smooth(double g) {
assert g > 0;
// smooth the function based such that it is g-Lipschitz
......
......@@ -512,6 +512,8 @@ public class GenEikMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
double desiredLen = getDesiredEdgeLength(p1, p2);
double ratio = len / desiredLen;
double absForce = f.apply(ratio);
//System.out.println(len);
//System.out.println(desiredLen);
VPoint force = p1p2.setMagnitude(absForce * desiredLen);
return force;
}
......
......@@ -356,6 +356,9 @@ public class GenRuppertsTriangulator<V extends IVertex, E extends IHalfEdge, F e
split(segment);
}
}
// we require this because split may other edges which are not direct neighbours encroached!
refineSimplex1D();
}
private void deEncrocheSgements() {
......@@ -490,7 +493,20 @@ public class GenRuppertsTriangulator<V extends IVertex, E extends IHalfEdge, F e
return p.distance(line.getVPoint1()) > GeometryUtils.DOUBLE_EPS && p.distance(line.getVPoint2()) > GeometryUtils.DOUBLE_EPS && diameterCircle.contains(p);
}
private boolean isEncroached(@NotNull final E segment) {
private boolean isEncroached(@NotNull final E segment) {
E seg = getMesh().isBoundary(segment) ? getMesh().getTwin(segment) : segment;
VPoint p1 = getMesh().toPoint(getMesh().getNext(seg));
if(isEncroached(seg, p1)) {
return true;
} else if(getMesh().isAtBoundary(seg)) {
return false;
} else {
VPoint p2 = getMesh().toPoint(getMesh().getNext(getMesh().getTwin(seg)));
return isEncroached(seg, p2);
}
}
/*private boolean isEncroached(@NotNull final E segment) {
E seg = getMesh().isBoundary(segment) ? getMesh().getTwin(segment) : segment;
VLine line = getMesh().toLine(seg);
VPoint midPoint = line.midPoint();
......@@ -510,7 +526,7 @@ public class GenRuppertsTriangulator<V extends IVertex, E extends IHalfEdge, F e
}
return false;
}
}*/
// TODO replace this!
private boolean isEncroachedExpensive(@NotNull final E segment) {
......
......@@ -20,6 +20,15 @@ public class PRuppertsTriangulator extends GenRuppertsTriangulator<PVertex, PHal
super(() -> new PMesh(), pslg, minAngle, circumRadiusFunc, true);
}
public PRuppertsTriangulator(
@NotNull final PSLG pslg,
@NotNull final Function<IPoint, Double> circumRadiusFunc,
final double minAngle,
final boolean createHoles) {
super(() -> new PMesh(), pslg, minAngle, circumRadiusFunc, createHoles);
}
public PRuppertsTriangulator(
@NotNull final PSLG pslg,
final double minAngle) {
......
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