Commit a927ca27 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen

add tests to solve the eikonal equation on the unstructured mesh for 2...

add tests to solve the eikonal equation on the unstructured mesh for 2 scenarios (bridge and filled chicken).
parent 9afa2262
Pipeline #167339 passed with stages
in 143 minutes and 52 seconds
# nVertices dimension nAttributes boundaryMarker
8 2 0 0
# vertexId x y
1 0.500000 0.500000
2 34.500000 0.500000
3 26.000000 21.000000
4 26.000000 41.000000
5 9.000000 41.000000
6 9.000000 21.000000
7 0.500000 59.500000
8 34.500000 59.500000
#
# nSegments boundaryMarker
8 0
# lineId vertexId1 vertexId2
1 2 8
2 8 7
3 7 1
4 1 2
5 3 4
6 4 5
7 5 6
8 6 3
#
# nHoles
1
# vertexId x y (of a vertex which lies inside the hole)
1 17.500000 31.000000
\ No newline at end of file
......@@ -15,11 +15,12 @@ public class Parameters {
public final static boolean uniform = false;
public final static String method = "Distmesh"; // "Distmesh" or "Density"
public final static double qualityMeasurement = 0.85;
public final static double qualityConvergence = 0.001;
public final static double MINIMUM = 0.25;
public final static double DENSITYWEIGHT = 2;
public final static int NPOINTS = 100000;
public final static int SAMPLENUMBER = 10;
public final static int SAMPLEDIVISION = 10;
public final static int SEGMENTDIVISION = 0;
public final static int MAX_NUMBER_OF_STEPS = 100;
public final static int MAX_NUMBER_OF_STEPS = 250;
}
......@@ -57,7 +57,7 @@ public class GenEikMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
private Collection<? extends IPoint> fixPoints;
private double scalingFactor;
private double deps;
private static final int MAX_STEPS = 250;
private static final int MAX_STEPS = Parameters.MAX_NUMBER_OF_STEPS;
private int nSteps;
private double initialEdgeLen;
private double maxMovement;
......@@ -301,6 +301,10 @@ public class GenEikMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
}
public Collection<V> getFixVertices() {
return refiner.getFixPoints();
}
public void initialize() {
initializeStep();
}
......@@ -343,8 +347,8 @@ public class GenEikMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
public boolean isFinished() {
synchronized (getMesh()) {
boolean converged = dQuality < 0.000000001;
return initializationFinished() && quality > Parameters.qualityMeasurement && converged || (maxMovement > 0 && maxMovement / initialEdgeLen < Parameters.DPTOL) || nSteps >= MAX_STEPS;
boolean converged = dQuality < Parameters.qualityConvergence;
return initializationFinished() && quality >= Parameters.qualityMeasurement && converged || (maxMovement > 0 && maxMovement / initialEdgeLen < Parameters.DPTOL) || nSteps >= MAX_STEPS;
}
}
......
......@@ -116,11 +116,42 @@ public class EikonalSolverFMMTriangulation<V extends IVertex, E extends IHalfEdg
F face = triangulation.locateFace(point.getX(), point.getY()).get();
if(!getMesh().isBoundary(face)) {
targetVertices.addAll(getMesh().getVertices(face));
for(F neighbourFace : getMesh().getFaceIt(face)) {
if(!getMesh().isBoundary(neighbourFace)) {
targetVertices.addAll(getMesh().getVertices(neighbourFace));
}
}
}
//initialFace(face, p -> point.distance(p));
}
}
/**
* Constructor for certain target points.
*
* @param timeCostFunction the time cost function t(x). Note F(x) = 1 / t(x).
* @param triangulation the triangulation the propagating wave moves on.
* @param targetVertices Points where the propagating wave starts i.e. points that are part of the target area.
*/
public EikonalSolverFMMTriangulation(@NotNull final ITimeCostFunction timeCostFunction,
@NotNull final IIncrementalTriangulation<V, E, F> triangulation,
@NotNull final Collection<V> targetVertices
) {
this.triangulation = triangulation;
this.calculationFinished = false;
this.timeCostFunction = timeCostFunction;
this.narrowBand = new PriorityQueue<>(pointComparator);
this.targetVertices = new HashSet<>();
this.distFunc = p -> IDistanceFunction.createToTargetPoints(targetVertices).apply(p);
for(V vertex : targetVertices) {
this.targetVertices.add(vertex);
for(V neighbouringVertices : getMesh().getAdjacentVertexIt(vertex)) {
this.targetVertices.add(neighbouringVertices);
}
}
}
/**
* Constructor for certain target shapes.
*
......
package org.vadere.simulator.models.potential.solver;
import org.jetbrains.annotations.NotNull;
import org.junit.Ignore;
import org.junit.Test;
import org.vadere.meshing.examples.MeshExamples;
......@@ -8,19 +9,31 @@ import org.vadere.meshing.mesh.gen.PFace;
import org.vadere.meshing.mesh.gen.PHalfEdge;
import org.vadere.meshing.mesh.gen.PMesh;
import org.vadere.meshing.mesh.gen.PVertex;
import org.vadere.meshing.mesh.impl.PMeshPanel;
import org.vadere.meshing.mesh.impl.PSLG;
import org.vadere.meshing.mesh.inter.IIncrementalTriangulation;
import org.vadere.meshing.mesh.triangulation.IEdgeLengthFunction;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.gen.GenEikMesh;
import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.utils.io.poly.MeshPSLGWriter;
import org.vadere.meshing.utils.io.poly.MeshPolyReader;
import org.vadere.meshing.utils.io.poly.MeshPolyWriter;
import org.vadere.meshing.utils.io.poly.PSLGGenerator;
import org.vadere.simulator.models.potential.solver.calculators.EikonalSolver;
import org.vadere.simulator.models.potential.solver.calculators.mesh.EikonalSolverFMMTriangulation;
import org.vadere.simulator.models.potential.solver.timecost.UnitTimeCostFunction;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.logging.Logger;
import org.vadere.util.math.IDistanceFunction;
import java.io.IOException;
import java.io.InputStream;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.stream.Collectors;
public class TestFMMEikMesh {
private static Logger log = Logger.getLogger(TestFMMEikMesh.class);
......@@ -57,4 +70,63 @@ public class TestFMMEikMesh {
System.out.println(mesh.toPythonTriangulation(v -> triangulation.getMesh().getDoubleData(v, "potential")));
}
@Ignore
@Test
public void testFilledChickenFMM() throws IOException {
testTriangulationFMM("/poly/filled_chicken.poly", new VPoint(2,2), 1.0);
}
@Ignore
@Test
public void testBridge() throws IOException {
testTriangulationFMM("/poly/bridge.poly", new VPoint(42.0, 46.0), 2.0);
}
private void testTriangulationFMM(@NotNull final String file, @NotNull final VPoint targetPoint, final double h0) throws IOException {
// 1. read the base PSLG-file
final InputStream inputStream = MeshExamples.class.getResourceAsStream(file);
PSLG pslg = PSLGGenerator.toPSLGtoVShapes(inputStream);
// 2. generate the unstructured mesh using EikMesh
Collection<VPoint> fixPoints = new ArrayList<>();
fixPoints.add(targetPoint);
var eikMesh = new GenEikMesh<>(
IDistanceFunction.create(pslg.getSegmentBound(), pslg.getHoles()),
p -> h0,
fixPoints,
h0,
pslg.getBoundingBox(),
pslg.getAllPolygons(),
() -> new PMesh());
PMeshPanel panel = new PMeshPanel(eikMesh.getMesh(), 600, 800);
panel.display();
// long for: eikMesh.generate();
while (!eikMesh.isFinished()) {
eikMesh.improve();
panel.repaint();
}
//eikMesh.generate();
// 3. solve the eikonal equation on the given mehsh
EikonalSolver solver = new EikonalSolverFMMTriangulation(
new UnitTimeCostFunction(),
eikMesh.getTriangulation(),
eikMesh.getFixVertices().stream().filter(v -> v.distance(targetPoint) < GeometryUtils.DOUBLE_EPS).collect(Collectors.toList()));
long ms = System.currentTimeMillis();
log.info("start FFM");
solver.initialize();
log.info("FFM finished");
log.info("time: " + (System.currentTimeMillis() - ms));
// 4. print the result to the console i.e. standard out
MeshPolyWriter<PVertex, PHalfEdge, PFace> meshPolyWriter = new MeshPolyWriter<>();
System.out.println(meshPolyWriter.to2DPoly(eikMesh.getMesh(), 1, i -> "potential", v -> false));
//System.out.println(eikMesh.getMesh().toPythonTriangulation(v -> eikMesh.getMesh().getDoubleData(v, "potential")));
}
}
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