classdef DiscretizedPDEModel
%DiscretizedPDEModel - Finite Element form of PDE
%This undocumented function may be removed in a future release.

%       Copyright 2015-2021 The MathWorks, Inc.

    properties
        % input parameters
        p, e, t, coefstruct;
        % Current values of global matrices
        K, A, F, Q, G, H, R;
        B, ud, Nu, Or;
        totalNumEqns, numConstrainedEqns;
        thePde;
        nrp;        
        % Flags indicating whether coefficient is dependent on u
        qu, gu, hu, ru, cu, au, fu;
        IsStationaryNonlinear;
    end
    
    
    methods(Access=public)
        function obj=DiscretizedPDEModel(thePde,p,e,t,coefstruct,u0,constCoeffsBCs)
            obj.qu=false; obj.gu=false; obj.hu=false; obj.ru=false; obj.cu=false;
            obj.au=false; obj.fu=false;
            obj.p = p;
            obj.e = e;
            obj.t = t;
            obj.thePde = thePde;
            obj.nrp = size(p,1);
            obj.coefstruct = coefstruct;
            tdummy = 0;
            obj = obj.initialDiscretization(u0,tdummy);
            
            if ~constCoeffsBCs
                obj = obj.checkSpatialCoefsForUDependence(u0,tdummy);
            end
            
            obj.IsStationaryNonlinear = (obj.qu || obj.gu || obj.hu || obj.ru ...
                        || obj.cu || obj.au || obj.fu);
        end
    end
    
    methods(Access=protected)
        
        function self=initialDiscretization(self,u0,tdummy)
            
            %checkFuncDepen determine whether coefficients are functions on u
            requiredMats =  'KAFQGHR';
            femat0 =  self.thePde.assembleSelectedFEMatrices(self.p, self.t, self.coefstruct, u0, tdummy, requiredMats, false);

            for matID = requiredMats
                self.(matID) = femat0.(matID);
            end
            
            [self.Nu,self.Or]=pdenullorth(self.H);
            if size(self.Or,2)==0
                ud0=zeros(size(self.K,2),1);
            else
                ud0=full(self.Or*((self.H*self.Or)\self.R));
            end
            
            [self.totalNumEqns, self.numConstrainedEqns]=size(self.Nu);
            
            self.B=self.Nu;
            self.ud=ud0;
        end
        
        function self = checkSpatialCoefsForUDependence(self,u0,tdummy)
            u1 = NaN*ones(size(u0,1), 1);
            femat =  self.thePde.assembleSelectedFEMatrices(self.p, self.t, self.coefstruct, u1, tdummy, 'KAFQGHR', false);
            % check if coefficients are functions of u
            self.cu = any(isnan(femat.K(:)));
            self.au = any(isnan(femat.A(:)));
            self.fu = any(isnan(femat.F(:)));
            self.qu = any(isnan(femat.Q(:)));
            self.gu = any(isnan(femat.G(:)));
            self.hu = any(isnan(femat.H(:)));
            self.ru = any(isnan(femat.R(:)));
        end
        
    end
end

