classdef GaussianIntegrationRule
%GAUSSIANINTEGRATIONRULE Implement Gaussian integration rules in 1D, 2D, and 3D
% Gaussian integration rules for selected orders are implemented for curves, 2D triangles,
% and 3D tetrahedrons.
% This undocumented function may be changed or removed in a future release.

%       Copyright 2015-2016 The MathWorks, Inc.
  
  properties
    points, wts;
  end
  
  properties(Constant)
      Curve=1, Triangle=2, Tetrahedron=3;
  end


  methods
    function obj = GaussianIntegrationRule(regionType, numPoints)
	  import pde.internal.elements.GaussianIntegrationRule;
      
      if(regionType == GaussianIntegrationRule.Triangle)
        obj=setTriPoints(obj, numPoints);
      elseif ( regionType == GaussianIntegrationRule.Curve)
        obj=setCurvePoints(obj, numPoints);
      else
        error(message('pde:pde3DBCImpl:GaussIntRegionType', regionType));
      end
      
    end
  end
  
  methods(Access=private)
    
    
    function self=setTriPoints(self, numPoints)
      jacFac = .5;
      if(numPoints==1)
        self.points=[1/3,1/3]';
        self.wts=jacFac;
      elseif(numPoints==3)
        self.points=[2/3 1/6; 1/6 2/3; 1/6 1/6]';
        self.wts=[1/3, 1/3, 1/3]*jacFac;
      else
        error(message('pde:pde3DBCImpl:GaussIntNumPts2D', numPoints));
      end
    end
    
    function self=setCurvePoints(self, numPoints)
      jacFac = .5; % le/2
      if(numPoints==2)
        self.points=[-1/sqrt(3); 1/sqrt(3)]';
        self.wts=[1, 1]*jacFac;
      elseif(numPoints==3)
        self.points=[-sqrt(0.6); 0; sqrt(0.6)]';
        self.wts=[ 5/9,    8/9,  5/9  ]*jacFac;
      else
        error(message('pde:pde3DBCImpl:GaussIntNumPts1D', numPoints));
      end
    end
    
    
  end
  
end

