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

physical tag for interior points fix #20

parent 57cd72bc
Pipeline #7967 passed with stages
in 2 minutes and 9 seconds
......@@ -66,6 +66,12 @@ domain.add_interior_curve([[-0.2,0.75],[0.1,0.6],[0.2,0.4],[0.1,0.4]],
domain.add_interior_curve([[0.1,0.4],[0.1,0.6],[0.1,0.8]],
"interior2",domain_srs,curve_type=CurveType.POLYLINE)
# %%
# Force several tagged mesh points
domain.add_interior_points([[0.8,0.13]],"forcedpoint1",domain_srs)
domain.add_interior_points([[0.8,0.8],[0.8,0.7]],"forcedpoint2",domain_srs)
# %%
# A callback which takes an numpy array of points and a projection as argument
# defines the mesh element size. In this case, a simple analytical function is
......
......@@ -9,10 +9,10 @@
import os
import urllib.request
import tarfile
if not os.path.isdir("data") :
urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz",
"data-test-1.tar.gz")
f = tarfile.open("data-test-1.tar.gz","r:*")
if not os.path.isdir("data2") :
urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-2.tar.gz",
"data-test-2.tar.gz")
f = tarfile.open("data-test-2.tar.gz","r:*")
f.extractall()
# %%
......@@ -37,20 +37,21 @@ domain = seamsh.geometry.Domain(domain_srs)
# In the shapefile, a field named "physical" defines the physical tag of each
# curve. If a re-projection is required, it will be done automatically.
domain.add_boundary_curves_shp("data/data_no_duplicate.shp","physical",CurveType.POLYLINE)
domain.add_boundary_curves_shp("data2/data_no_duplicate.shp","physical",CurveType.POLYLINE)
# %%
# Interior curves (and interior points) can be loaded in a similar way.
# In this case we do not assign any physical tag.
# Interior curves and interior points can be loaded in a similar way.
# Physical tags are optional
domain.add_interior_curves_shp("data/interior.shp",None,CurveType.STRICTPOLYLINE)
domain.add_interior_curves_shp("data2/interior.shp",None,CurveType.STRICTPOLYLINE)
domain.add_interior_points_shp("data2/interior_points.shp","physical")
# %%
# Seamsh provides helper classes to compute the element size field.
#
# field.Raster allows to load geotiff files (or any raster file supported by gdal).
bath_field = seamsh.field.Raster("data/medit.tiff")
bath_field = seamsh.field.Raster("data2/medit.tiff")
# %%
# field.Distance can be used to compute the distance from given (tagged) boundaries.
......
......@@ -99,6 +99,12 @@ class _Curve:
self.curve_type = curve_type
class _Point:
def __init__(self,x,tag) :
self.tag = tag
self.x = x
class Domain:
""" List the domain boundaries, forced mesh points and
forced mesh lines. """
......@@ -117,7 +123,8 @@ class Domain:
def _build_topology(self):
curvesiter = itertools.chain(self._curves, self._interior_curves)
allpoints = np.row_stack(list(itertools.chain(
(curve.points for curve in curvesiter), self._interior_points)))
(curve.points for curve in curvesiter),
(point.x for point in self._interior_points))))
self._points, unique_id, nbreakpoints = _generate_unique_points(
allpoints)
# assign points id in curves
......@@ -127,9 +134,10 @@ class Domain:
curve.pointsid = unique_id[cid:cid+n]
cid += n
# assign points id for interior points
self._interior_points_id = unique_id[cid:]
for point in self._interior_points :
point.pointid = unique_id[cid]
cid += 1
# break curves on repeated points
def split_curves(curves, nbk):
newcurves = []
for curve in curves:
......@@ -198,7 +206,7 @@ class Domain:
curve_type, interior, onlypoints)
return
if onlypoints:
self.add_interior_points(geometry.GetPoints(), projection)
self.add_interior_points(geometry.GetPoints(), tag, projection)
else:
assert(geometry.GetGeometryType() == 2)
if interior:
......@@ -229,16 +237,18 @@ class Domain:
self._add_geometry(i.geometry(), phys, layerproj, curve_type,
interior, points)
def add_interior_points(self, points: np.ndarray,
def add_interior_points(self, points: np.ndarray, physical_tag : str,
projection: osr.SpatialReference) -> None:
""" Add forced interior mesh points
Args:
x: the points [n,2]
physical_tag: the points physical tag
projection: the points coordinate system
"""
points = _ensure_valid_points(points, projection, self._projection)
self._interior_points.append(points)
x = _ensure_valid_points(points, projection, self._projection)
points = list(_Point(p,physical_tag) for p in x)
self._interior_points += points
def add_boundary_curve(self, points: np.ndarray, physical_tag: str,
projection: osr.SpatialReference,
......@@ -247,7 +257,7 @@ class Domain:
Args:
points: the curve control points [n,2]
physical_tag: the curve physical_tag
physical_tag: the curve physical tag
projection: the points coordinate system
curve_type: curve interpolation
"""
......@@ -256,7 +266,7 @@ class Domain:
self._curves.append(curve)
def add_interior_curve(self, points: np.ndarray, physical_tag: str,
projection: str,
projection: osr.SpatialReference,
curve_type: CurveType = CurveType.POLYLINE) -> None:
""" Adds a tagged curve inside the domain. The curve
is not part of the domain boundary and will be meshed on both sides.
......
......@@ -118,10 +118,13 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
tags = _gmsh_curve_geo(curve.curve_type, curve.pointsid)
embeded.extend(tags)
physicals.setdefault(curve.tag, []).extend(tags)
physicalspt = {}
for point in domain._interior_points:
physicalspt.setdefault(point.tag,[]).append(point.pointid+1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(1, embeded, 2, stag)
gmsh.model.mesh.embed(
0, [i+1 for i in domain._interior_points_id], 2, stag)
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):
......@@ -135,6 +138,9 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
for name, tags in physicals.items():
tag = gmsh.model.addPhysicalGroup(1, tags)
gmsh.model.setPhysicalName(1, tag, name)
for name, tags in physicalspt.items():
tag = gmsh.model.addPhysicalGroup(0, tags)
gmsh.model.setPhysicalName(0, tag, name)
tag = gmsh.model.addPhysicalGroup(2, [stag])
gmsh.model.setPhysicalName(2, stag, "domain")
# 1D mesh
......
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