Heat Equation

Ω=RegionIntersection[
Cuboid[{0,0,0},{3,2,1}]
];
MyPic=BoundaryDiscretizeRegion[Ω,
Axes->True,
Boxed->True,
PlotTheme->"Lines",
BoxRatios->Automatic,
AxesLabel->{"x","y","z"}
];
Manipulate[
Show[MyPic, Graphics3D[{Orange,Sphere[{x,y,z},r]}]],
{{x,1.5},r,3-r},
{{y,1},r,2-r},
{{z,0.5},r,1-r},
{{r,0.3},0,0.5}
]

({{Domain,  , Ω⊂(ℝ^3)}, {Subdomain,  , V⊂Ω}, {External Heat Source,  , Q:Ω×[0,∞)->ℝ}, {Thermal Energy Density,  , e:Ω×[0,∞)->ℝ}, {Heat Flux,  , ϕ:Ω×[0,∞)->(ℝ^3)}})

Conservation of Heat Energy for V is 
	(d/dt)((Underscript[∫∫∫, V])e dV)=-(Underscript[∯, ∂V]) ϕ.n dA+(Underscript[∫∫∫, V])Q dV
which is the same as 
	(Underscript[∫∫∫, V])((∂e)/(∂t)) dV=-(Underscript[∫∫∫, V])Div[ϕ] dV+(Underscript[∫∫∫, V])Q dV
which if it holds for every subdomain V and all the functions are smooth implies
	((∂e)/(∂t)) =-Div[ϕ] +Q
Additional assumptions and notation

({{Temperature,  , u:Ω×[0,∞)->ℝ}, {Constant Specific Heat,  , c∈(SuperPlus[ℝ])}, {Constant Mass Density,  , ρ∈(SuperPlus[ℝ])}, { ,  , e=ρ c u }, {Constant Thermal Conductivity,  , K∈(SuperPlus[ℝ])}, {Fourier "Law",  , ϕ=-K ∇u}})

Plug em in and you get the Heat Equation

      ρ c (Subscript[u, t]) =-K Div[∇u] +Q
      
 some folks write this as 
 	ρ c (Subscript[u, t]) =-K ∇.(∇u) +Q
 and some folks write this as 
 	ρ c (Subscript[u, t]) =-K Δu +Q
 they all mean the same thing.  Δu=Div[∇u] is the Laplacian

Other Conservation Laws

Lots of other things produce balance laws that look similar.  The key is that there is frequently a “flux” term that is turned into a volume integral using the divergence theorem. All such equations are called conservation laws. 
“Conservation of Energy”
“Conservation of Momentum”
“Conservation of Contaminants”
…
For right now one is enough.

Boundary Conditions

We need to decide what we are going to do on the boundaries of Ω. At each point x∈∂Ω on the boundary it is natural to think you can specify either u[x] or the heat flux ϕ[x].  If the point is pushed against an ice cube that is just on the edge of melting you are putting u=0.  If you insulate the point then you are making ϕ=0.  Putting the heat flux zero is the same as setting ((∂u)/(∂n))=0 because of the Fourier Law derivative as you go.

Specifying the u value is a Dirichlet boundary condition. 
Specifying the ∂u/∂n value is a Neumann Condition.
Other boundary conditions are possible.

Laplacian in Other Coordinate Systems

Mathematica knows lots of coordinate systems.  Here is the Laplacian for ones you already know.

{
Laplacian[u[x,y,z],{x,y,z},"Cartesian"] ,
Laplacian[u[r,θ,z],{r,θ,z}, "Cylindrical"],
Laplacian[u[r,θ,ϕ],{r,θ,ϕ}, "Spherical"]
}

There are lots of others.  In 2 dimensions we have Bipolar (one parameter), Cartesian, Confocal (two parameters), Elliptic (one paramater), Polar, and Planar Parabolic.

Clear[α,β,a]
TwoDCoords={
"Polar"->"Cartesian",
"PlanarParabolic"->"Cartesian",
{{"Elliptic",a}}->"Cartesian",
{{"Confocal",{α,β}},2}->"Cartesian",
{"Cartesian",2}->"Cartesian",
{{"Bipolar",a}}->"Cartesian"
};

TableForm[
Table[
CoordinateTransform[name,
CoordinateChartData[name[[1]],"StandardCoordinateNames"]
],
{name, TwoDCoords}],
TableHeadings->{
TwoDCoords[[All,1]],
{"x","y"}}
]

One obvious question is what sort of shapes are these good for.

{α,β}={2,1};
TabView[{
"Polar"->ParametricPlot[
{r Cos[θ],r Sin[θ]},
{r,0, 1},
{θ,0, 2π},
Mesh->Full],
"PlanarParabolic"->ParametricPlot[
{1/2 ((u^2)-(v^2)), u v},
{u,-1, 1},
{v,0, 1},
Mesh->Full],
"Elliptic a=1"->ParametricPlot[
{Cos[θ]Cosh[μ], Sin[θ] Sinh[μ]},
{θ,0, 2π},
{μ,0.2, 2},
Mesh->Full],
"bipolar a=1"->ParametricPlot[
a=1;
{ ((a Sinh[ξ])/(Cosh[ξ]-Cos[θ])),((a Sin[θ])/(Cosh[ξ]-Cos[θ]))},
{ξ,0.2,1.2},
{θ,0.3, 1.4},
Mesh->Full,
PlotRange->All],
"confocal α=2 β=1"->ParametricPlot[
{ (Sqrt[(((λ+α)(μ+α))/(α-β))]),(Sqrt[(((λ+β)(μ+β))/(β-α))])},
{λ,-0.9 α,-0.9β},
{μ,-0.9 β, 3},
Mesh->Full,
PlotRange->All]
}]

There are a lot more in 3D.

CoordinateChartData[{All,"Euclidean",3}]//TableForm

We are going to stick to the ones we are more comfortable with: Cartesian; Polar, Cylindrical and Spherical.

Analytical Capabilities of CAS

With an ODE you get arbitrary constants with a PDE you get arbitrary functions.

Here is a constant coefficient first order PDE.  The solution has a single arbitrary function.

u[x1,x2]/.DSolve[
3D[u[x1,x2],x1]+5D[u[x1,x2],x2]==x1 ,
u[x1,x2],
{x1,x2}][[1]]

If you add an initial condition it determines the function.

u[x1,x2]/.DSolve[{
3D[u[x1,x2],x1]+5D[u[x1,x2],x2]==x1 ,
u[x1,0]==Sin[x1]
},
u[x1,x2],
{x1,x2}][[1]]

It is almost always worth trying a CAS. This is what Mathematica looks like if it can not do anything worth while.

DSolve[{
 D[u[x,t],t]==D[u[x,t],{x,2}]
},
u,
{x,t}]

Mesh

Edge Disc then area

Mesh1=TriangulateMesh[
BoundaryDiscretizeGraphics[
Polygon[{{0,0},{2,0},{2,1},{0,1}}]
],
MaxCellMeasure->∞,
Axes->True
]

C1=MeshCoordinates[Mesh1];
P1=MeshCells[Mesh1,2];
n1= Length[C1];
MeshRegion[
Join[
C1,
C1/.{x_Real,y_Real}:>{x,-y},
C1/.{x_Real,y_Real}:>{-x,y},
C1/.{x_Real,y_Real}:>{-x,-y}
],
Join[
P1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+n1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+2n1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+3n1
] ,
Frame->True
]

Straight area

?DiscretizeRegion

MeshQualityGoal

Mesh1=DiscretizeRegion[
Polygon[{{0,0},{2,0},{2,1},{0,1}}],
MaxCellMeasure->∞,
MeshQualityGoal->"Minimal",
Axes->True
]

C1=MeshCoordinates[Mesh1];
P1=MeshCells[Mesh1,2];
n1= Length[C1];
MeshRegion[
Join[
C1,
C1/.{x_Real,y_Real}:>{x,-y},
C1/.{x_Real,y_Real}:>{-x,y},
C1/.{x_Real,y_Real}:>{-x,-y}
],
Join[
P1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+n1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+2n1,
P1/.{i1_Integer,i2_Integer,i3_Integer}:>{i1,i2,i3}+3n1
] ,
Frame->True
]