Inhomogeneous Heat Equation on Square Domain

This example shows how to solve the heat equation with a source term.

The basic heat equation with a unit source term is

`$\frac{\partial u}{\partial t}-\Delta u=1$`

This equation is solved on a square domain with a discontinuous initial condition and zero temperatures on the boundaries.

Create a square geometry centered at `x = 0` and `y = 0` with sides of length 2. Include a circle of radius 0.4 concentric with the square.

```R1 = [3;4;-1;1;1;-1;-1;-1;1;1]; C1 = [1;0;0;0.4]; C1 = [C1;zeros(length(R1) - length(C1),1)]; gd = [R1,C1]; sf = 'R1+C1'; ns = char('R1','C1')'; g = decsg(gd,sf,ns);```

Plot the geometry with edge and face labels.

```figure pdegplot(g,EdgeLabels="on",FaceLabels="on") axis([-1.1 1.1 -1.1 1.1]); axis equal```

Create an `femodel` object for transient thermal analysis and include the geometry.

```model = femodel(AnalysisType="thermalTransient", ... Geometry=g);```

Specify thermal properties of the material.

```model.MaterialProperties = ... materialProperties(ThermalConductivity=1,... MassDensity=1,... SpecificHeat=1);```

Specify internal heat source.

`model.FaceLoad = faceLoad(Heat=1);`

Set zero temperatures on all four outer edges of the square.

`model.EdgeBC(1:4) = edgeBC(Temperature=0);`

The discontinuous initial value is 1 inside the circle and zero outside. Specify zero initial temperature everywhere.

`model.FaceIC = faceIC(Temperature=0);`

Specify nonzero initial temperature inside the circle (Face 2).

`model.FaceIC(2) = faceIC(Temperature=1);`

Generate and plot a mesh.

```model = generateMesh(model); figure; pdemesh(model); axis equal```

Find the solution at 20 points in time between 0 and 0.1.

```nframes = 20; tlist = linspace(0,0.1,nframes); model.SolverOptions.ReportStatistics ='on'; result = solve(model,tlist);```
```97 successful steps 0 failed attempts 196 function evaluations 1 partial derivatives 21 LU decompositions 195 solutions of linear systems ```
`T = result.Temperature;`

Plot the solution.

```figure Tmax = max(max(T)); Tmin = min(min(T)); for j = 1:nframes pdeplot(result.Mesh,XYData=T(:,j),ZData=T(:,j)); axis([-1 1 -1 1 0 1]); Mv(j) = getframe; end```

To play the animation, use the `movie(Mv,1)` command.