We can make use of low level FEM programming tools for FE matrixes calculation but load vector let's calculate "by hands". For load vector we have:
$$b_i=\int_{0}^{L}g(t)\delta(x-a)\phi_i(x)dx=\begin{cases}g(t),& i=n_a\\0,& i\neq n_a\end{cases}$$
where $\phi_i(x)$-shape functions, $n_a$ - FE mesh node number which coordinate is $a$.
Input parameters of the model
Needs["NDSolve`FEM`"];
L = 0.2; (*dimension of calculation region*)
a = L/2; (*position of point heat source*)
rho = 2.7*10^3;
Cp = 1097.22;
k1 = 237; (*conductivity in 0=<x<=a*)
k2 = 237; (*conductivity in a<x<=L*)
rhocp = rho Cp;
g[t_] := 1 + Sin[\[Pi]*t] + Cos[\[Pi]*t] + Sin[2 \[Pi]*t] +
Cos[2 \[Pi]*t]
Mesh generation
h = L/100; (*spatial step size*)
Vert = Partition[Range[0., L, h], 1];
el = Table[{i, i + 1}, {i, 1, Length[Vert] - 1}];
Markers =
Table[m = el[[i]][[1]]; n = el[[i]][[2]];
If[ (Vert[[m]][[1]] + Vert[[n]][[1]])/2 < a, 0, 1], {i, 1,
Length[el]}];
mesh = ToElementMesh["Coordinates" -> Vert,
"MeshElements" -> {LineElement[el, Markers]}];
NodeSource = Floor[a/h] + 1;
FE matrices calculation
vd = NDSolve`VariableData[{"DependentVariables" -> {u},
"Space" -> {x}, "Time" -> t}];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];
initCoeffs =
InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-If[ElementMarker == 0, k1, k2]*
IdentityMatrix[1]}}, "LoadCoefficients" -> {{0.}},
"DampingCoefficients" -> {{rhocp}}];
methodData = InitializePDEMethodData[vd, sd];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
Solution of ODE system using Crank-Nicolson scheme
ProcDur = 4;
NumTimeStep = 400;
tau = ProcDur/NumTimeStep; (*temporal step size*)
theta = 0.5;
xlast = ConstantArray[0, {Length[Vert]}];
X = ConstantArray[0, {NumTimeStep + 1}];
X[[1]] = {0, xlast};
Do[
time = i*tau;
A = damping/tau + theta*stiffness;
b = (damping/tau - (1 - theta)*stiffness).xlast;
b[[NodeSource]] =
b[[NodeSource]] + theta*g[time] + (1 - theta)*g[time - tau];
x = LinearSolve[A, b, Method -> "Pardiso"];
xlast = x;
X[[i + 1]] = {time, x},
{i, 1, NumTimeStep}
]; // AbsoluteTiming
ClearAll[u]
vecu = Interpolation[X];
u[t_?NumericQ] := u[t] = ElementMeshInterpolation[{mesh}, vecu[t]]
Plot3D[u[t][x], {t, 0, ProcDur}, {x, 0, L}, PlotRange -> {0, 7*10^-5},
PlotPoints -> 100]
Let's examine integral heat balance
I1 = rhocp*NIntegrate[u[ProcDur][x], {x, 0, L}, AccuracyGoal -> 8]
I2 = Integrate[g[t], {t, 0, ProcDur}]
Abs[I1 - I2]/I1
3.99644
4
0.000890264
Heat balance is fulfilled with good accuracy.
Such approach also allows to set different PDE coefficient with the help of ElementMarker
when symmetry relatively $x=a$ fails.
