Csgmsh

.github/workflows/test.yaml

@@ -89,6 +89,9 @@ jobs:
       - name: Install Graphviz
         if: ${{ matrix.os == 'ubuntu-latest' }}
         run: sudo apt install -y graphviz
+      - name: Install LibGLU
+        if: ${{ matrix.os == 'ubuntu-latest' }}
+        run: sudo apt install -y libglu1
       - name: Install Nutils and dependencies
         id: install
         env:
@@ -97,7 +100,7 @@ jobs:
           python -um pip install --upgrade --upgrade-strategy eager wheel
           python -um pip install --upgrade --upgrade-strategy eager numpy$_numpy_version
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install "$_wheel[import_gmsh,export_mpl]"
+          python -um pip install "$_wheel[import_gmsh,export_mpl,cs]"
       - name: Install Scipy
         if: ${{ matrix.matrix-backend == 'scipy' }}
         run: python -um pip install --upgrade scipy
@@ -171,12 +174,15 @@ jobs:
         with:
           name: python-package
           path: dist/
+      - name: Install LibGLU
+        if: ${{ matrix.os == 'ubuntu' }}
+        run: sudo apt install -y libglu1
       - name: Install Nutils and dependencies
         id: install
         run: |
           python -um pip install --upgrade --upgrade-strategy eager wheel
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install "$_wheel[matrix_scipy,export_mpl]"
+          python -um pip install "$_wheel[matrix_scipy,export_mpl,cs]"
       - name: Test
         run: python -um unittest discover -b -q -t . -s examples
   test-sphinx:

examples/adaptivity.py

@@ -1,11 +1,44 @@
-from nutils import mesh, function, export, testing
+from nutils import mesh, function, export, testing, cs
 from nutils.solver import System
 from nutils.expression_v2 import Namespace
+from dataclasses import dataclass
+from typing import Union
 import numpy
 import treelog
 
 
-def main(etype: str = 'square',
+@dataclass
+class structured:
+
+    def generate(self):
+        topo, geom = mesh.rectilinear([2,2])
+        quadrant = topo[1:,:1].withboundary(corner='left,top')
+        return topo - quadrant, (geom - 1) / 2
+
+
+@dataclass
+class unstructured:
+
+    elemsize: float = .15
+
+    def generate(self):
+        square = cs.Rectangle(cs.Interval(-.5, .5), cs.Interval(-.5, .5))
+        quadrant = cs.Rectangle(cs.Interval(0, .5), cs.Interval(-.5, 0))
+        shapes = dict(dom=square-quadrant, corner=quadrant.boundary)
+        return mesh.csgmsh(shapes, elemsize=self.elemsize)
+
+
+@dataclass
+class multipatch:
+
+    def generate(self):
+        vertices = [0,-.5], [-.5,-.5], [0,0], [-.5,.5], [.5,0], [.5,.5]
+        patches = [0,1,2,3], [2,3,4,5]
+        topo, geom = mesh.multipatch(patchverts=vertices, patches=patches, nelems=1)
+        return topo.withboundary(corner='patch0-bottom,patch1-bottom'), geom
+
+
+def main(basemesh: Union[structured, unstructured, multipatch] = structured(),
          btype: str = 'h-std',
          degree: int = 2,
          nrefine: int = 5):
@@ -29,24 +62,22 @@ def main(etype: str = 'square',
 
     Parameters
     ----------
-    etype
-        Type of elements (square/triangle/mixed).
+    basemesh
+        Initial mesh: structured, unstructured, or multipatch.
     btype
         Type of basis function (h/th-std/spline), with availability depending
         on the configured element type.
     degree
         Polynomial degree.
     nrefine
         Number of refinement steps to perform.
+    unstructured
+        Generate triangulated domain if True, or structured domain if False.
     '''
 
-    domain, geom = mesh.unitsquare(2, etype)
-    geom -= .5 # shift domain center to origin
-
+    domain, geom = basemesh.generate()
     x, y = geom
     exact = (x**2 + y**2)**(1/3) * numpy.cos(numpy.arctan2(y+x, y-x) * (2/3))
-    selection = domain.select(exact, ischeme='gauss1')
-    domain = domain.subset(selection, newboundary='corner')
     linreg = LinearRegressor(bias=1)
 
     for irefine in treelog.iter.fraction('level', range(nrefine+1)):
@@ -136,8 +167,8 @@ def offset(self):
 
 class test(testing.TestCase):
 
-    def test_square_quadratic(self):
-        error, u = main(nrefine=2)
+    def test_structured(self):
+        error, u = main(nrefine=2, basemesh=structured())
         with self.subTest('degrees of freedom'):
             self.assertEqual(len(u), 149)
         with self.subTest('L2-error'):
@@ -152,32 +183,27 @@ def test_square_quadratic(self):
                 7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
                 bya+ZCPbWKRPpvgFaedebw==''')
 
-    def test_triangle_quadratic(self):
-        error, u = main(nrefine=2, etype='triangle')
+    def test_unstructured(self):
+        error, u = main(nrefine=2, basemesh=unstructured(elemsize=.5))
         with self.subTest('degrees of freedom'):
-            self.assertEqual(len(u), 98)
+            self.assertEqual(len(u), 97)
         with self.subTest('L2-error'):
-            self.assertAlmostEqual(error[0], 0.00138, places=5)
+            self.assertAlmostEqual(error[0], 0.00095, places=5)
         with self.subTest('H1-error'):
-            self.assertAlmostEqual(error[1], 0.05326, places=5)
-        with self.subTest('left-hand side'):
-            self.assertAlmostEqual64(u, '''
-                eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
-                bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
-                YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
+            self.assertAlmostEqual(error[1], 0.04006, places=5)
 
-    def test_mixed_linear(self):
-        error, u = main(nrefine=2, etype='mixed', degree=1)
+    def test_multipatch(self):
+        error, u = main(nrefine=2, basemesh=multipatch())
         with self.subTest('degrees of freedom'):
-            self.assertEqual(len(u), 34)
+            self.assertEqual(len(u), 93)
         with self.subTest('L2-error'):
-            self.assertAlmostEqual(error[0], 0.00450, places=5)
+            self.assertAlmostEqual(error[0], 0.00128, places=5)
         with self.subTest('H1-error'):
-            self.assertAlmostEqual(error[1], 0.11692, places=5)
+            self.assertAlmostEqual(error[1], 0.05662, places=5)
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
-                AQDBThbY''')
+                eNpjYHhnZGfCZNpqysAgYJxuYmK6C8gqM35ksspU2QyZhZBF6DhnUGuUZPzJON5ktclbEyWgiInRVuOJ
+                QNlE06WmD02FgPpajaabuJpuMuUy8zZrMdsMFCk24TKda/rN1MWsDsi/j1UNpjnodgEAAk420A==''')
 
 
 if __name__ == '__main__':

examples/platewithhole.py

@@ -1,4 +1,4 @@
-from nutils import mesh, function, export, testing
+from nutils import mesh, function, export, testing, cs
 from nutils.solver import System
 from nutils.expression_v2 import Namespace
 from dataclasses import dataclass
@@ -43,8 +43,8 @@ def generate(self, radius):
 
 
 @dataclass
-class NURBS:
-    '''Non-Uniform Radional B-Splines
+class IGA:
+    '''Isogeometric Analysis
 
     Generate a 1x2 structured topology, map it using quadratic NURBS to a
     square domain with circular cut-out, and refine several times before
@@ -81,7 +81,33 @@ def generate(self, radius):
         return topo.withboundary(hole='left', sym='top,bottom', far='right'), geom, nurbsbasis, 5
 
 
-def main(mode: Union[FCM, NURBS] = NURBS(),
+@dataclass
+class TRI:
+    '''Triangulation
+
+    Generate a triangulated mesh with second order geometry using gmsh.
+
+    Parameters
+    ----------
+    elemsize
+        Target element size.
+    degree
+        Polynomial degree for the basis.
+    '''
+
+    elemsize: int = .1
+    degree: int = 2
+
+    def generate(self, radius):
+        rect = cs.Rectangle()
+        circ = cs.Circle(radius=radius)
+        shapes = dict(dom=rect-circ, sym=rect.boundary['left|bottom'], far=rect.boundary['right|top'], hole=circ.boundary)
+        topo, geom = mesh.csgmsh(shapes, elemsize=self.elemsize, order=2)
+        basis = topo.basis('std', self.degree)
+        return topo, geom, basis, self.degree
+
+
+def main(mode: Union[FCM, IGA, TRI] = IGA(),
          radius: float = .5,
          traction: float = .1,
          poisson: float = .3):
@@ -96,15 +122,16 @@ def main(mode: Union[FCM, NURBS] = NURBS(),
 
     The script can be run in two modes: by specifying `mode=FCM`, the circular
     hole is cut out of a regular finite element mesh by means of the Finite Cell
-    Method; by specifying `mode=NURBS` a Non-Uniform Rational BSpline geometry is
+    Method; by specifying `mode=IGA` a Non-Uniform Rational BSpline geometry is
     created to map a regular domain onto the desired shape. Either mode supports
     sub-parameters which can be specified from the command-line by attaching them
     in curly braces (e.g. `FCM{nelems=20,degree=1}`).
 
     Parameters
     ----------
     mode
-        Discretization strategy: FCM (Finite Cell Method) or NURBS.
+        Discretization strategy: FCM (Finite Cell Method), IGA (Isogeometric
+        Analysis) or TRI (triangulation).
     radius
         Cut-out radius.
     traction
@@ -137,10 +164,10 @@ def main(mode: Union[FCM, NURBS] = NURBS(),
     log.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr))
 
     sqr = topo.boundary['sym'].integral('(u_i n_i)^2 dS' @ ns, degree=degree*2)
-    cons = System(sqr, trial='u').solve_constraints(droptol=1e-15)
+    cons = System(sqr, trial='u').solve_constraints(droptol=1e-10)
 
     sqr = topo.boundary['far'].integral('du_k du_k dS' @ ns, degree=20)
-    cons = System(sqr, trial='u').solve_constraints(droptol=1e-15, constrain=cons)
+    cons = System(sqr, trial='u').solve_constraints(droptol=1e-10, constrain=cons)
 
     res = topo.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*2)
     args = System(res, trial='u', test='v').solve(constrain=cons)
@@ -193,7 +220,7 @@ def test_mixed(self):
                 fyhZkYI=''')
 
     def test_nurbs0(self):
-        err, cons, args = main(mode=NURBS(nrefine=0))
+        err, cons, args = main(mode=IGA(nrefine=0))
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00200, places=5)
         with self.subTest('h1-error'):
@@ -206,7 +233,7 @@ def test_nurbs0(self):
                 eNpjYJh07qLhhnOTjb0vTDdmAAKVcy/1u85lGYforQDzFc6pGSedlzd+eP4ykA8AvkQRaA==''')
 
     def test_nurbs2(self):
-        err, cons, args = main(mode=NURBS(nrefine=2))
+        err, cons, args = main(mode=IGA(nrefine=2))
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00009, places=5)
         with self.subTest('h1-error'):
@@ -222,10 +249,17 @@ def test_nurbs2(self):
                 n2c23nZe3djqQqpx88XNxrOv7gOr0zwXZeBxztro/bmnRp7nVY1zgTjvvIXxSaBfnl3YYbzmygmgOgDU
                 Imlr''')
 
+    def test_tri(self):
+        err, cons, args = main(mode=TRI(elemsize=.5))
+        with self.subTest('l2-error'):
+            self.assertAlmostEqual(err[0], .00053, places=5)
+        with self.subTest('h1-error'):
+            self.assertAlmostEqual(err[1], .0131, places=4)
+
 
 if __name__ == '__main__':
     from nutils import cli
     cli.run(main)
 
 
-# example:tags=elasticity,FCM,NURBS
+# example:tags=elasticity,FCM,IGA

nutils/cs/__init__.py

@@ -0,0 +1,67 @@
+'''Constructive Solid Geometry helper for Gmsh'''
+
+
+from ._shape import (
+    Box,
+    Circle,
+    Cylinder,
+    Ellipse,
+    Entities,
+    Interval,
+    Line,
+    Path,
+    Point,
+    Rectangle,
+    Skeleton,
+    Sphere,
+    generate_mesh,
+)
+
+
+from ._field import (
+    AsField,
+    Ball,
+    LocalRefinement,
+    Distance,
+    Threshold,
+    Max,
+    Min,
+    set_background,
+    x,
+    y,
+    z,
+)
+
+
+def write(fname: str, entities: Entities, elemsize: AsField, order: int = 1) -> None:
+    'Create .msh file based on Constructive Solid Geometry description.'
+
+    import treelog, gmsh # type: ignore
+
+    gmsh.initialize()
+    try:
+        gmsh.option.setNumber('General.Terminal', 0)
+        gmsh.option.setNumber('Mesh.Binary', 1)
+        gmsh.option.setNumber('Mesh.ElementOrder', order)
+        gmsh.option.setNumber('Mesh.CharacteristicLengthExtendFromBoundary', 0)
+        gmsh.option.setNumber('Mesh.CharacteristicLengthFromPoints', 0)
+        gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 0)
+
+        gmsh.logger.start()
+        try:
+            #set_background(gmsh.model.mesh.field, elemsize)
+            generate_mesh(gmsh.model, entities, elemsize)
+            gmsh.write(fname)
+        finally:
+            with treelog.context('gmsh'):
+                for line in gmsh.logger.get():
+                    level, sep, msg = line.partition(': ')
+                    if level in ('Debug', 'Info', 'Warning', 'Error'):
+                        getattr(treelog, level.lower())(msg)
+            gmsh.logger.stop()
+
+    finally:
+        gmsh.finalize()
+
+
+# vim:sw=4:sts=4:et