Merge pull request #99 from Morpho-lang/meshgen

Meshgen module

.gitignore

@@ -59,3 +59,4 @@ dkms.conf
 
 test/FailedTests.txt
 *.png
+manual/src/manual.lyx~

benchmark/MeshGen/square.morpho

@@ -0,0 +1,10 @@
+// Square domain with a hole 
+import meshgen
+import plot
+
+var e0 = Domain(fn (x) ((x[0])^2+x[1]^2-1))
+var mg = MeshGen(e0, [-2..2:0.2, -2..2:0.2], quiet=true)
+var m = mg.build()
+
+print m
+

benchmark/benchmark.py

@@ -8,6 +8,7 @@
 from colored import stylize
 import subprocess
 
+"""
 languages = { "morpho": "morpho5",
               "m3": "morpho3",
               "fe": "evolver",
@@ -18,6 +19,9 @@
               "rb": "ruby",
               "dart": "dart"
             }
+"""
+
+languages = { "morpho": "morpho5" }
 
 # Gets the output generated
 def getoutput(filepath):

examples/meshgen/disk.morpho

@@ -0,0 +1,9 @@
+// Domain composed of a single disk
+import meshgen
+import plot
+
+var dom = fn (x) -(x[0]^2+x[1]^2-1)
+var mg = MeshGen(dom, [-1..1:0.2, -1..1:0.2], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/ellipse.morpho

@@ -0,0 +1,8 @@
+// Ellipse domain
+import meshgen
+import plot
+
+var e0 = Domain(fn (x) -((x[0]/2)^2+x[1]^2-1))
+var mg = MeshGen(e0, [-2..2:0.2, -1..1:0.2])
+var m = mg.build()
+Show(plotmesh(m, grade=1))

examples/meshgen/ellipsoidsection.morpho

@@ -0,0 +1,14 @@
+// 3D Ellipsoidal shell intersecting with a plane
+import meshgen
+import plot
+
+var dh = 0.2
+var e0 = Domain(fn (x) -((x[0]/2)^2+x[1]^2+x[2]^2-1))
+var e1 = Domain(fn (x) -((x[0]/2)^2+x[1]^2+x[2]^2-(0.5)^2))
+var e2 = HalfSpaceDomain(Matrix([0.25,0,0]), Matrix([1,0,0]))
+var dom = (e0.difference(e1)).intersection(e2)
+
+var mg = MeshGen(dom, [-2..2:dh, -1..1:dh, -1..1:dh], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/halfdisk.morpho

@@ -0,0 +1,11 @@
+// Domain composed of disk with half space removed
+import meshgen
+import plot
+
+var c = CircularDomain([0,0], 1)
+var hs = HalfSpaceDomain(Matrix([0,0]), Matrix([-1,0]))
+var dom = c.difference(hs) 
+var mg = MeshGen(dom, [-1..1:0.2, -1..1:0.2], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/overlappingdisks.morpho

@@ -0,0 +1,14 @@
+// Domain composed of overlapping disks
+import meshgen
+import plot
+
+// Create a complicated domain by composing circular disk domains
+var a = CircularDomain(Matrix([-0.5,0]), 1)
+var b = CircularDomain(Matrix([0.5,0]), 1)
+var c = CircularDomain(Matrix([0,0]), 0.3)
+var dom = a.union(b).difference(c)
+
+var mg = MeshGen(dom, [-2..2:0.1, -1..1:0.1], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/sphere.morpho

@@ -0,0 +1,10 @@
+// Sphere
+import meshgen
+import plot
+
+var dh = 0.2
+var dom = Domain(fn (x) -(x[0]^2+x[1]^2+x[2]^2-1))
+var mg = MeshGen(dom, [-1..1:dh, -1..1:dh, -1..1:dh], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/square.morpho

@@ -0,0 +1,11 @@
+// Square domain with a hole 
+import meshgen
+import plot
+
+var e0 = Domain(fn (x) ((x[0])^2+x[1]^2-1))
+var mg = MeshGen(e0, [-2..2:0.2, -2..2:0.2], quiet=false)
+//mg.maxiterations=10
+//mg.ttol=0.5
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/superellipse.morpho

@@ -0,0 +1,13 @@
+// Superellipse domain
+import meshgen
+import plot
+
+// Superellipse
+var e1 = Domain(fn (x) -((x[0]^4+x[1]^4)^(1/4)-1))
+var e2 = Domain(fn (x) -((x[0]^4+x[1]^4)^(1/4)-0.5))
+var dom = e1.difference(e2)
+
+var mg = MeshGen(dom, [-1..1:0.1, -1..1:0.1], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/superellipsoid.morpho

@@ -0,0 +1,11 @@
+// 3D Superellipsoid
+import meshgen
+import plot
+
+var dh = 0.2
+
+var dom = Domain(fn (x) -((x[0]^4+x[1]^4+x[2]^4)^(1/4)-1))
+var mg = MeshGen(dom, [-1..1:dh, -1..1:dh, -1..1:dh], quiet=false)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/meshgen/weighted.morpho

@@ -0,0 +1,10 @@
+// Circular domain with weight function
+import meshgen
+import plot
+
+// Weighted element size
+var hbar = fn (x) 1+x[0]/2 // larger elements as a function of x
+var mg = MeshGen(CircularDomain([0,0], 1), [-1..1:0.1, -1..1:0.1], quiet=false, weight=hbar)
+var m = mg.build()
+
+Show(plotmesh(m, grade=1))

examples/tactoid/tactoid.morpho

@@ -23,7 +23,7 @@ bnd.addgrade(0)
 var sigma=5.0*0.04 // Surface tension
 var W=3.0  // Anchoring
 
-var itermax = 8
+var itermax = 4
 
 // Set up the optimization problem
 var problem = OptimizationProblem(m)

morpho5/docs/index.rst

@@ -52,6 +52,7 @@ Morpho
    graphics
    implicitmesh
    kdtree
+   meshgen
    meshtools
    optimize
    plot

morpho5/docs/meshgen.md

@@ -0,0 +1,102 @@
+[comment]: # (Meshgen module help)
+[version]: # (0.5)
+
+# Meshgen
+[tagdelaunay]: # (delaunay)
+
+The `meshgen` module is used to create `Mesh` objects corresponding to a specified domain. It provides the `MeshGen` class to perform the meshing, which are created with the following arguments:
+
+    MeshGen(domain, boundingbox)
+
+Domains are specified by a scalar function that is positive in the region to be meshed and locally smooth. For example, to mesh the unit disk:
+
+    var dom = fn (x) -(x[0]^2+x[1]^2-1)
+
+A `MeshGen` object is then created and then used to build the `Mesh` like this:
+
+    var mg = MeshGen(dom, [-1..1:0.2, -1..1:0.2])
+    var m = mg.build()
+
+A bounding box for the mesh must be specified as a `List` of `Range` objects, one for each dimension. The increment on each `Range` gives an approximate scale for the size of elements generated.
+
+To facilitate convenient creation of domains, a `Domain` class is provided that provides set operations `union`, `intersection` and `difference`.
+
+`MeshGen` accepts a number of optional arguments:
+
+* `weight` A scalar weight function that controls mesh density.
+* `quiet` Set to `true` to suppress  `MeshGen` output.
+* `method` a list of options that controls the method used.
+
+Some method choices that are available include:
+
+* `"FixedStepSize"` Use a fixed step size in optimization.
+* `"StartGrid"` Start from a regular grid of points (the default).
+* `"StartRandom"` Start from a randomly generated collection of points.
+
+There are also a number of properties of a `MeshGen` object that can be set prior to calling `build` to control the operation of the mesh generation:
+
+* `stepsize`, `steplimit` Stepsize used internally by the `Optimizer`
+* `fscale` an internal "pressure"
+* `ttol` how far the vertices are allowed to move before retriangulation
+* `etol` energy tolerance for optimization problem
+
+`MeshGen` picks default values that cover a reasonable range of uses.
+
+[showsubtopics]: # (subtopics)
+
+## Domain
+[tagdomain]: # (domain)
+
+The `Domain` class is used to conveniently build a domain by composing simpler elements. 
+
+Create a `Domain` from a scalar function that is positive in the region of interest:
+
+    var dom = Domain(fn (x) -(x[0]^2+x[1]^2-1))
+
+You can pass it to `MeshGen` to specify the region to mesh: 
+
+    var mg = MeshGen(dom, [-1..1:0.2, -1..1:0.2])
+
+You can combine `Domain` objects using set operations `union`, `intersection` and `difference`: 
+
+    var a = CircularDomain(Matrix([-0.5,0]), 1)
+    var b = CircularDomain(Matrix([0.5,0]), 1)
+    var c = CircularDomain(Matrix([0,0]), 0.3)
+    var dom = a.union(b).difference(c)
+
+## CircularDomain
+[tagcirculardomain]: # (circulardomain)
+
+Conveniently constructs a `Domain` object correspondiong to a disk. Requires the position of the center and a radius as arguments. 
+
+Create a domain corresponding to the unit disk: 
+
+    var c = CircularDomain([0,0], 1)
+
+## HalfSpaceDomain
+[halfspacedomain]: # (halfspacedomain)
+
+Conveniently constructs a `Domain` object correspondiong to a half space defined by a plane at `x0` and a normal `n`:
+
+    var hs = HalfSpaceDomain(x0, n)
+
+Note `n` is an "outward" normal, so points into the *excluded* region.
+
+Half space corresponding to the allowed region `x<0`:
+
+    var hs = HalfSpaceDomain(Matrix([0,0,0]), Matrix([1,0,0]))
+
+Note that `HalfSpaceDomain`s cannot be meshed directly as they correspond to an infinite region. They are useful, however, for combining with other domains.
+
+Create half a disk by cutting a `HalfSpaceDomain` from a `CircularDomain`:
+
+    var c = CircularDomain([0,0], 1)
+    var hs = HalfSpaceDomain(Matrix([0,0]), Matrix([-1,0]))
+    var dom = c.difference(hs) 
+    var mg = MeshGen(dom, [-1..1:0.2, -1..1:0.2], quiet=false)
+    var m = mg.build()
+
+## MshGnDim
+[mshgndim]: # (mshgndim)
+
+The `MeshGen` module currently supports 2 and 3 dimensional meshes. Higher dimensional meshing will be available in a future release; please contact the developer if you are interested in this functionality.