Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Jonathan Lambrechts
seamsh
Commits
f437cd5c
Commit
f437cd5c
authored
Jul 29, 2020
by
Jonathan Lambrechts
Browse files
pep8
parent
1b57c63c
Changes
8
Hide whitespace changes
Inline
Side-by-side
examples/1-simple.py
View file @
f437cd5c
...
...
@@ -7,10 +7,10 @@
import
seamsh
from
seamsh.geometry
import
CurveType
import
numpy
as
np
from
osgeo
import
osr
from
osgeo
import
osr
# %%
# First, let's create a domain object and its associated projection.
# First,
let's create a domain object and its associated projection.
# Any projection supported by osgeo.osr can be used.
domain_srs
=
osr
.
SpatialReference
()
...
...
@@ -19,68 +19,73 @@ domain = seamsh.geometry.Domain(domain_srs)
# %%
# Boundary curves are created from a list of control points, an osr projection,
# a physical tag and a curve type. If the projection does not match the
domain's
# projection, the points are re-projected.
# a physical tag and a curve type. If the projection does not match the
#
domain's
projection,
the points are re-projected.
#
# A STRICTPOLYLINE curve linearly interpolates control points and forces a mesh
# point on each control point.
domain
.
add_boundary_curve
([[
0
,
0
],[
-
0.2
,
0.25
],[
-
0.2
,
0.75
],[
0
,
1
]],
"wall0"
,
domain_srs
,
curve_type
=
CurveType
.
STRICTPOLYLINE
)
domain
.
add_boundary_curve
([[
0
,
0
],
[
-
0.2
,
0.25
],
[
-
0.2
,
0.75
],
[
0
,
1
]],
"wall0"
,
domain_srs
,
curve_type
=
CurveType
.
STRICTPOLYLINE
)
# %%
# A POLYLINE is similar to a STRICTPOLYLINE but no mesh point is forced on the
# control points, the mesher can cut the corners.
# control points,
the mesher can cut the corners.
domain
.
add_boundary_curve
(
[[
0
,
1
],[
0.25
,
1.2
],[
0.75
,
1.2
],[
1
,
1
]]
,
"wall1"
,
domain_srs
,
curve_type
=
CurveType
.
POLYLINE
)
x
=
[[
0
,
1
],
[
0.25
,
1.2
],
[
0.75
,
1.2
],
[
1
,
1
]]
domain
.
add_boundary_curve
(
x
,
"wall1"
,
domain_srs
,
CurveType
.
POLYLINE
)
# %%
# A SPLINE uses a cubic spline interpolation through the control points
domain
.
add_boundary_curve
(
[[
0
,
0
],[
0.25
,
-
0.2
],[
0.75
,
-
0.2
],[
1
,
0
]]
,
"wall2"
,
domain_srs
,
curve_type
=
CurveType
.
SPLINE
)
x
=
[[
0
,
0
],
[
0.25
,
-
0.2
],
[
0.75
,
-
0.2
],
[
1
,
0
]]
domain
.
add_boundary_curve
(
x
,
"wall2"
,
domain_srs
,
CurveType
.
SPLINE
)
# %%
# A BSPLINE is smoother than a SPLINE but the control points are not on the
# curve.
domain
.
add_boundary_curve
(
[[
1
,
1
],[
1.2
,
0.75
],[
1.2
,
0.25
],[
1
,
0
]]
,
"wall3"
,
domain_srs
,
curve_type
=
CurveType
.
BSPLINE
)
x
=
[[
1
,
1
],
[
1.2
,
0.75
],
[
1.2
,
0.25
],
[
1
,
0
]]
domain
.
add_boundary_curve
(
x
,
"wall3"
,
domain_srs
,
CurveType
.
BSPLINE
)
# %%
# If the last point is the same as the first one, a periodic curve is
# If the last point is the same as the first one,
a periodic curve is
# created.
domain
.
add_boundary_curve
(
[[
0.4
,
0.6
],[
0.6
,
0.6
],[
0.6
,
0.4
],[
0.4
,
0.4
],[
0.4
,
0.6
]]
,
"island"
,
domain_srs
,
curve_type
=
CurveType
.
BSPLINE
)
x
=
[[
0.4
,
0.6
],
[
0.6
,
0.6
],
[
0.6
,
0.4
],
[
0.4
,
0.4
],
[
0.4
,
0.6
]]
domain
.
add_boundary_curve
(
x
,
"island"
,
domain_srs
,
CurveType
.
BSPLINE
)
# %%
# Interior curves are meshed on both sides, they do not define the domain
# Interior curves are meshed on both sides,
they do not define the domain
# boundaries but forces edge alignment and mesh points. When an interior curve
# touch a boundary or another interior curve, a control point should be
# touch a boundary or another interior curve,
a control point should be
# present in both curves at the intersection.
domain
.
add_interior_curve
(
[[
-
0.2
,
0.75
],[
0.1
,
0.6
],[
0.2
,
0.4
],[
0.1
,
0.4
]]
,
"interior1"
,
domain_srs
,
curve_type
=
CurveType
.
BSPLINE
)
domain
.
add_interior_curve
(
[[
0.1
,
0.4
],[
0.1
,
0.6
],[
0.1
,
0.8
]]
,
"interior2"
,
domain_srs
,
curve_type
=
CurveType
.
POLYLINE
)
x
=
[[
-
0.2
,
0.75
],
[
0.1
,
0.6
],
[
0.2
,
0.4
],
[
0.1
,
0.4
]]
domain
.
add_interior_curve
(
x
,
"interior1"
,
domain_srs
,
CurveType
.
BSPLINE
)
x
=
[[
0.1
,
0.4
],
[
0.1
,
0.6
],
[
0.1
,
0.8
]]
domain
.
add_interior_curve
(
x
,
"interior2"
,
domain_srs
,
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
)
x
=
[[
0.8
,
0.13
]]
domain
.
add_interior_points
(
x
,
"forcedpoint1"
,
domain_srs
)
x
=
[[
0.8
,
0.8
],
[
0.8
,
0.7
]]
domain
.
add_interior_points
(
x
,
"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
# defines the mesh element size. In this case,
a simple analytical function is
# used.
def
mesh_size
(
x
,
projection
)
:
return
(
0.005
+
(
1
+
np
.
sin
(
4
*
np
.
pi
*
x
[:,
1
]
*
np
.
pi
*
x
[:,
0
]))
*
0.01
)
def
mesh_size
(
x
,
projection
):
return
(
0.005
+
(
1
+
np
.
sin
(
4
*
np
.
pi
*
x
[:,
1
]
*
np
.
pi
*
x
[:,
0
]))
*
0.01
)
# %%
# Finally, the seamesh.gmsh.mesh function generates the mesh.
# Finally,
the seamesh.gmsh.mesh function generates the mesh.
seamsh
.
gmsh
.
mesh
(
domain
,
"basic_geometry.msh"
,
mesh_size
,
version
=
2.0
)
seamsh
.
gmsh
.
mesh
(
domain
,
"basic_geometry.msh"
,
mesh_size
,
version
=
2.0
)
examples/2-gis.py
View file @
f437cd5c
# %%
# GIS Data
# ========
# This example illustrates how the lirbary interacts with various GIS file
formats.
# directly from numpy arrays and osr spatial references.
#
#
For the sake of this example, we start by
download
ing
some data.
# This example illustrates how the lirbary interacts with various GIS file
#
formats
directly from numpy arrays and osr spatial references.
#
#
Let's
download some data.
import
os
import
urllib.request
import
tarfile
if
not
os
.
path
.
isdir
(
"data2"
)
:
url
lib
.
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:*"
)
if
not
os
.
path
.
isdir
(
"data2"
):
url
=
"ftp://braque.mema.ucl.ac.be/seamsh/data-test-2.tar.gz"
urllib
.
request
.
urlretrieve
(
url
,
"data-test-2.tar.gz"
)
f
=
tarfile
.
open
(
"data-test-2.tar.gz"
,
"r:*"
)
f
.
extractall
()
# %%
...
...
@@ -22,7 +22,7 @@ import seamsh
from
seamsh.geometry
import
CurveType
import
seamsh.geometry
import
numpy
as
np
from
osgeo
import
osr
from
osgeo
import
osr
# %%
# First, let's create a domain object and its associated projection.
...
...
@@ -37,47 +37,52 @@ 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
(
"data2/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.
# Physical tags are optional
domain
.
add_interior_curves_shp
(
"data2/interior.shp"
,
None
,
CurveType
.
STRICTPOLYLINE
)
domain
.
add_interior_points_shp
(
"data2/interior_points.shp"
,
"physical"
)
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
).
# field.Raster allows to load geotiff files (or any
gdal
raster file).
bath_field
=
seamsh
.
field
.
Raster
(
"data2/medit.tiff"
)
# %%
# field.Distance can be used to compute the distance from given (tagged) boundaries.
# A first field returns the distance from boundaries with physical tag "coast", or
# "island". The curves are sampled at regular intervals and the computed distance is
# actually the distance from the closest sampling point. The second argument sets
# the length of the interval between the sampling points.
# field.Distance can be used to compute the distance from given (tagged)
# boundaries. A first field returns the distance from boundaries with physical
# tag "coast", or "island". The curves are sampled at regular intervals and
# the computed distance is actually the distance from the closest sampling
# point. The second argument sets the length of the interval between the
# sampling points.
dist_coast
_field
=
seamsh
.
field
.
Distance
(
domain
,
100
,[
"coast"
,
"island"
])
dist_coast
=
seamsh
.
field
.
Distance
(
domain
,
100
,
[
"coast"
,
"island"
])
# %%
# A second distance field returns the distance from the island tagged "porquerolles",
# in the shp file.
# A second distance field returns the distance from the island tagged
# "porquerolles", in the shp file.
dist_porquerolles
=
seamsh
.
field
.
Distance
(
domain
,
20
,
[
"porquerolles"
])
dist_porquerolles_field
=
seamsh
.
field
.
Distance
(
domain
,
20
,[
"porquerolles"
])
# %%
# The mesh size is defined based on those fields, smaller elements are
prescribed
# around the Porquerolles island.
# The mesh size is defined based on those fields, smaller elements are
#
prescribed
around the Porquerolles island.
def
mesh_size
(
x
,
projection
)
:
s_coast
=
np
.
clip
((
dist_coast
_field
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles
_field
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
def
mesh_size
(
x
,
projection
):
s_coast
=
np
.
clip
((
dist_coast
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
return
s_dist
# %%
# Another option would be to take a mesh size proportional to the square root
# of the (clipped) bathymetry :
...
...
@@ -93,10 +98,11 @@ def mesh_size(x,projection) :
# special value "-" an interactive gmsh graphical window will pop up
# after each meshing step.
seamsh
.
gmsh
.
mesh
(
domain
,
"gis_mesh.msh"
,
mesh_size
,
intermediate_file_name
=
"debug"
)
seamsh
.
gmsh
.
mesh
(
domain
,
"gis_mesh.msh"
,
mesh_size
,
intermediate_file_name
=
"debug"
)
# %%
# The gmsh.convert_to_gis function can be used to convert a gmsh .msh file
# into a shape file or into a geo package file.
seamsh
.
gmsh
.
convert_to_gis
(
"gis_mesh.msh"
,
domain_srs
,
"gis_mesh.gpkg"
)
seamsh
.
gmsh
.
convert_to_gis
(
"gis_mesh.msh"
,
domain_srs
,
"gis_mesh.gpkg"
)
examples/3-coarsening.py
View file @
f437cd5c
...
...
@@ -9,9 +9,9 @@
# - Thanks to the fitting of the input data resolution to the final mesh
# element size, the quality of the elements near the boundary is improved.
# - The input curves do not have to form topologically valid closed loops.
# The only requirement is that the external boundary should be closed
# but the lines defining the boundaries can intersect each others.
#
# The only requirement is that the external boundary should be closed
# but the lines defining the boundaries can intersect each others.
#
# But :
#
# - Interior curves are not handled, they can be present as long as they
...
...
@@ -23,17 +23,17 @@
import
os
import
urllib.request
import
tarfile
if
not
os
.
path
.
isdir
(
"data"
)
:
url
lib
.
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
(
"data"
):
url
=
"ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz"
urllib
.
request
.
urlretrieve
(
url
,
"data-test-1.tar.gz"
)
f
=
tarfile
.
open
(
"data-test-1.tar.gz"
,
"r:*"
)
f
.
extractall
()
import
seamsh
from
seamsh.geometry
import
CurveType
import
seamsh.geometry
import
numpy
as
np
from
osgeo
import
osr
from
osgeo
import
osr
lonlat
=
osr
.
SpatialReference
()
lonlat
.
ImportFromProj4
(
"+proj=latlong +ellps=WGS84 +unit=degrees"
)
...
...
@@ -43,32 +43,36 @@ utm31 = osr.SpatialReference()
utm31
.
ImportFromProj4
(
"+proj=utm +ellps=WGS84 +zone=31"
)
domain
=
seamsh
.
geometry
.
Domain
(
utm31
)
domain
.
add_boundary_curves_shp
(
"data/data_no_duplicate.shp"
,
"physical"
,
CurveType
.
POLYLINE
)
domain
.
add_boundary_curves_shp
(
"data/data_no_duplicate.shp"
,
"physical"
,
CurveType
.
POLYLINE
)
# %%
# Clip the domain by a circle in UTM (zone32) coordinates
alphas
=
np
.
linspace
(
0
,
2
*
np
.
pi
,
200
)
circle_points
=
np
.
array
([[
272e3
,
4763e3
]])
+
50e3
*
np
.
column_stack
([
np
.
sin
(
alphas
),
np
.
cos
(
alphas
)])
domain
.
add_boundary_curve
(
circle_points
,
"open"
,
utm32
)
alphas
=
np
.
linspace
(
0
,
2
*
np
.
pi
,
200
)
circle_points
=
np
.
array
([[
272e3
,
4763e3
]])
\
+
50e3
*
np
.
column_stack
([
np
.
sin
(
alphas
),
np
.
cos
(
alphas
)])
domain
.
add_boundary_curve
(
circle_points
,
"open"
,
utm32
)
dist_coast
=
seamsh
.
field
.
Distance
(
domain
,
100
,
[
"coast"
,
"island"
])
dist_porquerolles
=
seamsh
.
field
.
Distance
(
domain
,
20
,
[
"porquerolles"
])
dist_coast_field
=
seamsh
.
field
.
Distance
(
domain
,
100
,[
"coast"
,
"island"
])
dist_porquerolles_field
=
seamsh
.
field
.
Distance
(
domain
,
20
,[
"porquerolles"
])
def
mesh_size
(
x
,
projection
)
:
s_coast
=
np
.
clip
((
dist_coast
_field
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles
_field
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
def
mesh_size
(
x
,
projection
):
s_coast
=
np
.
clip
((
dist_coast
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
return
s_dist
# %%
# The geometry.coarsen_boundaries function requires in input the coordinates
of one (any)
# point which is inside the domain.
# The geometry.coarsen_boundaries function requires in input the coordinates
#
of one (any)
point which is inside the domain.
coarse
=
seamsh
.
geometry
.
coarsen_boundaries
(
domain
,(
6.21
,
42.95
),
lonlat
,
mesh_size
)
coarse
=
seamsh
.
geometry
.
coarsen_boundaries
(
domain
,
(
6.21
,
42.95
),
lonlat
,
mesh_size
)
# %%
# The resulting coarsened domain can be meshed normally.
seamsh
.
gmsh
.
mesh
(
coarse
,
"coarse_boundary_mesh.msh"
,
mesh_size
)
seamsh
.
gmsh
.
mesh
(
coarse
,
"coarse_boundary_mesh.msh"
,
mesh_size
)
examples/4-advanced-gmsh.py
View file @
f437cd5c
...
...
@@ -2,19 +2,19 @@
# Advanced Gmsh / Quad Mesh
# =========================
# Seamsh uses gmsh internally. The whole gmsh python api is exposed through the
# seamsh.gmsh.gmsh module (the repetition is not a typo). This example is a
# variation on the previous one using the gmsh api to generate a quadrangular
# seamsh.gmsh.gmsh module (the repetition is not a typo). This example is a
# variation on the previous one using the gmsh api to generate a quadrangular
# instead of triangular mesh.
#
#
# We start as usual.
import
os
import
urllib.request
import
tarfile
if
not
os
.
path
.
isdir
(
"data"
)
:
url
lib
.
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
(
"data"
):
url
=
"ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz"
urllib
.
request
.
urlretrieve
(
url
,
"data-test-1.tar.gz"
)
f
=
tarfile
.
open
(
"data-test-1.tar.gz"
,
"r:*"
)
f
.
extractall
()
...
...
@@ -22,34 +22,37 @@ import seamsh
from
seamsh.geometry
import
CurveType
import
seamsh.geometry
import
numpy
as
np
from
osgeo
import
osr
from
osgeo
import
osr
domain_srs
=
osr
.
SpatialReference
()
domain_srs
.
ImportFromProj4
(
"+proj=utm +ellps=WGS84 +zone=31"
)
domain
=
seamsh
.
geometry
.
Domain
(
domain_srs
)
domain
.
add_boundary_curves_shp
(
"data/data_no_duplicate.shp"
,
"physical"
,
CurveType
.
POLYLINE
)
domain
.
add_boundary_curves_shp
(
"data/data_no_duplicate.shp"
,
"physical"
,
CurveType
.
POLYLINE
)
dist_coast
_field
=
seamsh
.
field
.
Distance
(
domain
,
100
,[
"coast"
,
"island"
])
dist_porquerolles
_field
=
seamsh
.
field
.
Distance
(
domain
,
20
,[
"porquerolles"
])
dist_coast
=
seamsh
.
field
.
Distance
(
domain
,
100
,
[
"coast"
,
"island"
])
dist_porquerolles
=
seamsh
.
field
.
Distance
(
domain
,
20
,
[
"porquerolles"
])
def
mesh_size
(
x
,
projection
)
:
s_coast
=
np
.
clip
((
dist_coast_field
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles_field
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
def
mesh_size
(
x
,
projection
):
s_coast
=
np
.
clip
((
dist_coast
(
x
,
projection
)
-
400
)
*
0.5
,
200
,
5000
)
s_porq
=
np
.
clip
((
dist_porquerolles
(
x
,
projection
)
-
200
)
*
0.5
,
50
,
5000
)
s_dist
=
np
.
minimum
(
s_coast
,
s_porq
)
return
s_dist
coarse
=
seamsh
.
geometry
.
coarsen_boundaries
(
domain
,(
8e5
,
4.68e6
),
domain_srs
,
mesh_size
,
20
)
coarse
=
seamsh
.
geometry
.
coarsen_boundaries
(
domain
,
(
8e5
,
4.68e6
),
domain_srs
,
mesh_size
)
# %%
# The gmsh options are set. The full list of gmsh mesh options can be found in
# `gmsh's documentation <https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options-list>`_.
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.RecombineAll"
,
1
)
;
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.RecombinationAlgorithm"
,
1
)
;
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.Algorithm"
,
8
)
;
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.RecombineAll"
,
1
)
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.RecombinationAlgorithm"
,
1
)
seamsh
.
gmsh
.
gmsh
.
option
.
setNumber
(
"Mesh.Algorithm"
,
8
)
# %%
# Eventually the mesh is generated.
seamsh
.
gmsh
.
mesh
(
coarse
,
"quad_mesh.msh"
,
mesh_size
)
seamsh
.
gmsh
.
mesh
(
coarse
,
"quad_mesh.msh"
,
mesh_size
)
seamsh/__init__.py
View file @
f437cd5c
# seamsh - Copyright (C) <2010-2020>
# <Universite catholique de Louvain (UCL), Belgium
#
#
# List of the contributors to the development of seamsh: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
Lesser General
# Public License as published by the Free Software Foundation,
either version
# 3 of the License, or (at your option) any later version.
#
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
#
Lesser General
Public License as published by the Free Software Foundation,
#
either version
3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING file). If not,
# along with this program (see COPYING file). If not,
# see <http://www.gnu.org/licenses/>.
from
.
import
field
...
...
seamsh/field.py
View file @
f437cd5c
# seamsh - Copyright (C) <2010-2020>
# <Universite catholique de Louvain (UCL), Belgium
#
#
# List of the contributors to the development of seamsh: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
Lesser General
# Public License as published by the Free Software Foundation,
either version
# 3 of the License, or (at your option) any later version.
#
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
#
Lesser General
Public License as published by the Free Software Foundation,
#
either version
3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING file). If not,
# along with this program (see COPYING file). If not,
# see <http://www.gnu.org/licenses/>.
from
.geometry
import
Domain
...
...
seamsh/geometry.py
View file @
f437cd5c
# seamsh - Copyright (C) <2010-2020>
# <Universite catholique de Louvain (UCL), Belgium
#
#
# List of the contributors to the development of seamsh: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
Lesser General
# Public License as published by the Free Software Foundation,
either version
# 3 of the License, or (at your option) any later version.
#
#
# This program (seamsh) is free software:
# you can redistribute it and/or modify it under the terms of the GNU
#
Lesser General
Public License as published by the Free Software Foundation,
#
either version
3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING file). If not,
# along with this program (see COPYING file). If not,
# see <http://www.gnu.org/licenses/>.
from
osgeo
import
ogr
,
osr
...
...
@@ -29,7 +29,6 @@ import ctypes as c
from
enum
import
Enum
import
typing
libdir
=
os
.
path
.
dirname
(
os
.
path
.
realpath
(
__file__
))
if
platform
.
system
()
==
"Windows"
:
libpath
=
os
.
path
.
join
(
libdir
,
"seamsh.dll"
)
...
...
@@ -100,7 +99,7 @@ class _Curve:
class
_Point
:
def
__init__
(
self
,
x
,
tag
)
:
def
__init__
(
self
,
x
,
tag
):
self
.
tag
=
tag
self
.
x
=
x
...
...
@@ -134,9 +133,10 @@ class Domain:
curve
.
pointsid
=
unique_id
[
cid
:
cid
+
n
]
cid
+=
n
# assign points id for interior points
for
point
in
self
.
_interior_points
:
for
point
in
self
.
_interior_points
:
point
.
pointid
=
unique_id
[
cid
]
cid
+=
1
# break curves on repeated points
def
split_curves
(
curves
,
nbk
):
newcurves
=
[]
...
...
@@ -237,7 +237,7 @@ class Domain:
self
.
_add_geometry
(
i
.
geometry
(),
phys
,
layerproj
,
curve_type
,
interior
,
points
)
def
add_interior_points
(
self
,
points
:
np
.
ndarray
,
physical_tag
:
str
,
def
add_interior_points
(
self
,
points
:
np
.
ndarray
,
physical_tag
:
str
,
projection
:
osr
.
SpatialReference
)
->
None
:
""" Add forced interior mesh points
...
...
@@ -247,7 +247,7 @@ class Domain:
projection: the points coordinate system
"""
x
=
_ensure_valid_points
(
points
,
projection
,
self
.
_projection
)
points
=
list
(
_Point
(
p
,
physical_tag
)
for
p
in
x
)
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
,
...
...
@@ -324,10 +324,8 @@ class Domain:
self
.
_add_shapefile
(
filename
,
physical_name_field
,
False
,
False
,
curve_type
)
from
.gmsh
import
_curve_sample
def
coarsen_boundaries
(
domain
:
Domain
,
x0
:
typing
.
Tuple
[
float
,
float
],
x0projection
:
osr
.
SpatialReference
,
mesh_size
:
MeshSizeCallback
)
->
Domain
:
...
...
@@ -341,34 +339,37 @@ def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
mesh_size: a function returning the desired mesh element size for given
coordinates
"""
x0
=
_ensure_valid_points
(
np
.
array
([
x0
]),
x0projection
,
domain
.
_projection
)[
0
]
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
:
def
mesh_size_half
(
x
,
p
):
return
mesh_size
(
x
,
p
)
/
2
for
curve
in
domain
.
_curves
:
cs
=
_curve_sample
(
curve
,
mesh_size_half
,
domain
.
_projection
)
sampled
.
append
(
cs
)
if
curve
.
tag
not
in
str2tag
:
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
))
tags
.
append
(
np
.
full
(
cs
.
shape
[
0
],
str2tag
[
curve
.
tag
],
dtype
=
np
.
int32
))
x
=
np
.
vstack
(
sampled
)
tags
=
np
.
concatenate
(
tags
)
x
,
unique_id
,
_
=
_generate_unique_points
(
x
)
utags
=
list
([
-
1
,
-
1
]
for
i
in
x
)
for
i
,
tag
in
zip
(
unique_id
,
tags
)
:
if
utags
[
i
][
0
]
==
-
1
:
utags
=
list
([
-
1
,
-
1
]
for
i
in
x
)
for
i
,
tag
in
zip
(
unique_id
,
tags
):
if
utags
[
i
][
0
]
==
-
1
:
utags
[
i
][
0
]
=
tag
found
=
False
j
=
i
while
j
!=
-
1
:
while
j
!=
-
1
:
found
|=
utags
[
j
][
0
]
==
tag
j
=
utags
[
j
][
1
]
if
not
found
:
if
not
found
:
utags
[
j
][
1
]
=
len
(
utags
)
utags
.
append
([
tag
,
-
1
])
utags
.
append
([
tag
,
-
1
])
tags
=
np
.
array
(
utags
).
reshape
([
-
1
])
x
=
np
.
copy
(
x
[:,
:
2
])
# avoid cocircular points
...
...
@@ -381,7 +382,7 @@ def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
r
=
c
.
c_void_p
(
tmp
.
ctypes
.
data
)
r
.
tmp
=
tmp
return
r
tri
=
Delaunay
(
x
)
tri
=
Delaunay
(
x
)
first
=
tri
.
find_simplex
(
x0
)
n_l
=
c
.
c_int
()
n_xo
=
c
.
c_int
()
...
...
@@ -401,8 +402,8 @@ def coarsen_boundaries(domain: Domain, x0: typing.Tuple[float, float],
lines
=
np
.
ctypeslib
.
frombuffer
(
linesbuf
.
contents
,
dtype
=
np
.
int32
)