55
$\begingroup$

Context

While studying manifold Learning I got interested in finding the eigenvectors of the Laplacian. (also in connection to this problem of solving the heat equation)

Following this and that amazing answer, I am interested in solving this Helmholtz equation in 3D

$ \triangledown^2 u(x,y,z) + k^2u(x,y,z) =0 \quad x,y,z \in \Omega\,, \quad u(x,y,z) = 0 \quad {\rm with}\quad x,y,z \in \partial\Omega $

where $\Omega =$ is some 3D boundary e.g. a ball, an ellipsoid, a regular 3D polygon etc.

I have played around with the 2D codes provided here to produce these first eigen modes of a snowflake (again beautiful code!):

Mathematica graphics

They look like this and are super-cool!

Imgur

but I would like to generalize their answer to 3D.

Question

How would one proceed in 3D, given that we have a 2D solution working?

Cheeky Attempt

I have modified slightly Mark McClure's code to make it 3D savvy, but I am no expert in this field

Needs["NDSolve`FEM`"]; helmholzSolve3D[g_, numEigenToCompute_Integer, opts : OptionsPattern[]] := Module[{u, x, y, z, t, pde, dirichletCondition, mesh, boundaryMesh, nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd, discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri, numEigen, res, eigenValues, eigenVectors, evIF},(*Discretize the region*) If[Head[g] === ImplicitRegion || Head[g] === ParametricRegion, mesh = ToElementMesh[DiscretizeRegion[g], opts], mesh = ToElementMesh[DiscretizeGraphics[g], opts]]; boundaryMesh = ToBoundaryMesh[mesh]; (*Set up the PDE and boundary condition*) pde = D[u[t, x, y, z], t] - Laplacian[u[t, x, y, z], {x, y, z}] + u[t, x, y, z] == 0; dirichletCondition = DirichletCondition[u[t, x, y, z] == 0, True]; (*Pre-process the equations to obtain the FiniteElementData in \ StateData*)nr = ToNumericalRegion[mesh]; {state} = NDSolve`ProcessEquations[{pde, dirichletCondition, u[0, x, y, z] == 0}, u, {t, 0, 1}, Element[{x, y, z}, nr]]; femdata = state["FiniteElementData"]; initBCs = femdata["BoundaryConditionData"]; methodData = femdata["FEMMethodData"]; initCoeffs = femdata["PDECoefficientData"]; (*Set up the solution*)vd = methodData["VariableData"]; sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}]; (*Discretize the PDE and boundary conditions*) discretePDE = DiscretizePDE[initCoeffs, methodData, sd]; discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]; (*Extract the relevant matrices and deploy the boundary conditions*) load = discretePDE["LoadVector"]; stiffness = discretePDE["StiffnessMatrix"]; damping = discretePDE["DampingMatrix"]; DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs]; (*Set the number of eigenvalues ignoring the Dirichlet positions*) pos = discreteBCs["DirichletMatrix"]["NonzeroPositions"][[All, 2]]; nDiri = Length[pos]; numEigen = numEigenToCompute + nDiri; (*Solve the eigensystem*) res = Eigensystem[{stiffness, damping}, -numEigen]; res = Reverse /@ res; eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]]; eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]]; evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors; (*Return the relevant information*){eigenValues, evIF, mesh}]

If I then define a 3D boundary 

 Ω = ImplicitRegion[0 <= x^2 + y^2 + z^2 <= 1, {x, y, z}]; RegionPlot3D[Ω, PlotStyle -> Opacity[0.5]]

Mathematica graphics

Naively this should give me the eigenmode:

{ev, if, mesh} = helmholzSolve3D[Ω, 1]; ev

but it actually crashes the kernel Mathematica (10.0.2).

Could anyone confirm this as a first step?

NB: Please do not loose sleep over this problem as it is mostly driven by curiosity :-)

PS: On the other hand I personally think this stuff is truly one of the best new useful features of Mathematica 10!

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ Regarding the spherically symmetric case: You're not going to get solutions that look like spherical harmonics, because solving the Dirichlet problem with FEM won't simultaneously solve the angular-momentum eigenvalue problem. Due to the degeneracies in a spherically symmetric system, the numerical solutions will be arbitrary linear combinations of spherical harmonics. You can see something analogous in the isotropic harmonic oscillator when solving with finite differences. One more reason not to use (Cartesian) FEM with spherical symmetry. $\endgroup$ Commented May 29, 2015 at 15:49
  • 1
    $\begingroup$ @Jens I was wondering about that. Thank you this is instructive for me. $\endgroup$ Commented May 29, 2015 at 16:53
  • 1
    $\begingroup$ @Jens: do you understand why the eigen-modes above do not quite rely on the symmetry of the boundary condition? Why .e.g for the second vector the two bumps are slightly tilted relative to the symmetric pattern of the snowflake? $\endgroup$ Commented May 29, 2015 at 17:01
  • 4
    $\begingroup$ @chris Absolutely. It's the same effect: the FEM solver has no knowledge about the discrete symmetry group of the boundary, which in this case is that of a hexagon, $D_{6h}$. Therefore, it cannot know how to label the eigenmodes by the irreducible representations of that group, and instead produces arbitrary linear combinations of them whenever there are degeneracies - which happens a lot here because $D_{6h}$ has two-dimensional irreducible representations (i.e., there are linearly independent, degenerate eigenfunctions related by symmetry operations). $\endgroup$ Commented May 29, 2015 at 17:25
  • 3
    $\begingroup$ So the clever way to approach such symmetric boundaries is to first reduce them to a fundamental domain with no remaining symmetries using the properties of the group. Essentially, the reflection axes then turn into boundaries with Dirichlet or Neumann conditions and you solve the wave problem only on a wedge-shaped segment of the hexagon. $\endgroup$ Commented May 29, 2015 at 17:29

2 Answers 2

Reset to default
43
+100
$\begingroup$

Version 11 has both symbolic and numeric eigensolvers, see here for an overview

Here is a slightly different way to do it. We write a function that converts any PDE (1D/2D/3D) into discretized system matices:

Needs["NDSolve`FEM`"] PDEtoMatrix[{pde_, Γ___}, u_, r__, o : OptionsPattern[NDSolve`ProcessEquations]] := Module[{ndstate, feData, sd, bcData, methodData, pdeData}, {ndstate} = NDSolve`ProcessEquations[Flatten[{pde, Γ}], u, Sequence @@ {r}, o]; sd = ndstate["SolutionData"][[1]]; feData = ndstate["FiniteElementData"]; pdeData = feData["PDECoefficientData"]; bcData = feData["BoundaryConditionData"]; methodData = feData["FEMMethodData"]; {DiscretizePDE[pdeData, methodData, sd], DiscretizeBoundaryConditions[bcData, methodData, sd], sd, methodData} ]

Example 1: An eigensolver is then something like this:

 {dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y], t] == Laplacian[u[t, x, y], {x, y}], u[0, x, y] == 0, DirichletCondition[u[t, x, y] == 0, True]}, u, {t, 0, 1}, {x, y} ∈ Rectangle[]]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; ContourPlot[#[x, y], {x, y} ∈ mesh, Frame -> False, ColorFunction -> ColorData["RedBlueTones"]] & /@ ifs

Mathematica graphics

This can be encapsulated as follows:

Helmholtz2D[bdry_, order_] := Module[{dPDE, dBC, sd, md, l, s, d, ifs, es, mesh, constraintMethod}, {dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y], t] == Laplacian[u[t, x, y], {x, y}], u[0, x, y] == 0, DirichletCondition[u[t, x, y] == 0, True]}, u, {t, 0, 1}, {x, y} ∈ bdry]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -order, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; {es, ifs, mesh} ]

Example 2: The the remaining problem in the question can then be solved like this:

RR = ImplicitRegion[ x^6 - 5 x^4 y z + 3 x^4 y^2 + 10 x^2 y^3 z + 3 x^2 y^4 - y^5 z + y^6 + z^6 <= 1, {{x, -1.25, 1.25}, {y, -1.25, 1.25}, {z, -1.25, 1.25}}]; mesh = ToElementMesh[RR, "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 31}] mesh["Wireframe"]

enter image description here

This creates a second order mesh with about 80T tets and 140T nodes. To discretize the PDE we use:

AbsoluteTiming[{dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y, z], t] == Laplacian[u[t, x, y, z], {x, y, z}], u[0, x, y, z] == 0, DirichletCondition[u[t, x, y, z] == 0, True]}, u, {t, 0, 1}, {x, y, z} ∈ mesh]; ] {6.24463, Null}

Get the eigenvalues and vectors:

l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; AbsoluteTiming[ First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] ] {13.484131`, {8.396796994677874`, 16.044484716974942`, 17.453692912770126`, 17.45703443132916`}}

Post process / visualize:

ifs = ElementMeshInterpolation[sd, #, "ExtrapolationHandler" -> {(Indeterminate &), "WarningMessage" -> False}] & /@ es[[2]];

Generate slices of the eigenfunctions in the region:

ctrs = Range @@ Join[mm = MinMax[ifs[[2]]["ValuesOnGrid"]], {Abs[Subtract @@ mm]/50}]; levels = Range[-1.25, 1.25, 0.25]; res = ContourPlot[ ifs[[2]][x, y, #], {x, -1.25, 1.25}, {y, -1.25, 1.25}, Frame -> False, ColorFunction -> ColorData["RedBlueTones"], Contours -> ctrs] & /@ levels; Show[{ RegionPlot3D[RR, PlotPoints -> 31, PlotStyle -> Directive[Opacity[0.25]]], Graphics3D[{Opacity[0.25], Flatten[MapThread[ Function[{gr, l}, Cases[gr, _GraphicsComplex] /. GraphicsComplex[coords_, rest__] :> GraphicsComplex[ Join[coords, ConstantArray[{l}, {Length[coords]}], 2], rest]], {res, levels}]]}] }, Boxed -> False, Background -> Gray]

enter image description here

Imgur

Example 3: As a self contained example, let us encapsulate the Helmholtz solver 

 Helmholtz3D[bdry_, order_] := Module[{dPDE, dBC, sd, md, l, s, d, ifs, es, mesh, constraintMethod}, {dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y, z], t] == Laplacian[u[t, x, y, z], {x, y, z}], u[0, x, y, z] == 0, DirichletCondition[u[t, x, y, z] == 0, True]}, u, {t, 0, 1}, {x, y, z} ∈ bdry]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; {es, ifs, mesh} ]

and consider

RR = ImplicitRegion[ x^4 + y^4 + z^4 < 1, {{x, -1, 1}, {y, -1, 1}, {z, -1, 1}}] {es, ifs, mesh} = Helmholtz3D[RR, nm=4]; mm = MinMax[ifs[[nm]]["ValuesOnGrid"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm]][Sequence @@ #]], Point[#]} &, mesh["Coordinates"]] // Graphics3D[#, Boxed -> False] &

Imgur

Example 4 Eigen modes on 3D Knot

Needs["NDSolve`FEM`"] f[t_] := With[{s = 3 t/2}, {(2 + Cos[s]) Cos[t], (2 + Cos[s]) Sin[t], Sin[s]} - {2, 0, 0}] v1[t_] := Cross[f'[t], {0, 0, 1}] // Normalize v2[t_] := Cross[f'[t], v1[t]] // Normalize g[t_, θ_] := f[t] + (Cos[θ] v1[t] + Sin[θ] v2[t])/2 gr = ParametricPlot3D[ Evaluate@g[t, θ], {t, 0, 4 Pi}, {θ, 0, 2 Pi}, Mesh -> None, MaxRecursion -> 4, Boxed -> False, Axes -> False]; tscale = 4; θscale = 0.5;(*scale roughly proportional to \ speeds*)dom = ToElementMesh[ FullRegion[2], {{0, tscale}, {0, θscale}},(*domain*) MaxCellMeasure -> {"Area" -> 0.001}]; coords = g[4 Pi #1/tscale, 2 Pi #2/θscale] & @@@ dom["Coordinates"];(*apply g*)bmesh2 = ToBoundaryMesh["Coordinates" -> coords, "BoundaryElements" -> dom["MeshElements"]]; emesh2 = ToElementMesh@bmesh2; RR = MeshRegion@emesh2

Mathematica graphics

{es, ifs, mesh} = Helmholtz3D[RR, nm = 4];

then

mm = MinMax[ifs[[nm]]["ValuesOnGrid"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm]][Sequence @@ #]], Point[#]} &, emesh2["Coordinates"]] // Graphics3D[#, Boxed -> False] &

Mathematica graphics

Improve this answer
$\endgroup$
25
  • $\begingroup$ @chris, added pictures and yes, this gives you an eigenmode. If it does not, have look if the ElementMeshInterpolation constructs the interpolating functions. $\endgroup$ Commented Jun 11, 2015 at 12:00
  • $\begingroup$ actually it was a matter of version of mathematica: your code works with 10.1 not with 10.0.2. $\endgroup$ Commented Jun 11, 2015 at 12:35
  • $\begingroup$ @chris, yes, I happen to have 10.1. $\endgroup$ Commented Jun 11, 2015 at 12:37
  • 1
    $\begingroup$ @chris I don't think it's a simple representation issue, and it's not just a diffusion equation. The Eigensystem step should yield the Helmholtz eigenvalues, but the boundary conditions aren't accounted for. If I understand correctly, the problem is simply that the result of DiscretizeBoundaryConditions is never actually used in setting up the matrices. What's missing is DeployBoundaryConditions. $\endgroup$ Commented Jun 11, 2015 at 22:01
  • 1
    $\begingroup$ @Jens that sounds like a more serious problem… so we are solving the PDE on the mesh but not satisfying the proper boundary condition at the edge of the mesh. It did strike me a bit when I look at the eigenmodes of the maps of France that somehow the borders seem to have little impact on their shape. $\endgroup$ Commented Jun 11, 2015 at 22:58
20
$\begingroup$

This slightly modified function

 Needs["NDSolve`FEM`"]; helmholzSolve3D[g_, numEigenToCompute_Integer, opts : OptionsPattern[]] := Module[{u, x, y, z, t, pde, dirichletCondition, mesh, boundaryMesh, nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd, discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri, numEigen, res, eigenValues, eigenVectors, evIF},(*Discretize the region*) If[Head[g] === ImplicitRegion || Head[g] === ParametricRegion, mesh = ToElementMesh[DiscretizeRegion[g,opts], opts], mesh = ToElementMesh[DiscretizeGraphics[g,opts], opts]]; boundaryMesh = ToBoundaryMesh[mesh]; (*Set up the PDE and boundary condition*) pde = D[u[t, x, y, z], t] - Laplacian[u[t, x, y, z], {x, y, z}] + u[t, x, y, z] == 0; dirichletCondition = DirichletCondition[u[t, x, y, z] == 0, True]; (*Pre-process the equations to obtain the FiniteElementData in \ StateData*)nr = ToNumericalRegion[mesh]; {state} = NDSolve`ProcessEquations[{pde, dirichletCondition, u[0, x, y, z] == 0}, u, {t, 0, 1}, Element[{x, y, z}, nr]]; femdata = state["FiniteElementData"]; initBCs = femdata["BoundaryConditionData"]; methodData = femdata["FEMMethodData"]; initCoeffs = femdata["PDECoefficientData"]; (*Set up the solution*)vd = methodData["VariableData"]; sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}]; (*Discretize the PDE and boundary conditions*) discretePDE = DiscretizePDE[initCoeffs, methodData, sd]; discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]; (*Extract the relevant matrices and deploy the boundary conditions*) load = discretePDE["LoadVector"]; stiffness = discretePDE["StiffnessMatrix"]; damping = discretePDE["DampingMatrix"]; DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs]; (*Set the number of eigenvalues ignoring the Dirichlet positions*) pos = discreteBCs["DirichletMatrix"]["NonzeroPositions"][[All, 2]]; nDiri = Length[pos]; numEigen = numEigenToCompute + nDiri; (*Solve the eigensystem*) res = Eigensystem[{stiffness, damping}, -numEigen]; res = Reverse /@ res; eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]]; eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]]; evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors; (*Return the relevant information*){eigenValues, evIF, mesh}]

works;

For a cuboid

 {ev, if, mesh} = helmholzSolve3D[Cuboid[], 4, MaxCellMeasure -> .05, AccuracyGoal -> 2]; ev

(* {38.2695,85.4791} *)

so that if we look at cross sections

Table[ContourPlot[if[[i]][x, y, 0.5], {x, 0, 1}, {y, 0, 1}], {i, 4}] // Partition[#, 2] & // GraphicsGrid

Imgur

or in 3D

Table[ContourPlot3D[if[[i]][x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Contours -> {-1/4, 1/4}],{i,4}]

Imgur

We can the boost up the resolution and look at higher order eigenmodes

{ev, if, mesh} = helmholzSolve3D[Cuboid[], 12, MaxCellMeasure -> 0.0025, "MaxBoundaryCellMeasure" -> 0.025]; Table[ ContourPlot[if[[i]][x, y, 0.5], {x, 0, 1}, {y, 0, 1}, Frame -> False, ColorFunction -> ColorData["RedBlueTones"]],

Imgur

Table[Image3D[ Table[if[[i]][x, y, z], {x, 0, 1, 0.025}, {y, 0, 1, 0.025}, {z, 0, 1, 0.025}], ColorFunction -> "RainbowOpacity"], {i, 1, 9}] // Partition[#, 3] & // GraphicsGrid

Imgur

For a ball

{ev, if, mesh} = helmholzSolve3D[Ball[], 12, MaxCellMeasure -> 0.025];

Then

 Table[ContourPlot[if[[i]][x, y, 0.], {x, -1, 1}, {y, -1, 1}, RegionFunction -> Function[{x, y}, Sqrt[x^2 + y^2] <= 1], Frame -> False, Axes -> False, ColorFunction -> ColorData["RedBlueTones"]], {i, 12}] // Partition[#, 2] & // GraphicsGrid

Mathematica graphics

Table[Image3D[ Table[If[x^2 + y^2 + z^2 < 1, if[[i]][x, y, z], 0], {x, -1, 1, 0.025}, {y, -1, 1, 0.025}, {z, -1, 1, 0.025}], ColorFunction -> "RainbowOpacity"], {i, 2, 11}] // Partition[#, 3] & // GraphicsGrid

Imgur

It also works for, say a cone 

 {ev, if, mesh} = helmholzSolve3D[Cone[], 4, MaxCellMeasure -> 0.25]

Ellipsoid

{ev, if, mesh} = helmholzSolve3D[Ellipsoid[{0, 0, 0}, {1, 2, 3}], 4, MaxCellMeasure -> 0.025] Table[ContourPlot[if[[i]][x, y, 0.1], {x, -1, 1}, {y, -2, 2}, RegionFunction -> Function[{x, y}, x^2/1 + y^2/2^2 < 1], Frame -> False, ColorFunction -> ColorData["RedBlueTones"], AspectRatio -> 1/2], {i, 1, 4}]

Mathematica graphics

Arbitrary boundary

It should work on say this cool boundary

RR = ImplicitRegion[ x^6 - 5 x^4 y z + 3 x^4 y^2 + 10 x^2 y^3 z + 3 x^2 y^4 - y^5 z + y^6 + z^6 <= 1, {x, y, z}];

Mathematica graphics

But unfortunately not with mathematica 10.0.2 because of this bug

If someone with 10.1 care to try this?

{ev, if, mesh} = helmholzSolve3D[RR, 4, MaxCellMeasure -> 0.25]

I get

MathKernel: /Jenkins/workspace/Component.TetGenLink.Linux-x86-64.10.0/Mathematica/Paclets/TetGenLink/CSource/TetGen/tetgen.cxx:19959: tetgenmesh::interresult tetgenmesh::scoutsubface(tetgenmesh::face*, tetgenmesh::triface*, int): Assertion `checksh.sh == dummysh' failed.

Update

with Mathematica 10.2 It just crashes the kernel.

Improve this answer
$\endgroup$
9
  • $\begingroup$ I think I get "TetGenTetrahedralize::reterr: Tetrahedralize returned an error, 3." Not really sure that everything loaded correctly, though I do get some nice images. $\endgroup$ Commented Jun 9, 2015 at 21:05
  • $\begingroup$ On 10.1, that is. $\endgroup$ Commented Jun 10, 2015 at 9:28
  • 1
    $\begingroup$ Very nice! To speed things up a bit you might want to use "Method" -> Arnoldi in Eigensystem. $\endgroup$ Commented Jun 10, 2015 at 19:13
  • $\begingroup$ @user21 it is much faster! you meant : Method -> "Arnoldi " I assume. $\endgroup$ Commented Jun 11, 2015 at 8:12
  • $\begingroup$ Ah, yes, of course. With this you should be able use finer meshes in 3D as well. $\endgroup$ Commented Jun 11, 2015 at 8:17

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.