PDE Model Boundary Conditions

Hi Everyone,

I have a question regarding the matlab PDE Model Boundary Conditions.

Currently, I am working with a 2D mesh with all egdes being adiabatic.

Files are available to be downloaded here:

https://drive.google.com/file/d/1EmpT8JeGhxN8ORLXpciN2HBSPIJK3ra_/view?usp=sharing

Thus,

        applyBoundaryCondition(model,'neumann','Edge',1,'q',0,'g',0);%Bottom Edge
        applyBoundaryCondition(model,'neumann','Edge',4,'q',0,'g',0);%Left Edge
        applyBoundaryCondition(model,'neumann','Edge',2,'q',0,'g',0);%Right Edge
        applyBoundaryCondition(model,'neumann','Edge',3,'q',0,'g',0);%Top Edge

However, I have a non-constant Boundary Condition on the Face of the mesh:

        applyBoundaryCondition(model,'neumann','face',1,'q',0,'g',@myfun);

The code is here:

load('dat.mat')
load('Q(x,y,t).mat')
%%
Interptest(Q,dat)
%%
function Interptest(Q,dat)
%Variables:
    for i=1:201
        t0_{i,1} = i;
    end
%
    tglob=0;
    tlist = [0;0.1;0.2;0.3];
    framerate = 100; %Frequency in Hertz
    x_width = 3.2e-2;
    y_height = 3.2e-2;
    xpix = 501;
    rat = x_width./xpix;
    freq = 0.01;
    ypix = xpix;
    x1 = linspace(0,x_width,xpix);
    y1 = linspace(0,y_height,ypix);
    tt_vec=linspace(0,2,201);
    [XX1,YY1] = meshgrid(x1,y1);
    [XX2,YY2] = meshgrid(linspace(0,x_width,xpix),linspace(0,y_height,ypix));
    [XXorig, YYorig] = meshgrid(linspace(0,x_width,size(Q{5,1}{1,1},2))',linspace(0,y_height,size(Q{5,1}{1,1},1))');
    [XX0, YY0] = meshgrid(linspace(0,x_width,size(dat.transftempK{1,1},2)),linspace(0,y_height,size(dat.transftempK{1, 1},1)));
     [XX3, YY3,tt_vec] = meshgrid(linspace(0,x_width,size(dat.transftempK{1,1},2)),linspace(0,y_height,size(dat.transftempK{1, 1},1)),linspace(0,2,201));
    Xvec0 = XX0(:);
    Yvec0 = YY0(:);
    XXorig_vec = XXorig(:);
    YYorig_vec = YYorig(:);
    
    
    Xvec1 = XX3(:);
    Yvec1 = YY3(:);
    ttvec1 = tt_vec(:);


%     tt_vec=linspace(0,2,201);
    Tu = 293.15;
    Tb = 77.15;
    Rho = 4.43*10^3;
    thickness = 1.5/1000;
    cvfunc = @(T) 413.5+0.4372.*T+1.443.*10.^-4.*T.^2;
    Sigma = 5.670.*10.^-8;
    lambdafunc = @(T) 2.369+1.316.*T.*10.^-2;
    tglob=0;

    for i=5:length(t0_)
        sprintf('Dataset %i of %i',i,length(t0_))
        model = createpde(3);
        geo = [3, 4, 0, x_width, x_width, 0, 0, 0, y_height, y_height]';
        sf = 'R1';
        ns = [82 49]';
        dl = decsg(geo,sf,ns);
        geometryFromEdges(model,dl);
        pdegplot(model,'Facelabels','on','EdgeLabels','on');
        hmax = 1e-3;
        mesh = generateMesh(model,'Hmax',hmax,'GeometricOrder','quadratic');

        m = 0;
        d = @(location,state) Rho.*cvfunc(state.u);
        c = @(location,state) lambdafunc(state.u);
        f = 0;
        a = 0;
%         x=@myfun;
        
        applyBoundaryCondition(model,'neumann','Edge',1,'q',0,'g',0);%Bottom Edge
        applyBoundaryCondition(model,'neumann','Edge',4,'q',0,'g',0);%Left Edge
        applyBoundaryCondition(model,'neumann','Edge',2,'q',0,'g',0);%Right Edge
        applyBoundaryCondition(model,'neumann','Edge',3,'q',0,'g',0);%Top Edge
        applyBoundaryCondition(model,'neumann','face',1,'q',0,'g',@myfun);

        specifyCoefficients(model,'m',m,...
            'd',d,...
            'c',c,...
            'a',a,...
            'f',f);
        model.SolverOptions.ReportStatistics='on';
        T0{i,1} = dat.transftempK{t0_{i,1},1};
        Tvec0 = T0{i,1}(:);
        T00 = griddata(Xvec0, Yvec0, Tvec0, model.Mesh.Nodes(1,:), model.Mesh.Nodes(2,:))';
        dummyResult = createPDEResults(model,[T00,T00],0:1*10^-12:1*10^-12,'time-dependent');
        setInitialConditions(model,dummyResult);
        %Interpolation unchanged here!

        tsel{i,1} = t0_{i,1}+tlist./freq;
        tsel{i,1} = tsel{i,1};
        disp('solve pde');
        results{i,1} = solvepde(model,tlist);
        for j=1:length(tlist)
            sprintf('Sim %i of %i \n',j,size(tlist,1))
            T_Sim_interp{i,1}{j,1} = rot90(flipud(griddata(results{i,1}.Mesh.Nodes(1,:)',results{i,1}.Mesh.Nodes(2,:)',results{i,1}.NodalSolution(:,j),XX1,YY1,'cubic')),-1);
        end
        

    end

function q = myfun(location,state)
    if(tglob~=state.time)
    tglob=state.time
    end
    q=0;
% q=interp3(Xvec1,Yvec1,ttvec1,state.u,location.x,location.y,state.time);
% disp(location.x);
% disp(location.y);
% disp(state.time);
end
end

Apparently, I have troubleshoot it and the myfun works when it is on the edges. It does not work when I place:

    applyBoundaryCondition(model,'neumann','Face',1,'q',0,'g',@myfun);

Did I make a mistake regarding createpde(1) or should I work with a 3D model instead?

I greatly appreciate if anyone could help regarding this.

Best regards,

Justin