From d46f4882129eb23d9ed9360fc3032f19adb43859 Mon Sep 17 00:00:00 2001 From: Benedikt Zoennchen Date: Sat, 13 Jun 2020 20:36:37 +0200 Subject: [PATCH] small performance improvements of the EikMesh algorithm. --- .../examples/MeshQuantityPrinting.java | 8 +- .../mesh/gen/IncrementalTriangulation.java | 2 +- .../org/vadere/meshing/mesh/inter/IMesh.java | 32 +++ .../meshing/mesh/inter/IPolyConnectivity.java | 37 +++- .../meshing/mesh/inter/ITriConnectivity.java | 2 +- .../improver/distmesh/Distmesh.java | 2 +- .../improver/distmesh/Parameters.java | 5 +- .../improver/eikmesh/gen/GenEikMesh.java | 95 ++++++--- .../eikmesh/gen/IEikMeshImprover.java | 2 +- .../mesh/triangulation/plots/RunTimeCPU.java | 2 +- .../plots/qualities/RunTimeCPU.java | 33 ++- .../gen/GenRegularRefinement.java | 6 +- .../GenUniformRefinementTriangulatorSFC.java | 20 +- .../meshing/utils/math/GeometryUtilsMesh.java | 56 +++++ .../vadere/simulator/examples/Curvature.java | 27 ++- .../calculators/mesh/AMeshEikonalSolver.java | 3 +- .../mesh/MeshEikonalSolverFMMIterative.java | 16 +- .../mesh/PredicateEdgeRefinement.java | 46 ++++ .../calculators/mesh/PredicateRefinement.java | 2 + .../TimeCostFunctionFactory.java | 2 +- ...TimeCostGaussianPedestrianDensityMesh.java | 197 ++++++++++++++++++ .../TimeCostPedestrianDensityMesh.java | 48 ++--- .../projects/migration/GeometryCleaner.java | 4 +- .../vadere/util/geometry/GeometryUtils.java | 18 +- 24 files changed, 555 insertions(+), 110 deletions(-) create mode 100644 VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateEdgeRefinement.java create mode 100644 VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostGaussianPedestrianDensityMesh.java diff --git a/VadereMeshing/src/org/vadere/meshing/examples/MeshQuantityPrinting.java b/VadereMeshing/src/org/vadere/meshing/examples/MeshQuantityPrinting.java index 7a6db671a..9081e764c 100644 --- a/VadereMeshing/src/org/vadere/meshing/examples/MeshQuantityPrinting.java +++ b/VadereMeshing/src/org/vadere/meshing/examples/MeshQuantityPrinting.java @@ -50,9 +50,9 @@ import java.util.stream.Collectors; public class MeshQuantityPrinting { public static void main(String... args) throws InterruptedException, IOException { - spaceFillingCurve2(); + //spaceFillingCurve2(); //uniformMeshDiscFunction(0.10); - //uniformMeshDiscFunctionDistMesh(0.05); + uniformMeshDiscFunctionDistMesh(0.05); //distMeshFail(0.05); //delaunyTri("/poly/a.poly"); } @@ -260,9 +260,9 @@ public class MeshQuantityPrinting { int iteration = 1; while (iteration < maxIteration+1) { distmesh.improve(); - if(iteration + 1 == maxIteration+1) { + /*if(iteration + 1 == maxIteration+1) { distmesh.reTriangulate(true); - } + }*/ bufferedWriterIllegalEdges.write(iteration + " " + distmesh.getNumberOfIllegalTriangles() + "\n"); bufferedWriterQualities1.write(printQualities(iteration, distmesh, f -> GeometryUtils.qualityOf(f))); bufferedWriterQualities2.write(printQualities(iteration, distmesh, f -> GeometryUtils.qualityLongestEdgeInCircle(f.p1, f.p2, f.p3))); diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/gen/IncrementalTriangulation.java b/VadereMeshing/src/org/vadere/meshing/mesh/gen/IncrementalTriangulation.java index 47b9cd36d..31e0d9cda 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/gen/IncrementalTriangulation.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/gen/IncrementalTriangulation.java @@ -321,7 +321,7 @@ public class IncrementalTriangulation removePredicate = face -> !polygon.contains(getMesh().toTriangle(face).midPoint()); + Predicate removePredicate = face -> !polygon.contains(getMesh().toMidpoint(face)); cdt.getTriangulation().shrinkBorder(removePredicate, true); IMesh holeMesh = incrementalTriangulation.getMesh(); diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/inter/IMesh.java b/VadereMeshing/src/org/vadere/meshing/mesh/inter/IMesh.java index a5e788b22..ccf10cdd9 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/inter/IMesh.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/inter/IMesh.java @@ -545,6 +545,9 @@ public interface IMesh * @return (optional) a boundary edge */ default Optional getBoundaryEdge(@NotNull final V vertex) { + if(isBoundary(getEdge(vertex))) { + return Optional.of(getEdge(vertex)); + } return streamEdges(vertex).filter(e -> isBoundary(e)).findAny(); } @@ -1360,6 +1363,24 @@ public interface IMesh return new VTriangle(new VPoint(vertices.get(0)), new VPoint(vertices.get(1)), new VPoint(vertices.get(2))); } + default VPoint toMidpoint(@NotNull final F face) { + assert getVertices(face).size() == 3 : "number of vertices of " + face + " is " + getVertices(face).size(); + E edge = getEdge(face); + V v1 = getVertex(edge); + V v2 = getVertex(getNext(edge)); + V v3 = getVertex(getPrev(edge)); + return GeometryUtils.getTriangleMidpoint(getX(v1), getY(v1), getX(v2), getY(v2), getX(v3), getY(v3)); + } + + default VPoint toCircumcenter(@NotNull final F face) { + assert getVertices(face).size() == 3 : "number of vertices of " + face + " is " + getVertices(face).size(); + E edge = getEdge(face); + V v1 = getVertex(edge); + V v2 = getVertex(getNext(edge)); + V v3 = getVertex(getPrev(edge)); + return GeometryUtils.getCircumcenter(getX(v1), getY(v1), getX(v2), getY(v2), getX(v3), getY(v3)); + } + /** * Returns the midpoint {@link VPoint} of a triangle defined by the face. * Assumption: The face represents a triangle, i.e. it has exactly 3 distinct points. This @@ -2362,6 +2383,17 @@ public interface IMesh } + default String toPythonValues(@NotNull final Function evalPoint) { + StringBuilder builder = new StringBuilder(); + List vertices = getVertices(); + for(V v : vertices) { + builder.append(evalPoint.apply(v) + ","); + } + builder.delete(builder.length()-1, builder.length()); + builder.append("\n"); + return builder.toString(); + } + /** * Constructs and returns a string which can be used to construct a matplotlib Triangulation * which is helpful to plot the mesh. diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/inter/IPolyConnectivity.java b/VadereMeshing/src/org/vadere/meshing/mesh/inter/IPolyConnectivity.java index 8ac4ec24f..3a129f164 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/inter/IPolyConnectivity.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/inter/IPolyConnectivity.java @@ -116,7 +116,7 @@ public interface IPolyConnectivity isAtBoundary(edge)).findAny().ifPresent(edge -> getMesh().setEdge(vertex, edge)); + getMesh().streamEdges(vertex).filter(edge -> getMesh().isBoundary(edge)).findAny().ifPresent(edge -> getMesh().setEdge(vertex, edge)); } /** @@ -632,14 +632,15 @@ public interface IPolyConnectivity createHole(@NotNull final F face, @NotNull final Predicate mergeCondition, final boolean deleteIsoletedVertices) { + default Optional createHole(@NotNull final F face, @NotNull final Predicate mergeCondition, final boolean deleteIsoletedVertices, final boolean vertexAdjust) { if(mergeCondition.test(face)) { getMesh().toHole(face); - shrinkBoundary(face, mergeCondition, deleteIsoletedVertices); + shrinkBoundary(face, mergeCondition, deleteIsoletedVertices, vertexAdjust); return Optional.of(face); } else { @@ -665,6 +666,10 @@ public interface IPolyConnectivity createHole(@NotNull final F face, @NotNull final Predicate mergeCondition, final boolean deleteIsoletedVertices) { + return createHole(face, mergeCondition, deleteIsoletedVertices, true); + } + /** * Shrinks the border as long as the removeCondition is satisfied i.e. a face will be removed if * it is at the border (during the shrinking process) and satisfies the condition. Like a virus this @@ -675,9 +680,14 @@ public interface IPolyConnectivity removeCondition, final boolean deleteIsolatedVertices, final boolean vertexAdjust) { + shrinkBoundary(getMesh().getBorder(), removeCondition, deleteIsolatedVertices, vertexAdjust); + } + default void shrinkBorder(final Predicate removeCondition, final boolean deleteIsolatedVertices) { - shrinkBoundary(getMesh().getBorder(), removeCondition, deleteIsolatedVertices); + shrinkBorder(removeCondition, deleteIsolatedVertices, true); } default void shrinkBoundary(final Predicate removeCondition, final boolean deleteIsolatedVertices) { @@ -697,8 +707,9 @@ public interface IPolyConnectivity removeCondition, final boolean deleteIsolatedVertices) { + default void shrinkBoundary(@NotNull final F boundary, final Predicate removeCondition, final boolean deleteIsolatedVertices, final boolean adjustVertices) { assert getMesh().isBoundary(boundary); List boundaryFaces = getMesh().getFaces(boundary); @@ -715,13 +726,17 @@ public interface IPolyConnectivity removeCondition, final boolean deleteIsolatedVertices) { + shrinkBoundary(boundary, removeCondition, deleteIsolatedVertices, true); + } + default void removeFacesAtBoundary(@NotNull final Predicate mergePredicate, @NotNull final Predicate errorPredicate) throws IllegalMeshException { mergeFaces(getMesh().getBorder(), mergePredicate, errorPredicate, true, 1); for(F face : getMesh().getHoles()) { @@ -1215,8 +1230,9 @@ public interface IPolyConnectivity neighbour.equals(boundary)).count() > 0; @@ -1357,7 +1373,9 @@ public interface IPolyConnectivity adjustVertex(v)); + if(adjustVertices) { + vertices.stream().filter(getMesh()::isAlive).forEach(v -> adjustVertex(v)); + } } } @@ -1371,6 +1389,9 @@ public interface IPolyConnectivityLinkedList. */ default LinkedList straightWalk2DGatherDirectional(@NotNull final F face, @NotNull final VPoint direction, @NotNull final Predicate additionalStopCondition) { - VPoint q = getMesh().toTriangle(face).midPoint(); + VPoint q = getMesh().toMidpoint(face); assert getMesh().toTriangle(face).contains(q); Predicate defaultStopCondion = e -> isRightOf(q.x, q.y, e); diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Distmesh.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Distmesh.java index aaba4d449..9c2014871 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Distmesh.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Distmesh.java @@ -112,7 +112,7 @@ public class Distmesh { } public void reTriangulate(boolean force) { - if(true || force || firstStep || maxMovementLen / initialEdgeLen > Parameters.TOL) { + if(force || firstStep || maxMovementLen / initialEdgeLen > Parameters.TOL) { firstStep = false; nTriangulations++; diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Parameters.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Parameters.java index c8fff387b..f4b46228e 100755 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Parameters.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/distmesh/Parameters.java @@ -17,7 +17,7 @@ public class Parameters { public final static boolean uniform = false; public final static String method = "Distmesh"; // "Distmesh" or "Density" public final static double qualityMeasurement = 0.95; - public final static double qualityConvergence = 0.001; + public final static double qualityConvergence = 0.01; public final static double MINIMUM = 0.25; public final static double DENSITYWEIGHT = 2; public final static int NPOINTS = 100000; @@ -25,5 +25,6 @@ public class Parameters { public final static int SAMPLEDIVISION = 10; public final static int SEGMENTDIVISION = 0; //TODO increase this - public final static int MAX_NUMBER_OF_STEPS = 100; + public final static int MAX_NUMBER_OF_STEPS = 200; + public final static int HIGHEST_LEGAL_TEST = 10; } diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/GenEikMesh.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/GenEikMesh.java index 8d68eb138..cfbc6c0ab 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/GenEikMesh.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/GenEikMesh.java @@ -3,12 +3,15 @@ package org.vadere.meshing.mesh.triangulation.improver.eikmesh.gen; import org.apache.commons.lang3.tuple.Pair; import org.jetbrains.annotations.NotNull; import org.vadere.meshing.mesh.IllegalMeshException; +import org.vadere.meshing.mesh.inter.IEdgeContainerBoolean; import org.vadere.meshing.mesh.inter.IFace; import org.vadere.meshing.mesh.inter.IHalfEdge; import org.vadere.meshing.mesh.inter.IMesh; import org.vadere.meshing.mesh.inter.IMeshSupplier; import org.vadere.meshing.mesh.inter.IIncrementalTriangulation; import org.vadere.meshing.mesh.inter.IVertex; +import org.vadere.meshing.mesh.inter.IVertexContainerBoolean; +import org.vadere.meshing.mesh.inter.IVertexContainerDouble; import org.vadere.meshing.mesh.iterators.EdgeIterator; import org.vadere.meshing.mesh.iterators.EdgeIteratorReverse; import org.vadere.meshing.mesh.triangulation.improver.IMeshImprover; @@ -129,6 +132,12 @@ public class GenEikMesh private static final String propVelocityY = "velocityY"; private static final String propAbsVelocity = "absVelocity"; + private final IVertexContainerBoolean fixpointC; + private final IEdgeContainerBoolean constraintC; + private final IVertexContainerDouble velocityXC; + private final IVertexContainerDouble velocityYC; + private final IVertexContainerDouble absVelocityC; + /** * Constructor to use EikMesh on an existing {@link org.vadere.meshing.mesh.inter.ITriangulation}, that is * EikMesh uses this triangulation as a bases. It will refineSimplex2D the triangulation by using a longest edge @@ -177,6 +186,12 @@ public class GenEikMesh } this.useSlidingLines = true; this.smoothBorder = false; + + this.fixpointC = triangulation.getMesh().getBooleanVertexContainer(propFixPoint); + this.constraintC = triangulation.getMesh().getBooleanEdgeContainer(propConstrained); + this.velocityXC = triangulation.getMesh().getDoubleVertexContainer(propVelocityX); + this.velocityYC = triangulation.getMesh().getDoubleVertexContainer(propVelocityY); + this.absVelocityC = triangulation.getMesh().getDoubleVertexContainer(propAbsVelocity); } /** @@ -244,6 +259,12 @@ public class GenEikMesh this.useSlidingLines = false; this.smoothBorder = true; + + this.fixpointC = triangulation.getMesh().getBooleanVertexContainer(propFixPoint); + this.constraintC = triangulation.getMesh().getBooleanEdgeContainer(propConstrained); + this.velocityXC = triangulation.getMesh().getDoubleVertexContainer(propVelocityX); + this.velocityYC = triangulation.getMesh().getDoubleVertexContainer(propVelocityY); + this.absVelocityC = triangulation.getMesh().getDoubleVertexContainer(propAbsVelocity); } /** @@ -288,6 +309,12 @@ public class GenEikMesh initialEdgeLen, distanceFunc, generateFixPoints()); + + this.fixpointC = refiner.getMesh().getBooleanVertexContainer(propFixPoint); + this.constraintC = refiner.getMesh().getBooleanEdgeContainer(propConstrained); + this.velocityXC = refiner.getMesh().getDoubleVertexContainer(propVelocityX); + this.velocityYC = refiner.getMesh().getDoubleVertexContainer(propVelocityY); + this.absVelocityC = refiner.getMesh().getDoubleVertexContainer(propAbsVelocity); } public GenEikMesh( @@ -337,6 +364,7 @@ public class GenEikMesh deps = 0.0001 * initialEdgeLen; computeFixPoints(); + refiner.getConstrains().forEach(e -> setConstraint(e, true)); if(useSlidingLines) { computeSlidingLines(); @@ -549,7 +577,7 @@ public class GenEikMesh private void computeForce(final V vertex) { // TODO: Get rid of IPoint IPoint p1 = getMesh().getPoint(vertex); - boolean isAtBoundary = getMesh().isAtBoundary(vertex); + boolean isAtBoundary = isBoundary(vertex); for(E edge : getMesh().getEdgeIt(vertex)) { @@ -765,6 +793,11 @@ public class GenEikMesh * @return true if the movement is legal, false otherwise */ private boolean isLegalMove(@NotNull final V vertex, double newX, double newY) { + // only test for early iterations! + if(nSteps > Parameters.HIGHEST_LEGAL_TEST) { + return true; + } + //if(vertex.distance(newX, newY) > GeometryUtils.DOUBLE_EPS) { // TODO: at the boundary vertices can still overtake each other. @@ -810,7 +843,7 @@ public class GenEikMesh */ private void updateFace(@NotNull F face) { if(canBreak(face) && isBreaking(face)) { - VPoint circumcenter = getMesh().toTriangle(face).getCircumcenter(); + VPoint circumcenter = getMesh().toCircumcenter(face); getTriangulation().splitTriangle(face, getMesh().createPoint(circumcenter.getX(), circumcenter.getY()), false); } } @@ -854,9 +887,9 @@ public class GenEikMesh newPosition = new VPoint(v1.getX(), v1.getY()); } else if(isFixPoint(v2)) { newPosition = new VPoint(v2.getX(), v2.getY()); - } else if(getMesh().isAtBoundary(v1) || (isSlidePoint(v1) && !isSlidePoint(v2))) { + } else if(isBoundary(v1) || (isSlidePoint(v1) && !isSlidePoint(v2))) { newPosition = new VPoint(v1.getX(), v1.getY()); - } else if(getMesh().isAtBoundary(v2) || (isSlidePoint(v2) && !isSlidePoint(v1))) { + } else if(isBoundary(v2) || (isSlidePoint(v2) && !isSlidePoint(v1))) { newPosition = new VPoint(v2.getX(), v2.getY()); } else { newPosition = new VPoint((v1.getX() + v2.getX()) * 0.5, (v1.getY() + v2.getY()) * 0.5); @@ -953,8 +986,9 @@ public class GenEikMesh if(faceToQuality(face) > 0.95) { E edge = getMesh().getEdge(face); if(edgeLengthFunc.apply(getMesh().toLine(edge).midPoint()) * 2.1 <= getMesh().toLine(edge).length()) { - VPoint circumcenter = getMesh().toTriangle(face).getCircumcenter(); - return getMesh().toTriangle(face).contains(circumcenter); + VPoint circumcenter = getMesh().toCircumcenter(face); + return triangulation.contains(circumcenter.getX(), circumcenter.getY(), face); + //getMesh().toTriangle(face).contains(circumcenter); } } return false; @@ -1041,7 +1075,7 @@ public class GenEikMesh */ private VPoint computeProjection(@NotNull final V vertex) { // we only project boundary vertices back - if(getMesh().isAtBoundary(vertex)) { + if(isBoundary(vertex)) { // TODO: get rid of VPoint VPoint position = getMesh().toPoint(vertex); @@ -1373,7 +1407,7 @@ public class GenEikMesh * This takes O(n) time where n is the number of removed faces which will be consumed. */ private void removeFacesAtBoundary() { - Predicate isOutside = f -> distanceFunc.apply(getMesh().toTriangle(f).midPoint()) > 0; + Predicate isOutside = f -> distanceFunc.apply(getMesh().toMidpoint(f)) > 0; Predicate isSeparated = f -> getMesh().isSeparated(f); //Predicate isInvalid = f -> !getTriangulation().isValid(f); Predicate isOfLowQuality = f -> faceToQuality(f) < Parameters.MIN_TRIANGLE_QUALITY && !isShortBoundaryEdge(f); @@ -1436,7 +1470,7 @@ public class GenEikMesh *

Shrinks the boundary such that there are no more triangles outside the boundary i.e. where the distance is positive.

*/ private void shrinkBoundary() { - Predicate removePredicate = face -> distanceFunc.apply(getTriangulation().getMesh().toTriangle(face).midPoint()) > 0; + Predicate removePredicate = face -> distanceFunc.apply(getTriangulation().getMesh().toMidpoint(face)) > 0; getTriangulation().shrinkBoundary(removePredicate, true); } @@ -1445,7 +1479,7 @@ public class GenEikMesh * Note the border is part of the whole boundary which is defined by the border and the holes.

*/ private void shrinkBorder() { - Predicate removePredicate = face -> distanceFunc.apply(getTriangulation().getMesh().toTriangle(face).midPoint()) > 0; + Predicate removePredicate = face -> distanceFunc.apply(getTriangulation().getMesh().toMidpoint(face)) > 0; getTriangulation().shrinkBorder(removePredicate, true); } @@ -1456,7 +1490,7 @@ public class GenEikMesh List faces = getTriangulation().getMesh().getFaces(); for(F face : faces) { if(!getTriangulation().getMesh().isDestroyed(face) && !getTriangulation().getMesh().isHole(face)) { - getTriangulation().createHole(face, f -> distanceFunc.apply(getTriangulation().getMesh().toTriangle(f).midPoint()) > 0, true); + getTriangulation().createHole(face, f -> distanceFunc.apply(getTriangulation().getMesh().toMidpoint(f)) > 0, true); } } } @@ -1542,38 +1576,41 @@ public class GenEikMesh } } + private boolean isBoundary(@NotNull final V vertex) { + return getMesh().isBoundary(getMesh().getEdge(vertex)); + } /* * The following methods are helper methods to quickly access properties saved on vertices, edges and faces */ private void setConstraint(E edge, boolean constraint) { - getMesh().setBooleanData(edge, propConstrained, constraint); + constraintC.setValue(edge, constraint); } private boolean isConstrained(E edge) { - return getMesh().getBooleanData(edge, propConstrained); + return constraintC.getValue(edge); } private void setFixPoint(V vertex, boolean fixPoint) { - getMesh().setBooleanData(vertex, propFixPoint, fixPoint); + fixpointC.setValue(vertex, fixPoint); } public boolean isFixPoint(V vertex) { - return getMesh().getBooleanData(vertex, propFixPoint); + return fixpointC.getValue(vertex); } private void setVelocity(V vertex, IPoint velocity) { - setVelocityX(vertex, velocity.getX()); - setVelocityY(vertex, velocity.getY()); + velocityXC.setValue(vertex, velocity.getX()); + velocityYC.setValue(vertex, velocity.getY()); } private void setVelocityX(V vertex, double velX) { - getMesh().setDoubleData(vertex, propVelocityX, velX); + velocityXC.setValue(vertex, velX); } private void setVelocityY(V vertex, double velY) { - getMesh().setDoubleData(vertex, propVelocityY, velY); + velocityYC.setValue(vertex, velY); } private void increaseVelocity(V vertex, IPoint dvelocity) { @@ -1582,21 +1619,21 @@ public class GenEikMesh } private double getVelocityX(V vertex) { - return getMesh().getDoubleData(vertex, propVelocityX); + return velocityXC.getValue(vertex); } private void increaseVelocityX(V vertex, double dVelX) { - double velX = getMesh().getDoubleData(vertex, propVelocityX); - getMesh().setDoubleData(vertex, propVelocityX, velX + dVelX); + double velX = velocityXC.getValue(vertex); + velocityXC.setValue(vertex, velX + dVelX); } private void increaseVelocityY(V vertex, double dVelY) { - double velY = getMesh().getDoubleData(vertex, propVelocityY); - getMesh().setDoubleData(vertex, propVelocityY, velY + dVelY); + double velY = velocityYC.getValue(vertex); + velocityYC.setValue(vertex, velY + dVelY); } private double getVelocityY(V vertex) { - return getMesh().getDoubleData(vertex, propVelocityY); + return velocityYC.getValue(vertex); } private VPoint getVelocity(V vertex) { @@ -1604,17 +1641,17 @@ public class GenEikMesh } private void setAbsVelocity(V vertex, double absVelocity) { - getMesh().setDoubleData(vertex, propAbsVelocity, absVelocity); + absVelocityC.setValue(vertex, absVelocity); } private void increaseAbsVelocity(V vertex, double dAbsVelocity) { - double absVel = getMesh().getDoubleData(vertex, propAbsVelocity); - getMesh().setDoubleData(vertex, propAbsVelocity, absVel + dAbsVelocity); + double absVel = absVelocityC.getValue(vertex); + absVelocityC.setValue(vertex, absVel + dAbsVelocity); } // setter to configure the algorithm strategy. private double getAbsVelocity(V vertex) { - return getMesh().getDoubleData(vertex, propAbsVelocity); + return absVelocityC.getValue(vertex); } public void setAllowEdgeSplits(final boolean allowEdgeSplits) { diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/IEikMeshImprover.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/IEikMeshImprover.java index 985a437f1..decfb111c 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/IEikMeshImprover.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/improver/eikmesh/gen/IEikMeshImprover.java @@ -26,7 +26,7 @@ public interface IEikMeshImprover outsidePredicate(@NotNull final IDistanceFunction distanceFunc) { - return f -> distanceFunc.apply(getMesh().toTriangle(f).midPoint()) > 0; + return f -> distanceFunc.apply(getMesh().toMidpoint(f)) > 0; } default Predicate lowQualityPredicate() { diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/RunTimeCPU.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/RunTimeCPU.java index 32c46d27b..ac0c700fc 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/RunTimeCPU.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/RunTimeCPU.java @@ -138,6 +138,6 @@ public class RunTimeCPU extends JFrame { } public static void main(String[] args) { - stepUniformRing(0.05, 0.05, 0.05); + stepUniformRing(0.05, 0.05, 0.05); } } diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/qualities/RunTimeCPU.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/qualities/RunTimeCPU.java index 21a74cfc7..0936e35f6 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/qualities/RunTimeCPU.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/plots/qualities/RunTimeCPU.java @@ -7,6 +7,7 @@ import org.vadere.meshing.mesh.gen.AMesh; import org.vadere.meshing.mesh.gen.AVertex; import org.vadere.meshing.mesh.inter.IMeshSupplier; import org.vadere.meshing.mesh.inter.IPointConstructor; +import org.vadere.meshing.mesh.triangulation.improver.distmesh.DistmeshPanel; import org.vadere.meshing.mesh.triangulation.improver.eikmesh.gen.GenEikMesh; import org.vadere.meshing.mesh.gen.MeshPanel; import org.vadere.meshing.mesh.triangulation.improver.eikmesh.EikMeshPoint; @@ -29,6 +30,10 @@ public class RunTimeCPU extends JFrame { private static final Logger log = Logger.getLogger(RunTimeCPU.class); + static { + log.setInfo(); + } + /** * Each geometry is contained this bounding box. */ @@ -79,7 +84,6 @@ public class RunTimeCPU extends JFrame { private static void stepAdaptiveRingEikMesh(double startLen, double endLen, double stepLen) { IMeshSupplier supplier = () -> new AMesh(); IDistanceFunction distanceFunc = p -> Math.abs(0.7 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 0.3; - IEdgeLengthFunction adaptiveEdgeLength = p -> 1.0 + Math.max(-distanceFunc.apply(p), 0) * 8.0; List obstacles = new ArrayList<>(); double initialEdgeLength = startLen; @@ -91,6 +95,8 @@ public class RunTimeCPU extends JFrame { while (initialEdgeLength >= minInitialEdgeLength) { initlialEdgeLengths.add(initialEdgeLength); + final double currentEdgeLen = initialEdgeLength; + IEdgeLengthFunction adaptiveEdgeLength = p -> currentEdgeLen + Math.max(-distanceFunc.apply(p), 0) * 0.4; GenEikMesh meshGenerator = new GenEikMesh<>( distanceFunc, adaptiveEdgeLength, @@ -98,6 +104,10 @@ public class RunTimeCPU extends JFrame { bbox, obstacles, supplier); + while (!meshGenerator.isInitialized()) { + meshGenerator.initialize(); + } + StopWatch overAllTime = new StopWatch(); int steps = 0; @@ -121,12 +131,12 @@ public class RunTimeCPU extends JFrame { nVertices.add(meshGenerator.getMesh().getVertices().size()); runTimes.add( overAllTime.getTime()); - /*PSMeshingPanel, AHalfEdge, AFace> distmeshPanel = new PSMeshingPanel(meshGenerator.getMesh(), f -> false, 1000, 800, bbox); + MeshPanel distmeshPanel = new MeshPanel<>(meshGenerator.getMesh(),1000, 800); JFrame frame = distmeshPanel.display(); frame.setVisible(true); frame.setTitle("uniformRing()"); frame.setDefaultCloseOperation(WindowConstants.DISPOSE_ON_CLOSE); - distmeshPanel.repaint();*/ + distmeshPanel.repaint(); initialEdgeLength = initialEdgeLength * 0.5; @@ -142,7 +152,7 @@ public class RunTimeCPU extends JFrame { private static void stepAdaptiveRingDistMesh(double startLen, double endLen, double stepLen) { IMeshSupplier supplier = () -> new AMesh(); IDistanceFunction distanceFunc = p -> Math.abs(0.7 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 0.3; - IEdgeLengthFunction adaptiveEdgeLength = p -> 1.0 + Math.max(-distanceFunc.apply(p), 0) * 8.0; + //IEdgeLengthFunction adaptiveEdgeLength = p -> 1.0 + Math.max(-distanceFunc.apply(p), 0) * 0.4; List obstacles = new ArrayList<>(); double initialEdgeLength = startLen; @@ -154,6 +164,8 @@ public class RunTimeCPU extends JFrame { while (initialEdgeLength >= minInitialEdgeLength) { initlialEdgeLengths.add(initialEdgeLength); + final double currentEdgeLen = initialEdgeLength; + IEdgeLengthFunction adaptiveEdgeLength = p -> currentEdgeLen + Math.max(-distanceFunc.apply(p), 0) * 0.1; Distmesh meshGenerator = new Distmesh(distanceFunc, adaptiveEdgeLength, initialEdgeLength, @@ -169,24 +181,24 @@ public class RunTimeCPU extends JFrame { meshGenerator.step(); overAllTime.suspend(); steps++; - } while (!meshGenerator.hasMaximalSteps()); + } while (steps <= 100); log.info("#vertices: " + meshGenerator.getPoints().size()); log.info("quality: " + meshGenerator.getQuality()); log.info("#step: " + steps); log.info("#tris: " + meshGenerator.getNumberOfReTriangulations()); log.info("overall time: " + overAllTime.getTime() + "[ms]"); - log.info("step avg time: " + overAllTime.getTime() / steps + "[ms]"); + log.info("step avg time: " + overAllTime.getTime() / steps + "[ms]\n"); nVertices.add(meshGenerator.getPoints().size()); runTimes.add( overAllTime.getTime()); - /*PSDistmeshPanel distmeshPanel = new PSDistmeshPanel(meshGenerator, 1000, 800, bbox, t -> false); + DistmeshPanel distmeshPanel = new DistmeshPanel(meshGenerator, 1000, 800, bbox, t -> false); JFrame frame = distmeshPanel.display(); frame.setVisible(true); frame.setTitle("uniformRing()"); frame.setDefaultCloseOperation(WindowConstants.DISPOSE_ON_CLOSE); - distmeshPanel.repaint();*/ + distmeshPanel.repaint(); initialEdgeLength = initialEdgeLength * 0.5; @@ -200,7 +212,8 @@ public class RunTimeCPU extends JFrame { } public static void main(String[] args) { - stepAdaptiveRingDistMesh(0.1, 0.005, 0.01); - //stepAdaptiveRingEikMesh(0.1, 0.005, 0.01); + //stepAdaptiveRingDistMesh(0.1, 0.005, 0.01); + stepAdaptiveRingEikMesh(0.1, 0.005, 0.01); + //stepAdaptiveRingEikMesh(0.005, 0.005, 0.01); } } diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenRegularRefinement.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenRegularRefinement.java index b8d16905b..920ab1de2 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenRegularRefinement.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenRegularRefinement.java @@ -96,10 +96,10 @@ public class GenRegularRefinement getLevel(e) == (maxLevel-1) && edgeRefinementPredicate.test(e); //this.edgeAddToRefine = e -> getLevel(e) == (maxLevel-1) && edgeRefinementPredicate.test(e); - VPoint p = new VPoint(5,5); - - this.edgeRefinementPredicate = e -> !getMesh().isBoundary(e) && getMesh().toTriangle(getMesh().getFace(e)).midPoint().distance(p) < 3.0 && (!isGreen(e) || getMesh().toLine(e).length() > 0.5); + //VPoint p = new VPoint(5,5); + //this.edgeRefinementPredicate = e -> !getMesh().isBoundary(e) && getMesh().toTriangle(getMesh().getFace(e)).midPoint().distance(p) < 3.0 && (!isGreen(e) || getMesh().toLine(e).length() > 0.5); + this.edgeRefinementPredicate = e -> getLevel(e) == (maxLevel-1) && (edgeRefinementPredicate.test(e)); this.finished = false; this.coarse = false; this.toRefine = new LinkedList<>(); diff --git a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenUniformRefinementTriangulatorSFC.java b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenUniformRefinementTriangulatorSFC.java index 89e918a10..d4eb161a7 100644 --- a/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenUniformRefinementTriangulatorSFC.java +++ b/VadereMeshing/src/org/vadere/meshing/mesh/triangulation/triangulator/gen/GenUniformRefinementTriangulatorSFC.java @@ -527,6 +527,7 @@ public class GenUniformRefinementTriangulatorSFC sierpinksyFaceOrder = sfc.asList().stream().map(node -> getMesh().getFace(node.getEdge())).collect(Collectors.toList()); shrinkBorder(); createHoles(); + adjustBoundaryVertexEdges(); //triangulation.smoothBorder(); sierpinksyFaceOrder.removeIf(face -> getMesh().isDestroyed(face) || getMesh().isHole(face)); @@ -617,8 +618,8 @@ public class GenUniformRefinementTriangulatorSFC */ private void shrinkBorder() { - Predicate removePredicate = face -> distFunc.apply(triangulation.getMesh().toTriangle(face).midPoint()) > 0; - triangulation.shrinkBorder(removePredicate, true); + Predicate removePredicate = face -> distFunc.apply(triangulation.getMesh().toMidpoint(face)) > 0; + triangulation.shrinkBorder(removePredicate, true, false); } @Override @@ -652,6 +653,7 @@ public class GenUniformRefinementTriangulatorSFC cdt = new GenConstrainSplitter<>(getTriangulation(), lines, GeometryUtils.DOUBLE_EPS, sfc); cdt.generate(false); + getTriangulation().setCanIllegalPredicate(e -> true); projections.putAll(cdt.getProjections()); constrains.addAll(cdt.getConstrains()); @@ -661,9 +663,21 @@ public class GenUniformRefinementTriangulatorSFC faces = triangulation.getMesh().getFaces(); for(F face : faces) { if(!triangulation.getMesh().isDestroyed(face) && !triangulation.getMesh().isHole(face)) { - triangulation.createHole(face, f -> distFunc.apply(triangulation.getMesh().toTriangle(f).midPoint()) > 0, true); + triangulation.createHole(face, f -> distFunc.apply(triangulation.getMesh().toMidpoint(f)) > 0, true, false); + } + } + } + + private void adjustBoundaryVertexEdges() { + for(F hole : getMesh().getHoles()) { + for(V v : getMesh().getVertexIt(hole)) { + triangulation.adjustVertex(v); } } + + for(V v : getMesh().getVertexIt(getMesh().getBorder())) { + triangulation.adjustVertex(v); + } } /*private void createHoles() { diff --git a/VadereMeshing/src/org/vadere/meshing/utils/math/GeometryUtilsMesh.java b/VadereMeshing/src/org/vadere/meshing/utils/math/GeometryUtilsMesh.java index 4cdf4590c..ab1312ed8 100644 --- a/VadereMeshing/src/org/vadere/meshing/utils/math/GeometryUtilsMesh.java +++ b/VadereMeshing/src/org/vadere/meshing/utils/math/GeometryUtilsMesh.java @@ -106,6 +106,62 @@ public class GeometryUtilsMesh { return new double[]{curvature, gaussianCurvature}; } + public static double[] curvature(@NotNull final IMesh mesh, + @NotNull final Function f, + @NotNull final Predicate valid, + @NotNull final E edge) { + double alpha = 0; + double beta = 0; + V v = mesh.getVertex(edge); + if(valid.test(v)) { + double vx = mesh.getX(v); + double vy = mesh.getY(v); + double vz = f.apply(v); + + if(!mesh.isAtBoundary(edge)) { + V vi = mesh.getTwinVertex(edge); + V vj = mesh.getTwinVertex(mesh.getTwin(mesh.getNext(edge))); + V vk = mesh.getTwinVertex(mesh.getPrev(mesh.getTwin(edge))); + + double vix = mesh.getX(vi); + double viy = mesh.getY(vi); + double viz = f.apply(vi); + + double vjx = mesh.getX(vj); + double vjy = mesh.getY(vj); + double vjz = f.apply(vj); + + if(valid.test(vi) && valid.test(vj) && valid.test(vk)) { + double eix = vix - vx; + double eiy = viy - vy; + double eiz = viz - vz; + + double ejx = vjx - vx; + double ejy = vjy - vy; + double ejz = vjz - vz; + + double ekx = mesh.getX(vk) - vx; + double eky = mesh.getY(vk) - vy; + double ekz = f.apply(vk) - vz; + + alpha = GeometryUtils.angle3D(eix, eiy, eiz, ejx, ejy, ejz); + double[] ni = new double[]{ 0, 0, 0 }; + double[] nk = new double[]{ 0, 0, 0 }; + + GeometryUtils.cross(new double[]{ eix, eiy, eiz}, new double[]{ ejx, ejy, ejz}, ni); + GeometryUtils.cross(new double[]{ ekx, eky, ekz }, new double[]{ eix, eiy, eiz}, nk); + GeometryUtils.norm3D(ni); + GeometryUtils.norm3D(nk); + + double betai = GeometryUtils.angle3D(nk[0], nk[1], nk[2], ni[0], ni[1], ni[2]); + beta = Math.sqrt(eix * eix + eiy * eiy) * betai; + } + } + } + + return new double[]{alpha, beta}; + } + /* public static double[] curvature( @NotNull final IMesh mesh, diff --git a/VadereSimulator/src/org/vadere/simulator/examples/Curvature.java b/VadereSimulator/src/org/vadere/simulator/examples/Curvature.java index d6b2daa1c..968d77522 100644 --- a/VadereSimulator/src/org/vadere/simulator/examples/Curvature.java +++ b/VadereSimulator/src/org/vadere/simulator/examples/Curvature.java @@ -24,12 +24,15 @@ import java.io.BufferedWriter; import java.io.File; import java.io.IOException; import java.util.Arrays; +import java.util.HashSet; +import java.util.List; +import java.util.Set; import java.util.function.Predicate; public class Curvature { public static void main(String ... args) throws IOException, InterruptedException { - adaptoveEikonalSolver(); + interativeEikonalSolver(); } public static void interativeEikonalSolver() throws InterruptedException, IOException { @@ -81,14 +84,24 @@ public class Curvature { panel.repaint(); var tri = improver.getTriangulation(); - for(int i = 0; i < 8; i++) { + for(int i = 0; i < 4; i++) { + + var curTri = tri; + List list = tri.getMesh().getBoundaryVertices(); + Set initialVertices = new HashSet<>(); + initialVertices.addAll(list); + list.stream().forEach(v -> curTri.getMesh().getAdjacentVertexIt(v).forEach(u -> initialVertices.add(u))); + + Set initialVertices2 = new HashSet<>(); + initialVertices.stream().forEach(v -> curTri.getMesh().getAdjacentVertexIt(v).forEach(u -> initialVertices2.add(u))); + initialVertices2.addAll(initialVertices); var solver = new MeshEikonalSolverFMM<>( p -> 1, //Arrays.asList(new VPoint(5, 5)), tri, - tri.getMesh().getBoundaryVertices(), - p -> -distanceFunction.apply(p)); + initialVertices2, + p -> Math.abs(distanceFunction.apply(p))); solver.solve(); System.out.println(solver.getPotential(5,5)); meshWriter.write(tri.getMesh().toPythonTriangulation(v -> solver.getPotential(v))); @@ -107,7 +120,7 @@ public class Curvature { } System.out.println(maxCurvature); - final double minEdgeLen = 0.025; + final double minEdgeLen = 0.1; final var pTri = tri; Predicate edgePredicate = e -> { VLine line = pTri.getMesh().toLine(e); @@ -121,7 +134,7 @@ public class Curvature { pTri.getTriPoints(face, x, y, z, "curvature"); double totalArea = GeometryUtils.areaOfPolygon(x, y); double curvature = InterpolationUtil.barycentricInterpolation(x, y, z, totalArea, p.getX(), p.getY()); - return len > minEdgeLen && curvature > 0.01; + return len > minEdgeLen && curvature > 0.25; }; var triangulation = new IncrementalTriangulation<>(tri.getMesh().clone(), e -> true); @@ -202,7 +215,7 @@ public class Curvature { //Arrays.asList(new VPoint(5, 5)), improver.getTriangulation(), improver.getTriangulation().getMesh().getBoundaryVertices(), - p -> -distanceFunction.apply(p)); + p -> Math.abs(distanceFunction.apply(p))); solver.solve(); meshWriter.write(solver.getTriangulation().getMesh().toPythonTriangulation(v -> solver.getPotential(v))); meshWriter.close(); diff --git a/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/AMeshEikonalSolver.java b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/AMeshEikonalSolver.java index 7e0411d9f..f83d25769 100644 --- a/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/AMeshEikonalSolver.java +++ b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/AMeshEikonalSolver.java @@ -228,7 +228,8 @@ public abstract class AMeshEikonalSolver refiner; + private @NotNull IDistanceFunction distFunc; + public MeshEikonalSolverFMMIterative(@NotNull final ITimeCostFunction timeCostFunction, @NotNull final IIncrementalTriangulation triangulation, @@ -53,6 +56,7 @@ public class MeshEikonalSolverFMMIterative(timeCostFunction, triangulation, targetVertices, distFunc)); + this.distFunc = distFunc; } public MeshEikonalSolverFMMIterative(@NotNull final MeshEikonalSolverFMM solver) { @@ -66,6 +70,8 @@ public class MeshEikonalSolverFMMIterative predicate = new PredicateRefinement<>(triangulation, clone, MIN_EDGE_LEN, MAX_CURVATURE); + //final Predicate predicate = new PredicateEdgeRefinement<>(solver, MIN_EDGE_LEN, MAX_CURVATURE); + //refiner = new GenRegularRefinement<>(triangulation, predicate, level); refiner = new GenRegularRefinement<>(clone, predicate, level); refiner.refine(); //refiner.coarse(); @@ -83,19 +89,25 @@ public class MeshEikonalSolverFMMIterative list = solver.getTriangulation().getMesh().getBoundaryVertices(); solver.solve(); + System.out.println(solver.getTriangulation().getMesh().toPythonTriangulation(v -> solver.getPotential(v))); level = 1; double lastMaxCurvature = 0; maxCurrentCurvature = Double.MAX_VALUE; - while(maxCurrentCurvature > MAX_CURVATURE && !hasConverged(lastMaxCurvature, maxCurrentCurvature) && level < MAX_ITERATIONS) { + while (level < MAX_ITERATIONS) { + //while(maxCurrentCurvature > MAX_CURVATURE && !hasConverged(lastMaxCurvature, maxCurrentCurvature) && level < MAX_ITERATIONS) { lastMaxCurvature = maxCurrentCurvature; //System.out.println(solver.getTriangulation().getMesh().toPythonTriangulation(v -> solver.getPotential(v))); maxCurrentCurvature = computeCurvature(solver.getTriangulation()); logger.debug("max curvature = " + maxCurrentCurvature); + solver.getTriangulation().removeTriEventListener(solver); IIncrementalTriangulation refinedTriangulation = refine(solver.getTriangulation()); solver = new MeshEikonalSolverFMM<>(solver, refinedTriangulation, refinedTriangulation.getMesh().getBoundaryVertices()); solver.solve(); + + System.out.println(solver.getTriangulation().getMesh().toPythonTriangulation(v -> solver.getPotential(v))); level++; } diff --git a/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateEdgeRefinement.java b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateEdgeRefinement.java new file mode 100644 index 000000000..14af5bb13 --- /dev/null +++ b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateEdgeRefinement.java @@ -0,0 +1,46 @@ +package org.vadere.simulator.models.potential.solver.calculators.mesh; + +import org.jetbrains.annotations.NotNull; +import org.vadere.meshing.mesh.inter.IFace; +import org.vadere.meshing.mesh.inter.IHalfEdge; +import org.vadere.meshing.mesh.inter.IVertex; +import org.vadere.meshing.utils.math.GeometryUtilsMesh; +import org.vadere.util.geometry.shapes.VLine; +import java.util.function.Predicate; + +public class PredicateEdgeRefinement implements Predicate { + private final MeshEikonalSolverFMM solver; + private double minEdgeLen; + private double maxCurvature; + + public PredicateEdgeRefinement( + @NotNull final MeshEikonalSolverFMM solver, + final double minEdgeLen, + final double maxCurvature + ){ + this.solver = solver; + this.minEdgeLen = minEdgeLen; + this.maxCurvature = maxCurvature; + } + + @Override + public boolean test(E e) { + V vertex = solver.getMesh().getVertex(e); + V vertexTwin = solver.getMesh().getTwinVertex(e); + VLine line = solver.getMesh().toLine(e); + double len = line.length(); + double maxCurvature = 0.0; + if(solver.isBurned(solver.getMesh().getVertex(e)) && solver.isBurned(solver.getMesh().getTwinVertex(e))) { + for (E edge : solver.getMesh().getEdgeIt(vertex)) { + double[] curvatures = GeometryUtilsMesh.curvature(solver.getMesh(), v -> solver.getPotential(v), v -> solver.isBurned(v), edge); + maxCurvature = Math.max(maxCurvature, curvatures[1]); + } + + for (E edge : solver.getMesh().getEdgeIt(vertexTwin)) { + double[] curvatures = GeometryUtilsMesh.curvature(solver.getMesh(), v -> solver.getPotential(v), v -> solver.isBurned(v), edge); + maxCurvature = Math.max(maxCurvature, curvatures[1]); + } + } + return len > minEdgeLen && maxCurvature > this.maxCurvature; + } +} diff --git a/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateRefinement.java b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateRefinement.java index 5b480cdcb..985797a60 100644 --- a/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateRefinement.java +++ b/VadereSimulator/src/org/vadere/simulator/models/potential/solver/calculators/mesh/PredicateRefinement.java @@ -5,6 +5,7 @@ import org.vadere.meshing.mesh.inter.IFace; import org.vadere.meshing.mesh.inter.IHalfEdge; import org.vadere.meshing.mesh.inter.IIncrementalTriangulation; import org.vadere.meshing.mesh.inter.IVertex; +import org.vadere.meshing.utils.math.GeometryUtilsMesh; import org.vadere.util.geometry.GeometryUtils; import org.vadere.util.geometry.shapes.VLine; import org.vadere.util.geometry.shapes.VPoint; @@ -33,6 +34,7 @@ public class PredicateRefinement -topography.distanceToObstacle(p) ); - TimeCostPedestrianDensityMesh timeCostPedestrianDensityMesh = new TimeCostPedestrianDensityMesh<>( + TimeCostGaussianPedestrianDensityMesh timeCostPedestrianDensityMesh = new TimeCostGaussianPedestrianDensityMesh<>( timeCostObstacleDensity, triangulation, loadingStrategy, diff --git a/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostGaussianPedestrianDensityMesh.java b/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostGaussianPedestrianDensityMesh.java new file mode 100644 index 000000000..6b53f3be2 --- /dev/null +++ b/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostGaussianPedestrianDensityMesh.java @@ -0,0 +1,197 @@ +package org.vadere.simulator.models.potential.timeCostFunction; + +import org.jetbrains.annotations.NotNull; +import org.jetbrains.annotations.Nullable; +import org.vadere.meshing.mesh.inter.IFace; +import org.vadere.meshing.mesh.inter.IHalfEdge; +import org.vadere.meshing.mesh.inter.IIncrementalTriangulation; +import org.vadere.meshing.mesh.inter.IVertex; +import org.vadere.meshing.mesh.inter.IVertexContainerDouble; +import org.vadere.meshing.mesh.triangulation.triangulator.gen.GenRegularRefinement; +import org.vadere.meshing.utils.io.IOUtils; +import org.vadere.meshing.utils.math.GeometryUtilsMesh; +import org.vadere.simulator.models.potential.solver.timecost.ITimeCostFunctionMesh; +import org.vadere.simulator.models.potential.timeCostFunction.loading.IPedestrianLoadingStrategy; +import org.vadere.state.attributes.scenario.AttributesAgent; +import org.vadere.state.scenario.Pedestrian; +import org.vadere.state.scenario.Topography; +import org.vadere.util.geometry.GeometryUtils; +import org.vadere.util.geometry.shapes.IPoint; +import org.vadere.util.geometry.shapes.VTriangle; + +import java.io.BufferedWriter; +import java.io.File; +import java.io.IOException; +import java.util.Optional; +import java.util.Set; +import java.util.function.Predicate; + +/** + * A TimeCostFunction which reduces the travelling speed decreases with the obstacle density. + * + * @param + * @param + * @param + */ +public class TimeCostGaussianPedestrianDensityMesh implements ITimeCostFunctionMesh { + + public static final String nameAgentDensity = "agent_density"; + + private final ITimeCostFunctionMesh timeCostFunction; + private final Topography topography; + private final IIncrementalTriangulation triangulation; + private final IVertexContainerDouble densities; + private final IPedestrianLoadingStrategy loadingStrategy; + private boolean updated; + private GenRegularRefinement refiner; + + private final double R = 0.7; + private final int influenceRadius = 7; + private final double a; + private final double Sp; + private final double c; + + // to debug + //private MeshPanel debugPanel; + private int step; + private BufferedWriter meshWriter; + + public TimeCostGaussianPedestrianDensityMesh( + @NotNull final ITimeCostFunctionMesh timeCostFunction, + @NotNull final IIncrementalTriangulation triangulation, + @NotNull final IPedestrianLoadingStrategy loadingStrategy, + final AttributesAgent attributesAgent, + final Topography topography){ + this.timeCostFunction = timeCostFunction; + this.loadingStrategy = loadingStrategy; + this.triangulation = triangulation; + this.topography = topography; + this.densities = triangulation.getMesh().getDoubleVertexContainer(nameAgentDensity); + this.updated = false; + this.refiner = new GenRegularRefinement<>(triangulation, e -> false); + this.refiner.setCoarsePredicate(v -> coarse(v)); + this.refiner.setEdgeRefinementPredicate(e -> refine(e)); + this.step = 0; + //TODO duplicated code + double dia = attributesAgent.getRadius() * 2.0; + Sp = (dia * dia * Math.sqrt(3)) * 0.5; + a = -1 / (2 * R * R); + c = 2 * Math.PI * R * R; + /*try { + this.meshWriter = IOUtils.getWriter("floorField_densities.txt", new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/")); + } catch (IOException e) { + e.printStackTrace(); + }*/ + + //debugPanel = new MeshPanel<>(triangulation.getMesh(), 1000, 1000); + //debugPanel.display("debug"); + } + + private void refineMesh() { + long ms = System.currentTimeMillis(); + refiner.coarse(); + //debugPanel.paintImmediately(0, 0, debugPanel.getWidth(), debugPanel.getHeight()); + refiner.refine(); + long runTime = System.currentTimeMillis() - ms; + refiner.getMesh().garbageCollection(); + //debugPanel.paintImmediately(0, 0, debugPanel.getWidth(), debugPanel.getHeight()); + //System.out.println("runTime refinement = " + runTime); + } + + private boolean coarse(@NotNull final V vertex) { + //return true; + for(Pedestrian pedestrian : topography.getPedestrianDynamicElements().getElements()) { + if(pedestrian.getPosition().distanceSq(triangulation.getMesh().toPoint(vertex)) > influenceRadius * influenceRadius) { + return true; + } + } + return false; + } + + private boolean refine(@NotNull final E e) { + //return refiner.getLevel(e) < 2; + if(!triangulation.getMesh().isBoundary(e)) { + VTriangle triangle = triangulation.getMesh().toTriangle(triangulation.getMesh().getFace(e)); + if(/*!refiner.isGreen(e) || */triangulation.getMesh().toLine(e).length() > 1.0) { + for(Pedestrian pedestrian : topography.getPedestrianDynamicElements().getElements()) { + if(pedestrian.getPosition().distanceSq(triangle.midPoint()) < influenceRadius * influenceRadius) { + return true; + } + } + } + } + return false; + } + + private double density(final double x1, final double y1, final double x2, final double y2, @Nullable final Pedestrian ped) { + double dist = GeometryUtils.lengthSq(x1 - x2, y1 - y2); + double density = loadingStrategy.calculateLoading(ped) * (Sp / (c)) * Math.exp(a * dist); + return density; + } + + @Override + public double costAt(@NotNull final V v, @Nullable final Object caller) { + return timeCostFunction.costAt(v) + densities.getValue(v); + } + + @Override + public double costAt(V v) { + return costAt(v, null); + } + + @Override + public double costAt(@NotNull final IPoint p) { + F face = triangulation.locate(p).get(); + double cost = 0; + if (!triangulation.getMesh().isBoundary(face)) { + cost = GeometryUtilsMesh.barycentricInterpolation(face, triangulation.getMesh(), + v -> densities.getValue(v), p.getX(), p.getY()); + } + return timeCostFunction.costAt(p) + cost; + } + + @Override + public boolean needsUpdate() { + return true; + } + + @Override + public void update() { + timeCostFunction.update(); + if(step % 10 == 0) { + //refineMesh(); + } + step++; + + //long ms = System.currentTimeMillis(); + + + var mesh = triangulation.getMesh(); + densities.reset(); + + for (Pedestrian element : topography.getPedestrianDynamicElements().getElements()) { + Optional optional = triangulation.locateFace(element.getPosition(), element); + assert optional.isPresent(); + + if (optional.isPresent()) { + F pedFace = optional.get(); + Predicate predicate = v -> GeometryUtils.lengthSq(mesh.getX(v) - element.getPosition().x, + mesh.getY(v) - element.getPosition().y) < influenceRadius * influenceRadius; + Set closeVertices = triangulation.getVertices(element.getPosition().getX(), element.getPosition().getY(), pedFace, predicate); + for (V v : closeVertices) { + double density = densities.getValue(v) + density(element.getPosition().x, element.getPosition().y, mesh.getX(v), mesh.getY(v), element); + densities.setValue(v, density); + } + } + } + + /*try { + this.meshWriter.write(triangulation.getMesh().toPythonTriangulation(v -> densities.getValue(v) + triangulation.getMesh().getDoubleData(v, TimeCostObstacleDensityMesh.nameObstacleDensity))); + this.meshWriter.flush(); + } catch (IOException e) { + e.printStackTrace(); + }*/ + //long runTime = System.currentTimeMillis() - ms; + //System.out.println("runTime of density computation = " + runTime); + } +} diff --git a/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostPedestrianDensityMesh.java b/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostPedestrianDensityMesh.java index 6040ad61c..962c681d9 100644 --- a/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostPedestrianDensityMesh.java +++ b/VadereSimulator/src/org/vadere/simulator/models/potential/timeCostFunction/TimeCostPedestrianDensityMesh.java @@ -2,18 +2,14 @@ package org.vadere.simulator.models.potential.timeCostFunction; import org.jetbrains.annotations.NotNull; import org.jetbrains.annotations.Nullable; -import org.vadere.meshing.mesh.gen.MeshPanel; -import org.vadere.meshing.mesh.gen.MeshRenderer; -import org.vadere.meshing.mesh.gen.PVertex; import org.vadere.meshing.mesh.inter.IFace; import org.vadere.meshing.mesh.inter.IHalfEdge; import org.vadere.meshing.mesh.inter.IIncrementalTriangulation; import org.vadere.meshing.mesh.inter.IVertex; import org.vadere.meshing.mesh.inter.IVertexContainerDouble; import org.vadere.meshing.mesh.triangulation.triangulator.gen.GenRegularRefinement; -import org.vadere.meshing.utils.color.Colors; +import org.vadere.meshing.utils.io.IOUtils; import org.vadere.meshing.utils.math.GeometryUtilsMesh; -import org.vadere.simulator.models.potential.solver.timecost.ITimeCostFunction; import org.vadere.simulator.models.potential.solver.timecost.ITimeCostFunctionMesh; import org.vadere.simulator.models.potential.timeCostFunction.loading.IPedestrianLoadingStrategy; import org.vadere.state.attributes.scenario.AttributesAgent; @@ -23,7 +19,9 @@ import org.vadere.util.geometry.GeometryUtils; import org.vadere.util.geometry.shapes.IPoint; import org.vadere.util.geometry.shapes.VTriangle; -import java.awt.*; +import java.io.BufferedWriter; +import java.io.File; +import java.io.IOException; import java.util.Optional; import java.util.Set; import java.util.function.Predicate; @@ -37,7 +35,7 @@ import java.util.function.Predicate; */ public class TimeCostPedestrianDensityMesh implements ITimeCostFunctionMesh { - public static final String nameObstacleDensity = "agent_density"; + public static final String nameAgentDensity = "agent_density"; private final ITimeCostFunctionMesh timeCostFunction; private final Topography topography; @@ -47,15 +45,12 @@ public class TimeCostPedestrianDensityMesh refiner; - private final double R = 0.7; - private final int influenceRadius = 5; - private final double a; - private final double Sp; - private final double c; + private final double influenceRadius = 7.0; // to debug //private MeshPanel debugPanel; private int step; + private BufferedWriter meshWriter; public TimeCostPedestrianDensityMesh( @NotNull final ITimeCostFunctionMesh timeCostFunction, @@ -67,17 +62,20 @@ public class TimeCostPedestrianDensityMesh(triangulation, e -> false); this.refiner.setCoarsePredicate(v -> coarse(v)); this.refiner.setEdgeRefinementPredicate(e -> refine(e)); this.step = 0; //TODO duplicated code - double dia = attributesAgent.getRadius() * 2.0; - Sp = (dia * dia * Math.sqrt(3)) * 0.5; - a = -1 / (2 * R * R); - c = 2 * Math.PI * R * R; + + try { + this.meshWriter = IOUtils.getWriter("floorField_densities.txt", new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/")); + } catch (IOException e) { + e.printStackTrace(); + } + //debugPanel = new MeshPanel<>(triangulation.getMesh(), 1000, 1000); //debugPanel.display("debug"); } @@ -107,7 +105,7 @@ public class TimeCostPedestrianDensityMesh 5) { + if(/*!refiner.isGreen(e) || */triangulation.getMesh().toLine(e).length() > 1.0) { for(Pedestrian pedestrian : topography.getPedestrianDynamicElements().getElements()) { if(pedestrian.getPosition().distanceSq(triangle.midPoint()) < influenceRadius * influenceRadius) { return true; @@ -118,12 +116,6 @@ public class TimeCostPedestrianDensityMesh closeVertices = triangulation.getVertices(element.getPosition().getX(), element.getPosition().getY(), pedFace, predicate); for (V v : closeVertices) { - double density = densities.getValue(v) + density(element.getPosition().x, element.getPosition().y, mesh.getX(v), mesh.getY(v), element); + double density = densities.getValue(v) + loadingStrategy.calculateLoading(element) / (influenceRadius * influenceRadius * Math.PI); densities.setValue(v, density); } } } + try { + this.meshWriter.write(triangulation.getMesh().toPythonValues(v -> densities.getValue(v) + triangulation.getMesh().getDoubleData(v, TimeCostObstacleDensityMesh.nameObstacleDensity))); + this.meshWriter.flush(); + } catch (IOException e) { + e.printStackTrace(); + } //long runTime = System.currentTimeMillis() - ms; //System.out.println("runTime of density computation = " + runTime); } diff --git a/VadereSimulator/src/org/vadere/simulator/projects/migration/GeometryCleaner.java b/VadereSimulator/src/org/vadere/simulator/projects/migration/GeometryCleaner.java index 0492d5160..027d1cb3d 100644 --- a/VadereSimulator/src/org/vadere/simulator/projects/migration/GeometryCleaner.java +++ b/VadereSimulator/src/org/vadere/simulator/projects/migration/GeometryCleaner.java @@ -139,7 +139,7 @@ public class GeometryCleaner { IDistanceFunction distanceFunction = IDistanceFunction.create(bound, polygons); // 3. compute the segment-bounding simple polygon by using the distance function - Predicate removePredicate = face -> distanceFunction.apply(triangulation.getMesh().toTriangle(face).midPoint()) > 0; + Predicate removePredicate = face -> distanceFunction.apply(triangulation.getMesh().toMidpoint(face)) > 0; triangulation.shrinkBorder(removePredicate, true); VPolygon boundingPolygon = GeometryUtils.toPolygon(triangulation.getMesh().getPoints(triangulation.getMesh().getBorder())); @@ -147,7 +147,7 @@ public class GeometryCleaner { List faces = triangulation.getMesh().getFaces(); for(PFace face : faces) { if(!triangulation.getMesh().isBorder(face) && !triangulation.getMesh().isDestroyed(face) && !triangulation.getMesh().isHole(face)) { - triangulation.createHole(face, f -> distanceFunction.apply(triangulation.getMesh().toTriangle(f).midPoint()) > 0, true); + triangulation.createHole(face, f -> distanceFunction.apply(triangulation.getMesh().toMidpoint(f)) > 0, true); } } diff --git a/VadereUtils/src/org/vadere/util/geometry/GeometryUtils.java b/VadereUtils/src/org/vadere/util/geometry/GeometryUtils.java index 116e745e6..a594ad7fb 100644 --- a/VadereUtils/src/org/vadere/util/geometry/GeometryUtils.java +++ b/VadereUtils/src/org/vadere/util/geometry/GeometryUtils.java @@ -226,16 +226,18 @@ public class GeometryUtils { * @return the circumcenter of a triangle */ public static VPoint getCircumcenter(@NotNull final IPoint p1, @NotNull final IPoint p2, @NotNull final IPoint p3) { - assert !p1.equals(p2) && !p1.equals(p3) && !p2.equals(p3); + return getCircumcenter(p1.getX(), p1.getY(), p2.getX(), p2.getY(), p3.getX(), p3.getY()); + } - double d = 2 * (p1.getX() * (p2.getY() - p3.getY()) + p2.getX() * (p3.getY() - p1.getY()) + p3.getX() * (p1.getY() - p2.getY())); - double x = ((p1.getX() * p1.getX() + p1.getY() * p1.getY()) * (p2.getY() - p3.getY()) - + (p2.getX() * p2.getX() + p2.getY() * p2.getY()) * (p3.getY() - p1.getY()) - + (p3.getX() * p3.getX() + p3.getY() * p3.getY()) * (p1.getY() - p2.getY())) / d; - double y = ((p1.getX() * p1.getX() + p1.getY() * p1.getY()) * (p3.getX() - p2.getX()) - + (p2.getX() * p2.getX() + p2.getY() * p2.getY()) * (p1.getX() - p3.getX()) - + (p3.getX() * p3.getX() + p3.getY() * p3.getY()) * (p2.getX() - p1.getX())) / d; + public static VPoint getCircumcenter(final double x1, final double y1, final double x2, final double y2, final double x3, final double y3) { + double d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)); + double x = ((x1 * x1 + y1 * y1) * (y2 - y3) + + (x2 * x2 + y2 * y2) * (y3 - y1) + + (x3 * x3 + y3 * y3) * (y1 - y2)) / d; + double y = ((x1 * x1 + y1 * y1) * (x3 - x2) + + (x2 * x2 + y2 * y2) * (x1 - x3) + + (x3 * x3 + y3 * y3) * (x2 - x1)) / d; return new VPoint(x,y); } -- 2.26.2