Commit ee67ab99 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Merge branch 'master' into distance_from_point

parents 655b47a6 c4e47416
Pipeline #8388 passed with stages
in 2 minutes and 36 seconds
......@@ -25,7 +25,9 @@ from scipy.spatial import cKDTree
import logging
from typing import List
import itertools
from .gmsh import _curve_sample_fix
from .gmsh import _curve_sample
import sys
import time
def _ensure_valid_points(x, pfrom, pto):
......@@ -35,6 +37,10 @@ def _ensure_valid_points(x, pfrom, pto):
x = np.asarray(x)
return x
def _lineup() :
if sys.stdout.isatty() :
print("\033[F\033[K",end="") # Cursor up one line and clear
class Distance:
"""Callable evaluating the distance to a set of discretized curves.
......@@ -53,12 +59,18 @@ class Distance:
if None, all curves are taken into account.
"""
points = []
for curve in itertools.chain(domain._curves, domain._interior_curves):
print("")
tic = time.time()
for icurve,curve in enumerate(itertools.chain(domain._curves, domain._interior_curves)):
if (tags is None) or (curve.tag in tags):
points.append(_curve_sample_fix(curve, sampling))
_lineup()
print("sampling curve {} for distance computation".format(icurve))
points.append(_curve_sample(curve, lambda x,proj : np.full([x.shape[0]],sampling),
None))
for point in itertools.chain(domain._interior_points):
if (tags is None) or (point.tag in tags):
points.append(point.x)
print("time : ",time.time()-tic)
points = np.vstack(points)
self._tree = cKDTree(points)
self._projection = domain._projection
......
......@@ -28,6 +28,7 @@ import itertools
import ctypes as c
from enum import Enum
import typing
import sys
libdir = os.path.dirname(os.path.realpath(__file__))
if platform.system() == "Windows":
......@@ -40,10 +41,15 @@ else:
lib = c.CDLL(libpath)
_transform_cache = {}
def _ensure_valid_points(x, pfrom, pto):
x = np.asarray(x)[:, :2]
if not pfrom.IsSame(pto):
x = osr.CoordinateTransformation(pfrom, pto).TransformPoints(x)
trans = _transform_cache.get((pfrom,pto),None)
if trans is None :
trans = osr.CoordinateTransformation(pfrom, pto)
_transform_cache[(pfrom,pto)] = trans
x = trans.TransformPoints(x)
x = np.asarray(x)[:, :2]
return x
......@@ -86,7 +92,9 @@ def _generate_unique_points(x):
cid += 1
xo = np.ndarray([cid, 2])
xo[unique_id, :] = x
return xo, unique_id, nbreakpoints
breakpoints = np.zeros([xo.shape[0]],dtype=bool)
breakpoints[:nbreakpoints] = True
return xo, unique_id, breakpoints
class _Curve:
......@@ -124,25 +132,29 @@ class Domain:
allpoints = np.row_stack(list(itertools.chain(
(curve.points for curve in curvesiter),
(point.x for point in self._interior_points))))
self._points, unique_id, nbreakpoints = _generate_unique_points(
self._points, unique_id, breakpoints = _generate_unique_points(
allpoints)
# assign points id in curves
cid = 0
touched_interior = np.zeros((self._points.shape[0],),bool)
for curve in itertools.chain(self._curves, self._interior_curves):
n = curve.points.shape[0]
curve.pointsid = unique_id[cid:cid+n]
touched_interior[curve.pointsid] = True
cid += n
# assign points id for interior points
for point in self._interior_points:
point.pointid = unique_id[cid]
touched_interior[point.pointid] = True
cid += 1
breakpoints = np.logical_and(touched_interior,breakpoints)
# break curves on repeated points
def split_curves(curves, nbk, loops=[]):
def split_curves(curves, breakpoints, loops=[]):
newcurves = []
newmap = []
for curve in curves:
breaks = np.where(curve.pointsid[1:-1] < nbk)[0]
breaks = np.where(breakpoints[curve.pointsid[1:-1]])[0]
if len(breaks) == 0:
newcurves.append(curve)
newmap.append([len(newcurves)-1])
......@@ -202,9 +214,8 @@ class Domain:
assert(i != -1)
except AssertionError:
raise ValueError("Invalid topology")
self._curves = split_curves(self._curves, nbreakpoints, self._curveloops)
self._interior_curves = split_curves(
self._interior_curves, nbreakpoints)
self._curves = split_curves(self._curves, breakpoints, self._curveloops)
self._interior_curves = split_curves(self._interior_curves, breakpoints)
loopbboxarea = np.zeros([len(self._curveloops)])
for i, loop in enumerate(self._curveloops):
lpts = np.row_stack([self._curves[j].points for j, o in loop])
......@@ -232,6 +243,7 @@ class Domain:
def _add_shapefile(self, filename, physical_name_field,
interior, points, curve_type):
print("importing features from '{}'\n".format(filename))
if filename[-5:] == ".gpkg" :
driver = ogr.GetDriverByName('GPKG')
else :
......@@ -240,6 +252,8 @@ class Domain:
layer = data.GetLayer()
layerdef = layer.GetLayerDefn()
physfield = None
count = 0
tic = time.time()
if physical_name_field is not None:
for i in range(layerdef.GetFieldCount()):
field_name = layerdef.GetFieldDefn(i).GetName()
......@@ -250,9 +264,13 @@ class Domain:
"' not found in shapefile")
layerproj = layer.GetSpatialRef()
for i in layer:
phys = i.GetField(physfield) if physfield else "boundary"
phys = i.GetField(physfield) if not (physfield is None) else "boundary"
self._add_geometry(i.geometry(), phys, layerproj, curve_type,
interior, points)
_lineup()
print("{} features imported".format(count))
count += 1
print("time : ",time.time()-tic)
def add_interior_points(self, points: np.ndarray, physical_tag: str,
projection: osr.SpatialReference) -> None:
......@@ -343,6 +361,12 @@ class Domain:
from .gmsh import _curve_sample
def _lineup() :
if sys.stdout.isatty() :
print("\033[F\033[K",end="") # Cursor up one line and clear
import time
def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
x0projection: osr.SpatialReference,
mesh_size: MeshSizeCallback) -> Domain:
......@@ -364,15 +388,21 @@ def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
str2tag = {}
def mesh_size_half(x, p):
return mesh_size(x, p)*0.25
for curve in domain._curves:
return mesh_size(x, p)*0.5
ncurves = len(domain._curves)
print("")
tic = time.time()
for icurve, curve in enumerate(domain._curves):
_lineup()
print("sampling curve {}/{}".format(icurve,ncurves))
cs = _curve_sample(curve, mesh_size_half, domain._projection)
sampled.append(cs)
if curve.tag not in str2tag:
str2tag[curve.tag] = maxtag
maxtag += 1
tags.append(np.full(cs.shape[0], str2tag[curve.tag], dtype=np.int32))
print("time : ", time.time()-tic)
x = np.vstack(sampled)
tags = np.concatenate(tags)
x, unique_id, _ = _generate_unique_points(x)
......
......@@ -30,8 +30,14 @@ import atexit
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("General.Verbosity",2)
GmshOptionValue = typing.Union[float, str, typing.Tuple[float, float, float]]
import sys
def _lineup() :
if sys.stdout.isatty() :
print("\033[F\033[K",end="") # Cursor up one line and clear
def _gmsh_curve_geo(curve_type: CurveType, pointsid):
pts = list(i+1 for i in pointsid)
......@@ -48,7 +54,29 @@ def _gmsh_curve_geo(curve_type: CurveType, pointsid):
raise ValueError("unknown curve_type '%s'" % curve_type)
def _curve_sample_fix(curve, lc):
def _curve_sample_gmsh_tag(tag, lc, projection):
bounds = gmsh.model.getParametrizationBounds(1,tag)
xi = np.linspace(bounds[0][0],bounds[1][0],5)
count = 0
x = gmsh.model.getValue(1,tag,xi).reshape([-1,3])
size = lc(x,projection)
while True :
count += 1
dist = np.linalg.norm(x[1:,:]-x[:-1,:],axis=1)
target = np.minimum(size[1:],size[:-1])
if not np.any(target<dist) :
break
exi = (xi[1:][target<dist]+xi[:-1][target<dist])/2
ex = gmsh.model.getValue(1,tag,exi).reshape([-1,3])
esize = lc(ex,projection)
xi = np.hstack([xi,exi])
s = np.argsort(xi)
xi = xi[s]
x = np.row_stack([x,ex])[s,:]
size = np.hstack([size,esize])[s]
return x[:,:2], xi, size
def _curve_sample_full(curve, lc, projection):
gmsh.model.add(str(uuid.uuid4()))
tags = list([gmsh.model.geo.addPoint(*x, 0) -
1 for x in curve.points[:-1, :]])
......@@ -58,48 +86,30 @@ def _curve_sample_fix(curve, lc):
tags.append(gmsh.model.geo.addPoint(*curve.points[-1, :], 0)-1)
ltag = _gmsh_curve_geo(curve.curve_type, tags)
gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(1)
nodes = gmsh.model.mesh.getNodes(1, ltag[0], includeBoundary=True)[
1].reshape([-1, 3])
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 1e22)
bounds = gmsh.model.getParametrizationBounds(1,ltag[0])
xi = np.linspace(bounds[0][0],bounds[1][0],5)
count = 0
x = gmsh.model.getValue(1,ltag[0],xi).reshape([-1,3])
size = lc(x,projection)
while True :
count += 1
dist = np.linalg.norm(x[1:,:]-x[:-1,:],axis=1)
target = np.minimum(size[1:],size[:-1])
if not np.any(target<dist) :
break
exi = (xi[1:][target<dist]+xi[:-1][target<dist])/2
ex = gmsh.model.getValue(1,ltag[0],exi).reshape([-1,3])
esize = lc(ex,projection)
xi = np.hstack([xi,exi])
s = np.argsort(xi)
xi = xi[s]
x = np.row_stack([x,ex])[s,:]
size = np.hstack([size,esize])[s]
gmsh.model.remove()
return nodes[:, :2]
return x[:,:2], xi, size
def _curve_sample(curve, lc, projection):
gmsh.option.setNumber("Mesh.CharacteristicLengthFromPoints", 0)
gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-3)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromParametricPoints",
1)
gmsh.model.add(str(uuid.uuid4()))
tags = list([gmsh.model.geo.addPoint(*x, 0) -
1 for x in curve.points[:-1, :]])
if np.linalg.norm(curve.points[0, :]-curve.points[-1, :]) < 1e-8:
tags.append(tags[0])
else:
tags.append(gmsh.model.geo.addPoint(*curve.points[-1, :], 0)-1)
ltag = _gmsh_curve_geo(curve.curve_type, tags)
gmsh.model.geo.synchronize()
nadapt1d = 3
for i in range(nadapt1d):
gmsh.model.mesh.generate(1)
for (dim, tag) in gmsh.model.getEntities(1):
_, x, u = gmsh.model.mesh.getNodes(dim, tag, True)
size = lc(x.reshape([-1, 3]), projection)
gmsh.model.mesh.setSizeAtParametricPoints(dim, tag, u, size)
gmsh.model.mesh.clear()
gmsh.model.mesh.generate(1)
nodes = gmsh.model.mesh.getNodes(1, ltag[0], includeBoundary=True)[
1].reshape([-1, 3])
gmsh.model.remove()
return nodes[:, :2]
return _curve_sample_full(curve, lc, projection)[0]
def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
......@@ -119,11 +129,10 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
domain._build_topology()
nadapt = 3
nadapt1d = 4
for curve in domain._curves:
curve.mesh_size = mesh_size(curve.points, domain._projection)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromPoints", 0)
gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-3)
gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-5)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 1)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
......@@ -156,16 +165,6 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
gmsh.model.mesh.embed(1, embeded, 2, stag)
gmsh.model.mesh.embed(
0, [p.pointid+1 for p in domain._interior_points], 2, stag)
it = 0
for cl in domain._curveloops:
for i, (l, o) in enumerate(cl):
curve = domain._curves[l]
sizes = curve.mesh_size if o else np.flip(curve.mesh_size)
u0, u1 = gmsh.model.getParametrizationBounds(1, abs(ctags[it]))
u = np.linspace(u0, u1, curve.points.shape[0])
gmsh.model.mesh.setSizeAtParametricPoints(
1, ctags[it], u, sizes)
it += 1
for name, tags in physicals.items():
tag = gmsh.model.addPhysicalGroup(1, tags)
gmsh.model.setPhysicalName(1, tag, name)
......@@ -175,16 +174,12 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
tag = gmsh.model.addPhysicalGroup(2, [stag])
gmsh.model.setPhysicalName(2, stag, "domain")
# 1D mesh
for i in range(nadapt1d):
gmsh.model.mesh.generate(1)
if intermediate_file_name is not None \
and intermediate_file_name != "-":
gmsh.write(intermediate_file_name+"_1d_"+str(i)+".msh")
for (dim, tag) in gmsh.model.getEntities(1):
_, x, u = gmsh.model.mesh.getNodes(dim, tag, True)
size = mesh_size(x.reshape([-1, 3]), domain._projection)
gmsh.model.mesh.setSizeAtParametricPoints(dim, tag, u, size)
gmsh.model.mesh.clear()
bounds = gmsh.model.getParametrizationBounds(1,2)
for icurve,(dim, tag) in enumerate(gmsh.model.getEntities(1)) :
_lineup()
print("sample curve for mesh_size",icurve,dim,tag);
_, xi, size = _curve_sample_gmsh_tag(tag, lambda x,p : mesh_size(x,p)/2, domain._projection)
gmsh.model.mesh.setSizeAtParametricPoints(dim, tag, xi, size*2)
gmsh.model.mesh.generate(1)
if intermediate_file_name is not None:
if intermediate_file_name == "-":
......
......@@ -497,7 +497,7 @@ static int *color_triangles(const double *x, int n_tri, const int *tri, const in
const double *p0 = x+2*i0;
const double *p1 = x+2*i1;
double d2 = (p1[0]-p0[0])*(p1[0]-p0[0])+(p1[1]-p0[1])*(p1[1]-p0[1]);
double lmin = fmin(length[i0],length[i1])*1.2;//*1.1;
double lmin = fmin(length[i0],length[i1]);
if(d2 > lmin*lmin) {
tri_color[n] = 1;
stack[++stack_p] = n;
......
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