classdef surfaceToSurfaceSettings

% Copyright 2023 The MathWorks, Inc.    

    properties (SetAccess=private)
        ViewFactors {mustBeNumeric,mustBeFinite,mustBeNonnegative}
        ViewFactorMethod
    end

    properties
        Enclosures % Public to allow edits on PerfectEnclosure.
    end

    properties (Dependent, SetAccess=private)
        ParticipatingEnclosures
    end
    
    properties(Hidden,SetAccess=private)
        Area
        NumRays
        NumRaySources
        DefinedEnclosures
        ParticipatingEnclosures_I
        ElementsToNodesProjection
        NodesToElementProjection
        Mesh
    end

    methods
        function self = surfaceToSurfaceSettings(options)
            arguments
                options.Enclosures (1,:) dictionary
                options.ParticipatingEnclosures (1,:) string {mustBeNonempty}
                options.Mesh (1,1) pde.FEMesh
                options.ViewFactorMethod string {mustBeScalarOrEmpty,mustBeMember(options.ViewFactorMethod,...
                                                   ["areaintegral","montecarlo"])} = "areaintegral"
                options.NumRays {mustBeNumeric,mustBeInteger,mustBePositive,mustBeNonsparse}
                options.NumRaySources {mustBeNumeric,mustBeInteger,mustBePositive,mustBeNonsparse}
            end

            for k=string(fields(options))'
                self.(k) = options.(k);
            end
            enclosureNames = keys(self.Enclosures);
            self.DefinedEnclosures = enclosureNames';
            allBoundaryFacets = [self.Enclosures(enclosureNames).BoundaryFacet];
            [F, Area,B] = generateViewFactors(self,[self.Enclosures(enclosureNames).Face]);
            self.ViewFactors = F;
            self.Area = Area;
            self.ElementsToNodesProjection = B;

            numFacets = size(allBoundaryFacets,2);
            rowIdx = repmat((1:numFacets)',1,size(allBoundaryFacets,1));
            colIdx = allBoundaryFacets';
            val = 1/size(allBoundaryFacets,1);
            neqn = size(self.Mesh.Nodes,2);
            self.NodesToElementProjection = sparse(rowIdx,colIdx,val,numFacets,neqn);
        end

        function val = get.ParticipatingEnclosures(s2sObj)
            val = s2sObj.ParticipatingEnclosures_I;
        end

        function s2sObj = set.ParticipatingEnclosures(s2sObj,val)
            s2sObj.ParticipatingEnclosures_I = val;
        end

        function obj = set.Enclosures(obj,newEnc)
            if ~isa(newEnc,"dictionary") || ~isConfigured(newEnc) 
                error(message("pde:femodel:InvalidEnclosuresDictionary"));
            end

            [keyType,valueType] = types(newEnc);
            if keyType ~="string" || valueType ~= "enclosureDefinition"
                error(message("pde:femodel:InvalidEnclosuresDictionary"));
            end

            obj.Enclosures = newEnc;
        end
    end

    methods(Hidden)

        function [F, Area,B] = generateViewFactors(self, faceIDs)
            allBoundaryFacets = self.Mesh.boundaryFacets("Face",faceIDs);
            NodesMat = self.Mesh.Nodes;

            orig_state = warning;
            restoreWarningState = onCleanup(@() warning(orig_state));
            warning("off","MATLAB:triangulation:PtsNotInTriWarnId");
            allTriCompute = triangulation(allBoundaryFacets(1:3,:)',NodesMat');

            numTri = size(allTriCompute.ConnectivityList,1);
            F = zeros(numTri);

            boundaryDim = size(NodesMat,1)-1;
            numFacetNodes = size(allBoundaryFacets,1);

            if(numFacetNodes == 3)
                sfCalc = @pde.internal.elements.shapeTri3;
                dsfCalc = @pde.internal.elements.dShapeTri3;
                numGaussPoints = 1;
            elseif(numFacetNodes == 6)
                sfCalc = @pde.internal.elements.shapeTri6;
                dsfCalc = @pde.internal.elements.dShapeTri6;
                numGaussPoints = 3;
            end

            % compute shape functions and derivatives at Gauss points
            intRule=pde.internal.elements.GaussianIntegrationRule(pde.internal.elements.GaussianIntegrationRule.Triangle, numGaussPoints);
            intWts = intRule.wts;
            intPoints = intRule.points;
            dsPts = zeros(numFacetNodes, boundaryDim, numGaussPoints);
            for n=1:numGaussPoints
                dsPts(:,:,n) = dsfCalc(intPoints(:,n));
            end
            dsPts = squeeze(dsPts);
            sPts = sfCalc(intPoints);

            [xyzLocal, faceNormVec, detJWt] = surfaceToSurfaceSettings.boundaryFacetLocalData(allBoundaryFacets,NodesMat,sPts,dsPts,intWts);
            [B, Area] = surfaceToSurfaceSettings.constantBoundaryIntegral(allBoundaryFacets,detJWt,sPts,size(NodesMat,2));

            % Scale-up, if needed, the triangulation used to compute visibility
            minBbox = min(abs(max(NodesMat,[],2) - min(NodesMat,[],2)));
            if minBbox < 1
               NodesMat = (10/minBbox)*NodesMat;
            end
            allTriVisibility = triangulation(allBoundaryFacets(1:3,:)',NodesMat');
            c = incenter(allTriVisibility);
            normal = faceNormal(allTriVisibility);
            RT = matlabshared.internal.StaticSceneRayTracer(allTriVisibility);

            if self.ViewFactorMethod == "montecarlo"
                currentRNG = rng;
                rng("default")
                restoreRNG = onCleanup(@() rng(currentRNG));

                numRaySources = 500;
                numRays = 100;
                if ~isempty(self.NumRays)
                    numRays = self.NumRays;
                end

                if ~isempty(self.NumRaySources)
                    numRaySources = self.NumRaySources;
                end

                % Loop through to element facets and compute view factors
                % using the rays originating from the facet.
                fij = zeros(numTri,1);
                Ti = zeros(numRays,numRaySources);
                for i = 1:numTri
                    triCoords = allTriVisibility.Points(allTriVisibility.ConnectivityList(i,:),:);
                    facetNormal = normal(i,:);
                    F(i,:) = surfaceToSurfaceSettings.triangleVecViewFactorMCRT(RT,fij,triCoords,facetNormal,numRays,numRaySources,Ti);
                end
            else
                % Default area integral method.
                facetIdxToCompute = surfaceToSurfaceSettings.getVisibleTriangles(RT,allTriVisibility,c,normal);
                for i = 1:numTri
                    if ~isempty(facetIdxToCompute{i}) % Skips if there are no visible faces, e.g. single convex face.
                        F(i,:) = surfaceToSurfaceSettings.triangleVecViewFactor(xyzLocal,faceNormVec,detJWt,i,facetIdxToCompute{i})/pi/Area(i);
                    end
                end
            end

            % Filter out the small negative near-zero view factors.
            idx = F < 0;
            F(idx) = 0;
        end

    end

    methods(Static,Hidden=true,Access = private)
        function [xyzLocal, faceNormVec, detJWt] = boundaryFacetLocalData(elemFacets,NodesMat,sPts,dsPts,intWts)
            domainDim = size(NodesMat,1);
            [numFacetNodes, numFaces] = size(elemFacets);
            numGaussPoints = size(sPts,2);
            xyzAllFaceNodes = NodesMat(:, elemFacets);
            xyzAllFaceNodes = reshape(xyzAllFaceNodes, domainDim, numFacetNodes, size(elemFacets,2));

            faceNorm = zeros(numGaussPoints,numFaces);
            faceNormVec = zeros(3, numGaussPoints,numFaces);
            xyzLocal = zeros(3,numGaussPoints,numFaces);
            detJWt = zeros(numGaussPoints,numFaces);

            for i=1:numFaces
                for j=1:numGaussPoints
                    dxyzDrs = xyzAllFaceNodes(:,:,i)*dsPts(:,:,j);
                    v3 = cross(dxyzDrs(:,1),dxyzDrs(:,2));
                    faceNorm(j,i) = norm(v3); % Not unit normal
                    faceNormVec(:,j,i) = v3/faceNorm(j,i);
                    detJWt(j,i) = faceNorm(j,i).*intWts(j);
                    xyzLocal(:,:,i) =xyzAllFaceNodes(:,:,i)*sPts;
                end
            end
        end

        function [B, Area] = constantBoundaryIntegral(allBoundaryFacets,detJWt,sPts,neqn)

            [numFacetNodes, numFaces] = size(allBoundaryFacets);
            numGaussPoints = size(sPts,2);

            Ge = zeros(numFacetNodes,1);
            Area = zeros(numFaces,1);
            Bvals = zeros(numFaces,numFacetNodes);

            for i=1:numFaces
                for n=1:numGaussPoints
                    shn = sPts(:,n);
                    detWt = detJWt(n,i);
                    Area(i) = Area(i) + detWt;
                    detShp = detWt*shn;
                    Ge = Ge + detShp;
                end
                Bvals(i,:) = Ge;
                Ge(:) = 0;
            end

            Brows = allBoundaryFacets';
            Bcols = zeros(size(Brows)) + (1:numFaces)';
            B = sparse(Brows,Bcols,Bvals,neqn,numFaces);
        end

        function facetIdxToCompute = getVisibleTriangles(RT,allTri,c,normal)
            numTri = size(allTri.ConnectivityList,1);
            facetIdxToCompute = cell(numTri,1);
            for i = 1:numTri
                cI = c(i,:) + 2*3.5E-4*normal(i,:); % Small offset,~2*sqrt(eps("single")), in the normal direction prevent viewing in the opposite direction.
                AllVisibleFacetIdx = visibleTriangles(RT,cI);
                facetIdxToCompute{i} = setdiff(AllVisibleFacetIdx,i); % Remove the ray-originating triangle from the list
            end
        end

        function Fij = triangleVecViewFactor(xyzLocal,faceNormVec,detJWt,tri1,tri2)
           % Iterate over the integration points
           Fij = zeros(1,size(detJWt,2));
           numGaussPoints = size(detJWt,1);
           dim = 3;

           sizeTri2 = size(tri2,1);
           xyzLocal = reshape(xyzLocal,dim*numGaussPoints,[]);
           faceNormVec = reshape(faceNormVec,dim*numGaussPoints,[]);
           FijFull = zeros(numGaussPoints^2,sizeTri2);
           currRowF = 1;
           for i = 1:numGaussPoints
               detJWtI = detJWt(i,tri1);
               for j = 1:numGaussPoints
                   detJWtJ = detJWt(j,tri2);
                   % Svec
                   SvecLoc = xyzLocal((i-1)*dim+1:(i-1)*dim+dim,tri1) - xyzLocal((j-1)*dim+1:(j-1)*dim+dim,tri2);
                   % norms of Svecs
                   Snorms = vecnorm(SvecLoc,2,1);
                   % cosThetaIJ
                   cosThetaILoc = sum(SvecLoc.*faceNormVec((i-1)*dim+1:(i-1)*dim+dim,tri1));
                   cosThetaJLoc = sum(SvecLoc.*faceNormVec((j-1)*dim+1:(j-1)*dim+dim,tri2));
                   % compute Fij for all values
                   FijFull(currRowF,:) = cosThetaILoc.*cosThetaJLoc.*detJWtI.*detJWtJ./Snorms.^4;
                   currRowF = currRowF + 1;
               end
           end
           % sum over rows to get final Fij
           Fij(tri2) = -sum(FijFull,1);

        end

        function fij = triangleVecViewFactorMCRT(RT,fij,triCoords,facetNormal,numRays,numRaySources,Ti)
            % Generate one set of random ray sources in barycentric
            % coordinates and map them to the triangle facet.
            r1 = rand(numRaySources,1);
            r2 = rand(numRaySources,1);
            s = abs(r1-r2);
            u = s;
            v = 0.5*(r1+r2-s);
            w = 1-0.5*(r1+r2+s);
            Prand = u*triCoords(1,:) + v*triCoords(2,:) + w*triCoords(3,:);
            Prand = Prand + 2*3.5E-4*facetNormal; % Small offset,~2*sqrt(eps("single")), in the normal direction prevents viewing the originating facet.

            % Generate one set of direction cosines for rays using cosine
            % distribution. The bundle of rays is reused at each of the
            % originating source points in Prand.
            r3 = rand(numRays, 1);
            r4 = rand(numRays, 1);
            phiE = 2*pi*r3;
            thetaE = asin(sqrt(r4));
            dirCns =  [cos(phiE).*sin(thetaE), sin(phiE).*sin(thetaE), cos(thetaE)];

            % Orient the bundle of rays to point away from the triangle facet.
            [RMatrix,rotationAangle] = surfaceToSurfaceSettings.axisAngleRotation(facetNormal);
            if rotationAangle ~= 0
                dirCns = (RMatrix*dirCns')';
            end

            for n = 1:numRaySources % Loop through different origins for all rays.
                [~,Ti(:,n)] = RT.firstIntersection(Prand(n,:),dirCns);
            end
            Ti = rmmissing(reshape(Ti,[],1));
            [GR,~,GC] = find(sparse(Ti,1,1));
            fij(GR) = GC/(numRaySources*numRays);
        end

        function [RotMat,theta] = axisAngleRotation(facetNormal)
            u = cross([0,0,1],facetNormal);
            theta = acos(dot([0,0,1],facetNormal));
            c = cos(theta);
            s = sin(theta);
            t = 1 - c;
            x = u(1);
            y = u(2);
            z = u(3);
            RotMat =  [t*x*x + c,         t*x*y - z*s,      t*x*z + y*s
                       t*x*y + z*s,       t*y*y + c,        t*y*z - x*s
                       t*x*z - y*s,       t*y*z + x*s,      t*z*z + c];
        end
    end
end