Hot answers tagged

52

Background Details about multigrid solvers can be found in this pretty neat script by Volker John. That's basically the source from which I drew the information to implement the V-cycle solver below. In a nutshell, a multigrid solver builds on two ingredients: A hierarchy of linear systems (with so-called prolongation operators mapping between them) and a ...


40

We need quite a bit of preparation. In the first place we need methods to compute cell adjacency matrices from here. I copied the code for completeness. CellAdjacencyMatrix[R_MeshRegion, d_, 0] := If[MeshCellCount[R, d] > 0, Unitize[R["ConnectivityMatrix"[d, 0]]], {} ]; CellAdjacencyMatrix[R_MeshRegion, 0, d_] := If[MeshCellCount[R, d] > 0, ...


26

You are combining the images in the form Show[Graphics[simplePrimitives], complicatedRegionPlot] The options in the resulting figure are inherited from the first term, namely Graphics[simplePrimitives]. This does not include the "TransparentPolygonMesh" -> True generated by RegionPlot. You see the mesh as a result. If you combine things as follows: ...


25

As announced before, here my take on the mean curvature flow for surfaces. The code is rather lengthy and I tried to recycle as much as possible from this post about finding minimal surfaces (solving Plateau's problem). Please find the code at the end of this post. Background Mean curvature flow is the $L^2$-gradient flow of the area functional on the ...


23

It seems to me that the logo has three semitransparent layers of triangle meshes. One can start with discretized sphere reg = DiscretizeGraphics[Sphere[], MaxCellMeasure -> {"Length" -> 0.8}] Or with Simon's Geodesate. Then the function for disks in 3D is helpful disk[pos_, {nx_, ny_, nz_}, r_, n_: 16] := With[{θ = ArcTan[Sqrt[nx^2 + ny^2], nz], ...


21

Here's a possible approach. First use TetGen to tetrahedralize the data: Needs["TetGenLink`"] {pts, tetrahedra} = TetGenDelaunay[data3D]; Next define a function to compute the radius of the circumsphere of a tetrahedron (formula from Wikipedia) csr[{aa_, bb_, cc_, dd_}] := With[{a = aa - dd, b = bb - dd, c = cc - dd}, Norm[a.a Cross[b, c] + b.b Cross[...


20

There is a neat (and widely known) trick to produce n approximately evenly distributed points in a disk, or on any surface of revolution. Points are placed on a spiral, each time turning by the "golden angle". The reference is: Vogel, H (1979). "A better way to construct the sunflower head". Mathematical Biosciences. 44 (44): 179–189.4. We can adapt the ...


20

Based on the tutorial here QuadElement meshes behave exactly the same as TriangleElement meshes, with the exception that, for linear quad elements, four incidents per element are needed, and, for quadratic elements, eight incidents per element are needed. here is another possible answer: Needs["NDSolve`FEM`"] rectangleA = {{0, 4}, {8, 8}}; rectangleB = ...


20

Lets start by making a symbolic representation of maze, with BooleanRegion using graphic primitives. That will be our basis for discretization by any method. Number of vertices of graph should be small to make experimentation easier. (* "IGraphM" package downloaded from https://github.com/szhorvat/IGraphM *) Get["IGraphM`"] n = 3; (* Number of vertices per ...


19

First, you can generate your random points like so: SeedRandom[1]; pts = RandomReal[{0, 12}, {100, 2}]; The DelaunayTriangulation command returns an adjacency list representation of the triangulation. Needs["ComputationalGeometry`"]; dt = DelaunayTriangulation[pts]; dt // Column This says that the first point should be connected to the 2nd, the 24th, etc....


19

first part..i had lying around.. poly = Random[Real, {1, 2}] {Cos[#], Sin[#]} & /@ Sort[Table[Random[Real, {0, 2 Pi}], {5}]] isLeft[P2_, {P0_, P1_}] := -Sign@Det@{P2 - P0, P1 - P0}; pinpoly[p_, poly_] := Module[{ed},(*winding rule*) ed = Partition[Append[poly, poly[[1]]], {2}, 1]; Count[ed,pr_ /; (pr[[1, 2]] <= p[[2]] < pr[[2, 2]] &&...


19

The blue line occurs at the edge of the function, where ϕ wraps from 2π to 0. We can get rid of it by adding BoundaryStyle -> None: SphericalPlot3D[ Abs[.5 + Sin[2 ϕ]/2] Sin[θ] + Abs[.5 + Sin[2 (ϕ + π/2)]/2] Sin[θ], {θ, 0, π}, {ϕ, 0, 2 π}, PlotStyle -> {Opacity[0.3], Yellow}, BoxRatios -> {1, 1, 1/2}, MeshFunctions -> {#3 &}, ...


19

Catmull-Clark Subdivision Indeed, I have some code for Catmull-Clark subdivision and I planned to post it here for quite some time. This seems to be a good opportunity. The code is optimized for performance, so it involves a lot of CompiledFunction and SparseArray hacks. I am sorry if you find it somewhat unidiomatic. CatmullClarkSubdivisionMatrix creates ...


18

Table[drawtriangulation[mesh @@ example, First@example, AspectRatio -> Automatic], {example, {circle, circle34, ellipseeye}}] // GraphicsRow Calculating specifications for these examples: (* distance function, bounding box, fixed points, number of initial points, max iterations, min triangle quality *) circle = {Sqrt[#1^2 + #2^2] - 1. &, {{-...


17

You could use the (undocumented) option Method -> {"TransparentPolygonMesh" -> True} for this, e.g. Show[Graphics[Point[{p1, p2}]], RegionPlot[{d[{x, y}, p1, M1] < d[{x, y}, p2, M2], d[{x, y}, p1, M1] > d[{x, y}, p2, M2]}, {x, -4, 4}, {y, -4, 4}], Method -> {"TransparentPolygonMesh" -> True}] which produce


17

Here is a different approach. The idea is to first generate a second order triangle mesh. Next a center coordinate is added in every triangle. Then we split every triangle into three first order quads making use of the newly added center coordinate and the initial second order triangle incidents. So we add a node in the center of the element and create the ...


17

Using the Geodesics in Heat Algorithm implemented here, we can calculate the distances of all vertices on the surface to a given vertex. By repeating this algorithm on a selected subset of vertices on the surface, we can calculate readily how close all the other vertices on these surfaces are, and label them according to which selected vertex they are ...


16

There is good news and bad news. So right now there is no option to make MeshRegion use Hex elements. However, you can use ElementMesh to do that: Needs["NDSolve`FEM`"] em = ToElementMesh[Cuboid[{0, 0, 0}, {1, 1, 1}], MaxCellMeasure -> 1/4^3] (* ElementMesh[{{0., 1.}, {0., 1.}, {0., 1.}}, {HexahedronElement[ "<" 64 ">"]}] *) The bad news is ...


16

Update: With the function top defined in the original post you can replicate all the cool things you see in rm-rf's answer in the linked Q/A. For example, with a slight modification of gr1, i.e., Graphics3D[hexTile[20, 20] /. Polygon[l_] :> {Directive[Orange, Opacity[0.8], Specularity[White, 30]], Polygon[l], Polygon[{Pi/5, 0} + {-1, 1} # & /...


16

Method 1: Construct mesh elements manually We can triangulate a periodic quad-lattice on the surface: n = {180, 20}; (* number of points in each direction *) pts = Table[ g[4. Pi/n[[1]] t, 2. Pi/ n[[2]] θ], {t, n[[1]]}, {θ, n[[2]]}]; idcs = {{{1, 2, 4}, {1, 4, 3}}, {{1, 2, 3}, {2, 4, 3}}}; (* for a diamond pattern *) tri = 1 + Array[Function[quad, ...


16

What you want is the 2D alpha shape to try to get close to the outline you seek. Of course, it's no longer a Delaunay triangulation since you're deleting certain polygons from the DelaunayMesh. We'll adopt my answer from this post. Here it is for a 2D point set: alphaShapes2D[points_, crit_] := Module[{alphacriteria, del = Quiet @ DelaunayMesh @ points, ...


16

That's not possible in general in V11.0. The only primitive from which a quad mesh can be generated is ToElementMesh[Rectangle[]] or ToElementMesh[FullRegion[2], {{0, 1}, {0, 1}}] So to get a quad mesh for a more complicated shape you'd need to program a bit and use ToElementMesh["Coordinates"->..., "MeshElements"->{QuadElement[....]}] As ...


16

this is a bit of a hack but it works for this example: generate separate rectangle meshes -- note the connected edge nodes must be exactly coincident. << NDSolve`FEM` m1 = ToElementMesh[Rectangle[{-10, 0}, {0, 10}], "MeshElementType" -> QuadElement, MaxCellMeasure -> 4, "MeshOrder" -> 1]; m2 = ToElementMesh[Rectangle[{-20, 10}, {0, ...


15

Here's my go at it. This tells you if two line segments intersect (unless they lie on the same line, in which case it fails horribly): ClearAll[segmentsIntersect]; segmentsIntersect[{a_, b_}, {p_, q_}] := Module[{s, t, soln}, soln = NSolve[a + t (b - a) == p + s (q - p), {s, t}]; If[Length@soln == 0, False, (0 <= s <= 1 && 0 <= t ...


15

It might be easier to use TriangularSurfacePlot3D to find the Delaunay triangulation of the points. For example, Needs["ComputationalGeometry`"]; triangles[points_] := Module[{pl}, pl = TriangularSurfacePlot[ArrayPad[points, {{0, 0}, {0, 1}}]]; Cases[pl, Polygon[a_] :> Flatten[(Position[points, #[[{1, 2}]]] & /@ a)], Infinity]] Graphics[...


15

You have to create your own mesh and you have to convert your u and v to mesh interpolations. (In the example in the documentation, NDSolveValue does this itself in constructing uif, vif.) Example: Needs["NDSolve`FEM`"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}]; u = Function[{x, y}, x (y - 0.5)/25]; v = Function[{x, y}, -x^2/50]; uif = ...


15

I discarded my previous approach to generate cubes, then fuse them together, since it seems to do a lot of wasted work. Instead, I propose here my version of a cartesian mesher. The approach is conceptually the same as the one delineated by Zviovich, but I wasn't entirely satisfied with his results, as it seems to me that his process still leads to ...


15

There are no surface element shape functions. There are, however, the normal shape functions. Load the package: Needs["NDSolve`FEM`"] This gives you the shape functions for implemented elements (see documentation) elementOrder = 1; ElementShapeFunction[TriangleElement, elementOrder][r, s] {1 - r - s, r, s} ElementShapeFunction[TriangleElement, 2][r, s] {...


15

I think, I found a general and even faster way, but I haven't tested it for $1$- or $3$-dimensional MeshRegions. The following function computes the cell-vertex-adjacency matrix A. Two $d$-dimensional cells ($d>0$) are adjacent if they share at least $d$ common points. We can find these pairs by looking for entries $\geq d$ in A.Transpose[A]. ToPack = ...


Only top voted, non community-wiki answers of a minimum length are eligible