Authentication method changed. UCLouvain users, please use the "UCLouvain SSO" button to connect on the website and use ssh + ssh key (https://git.immc.ucl.ac.be/-/profile/keys) or https + personnal access token (https://git.immc.ucl.ac.be/-/profile/personal_access_tokens) in your git clients.

Commit 1b57c63c authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

automatic curve sampling in uref fix #4

parent f0338120
Pipeline #7968 passed with stages
in 2 minutes and 11 seconds
......@@ -63,11 +63,9 @@ def mesh_size(x,projection) :
# %%
# The geometry.coarsen_boundaries function requires in input the coordinates of one (any)
# point which is inside the domain. The last parameters should be at least two times smaller
# than the smallest value returned by the mesh_size callback, otherwise the resulting
# geometry will be invalid.
# point which is inside the domain.
coarse = seamsh.geometry.coarsen_boundaries(domain,(6.21,42.95),lonlat,mesh_size,20)
coarse = seamsh.geometry.coarsen_boundaries(domain,(6.21,42.95),lonlat,mesh_size)
# %%
# The resulting coarsened domain can be meshed normally.
......
......@@ -25,7 +25,7 @@ from scipy.spatial import cKDTree
import logging
from typing import List
import itertools
from .gmsh import _curve_sample
from .gmsh import _curve_sample_fix
def _ensure_valid_points(x, pfrom, pto):
......@@ -55,7 +55,7 @@ class Distance:
points = []
for curve in itertools.chain(domain._curves, domain._interior_curves):
if (tags is None) or (curve.tag in tags):
points.append(_curve_sample(curve, sampling))
points.append(_curve_sample_fix(curve, sampling))
points = np.vstack(points)
self._tree = cKDTree(points)
self._projection = domain._projection
......
......@@ -330,8 +330,7 @@ from .gmsh import _curve_sample
def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
x0projection: osr.SpatialReference,
mesh_size: MeshSizeCallback,
sampling: float) -> Domain:
mesh_size: MeshSizeCallback) -> Domain:
""" Creates a new Domain with the same projection and coarsened
boundaries.
......@@ -341,18 +340,15 @@ def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
x0projection: the coordinates system of x0.
mesh_size: a function returning the desired mesh element size for given
coordinates
sampling: should be (at least) two times smaller than the smallest
target mesh element size, otherwise the coarsening algorithm will
fail. This parameter will be removed (and determined automatically)
in the future.
"""
x0 = _ensure_valid_points(np.array([x0]),x0projection,domain._projection)[0]
sampled = []
tags = []
maxtag = 1
str2tag = {}
mesh_size_half = lambda x,p : mesh_size(x,p)/2
for curve in domain._curves :
cs = _curve_sample(curve, sampling)
cs = _curve_sample(curve, mesh_size_half, domain._projection)
sampled.append(cs)
if curve.tag not in str2tag :
str2tag[curve.tag] = maxtag
......
......@@ -48,7 +48,7 @@ def _gmsh_curve_geo(curve_type: CurveType, pointsid):
raise ValueError("unknown curve_type '%s'" % curve_type)
def _curve_sample(curve, lc):
def _curve_sample_fix(curve, lc):
gmsh.model.add(str(uuid.uuid4()))
tags = list([gmsh.model.geo.addPoint(*x, 0) -
1 for x in curve.points[:-1, :]])
......@@ -70,6 +70,38 @@ def _curve_sample(curve, lc):
return nodes[:, :2]
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]
def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
version: float = 4.0, intermediate_file_name: str = None) -> None:
""" Calls gmsh to generate a mesh from a geometry and a mesh size callback
......
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