Commit d46f4882 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen

small performance improvements of the EikMesh algorithm.

parent 1b2630ed
Pipeline #271855 failed with stages
in 2 minutes and 10 seconds
......@@ -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)));
......
......@@ -321,7 +321,7 @@ public class IncrementalTriangulation<V extends IVertex, E extends IHalfEdge, F
// remove all faces outside the hole
VPolygon polygon = getMesh().toPolygon(hole);
Predicate<F> removePredicate = face -> !polygon.contains(getMesh().toTriangle(face).midPoint());
Predicate<F> removePredicate = face -> !polygon.contains(getMesh().toMidpoint(face));
cdt.getTriangulation().shrinkBorder(removePredicate, true);
IMesh<V, E, F> holeMesh = incrementalTriangulation.getMesh();
......
......@@ -545,6 +545,9 @@ public interface IMesh<V extends IVertex, E extends IHalfEdge, F extends IFace>
* @return (optional) a boundary edge
*/
default Optional<E> 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<V extends IVertex, E extends IHalfEdge, F extends IFace>
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<V extends IVertex, E extends IHalfEdge, F extends IFace>
}
default String toPythonValues(@NotNull final Function<V, Double> evalPoint) {
StringBuilder builder = new StringBuilder();
List<V> 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.
......
......@@ -116,7 +116,7 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
* @param vertex v
*/
default void adjustVertex(@NotNull final V vertex) {
getMesh().streamEdges(vertex).filter(edge -> 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<V extends IVertex, E extends IHalfEdge, F ext
* @param face they face which will be transformed into a hole
* @param mergeCondition the merge condition
* @param deleteIsoletedVertices if true isolated vertices, i.e. vertices without any edges, will be removed from the mesh
* @param vertexAdjust true means that boundary vertices will get their boundary edge as edge, false means there is no guarantee that this adjustment is made
* @return (optional) the hole or face itself it the face does not fulfill the merge condition
* or empty if due to the creation of the hole all faces will be removed!
*/
default Optional<F> createHole(@NotNull final F face, @NotNull final Predicate<F> mergeCondition, final boolean deleteIsoletedVertices) {
default Optional<F> createHole(@NotNull final F face, @NotNull final Predicate<F> 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<V extends IVertex, E extends IHalfEdge, F ext
}*/
}
default Optional<F> createHole(@NotNull final F face, @NotNull final Predicate<F> 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<V extends IVertex, E extends IHalfEdge, F ext
*
* @param removeCondition the remove condition
* @param deleteIsolatedVertices true then isolated vertices (they are not connected to an edge) will be removed.
* @param vertexAdjust true means that boundary vertices will get their boundary edge as edge, false means there is no guarantee that this adjustment is made
*/
default void shrinkBorder(final Predicate<F> removeCondition, final boolean deleteIsolatedVertices, final boolean vertexAdjust) {
shrinkBoundary(getMesh().getBorder(), removeCondition, deleteIsolatedVertices, vertexAdjust);
}
default void shrinkBorder(final Predicate<F> removeCondition, final boolean deleteIsolatedVertices) {
shrinkBoundary(getMesh().getBorder(), removeCondition, deleteIsolatedVertices);
shrinkBorder(removeCondition, deleteIsolatedVertices, true);
}
default void shrinkBoundary(final Predicate<F> removeCondition, final boolean deleteIsolatedVertices) {
......@@ -697,8 +707,9 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
*
* @param removeCondition the remove condition
* @param deleteIsolatedVertices true then isolated vertices (they are not connected to an edge) will be removed.
* @param adjustVertices true means that boundary vertices will get their boundary edge as edge, false means there is no guarantee that this adjustment is made
*/
default void shrinkBoundary(@NotNull final F boundary, final Predicate<F> removeCondition, final boolean deleteIsolatedVertices) {
default void shrinkBoundary(@NotNull final F boundary, final Predicate<F> removeCondition, final boolean deleteIsolatedVertices, final boolean adjustVertices) {
assert getMesh().isBoundary(boundary);
List<F> boundaryFaces = getMesh().getFaces(boundary);
......@@ -715,13 +726,17 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
nextNeighbouringFaces.add(face);
}
}
removeFaceAtBoundary(neighbouringFace, boundary, deleteIsolatedVertices);
removeFaceAtBoundary(neighbouringFace, boundary, deleteIsolatedVertices, adjustVertices);
}
}
neighbouringFaces = nextNeighbouringFaces;
} while (!neighbouringFaces.isEmpty());
}
default void shrinkBoundary(@NotNull final F boundary, final Predicate<F> removeCondition, final boolean deleteIsolatedVertices) {
shrinkBoundary(boundary, removeCondition, deleteIsolatedVertices, true);
}
default void removeFacesAtBoundary(@NotNull final Predicate<F> mergePredicate, @NotNull final Predicate<F> errorPredicate) throws IllegalMeshException {
mergeFaces(getMesh().getBorder(), mergePredicate, errorPredicate, true, 1);
for(F face : getMesh().getHoles()) {
......@@ -1215,8 +1230,9 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
* @param face the face that will be removed from the mesh
* @param boundary the boundary which has to be a neighbouring boundary of the face
* @param deleteIsolatedVertices true means that all vertices with degree smaller equals 1 will be removed as well
* @param adjustVertices true means that boundary vertices will get their boundary edge as edge, false means there is no guarantee that this adjustment is made
*/
default void removeFaceAtBoundary(@NotNull final F face, @NotNull final F boundary, final boolean deleteIsolatedVertices) {
default void removeFaceAtBoundary(@NotNull final F face, @NotNull final F boundary, final boolean deleteIsolatedVertices, final boolean adjustVertices) {
if(!getMesh().isDestroyed(face)) {
assert getMesh().streamFaces(face).filter(neighbour -> neighbour.equals(boundary)).count() > 0;
......@@ -1357,9 +1373,11 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
getMesh().destroyEdge(h1);
// adjust vertices such that we speed up the querry isBoundary(vertex).
if(adjustVertices) {
vertices.stream().filter(getMesh()::isAlive).forEach(v -> adjustVertex(v));
}
}
}
if(nEdges > 0) {
getMesh().destroyFace(face);
......@@ -1371,6 +1389,9 @@ public interface IPolyConnectivity<V extends IVertex, E extends IHalfEdge, F ext
}
}
default void removeFaceAtBoundary(@NotNull final F face, @NotNull final F boundary, final boolean deleteIsolatedVertices) {
removeFaceAtBoundary(face, boundary, deleteIsolatedVertices);
}
/**
* Removes a face from the mesh by removing all boundary edges of the face.
......
......@@ -1953,7 +1953,7 @@ public interface ITriConnectivity<V extends IVertex, E extends IHalfEdge, F exte
* @return all visited faces in a first visited first in ordered queue, i.e. <tt>LinkedList</tt>.
*/
default LinkedList<E> straightWalk2DGatherDirectional(@NotNull final F face, @NotNull final VPoint direction, @NotNull final Predicate<E> additionalStopCondition) {
VPoint q = getMesh().toTriangle(face).midPoint();
VPoint q = getMesh().toMidpoint(face);
assert getMesh().toTriangle(face).contains(q);
Predicate<E> defaultStopCondion = e -> isRightOf(q.x, q.y, e);
......
......@@ -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++;
......
......@@ -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;
}
......@@ -26,7 +26,7 @@ public interface IEikMeshImprover<V extends IVertex, E extends IHalfEdge, F exte
}
default Predicate<F> outsidePredicate(@NotNull final IDistanceFunction distanceFunc) {
return f -> distanceFunc.apply(getMesh().toTriangle(f).midPoint()) > 0;
return f -> distanceFunc.apply(getMesh().toMidpoint(f)) > 0;
}
default Predicate<F> lowQualityPredicate() {
......
......@@ -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<AVertex, AHalfEdge, AFace> 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<VShape> 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<AVertex, AHalfEdge, AFace> 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<MeshPoint, AVertex<MeshPoint>, AHalfEdge<MeshPoint>, AFace<MeshPoint>> distmeshPanel = new PSMeshingPanel(meshGenerator.getMesh(), f -> false, 1000, 800, bbox);
MeshPanel<AVertex, AHalfEdge, AFace> 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<AVertex, AHalfEdge, AFace> 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<VShape> 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);
}
}
......@@ -96,10 +96,10 @@ public class GenRegularRefinement<V extends IVertex, E extends IHalfEdge, F exte
//this.edgeRefinementPredicate = e -> 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<>();
......
......@@ -527,6 +527,7 @@ public class GenUniformRefinementTriangulatorSFC<V extends IVertex, E extends IH
List<F> 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<V extends IVertex, E extends IH
* Note the border is part of the whole boundary which is defined by the border and the holes.</p>
*/
private void shrinkBorder() {
Predicate<F> removePredicate = face -> distFunc.apply(triangulation.getMesh().toTriangle(face).midPoint()) > 0;
triangulation.shrinkBorder(removePredicate, true);
Predicate<F> removePredicate = face -> distFunc.apply(triangulation.getMesh().toMidpoint(face)) > 0;
triangulation.shrinkBorder(removePredicate, true, false);
}
@Override
......@@ -652,6 +653,7 @@ public class GenUniformRefinementTriangulatorSFC<V extends IVertex, E extends IH
GenConstrainSplitter<V, E, F> 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<V extends IVertex, E extends IH
List<F> 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() {
......
......@@ -106,6 +106,62 @@ public class GeometryUtilsMesh {
return new double[]{curvature, gaussianCurvature};
}
public static <V extends IVertex, E extends IHalfEdge, F extends IFace> double[] curvature(@NotNull final IMesh<V, E, F> mesh,
@NotNull final Function<V, Double> f,
@NotNull final Predicate<V> 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 <V extends IVertex, E extends IHalfEdge, F extends IFace> double[] curvature(
@NotNull final IMesh<V, E, F> mesh,
......
......@@ -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<PVertex> list = tri.getMesh().getBoundaryVertices();
Set<PVertex> initialVertices = new HashSet<>();
initialVertices.addAll(list);
list.stream().forEach(v -> curTri.getMesh().getAdjacentVertexIt(v).forEach(u -> initialVertices.add(u)));
Set<PVertex> 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<PHalfEdge> 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();
......
......@@ -228,7 +228,8 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
// this might happen for vertices which are close at the initialVertex!
if(list.isEmpty()) {
//logger.warn("could not find virtual support for non-acute triangle.");
potential = computePotential(edge, next, prev, -1.0);
//potential = computePotential(edge, next, prev, -1.0);
//potential = computePotential(edge, next, prev, getCosPhi(edge));
} else {
DoubleArrayList cosPhis = getVirtualSupportCosPhi(edge);
......
......@@ -17,6 +17,7 @@ import org.vadere.util.logging.Logger;
import org.vadere.util.math.IDistanceFunction;
import java.util.Collection;
import java.util.List;
import java.util.function.Function;
import java.util.function.Predicate;
......@@ -36,7 +37,7 @@ public class MeshEikonalSolverFMMIterative<V extends IVertex, E extends IHalfEdg
private static final double MIN_CURVATURE = 0.1;
private static final double MAX_CURVATURE = 0.1;
private static final double MAX_CURVATURE = 0.3;
private static final double MIN_EDGE_LEN = 0.1;
......@@ -46,6 +47,8 @@ public class MeshEikonalSolverFMMIterative<V extends IVertex, E extends IHalfEdg
private GenRegularRefinement<V, E, F> refiner;
private @NotNull IDistanceFunction distFunc;
public MeshEikonalSolverFMMIterative(@NotNull final ITimeCostFunction timeCostFunction,
@NotNull final IIncrementalTriangulation<V, E, F> triangulation,
......@@ -53,6 +56,7 @@ public class MeshEikonalSolverFMMIterative<V extends IVertex, E extends IHalfEdg
@NotNull final IDistanceFunction distFunc
) {
this(new MeshEikonalSolverFMM<>(timeCostFunction, triangulation, targetVertices, distFunc));
this.distFunc = distFunc;
}
public MeshEikonalSolverFMMIterative(@NotNull final MeshEikonalSolverFMM<V, E, F> solver) {
......@@ -66,6 +70,8 @@ public class MeshEikonalSolverFMMIterative<V extends IVertex, E extends IHalfEdg
// (2) refine if necessary
final Predicate<E> predicate = new PredicateRefinement<>(triangulation, clone, MIN_EDGE_LEN, MAX_CURVATURE);
//final Predicate<E> 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<V extends IVertex, E extends IHalfEdg
public void solve() {
if (!calculationFinished) {
List<V> 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());