Hi everyone, I am trying to implement a thermo-mechanical coupled constitutive model of bulk metallic glass into ABAQUS by VUMAT. The rate of temperature change is a fuction of current temperature and rate of plastic strain in the material with mathematic expression

dT/dt = A Laplacian(T) + B dP/dt

Where T is temperature, Laplacian(T) = d^2(T)/(dx)^2 + d^2(T)/(dy)^2, A and B are material coefficient, and dP/dt is plastic strain rate. In other words, the temperature change of the metallic glass are caused by two mechanisms, one is diffusion of heat (Laplacian), and another is heat generation due to plastic dissipation. In order to implement it into VUMAT, I convert the PDE above into forward Euler incremental equation as

dT= A(dt) Laplacian(T) + B dP

So far I manage to determine value of A, B and dP. If I could find the value of Laplacian(T) on each material point, I can solve the 2nd equation to determine value of dT as dt will be provided by ABAQUS. Then, I can use the dT value to update temperature at each material point by equation

TempNew = TempOld + dt

However, my problem is I dont know how to implement the Laplacian into VUMAT. It is because by using finite difference method, Laplacian of a function takes the form as

Laplacian ( f(x,y) ) = ( f(x-h,y) + f(x+h,y) + f(x,y-h) + f(x,y+h) - 4f(x,y) ) / h^2

which means in order to determine value of Laplacian(T) in a material point, I need to call back temperature value of 4 nearby material points, and the h value, which is the distance between the material point I want to study and its neighbor material point. Is anyone has experience in solving this kind of problems? Or my problem can simply be solved by using some built in function in ABAQUS? Any opinion or suggestion on how to dealing with Laplacian term using VUMAT are welcomed, and that would be great if can provide the example solution for this kind of problem.

Thank you.