From b449725a3299696919051ef567b67fa0ecf35275 Mon Sep 17 00:00:00 2001 From: softmattertheory Date: Fri, 25 Feb 2022 08:13:52 -0500 Subject: [PATCH 1/2] Meshgen module added Including documentation and examples --- .gitignore | 1 + examples/meshgen/disk.morpho | 9 + examples/meshgen/ellipse.morpho | 8 + examples/meshgen/ellipsoidsection.morpho | 14 + examples/meshgen/halfdisk.morpho | 11 + examples/meshgen/overlappingdisks.morpho | 14 + examples/meshgen/sphere.morpho | 10 + examples/meshgen/square.morpho | 11 + examples/meshgen/superellipse.morpho | 13 + examples/meshgen/superellipsoid.morpho | 11 + examples/meshgen/weighted.morpho | 10 + morpho5/docs/index.rst | 1 + morpho5/docs/meshgen.md | 96 ++++++ morpho5/modules/meshgen.morpho | 382 +++++++++++++++++++++++ morpho5/modules/plot.morpho | 27 +- 15 files changed, 613 insertions(+), 5 deletions(-) create mode 100644 examples/meshgen/disk.morpho create mode 100644 examples/meshgen/ellipse.morpho create mode 100644 examples/meshgen/ellipsoidsection.morpho create mode 100644 examples/meshgen/halfdisk.morpho create mode 100644 examples/meshgen/overlappingdisks.morpho create mode 100644 examples/meshgen/sphere.morpho create mode 100644 examples/meshgen/square.morpho create mode 100644 examples/meshgen/superellipse.morpho create mode 100644 examples/meshgen/superellipsoid.morpho create mode 100644 examples/meshgen/weighted.morpho create mode 100644 morpho5/docs/meshgen.md create mode 100644 morpho5/modules/meshgen.morpho diff --git a/.gitignore b/.gitignore index cf528893..9b40466a 100644 --- a/.gitignore +++ b/.gitignore @@ -59,3 +59,4 @@ dkms.conf test/FailedTests.txt *.png +manual/src/manual.lyx~ diff --git a/examples/meshgen/disk.morpho b/examples/meshgen/disk.morpho new file mode 100644 index 00000000..598b7bf5 --- /dev/null +++ b/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)) diff --git a/examples/meshgen/ellipse.morpho b/examples/meshgen/ellipse.morpho new file mode 100644 index 00000000..719eba87 --- /dev/null +++ b/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], quiet=false) +var m = mg.build() +Show(plotmesh(m, grade=1)) diff --git a/examples/meshgen/ellipsoidsection.morpho b/examples/meshgen/ellipsoidsection.morpho new file mode 100644 index 00000000..7339bd07 --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/examples/meshgen/halfdisk.morpho b/examples/meshgen/halfdisk.morpho new file mode 100644 index 00000000..3829b182 --- /dev/null +++ b/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)) diff --git a/examples/meshgen/overlappingdisks.morpho b/examples/meshgen/overlappingdisks.morpho new file mode 100644 index 00000000..2cf6419e --- /dev/null +++ b/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)) diff --git a/examples/meshgen/sphere.morpho b/examples/meshgen/sphere.morpho new file mode 100644 index 00000000..74faee38 --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/examples/meshgen/square.morpho b/examples/meshgen/square.morpho new file mode 100644 index 00000000..31312d5e --- /dev/null +++ b/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)) diff --git a/examples/meshgen/superellipse.morpho b/examples/meshgen/superellipse.morpho new file mode 100644 index 00000000..f53ff862 --- /dev/null +++ b/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)) diff --git a/examples/meshgen/superellipsoid.morpho b/examples/meshgen/superellipsoid.morpho new file mode 100644 index 00000000..185aacff --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/examples/meshgen/weighted.morpho b/examples/meshgen/weighted.morpho new file mode 100644 index 00000000..4a8af3c5 --- /dev/null +++ b/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)) diff --git a/morpho5/docs/index.rst b/morpho5/docs/index.rst index 33072688..a46cc35e 100644 --- a/morpho5/docs/index.rst +++ b/morpho5/docs/index.rst @@ -52,6 +52,7 @@ Morpho graphics implicitmesh kdtree + meshgen meshtools optimize povray diff --git a/morpho5/docs/meshgen.md b/morpho5/docs/meshgen.md new file mode 100644 index 00000000..28041f3b --- /dev/null +++ b/morpho5/docs/meshgen.md @@ -0,0 +1,96 @@ +[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 + +var a = CircularDomain([-0.5,0], 1) +var b = CircularDomain([0.5,0], 1) + +var dom = a.union(b).difference(c) + +You can combine `Domain` objects using set operations `union`, `intersection` and `difference`: + + var a = (s1.intersection(s2)).union(s3) + +## 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. \ No newline at end of file diff --git a/morpho5/modules/meshgen.morpho b/morpho5/modules/meshgen.morpho new file mode 100644 index 00000000..2e57992f --- /dev/null +++ b/morpho5/modules/meshgen.morpho @@ -0,0 +1,382 @@ +/* Meshgen - a simple morpho mesh generator + * + * Inspired by the distmesh algorithm presented in: + * A Simple Mesh Generator in MATLAB, Per-Olof Persson and Gilbert Strang, SIAM Rev., 46(2), 329–345 (2004) + * Reimplemented using Morpho's optimization capabilities */ + +import functionals +import optimize +import meshtools +import delaunay + +/* ******************************************************** + * The Domain of interest is specified by a generalization + * of a signed distance function: wherever the function is + * positive will be meshed. Domains are constructed with + * the function of interest and can be composed by set + * operations intersection, union, difference. + * ******************************************************** */ + +class Domain { + init (func) { + self.func = func + } + + intersection(dom) { // + fn fint (x) { + return min(self.func(x), dom.func(x)) + } + return Domain(fint) + } + + difference(dom) { + fn fdiff (x) { + return min(self.func(x), -dom.func(x)) + } + return Domain(fdiff) + } + + union(dom) { + fn funion (x) { + return max(self.func(x), dom.func(x)) + } + return Domain(funion) + } +} + +fn CircularDomain(x0, r) { + var xx = x0 + if (islist(x0)) xx = Matrix(x0) + + fn circle (x) { + var dx = x - xx + return (r^2-dx.inner(dx)) + } + + return Domain(circle) +} + +// A Half Space defined by a plane at x0 and a normal n. +// n is an "outward" normal, so points into the excluded region +fn HalfSpaceDomain(x0, n) { + var xx0 = x0 + if (islist(xx0)) xx0=Matrix(x0) + var nn = n + if (islist(nn)) xx0=Matrix(nn) + + return Domain(fn (x) -(x-x0).inner(n)) +} + +/* ******************************************************** + * Functionals to control mesh quality + * ******************************************************** */ + +/* ------------------------------------------------------- + * One sided hookean elasticity: + * Energy F = 1/2 ((L-L0)/L0)^2 for Llen0) return 0 + return ((len-len0)/len0)^2/2 + } + + gradientfn(mesh, vert, id, el, grd) { + var x0 = vert.column(el[0]) + var x1 = vert.column(el[1]) + var dx = x0-x1 + + var len = dx.norm() + var len0 = self.ref[self.grade, id] + if (len>len0) return + var g = dx*(len-len0)/(len0^2*len) + + grd.setcolumn(el[0], g+grd.column(el[0])) + grd.setcolumn(el[1], -g+grd.column(el[1])) + } +} + +/* ------------------------------------------------------- + * Evaluates a scalar potential at the midpoint of + * an element. Will be used to detect elements inside + * forbidden areas. + * ------------------------------------------------------- */ +class ScalarPotentialMidpoint is Functional { + init(hbar) { + self.func = hbar + super.init(1) + } + + integrandfn(mesh, vert, id, el) { + var dim = el.count() + var x = vert.column(el[0]) + for (var i=1; iself.ttol) return true + } + } + + if (dx.norm()>self.ttol) return true // old criterion + return super.hasconverged() + } + + didconverge() { // Checks if we converged according to superclass + return super.hasconverged() + } +} + +/* ******************************************************** + * Mesh generator + * ******************************************************** */ + +var MshGnDim = Error("MshGnDim", "MeshGen only supports 2 or 3 dimensions") + +class MeshGen { + init (f, bbox, weight=nil, quiet=false, method=nil) { + self.func = f // The function to call + if (!iscallable(f) && isobject(f) && f.has("func")) self.func = f.func // You can give a Domain object directly + + self.bbox = bbox // Bounding box + self.weight = weight // Weight function + + fn defaultweight (x) { return 1 } // Provide a default weight function + if (isnil(weight)) self.weight = defaultweight + + self.dim = bbox.count() // Dimensionality of problem + + self.quiet = quiet // Whether to report output + + var sep = [] + for (range in bbox) sep.append(range[1]-range[0]) // Get stepsize + self.h0 = min(sep) // Element separation + + self.stepsize = self.h0/5 // Stepsize for optimizer + self.steplimit = self.h0/2 // Steplm for optimizer + if (self.dim>2) { self.stepsize/=2; self.steplimit/=2 } + + self.mesh = nil // Mesh + self.problem = nil // Optimization problem + + self.fscale = 1.2 // Internal pressure + if (self.dim>2) self.fscale = 1.1 + + self.etol = 1e-3 // Energy tolerance for optimization problem (it's deliberately loose) + self.ttol = 0.5 // How much motion of the vertices before retriangulation + self.xtol = 1e-12 // Tolerance for point comparisons + + self.method = [] // Method selection + if (isstring(method)) self.method.append(method) + if (islist(method)) for (m in method) self.method.append(m) + + self.maxiterations = 100 + } + + initialgridmesh() { // Create the initial mesh + if (!self.quiet) print "Creating initial mesh on regular grid" + var f = self.func + var vert = [] + if (self.dim==2) { + for (v in self.bbox[1]) { + for (u in self.bbox[0]) { + if (f(Matrix([u,v]))>0) vert.append([u,v]) + } + } + } else if (self.dim==3) { + for (w in self.bbox[2]) { + for (v in self.bbox[1]) { + for (u in self.bbox[0]) { + if (f(Matrix([u,v,w]))>0) vert.append([u,v,w]) + } + } + } + } else { + MshGnDim.throw() + } + self.retriangulate(vert) + } + + initialrandommesh() { // Create an initial mesh from random points + var vert = [] + var bnds = [] + var vol = 1 // compute the generalized volume of the element + for (b in self.bbox) { + var sep = b[b.enumerate(-1)-1]-b[0] // range of this dimension + bnds.append([b[0], sep]) + vol *= sep + } + + var n = ceil(vol/(self.h0)^(self.dim)) + if (!self.quiet) print "Creating initial mesh with ${n} random points" + + for (var nsuccess=0; nsuccess0) { + vert.append(pt) + nsuccess+=1 + } + } + + self.retriangulate(vert) + } + + selectbboxpts() { + var bnds = [], xtol = self.xtol + for (range,k in self.bbox) { + bnds.append(bounds(range)) + } + + fn selectOnBbox2D(x,y) { + return ( abs(x-bnds[0][0])0) { + mb.addface(t) + } + } + + self.mesh=mb.build() + self.mesh.addgrade(1) + } + + vertices() { // Extract a list of vertices from the current mesh + var v = self.mesh.vertexmatrix() + var vlist = [] + for (id in 0...v.dimensions()[1]) vlist.append(v.column(id)) + return vlist + } + + setupproblem() { // Setup optimization problem + self.problem = OptimizationProblem(self.mesh) + // Build reference lengths + var lengths = Length().integrand(self.mesh) + + // Evaluate the weight function at the midpoint + var hbar = ScalarPotentialMidpoint(self.weight).integrand(self.mesh) + + // Calculate preferred lengths + var len0 = Field(self.mesh, grade=1) + + fn lnorm(l, n) { // n-norm of a matrix + var sum = 0 + for (a in l) sum+=a^n + return sum^(1/n) + } + + //var lscale=self.fscale*(lengths.inner(lengths)/hbar.inner(hbar))^(1/2) + var lscale=self.fscale*(lnorm(lengths, self.dim)/lnorm(hbar, self.dim)) + for (id in 0...self.mesh.count(1)) len0[id]=lscale*hbar[id] + + // Impose elasticity + var lel = OneSidedHookeElasticity(len0) + self.problem.addenergy(lel) + + // Use a one sided level set constraint to keep things in the domain + var ls + if (self.dim==2) { // [TODO]: Rewrite this with varargs when available + ls = ScalarPotential(fn (x,y) self.func(Matrix([x,y]))) + } else if (self.dim==3) { + ls = ScalarPotential(fn (x,y,z) self.func(Matrix([x,y,z]))) + } + self.problem.addlocalconstraint(ls, onesided=true) + } + + optimize() { // Perform optimization + var opt = MGShapeOptimizer(self.problem, self.mesh, ttol=self.ttol*self.h0, localcheck=self.method.ismember("LocalCheck")) + opt.fix(self.selectbboxpts()) + opt.stepsize = self.stepsize + opt.steplimit = self.steplimit + opt.etol = self.etol // Use a fairly loose convergence criterion + opt.quiet = self.quiet + if (self.method.ismember("FixedStepSize")) { + opt.relax(1) + } else { + opt.linesearch(1) + opt.conjugategradient(100) + } + + return opt.didconverge() // Returns true if we converged + } + + build(outputdim=3) { // Build the mesh + if (self.method.ismember("StartRandom")) { + self.initialrandommesh() + } else self.initialgridmesh() + + for (iter in 0...self.maxiterations) { + if (!self.quiet) print "Mesh generator iteration ${iter}." + self.setupproblem() + + if (self.optimize()) break // Converged + + self.retriangulate(self.vertices()) + } + + //if (self.dim3) dim=3 + + for (i in 0...nv) { + var pt = Matrix(3) + var v = vert.column(i) + for (j in 0...dim) pt[j]=v[j] + out.append(pt) + } + + return out +} + /** Plots elements of a mesh * @param[in] mesh - mesh to plot * @param[in] selection - selection @@ -62,10 +79,10 @@ fn plotmesh(mesh, selection=nil, grade=nil, color=nil, filter=nil, transmit=nil if (islist(grade)) for (q in grade) glist.append(q) if (isnumber(grade)) glist.append(grade) - var vert = mesh.vertexmatrix() + var vert = vertto3d(mesh.vertexmatrix()) var flat = nil - var bb = bbox(vert) + var bb = bbox(mesh.vertexmatrix()) var size[bb.count()] for (x, i in bb) { size[i]=x[1]-x[0] @@ -78,7 +95,7 @@ fn plotmesh(mesh, selection=nil, grade=nil, color=nil, filter=nil, transmit=nil for (i in 0...np) { if (isselection(selection) && !selection.isselected(0, i)) continue // Skip unselected components - g.display(Sphere(vert.column(i), 0.025, color=plotcolorforelement(vcol, i), filter=filter, transmit=transmit)) + g.display(Sphere(vert[i], 0.025, color=plotcolorforelement(vcol, i), filter=filter, transmit=transmit)) } } @@ -90,7 +107,7 @@ fn plotmesh(mesh, selection=nil, grade=nil, color=nil, filter=nil, transmit=nil for (i in 0...nl) { if (isselection(selection) && !selection.isselected(1, i)) continue // Skip unselected components var el = lines.rowindices(i) - g.display(Cylinder(vert.column(el[0]), vert.column(el[1]), + g.display(Cylinder(vert[el[0]], vert[el[1]], aspectratio=0.05, color=plotcolorforelement(lcol, i), filter=filter, transmit=transmit)) } } @@ -121,7 +138,7 @@ fn plotmesh(mesh, selection=nil, grade=nil, color=nil, filter=nil, transmit=nil var el = faces.rowindices(i) var v[3] - for (p, i in el) v[i]=vert.column(p) + for (p, i in el) v[i]=vert[p] var normal if (flat) { normal=Matrix(bb.count()) From c4e91141201a0690d1ffbbfa059034628bdc3899 Mon Sep 17 00:00:00 2001 From: softmattertheory Date: Fri, 25 Feb 2022 20:14:29 -0500 Subject: [PATCH 2/2] Update to meshgen Update meshgen to improve documentation, improve benchmarking and add in a unit test --- benchmark/MeshGen/square.morpho | 10 ++++++++++ benchmark/benchmark.py | 4 ++++ examples/meshgen/ellipse.morpho | 2 +- examples/tactoid/tactoid.morpho | 2 +- morpho5/docs/meshgen.md | 18 ++++++++++++------ test/{packages => modules}/constants.morpho | 0 test/modules/meshgen/disk.morpho | 10 ++++++++++ .../meshtools/meshrefine.morpho | 0 .../meshtools/polyhedron.morpho | 0 .../meshtools/tetrahedron.mesh | 0 test/{packages => modules}/optimize/cg.morpho | 0 11 files changed, 38 insertions(+), 8 deletions(-) create mode 100644 benchmark/MeshGen/square.morpho rename test/{packages => modules}/constants.morpho (100%) create mode 100644 test/modules/meshgen/disk.morpho rename test/{packages => modules}/meshtools/meshrefine.morpho (100%) rename test/{packages => modules}/meshtools/polyhedron.morpho (100%) rename test/{packages => modules}/meshtools/tetrahedron.mesh (100%) rename test/{packages => modules}/optimize/cg.morpho (100%) diff --git a/benchmark/MeshGen/square.morpho b/benchmark/MeshGen/square.morpho new file mode 100644 index 00000000..268d2ff0 --- /dev/null +++ b/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 + diff --git a/benchmark/benchmark.py b/benchmark/benchmark.py index c84040ac..3178a912 100644 --- a/benchmark/benchmark.py +++ b/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): diff --git a/examples/meshgen/ellipse.morpho b/examples/meshgen/ellipse.morpho index 719eba87..e1d8e706 100644 --- a/examples/meshgen/ellipse.morpho +++ b/examples/meshgen/ellipse.morpho @@ -3,6 +3,6 @@ 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], quiet=false) +var mg = MeshGen(e0, [-2..2:0.2, -1..1:0.2]) var m = mg.build() Show(plotmesh(m, grade=1)) diff --git a/examples/tactoid/tactoid.morpho b/examples/tactoid/tactoid.morpho index 8dd20a83..02e67bd3 100644 --- a/examples/tactoid/tactoid.morpho +++ b/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) diff --git a/morpho5/docs/meshgen.md b/morpho5/docs/meshgen.md index 28041f3b..d69860c5 100644 --- a/morpho5/docs/meshgen.md +++ b/morpho5/docs/meshgen.md @@ -44,19 +44,25 @@ There are also a number of properties of a `MeshGen` object that can be set prio [showsubtopics]: # (subtopics) -## Domain +## Domain [tagdomain]: # (domain) -The `Domain` class is used to conveniently build a +The `Domain` class is used to conveniently build a domain by composing simpler elements. -var a = CircularDomain([-0.5,0], 1) -var b = CircularDomain([0.5,0], 1) +Create a `Domain` from a scalar function that is positive in the region of interest: -var dom = a.union(b).difference(c) + 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 = (s1.intersection(s2)).union(s3) + 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) diff --git a/test/packages/constants.morpho b/test/modules/constants.morpho similarity index 100% rename from test/packages/constants.morpho rename to test/modules/constants.morpho diff --git a/test/modules/meshgen/disk.morpho b/test/modules/meshgen/disk.morpho new file mode 100644 index 00000000..8f64da01 --- /dev/null +++ b/test/modules/meshgen/disk.morpho @@ -0,0 +1,10 @@ +// 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=true) +var m = mg.build() + +print ismesh(m) +// expect: true diff --git a/test/packages/meshtools/meshrefine.morpho b/test/modules/meshtools/meshrefine.morpho similarity index 100% rename from test/packages/meshtools/meshrefine.morpho rename to test/modules/meshtools/meshrefine.morpho diff --git a/test/packages/meshtools/polyhedron.morpho b/test/modules/meshtools/polyhedron.morpho similarity index 100% rename from test/packages/meshtools/polyhedron.morpho rename to test/modules/meshtools/polyhedron.morpho diff --git a/test/packages/meshtools/tetrahedron.mesh b/test/modules/meshtools/tetrahedron.mesh similarity index 100% rename from test/packages/meshtools/tetrahedron.mesh rename to test/modules/meshtools/tetrahedron.mesh diff --git a/test/packages/optimize/cg.morpho b/test/modules/optimize/cg.morpho similarity index 100% rename from test/packages/optimize/cg.morpho rename to test/modules/optimize/cg.morpho