Commit bb1a9437 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

small bug fix for the mesh, tex mesh generator

parent 9b5f4ae8
......@@ -114,6 +114,7 @@ public class IncrementalTriangulation<P extends IPoint, V extends IVertex<P>, E
p1 = mesh.insertVertex(bound.getX() + 2 * max + gap, bound.getY() - gap);
p2 = mesh.insertVertex(bound.getX() + (max+2*gap)/2, bound.getY() + 2 * max + gap);
// counter clockwise!
superTriangle = mesh.createFace(p0, p1, p2);
borderFace = mesh.getTwinFace(mesh.getEdge(superTriangle));
......
......@@ -19,6 +19,7 @@ import java.util.stream.Stream;
public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P>, PFace<P>> {
private List<PFace<P>> faces;
//private List<PFace<P>> borderFaces;
private PFace<P> boundary;
private List<PHalfEdge<P>> edges;
private IPointConstructor<P> pointConstructor;
......@@ -26,6 +27,7 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
public PMesh(final IPointConstructor<P> pointConstructor) {
this.faces = new ArrayList<>();
//this.borderFaces = new ArrayList<>();
this.edges = new ArrayList<>();
this.vertices = new HashSet<>();
this.boundary = new PFace<>(true);
......@@ -290,4 +292,9 @@ public class PMesh<P extends IPoint> implements IMesh<P, PVertex<P>, PHalfEdge<P
public List<PFace<P>> getFaces() {
return streamFaces().filter(face -> !face.isBorder()).filter(face -> !face.isDestroyed()).collect(Collectors.toList());
}
@Override
public List<PHalfEdge<P>> getBoundaryEdges() {
return streamEdges().filter(edge -> edge.isBoundary()).filter(edge -> !edge.isDestroyed()).collect(Collectors.toList());
}
}
......@@ -22,6 +22,7 @@ import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.Random;
import java.util.Set;
/**
......@@ -40,6 +41,7 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
private IMesh<P, V, E, F> mesh;
private static final Logger logger = LogManager.getLogger(UniformRefinementTriangulation.class);
private final IDistanceFunction distFunc;
private final static Random random = new Random();
public UniformRefinementTriangulation(
final ITriangulation<P, V, E, F> triangulation,
......@@ -64,7 +66,7 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
for(E edge : mesh.getEdgeIt(mesh.getBoundary())) {
if(!isCompleted(edge) && !points.contains(mesh.getVertex(edge))) {
toRefineEdges.add(edge);
toRefineEdges.add(mesh.getTwin(edge));
}
}
......@@ -142,17 +144,20 @@ public class UniformRefinementTriangulation<P extends IPoint, V extends IVertex<
return shape.intersects(line) || shape.contains(line.getP1()) || shape.contains(line.getP2());
}
private boolean isCompleted(final E edge) {
private boolean isCompleted(E edge) {
if(mesh.isBoundary(edge)){
edge = mesh.getTwin(edge);
}
F face = mesh.getFace(edge);
F twin = mesh.getTwinFace(edge);
VTriangle triangle = mesh.toTriangle(face);
VTriangle twinTriangle = mesh.toTriangle(twin);
VLine line = mesh.toLine(edge);
return line.length() <= lenFunc.apply(line.midPoint())
|| (!triangle.intersect(bbox) && !twinTriangle.intersect(bbox))
|| boundary.stream().anyMatch(shape -> shape.contains(triangle.getBounds2D()) || shape.contains(twinTriangle.getBounds2D()));
return (line.length() <= lenFunc.apply(line.midPoint()) && random.nextDouble() < 0.96)
|| (!triangle.intersect(bbox) && (mesh.isBoundary(twin) || !mesh.toTriangle(twin).intersect(bbox)))
|| boundary.stream().anyMatch(shape -> shape.contains(triangle.getBounds2D()) || (!mesh.isBoundary(twin) && shape.contains(mesh.toTriangle(twin).getBounds2D())));
}
private Collection<E> refine(final E edge) {
......
......@@ -29,6 +29,7 @@ import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.Set;
import java.util.stream.Collectors;
import java.util.stream.Stream;
import java.util.stream.StreamSupport;
......@@ -54,6 +55,10 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
return new VLine(new VPoint(getVertex(getPrev(halfEdge))), new VPoint(getVertex(halfEdge)));
}
default VPoint toPoint(@NotNull V vertex) {
return new VPoint(new VPoint(vertex));
}
E getEdge(@NotNull V vertex);
E getEdge(@NotNull F face);
......@@ -86,6 +91,27 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
* @return
*/
boolean isBoundary(@NotNull F face);
default boolean isNeighbourBoundary(@NotNull F face){
for(F neighbourFace : getFaceIt(face)) {
if(isBoundary(neighbourFace)) {
return true;
}
}
return false;
}
default Optional<E> getLinkToBoundary(@NotNull F face){
for(E edge : getEdgeIt(face)) {
if(isBoundary(getTwin(edge))) {
return Optional.of(edge);
}
}
return Optional.empty();
}
boolean isBoundary(@NotNull E halfEdge);
boolean isDestroyed(@NotNull F face);
......@@ -181,6 +207,8 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
List<F> getFaces();
List<E> getBoundaryEdges();
Stream<F> streamFaces();
Stream<E> streamEdges();
......@@ -416,6 +444,18 @@ public interface IMesh<P extends IPoint, V extends IVertex<P>, E extends IHalfEd
return edges;
}
default Set<VLine> getLines() {
return streamEdges().map(edge -> toLine(edge)).collect(Collectors.toSet());
}
default Set<VPoint> getUniquePoints() {
return streamVertices().map(vertex -> toPoint(vertex)).collect(Collectors.toSet());
}
default Collection<P> getPoints() {
return streamVertices().map(vertex -> getPoint(vertex)).collect(Collectors.toList());
}
Iterable<E> getEdgeIt();
......
......@@ -269,15 +269,82 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
return remFace;
}
/**
* Removes a face from the mesh by removing all boundary edges of the face.
* If there are no boundary edges the face will be converted to be a part of the boundary
* itself i.e. a hole.
* Mesh changing method.
*
* @param face the face that will be removed from the mesh
* @param deleteIsolatedVertices true means that all vertices with degree <= 1 will be removed as well
*/
default void makeHole(final V vertex) {
List<F> toDeleteFaces = getMesh().getFaces(vertex);
for(F face : toDeleteFaces) {
removeFace(face, true);
}
}
/*default void fill_hole (final V v, final List<E> deletedEdges)
{
// uses the fact that the hole is starshaped
// with repect to v->point()
typedef std::list<Edge> Hole;
Face_handle ff, fn;
int ii , in;
Vertex_handle v0, v1, v2;
Bounded_side side;
//stack algorithm to create faces
// create face v0,v1,v2
//if v0,v1,v2 are finite vertices
// and form a left_turn
// and triangle v0v1v2 does not contain v->point()
if( hole.size() != 3) {
typename Hole::iterator hit = hole.begin();
typename Hole::iterator next= hit;
while( hit != hole.end() && hole.size() != 3) {
ff = (*hit).first;
ii = (*hit).second;
v0 = ff->vertex(cw(ii));
v1 = ff->vertex(ccw(ii));
if( !is_infinite(v0) && !is_infinite(v1)) {
next=hit; next++;
if(next == hole.end()) next=hole.begin();
fn = (*next).first;
in = (*next).second;
v2 = fn->vertex(ccw(in));
if ( !is_infinite(v2) &&
orientation(v0->point(), v1->point(), v2->point()) == LEFT_TURN ) {
side = bounded_side(v0->point(),
v1->point(),
v2->point(),
v->point());
if( side == ON_UNBOUNDED_SIDE ||
(side == ON_BOUNDARY && orientation(v0->point(),
v->point(),
v2->point()) == COLLINEAR &&
collinear_between(v0->point(),v->point(),v2->point()) ))
{
//create face
Face_handle newf = create_face(ff,ii,fn,in);
typename Hole::iterator tempo=hit;
hit = hole.insert(hit,Edge(newf,1)); //push newf
hole.erase(tempo); //erase ff
hole.erase(next); //erase fn
if (hit != hole.begin() ) --hit;
continue;
}
}
}
++hit;
}
}*/
/**
* Removes a face from the mesh by removing all boundary edges of the face.
* If there are no boundary edges the face will be converted to be a part of the boundary
* itself i.e. a hole.
* Mesh changing method.
*
* @param face the face that will be removed from the mesh
* @param deleteIsolatedVertices true means that all vertices with degree <= 1 will be removed as well
*/
default void removeFace(@NotNull final F face, final boolean deleteIsolatedVertices) {
assert !getMesh().isDestroyed(face);
......@@ -286,6 +353,18 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
F boundary = getMesh().createFace(true);
/*for(E edge : getMesh().getEdgeIt(face)) {
if(getMesh().isBoundary(getMesh().getTwin(edge))) {
boundary = getMesh().getTwinFace(edge);
break;
}
}
// no boundary around
if(boundary == null) {
boundary = getMesh().createFace(true);
}*/
for(E edge : getMesh().getEdgeIt(face)) {
getMesh().setFace(edge, boundary);
if(getMesh().isBoundary(getMesh().getTwin(edge))) {
......
......@@ -9,6 +9,7 @@ import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VTriangle;
import java.util.ArrayList;
import java.util.List;
import java.util.Optional;
import java.util.Random;
......@@ -476,15 +477,23 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
@Override
default Optional<F> locateFace(final double x, final double y) {
Optional<F> optFace;
if(getMesh().getNumberOfFaces() > 1) {
return locateFace(x, y, getMesh().getFace());
optFace = locateFace(x, y, getMesh().getFace());
}
else if(getMesh().getNumberOfFaces() == 1) {
return Optional.of(getMesh().getFace());
optFace = Optional.of(getMesh().getFace());
}
else {
optFace = Optional.empty();
}
if(optFace.isPresent() && getMesh().isBoundary(optFace.get())) {
return Optional.empty();
}
else {
return optFace;
}
}
/*default Optional<P> locateVertex(double x, double y, F startFace) {
......@@ -592,6 +601,31 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
return face;
}
default F walkThroughHole(double x1, double y1, @NotNull F borderFace) {
assert getMesh().isBoundary(borderFace);
double minDistance = Double.MAX_VALUE;
F choosenFace = borderFace;
for(E edge : getMesh().getEdgeIt(borderFace)) {
V v1 = getMesh().getVertex(getMesh().getPrev(edge));
V v2 = getMesh().getVertex(edge);
// get out of the hole/border. Note: non-border faces are ccw oriented but border-faces are cw oriented!
if (GeometryUtils.isRightOf(v2, v1, x1, y1)) {
double distance = GeometryUtils.distanceToLineSegment(v1, v2, x1, y1);
if(distance < minDistance) {
minDistance = distance;
choosenFace = getMesh().getTwinFace(edge);
}
}
}
return choosenFace;
}
default F marchRandom2D(double x1, double y1, F startFace) {
boolean first = true;
F face = startFace;
......@@ -599,6 +633,19 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
int count = 0;
while (true) {
if(getMesh().isBoundary(face)) {
F tmpFace = walkThroughHole(x1, y1, face);
if(getMesh().isBoundary(tmpFace)) {
assert tmpFace == face;
return tmpFace;
}
else {
face = tmpFace;
}
}
count++;
boolean goLeft = random.nextBoolean();
//boolean goLeft = true;
......
......@@ -4,6 +4,8 @@ import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.ConstantLineIterator;
import org.vadere.util.geometry.mesh.gen.PFace;
import org.vadere.util.geometry.mesh.gen.PHalfEdge;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.mesh.inter.IHalfEdge;
import org.vadere.util.geometry.mesh.inter.IMesh;
......@@ -22,6 +24,7 @@ import org.vadere.util.potential.PathFindingTag;
import org.vadere.util.potential.timecost.ITimeCostFunction;
import org.vadere.util.geometry.mesh.iterators.FaceIterator;
import org.vadere.util.triangulation.IPointConstructor;
import org.vadere.util.triangulation.adaptive.MeshPoint;
import java.util.ArrayList;
import java.util.Arrays;
......@@ -102,6 +105,31 @@ public class EikonalSolverFMMTriangulation<P extends PotentialPoint, V extends I
}
}
public EikonalSolverFMMTriangulation(final ITimeCostFunction timeCostFunction,
final ITriangulation<P, V, E, F> triangulation,
final Collection<E> targetEdges
) {
this.triangulation = triangulation;
this.calculationFinished = false;
this.timeCostFunction = timeCostFunction;
this.targetAreas = new ArrayList<>();
this.narrowBand = new PriorityQueue<>(pointComparator);
this.mesh = triangulation.getMesh();
for(E halfEdge : targetEdges) {
P potentialPoint = mesh.getPoint(halfEdge);
double distance = 0.0;
if(potentialPoint.getPathFindingTag() != PathFindingTag.Undefined) {
narrowBand.remove(new FFMHalfEdge(halfEdge));
}
potentialPoint.setPotential(Math.min(potentialPoint.getPotential(), distance * timeCostFunction.costAt(potentialPoint)));
potentialPoint.setPathFindingTag(PathFindingTag.Reached);
narrowBand.add(new FFMHalfEdge(halfEdge));
}
}
private void initializeTargetAreas() {
for(VRectangle rectangle : targetAreas) {
......@@ -276,7 +304,7 @@ public class EikonalSolverFMMTriangulation<P extends PotentialPoint, V extends I
private double computeValue(final P point, final F face) {
// check whether the triangle does contain useful data
List<P> points = mesh.getPoints(face);
E halfEdge = mesh.getEdges(face).stream().filter(p -> mesh.getVertex(p).equals(point)).findAny().get();
//E halfEdge = mesh.getEdges(face).stream().filter(p -> mesh.getPoint(p).equals(point)).findAny().get();
points.removeIf(p -> p.equals(point));
assert points.size() == 2;
......
package org.vadere.util.tex;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.mesh.inter.IHalfEdge;
import org.vadere.util.geometry.mesh.inter.IMesh;
import org.vadere.util.geometry.mesh.inter.IVertex;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPoint;
public class TexGraphGenerator {
public static <P extends IPoint, V extends IVertex<P>, E extends IHalfEdge<P>, F extends IFace<P>> String meshToGraph(final IMesh<P, V, E, F> mesh){
StringBuilder builder = new StringBuilder();
builder.append("\\begin{tikzpicture}[scale=1.0]\n");
for(VPoint point : mesh.getUniquePoints()) {
//builder.append("\\draw[fill=black] ("+point.getX()+","+point.getY()+") circle (3pt); \n");
}
builder.append("\\draw ");
for(VLine line : mesh.getLines()) {
builder.append("("+line.getX1()+","+line.getY1()+") -- ("+line.getX2()+","+line.getY2()+")\n");
}
builder.append(";\n");
builder.append("\\end{tikzpicture}");
return builder.toString();
}
}
/*
%% vertices
\draw[fill=black] (0,0) circle (3pt);
\draw[fill=black] (4,0) circle (3pt);
\draw[fill=black] (2,1) circle (3pt);
\draw[fill=black] (2,3) circle (3pt);
%% vertex labels
\node at (-0.5,0) {1};
\node at (4.5,0) {2};
\node at (2.5,1.2) {3};
\node at (2,3.3) {4};
\begin{tikzpicture}[thick,scale=0.8]
% The following path utilizes several useful tricks and features:
% 1) The foreach statement is put inside a path, so all the edges
% will in fact be a the same path.
% 2) The node construct is used to draw the nodes. Nodes are special
% in the way that they are drawn *after* the path is drawn. This
% is very useful in this case because the nodes will be drawn on
% top of the path and therefore hide all edge joins.
% 3) Simple arithmetics can be used when specifying coordinates.
\draw \foreach \x in {0,36,...,324}
{
(\x:2) node {} -- (\x+108:2)
(\x-10:3) node {} -- (\x+5:4)
(\x-10:3) -- (\x+36:2)
(\x-10:3) --(\x+170:3)
(\x+5:4) node {} -- (\x+41:4)
};
\end{tikzpicture}
*/
\ No newline at end of file
......@@ -29,7 +29,6 @@ import java.util.Set;
import java.util.function.BiFunction;
import java.util.function.Function;
import java.util.stream.Collectors;
import java.util.stream.DoubleStream;
/**
* @author Benedikt Zoennchen
......@@ -37,6 +36,7 @@ import java.util.stream.DoubleStream;
public class PSMeshing {
private BiFunction<Double, Double, Double> f = (l, l_0) -> Math.max(l_0 - l, 0.0);
private Function<Integer, Integer> funcNextTri = i -> (i*i)+1;
private static final Logger log = LogManager.getLogger(PSMeshing.class);
......@@ -49,26 +49,34 @@ public class PSMeshing {
private double scalingFactor;
private PriorityQueue<Pair<PFace<MeshPoint>, Double>> heap;
private double initialEdgeLen = 10.0;
private double initialEdgeLen = 0.4;
private double deps;
private double sumOfqDesiredEdgeLength;
private double sumOfqLengths;
private final int retrangulationRate = 10;
private final int retrangulationRate = 1;
private final int resplitRate = 1;
private int retrangulationCounter = 0;
private int nextTriangualtion = 1;
private double maxMovement = 0;
private int stepCounter = 0;
private boolean initialized = false;
public PSMeshing(
final IDistanceFunction distanceFunc,
final IEdgeLengthFunction edgeLengthFunc,
final double initialEdgeLen,
final VRectangle bound,
final Collection<? extends VShape> obstacleShapes) {
this.bound = bound;
this.heap = new PriorityQueue<>(new FaceComparator());
this.heap = new PriorityQueue<>(new FacePairComparator());
this.distanceFunc = distanceFunc;
this.edgeLengthFunc = edgeLengthFunc;
this.initialEdgeLen = initialEdgeLen;
this.obstacleShapes = obstacleShapes;
this.deps = 1.4901e-8 * initialEdgeLen;
this.triangulation = ITriangulation.createPTriangulation(IPointLocator.Type.DELAUNAY_HIERARCHY, bound, (x, y) -> new MeshPoint(x, y, false));
this.nextTriangualtion = 1;
}
/**
......@@ -76,20 +84,17 @@ public class PSMeshing {
*/
public void initialize() {
log.info("##### (start) compute a uniform refined triangulation #####");
UniformRefinementTriangulation uniformRefinementTriangulation = new UniformRefinementTriangulation(triangulation, bound, obstacleShapes, p -> 0.7, distanceFunc);
UniformRefinementTriangulation uniformRefinementTriangulation = new UniformRefinementTriangulation(triangulation, bound, obstacleShapes, p -> edgeLengthFunc.apply(p) * initialEdgeLen, distanceFunc);
uniformRefinementTriangulation.compute();
//computeScalingFactor();
for(PFace<MeshPoint> face : getMesh().getFaces()) {
heap.add(Pair.of(face, faceToQuality(face)));
}
initialized = true;
log.info("##### (end) compute a uniform refined triangulation #####");
}
public void improve() {
public ITriangulation<MeshPoint, PVertex<MeshPoint>, PHalfEdge<MeshPoint>, PFace<MeshPoint>> getTriangulation() {
return triangulation;
}
/*public void improve() {
computeScalingFactor();
int steps = heap.size();
......@@ -132,14 +137,41 @@ public class PSMeshing {
}
//retriangulate();
}
}*/
public void execute() {
if(!initialized) {
initialize();
}
double quality = getQuality();
while (quality < Parameters.qualityMeasurement) {
step();
quality = getQuality();
log.info("quality = " + quality);
}
public double step() {
computeScalingFactor();
computeForces();