Hot answers tagged

57

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 ...


41

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, ...


36

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 ...


24

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[...


24

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], ...


24

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 ...


22

Here's an approach that unions the primitives as rasters, meshes, and smooths. Data from OP: showmaze = Uncompress[FromCharacterCode @@ ImageData[Import["https://i.stack.imgur.com/XVJcP.png"], "Byte"]]; Primitives with an extended inlet and outlet: prims = CapsuleShape @@@ Cases[showmaze, _Cylinder, ∞]; prims = prims /. {{5., 5., 5.} -> {5.5, 5., 5.}, ...


22

There are many ways to do this, modifying, improving my method or doing a completely different thing. My goal here is to show a very basic idea that should give you a start. LocatorPane and Manipulate give means of interactive addition/deletion and dragging of points in 2D plane. The problem is how to add an edge -- there has to be interaction between 2 ...


20

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]] &&...


20

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 &}, ...


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 = ...


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....


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. &, {{-...


18

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} # & /...


18

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 ...


18

3D Example The problem with direct solvers is that starting in 3 dimensions, their performance for dealing with matrices stemming from PDEs drops rapidly. This is why I wanted to show at least one 3-dimensional example. As there is no immediate analogon for Loop subdivision of tetrahedral meshes, I use hexahedral meshes instead. Preparation Here are some ...


18

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 ...


17

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 ...


17

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, ...


17

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, ...


17

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 ...


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 ...


17

Here is what I was suggesting in comments: SeedRandom[] relaxed = Nest[ PropertyValue[{VoronoiMesh[#, {{-1, 1}, {-1, 1}}], {2, All}}, MeshCellCentroid] &, RandomReal[{-1, 1}, {45, 2}], 500 ]; mesh = VoronoiMesh[relaxed, {{-1, 1}, {-1, 1}}, MeshCellStyle -> {1 -> White}]; Then extract the cell primitives corresponding to the interior ...


16

Here's my take on Lloyd's algorithm. The code I present here should not be too hard to encapsulate into a routine; I have only elected to present it in this way so that I can animate the progress of the relaxation method: BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *) pts = RandomReal[{-1, 1}, {50, 2}]]; pl = With[{...


16

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] {...


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, ...


16

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