<commit_before>"""
A class hierarchy for numerical integration on reference domains in 1,2, and 3D.
"""
from EGeometry import *
from Profiling import logEvent
from math import *

class Q_base:
    """
    The base class for quadrature methods.
    """
    def __init__(self,order=1):
        self.points=()
        self.weights=()
        self.pointsAll=(())
        self.weightsAll=(())
        self.order=order
    def setOrder(self,k):
        self.order=k
        if self.order > len(self.pointsAll):
            logEvent("WARNING Q_base requested order=%d > max allowed=%d setting order to max" % (self.order,
                                                                                               len(self.pointsAll)))
            self.order= len(self.pointsAll)
        self.points=self.pointsAll[self.order-1]
        self.weights=self.weightsAll[self.order-1]

class GaussPoint(Q_base):
    """
    A dummy class for integrating the boundary of the unit interval (i.e. a point).
    """
    def __init__(self,order=1):
        Q_base.__init__(self,order)
        self.pointsAll=(
            (EVec(0.0),),)
        self.weightsAll=(
            (1.0,),)
        self.setOrder(order)
    def setOrder(self,k):
        self.order = k
        self.points = self.pointsAll[0]
        self.weights = self.weightsAll[0]

LobattoPoint = GaussPoint

class GaussEdge(Q_base):
    """
    Gaussian Quadrature on the unit interval.
    """
    def __init__(self,order=1):
        #mwf for convenience, see Guern and Ermond
        a1 = 0.5*sqrt((15.0+2.0*sqrt(30))/35.0)
        a2 = 0.5*sqrt((15.0-2.0*sqrt(30))/35.0)
        w1 = 0.25 - sqrt(5./6.)/12.
        w2 = 0.25 + sqrt(5./6.)/12.
        Q_base.__init__(self,order)
        self.pointsAll=(
            (EVec(0.5),),
            (EVec((sqrt(3.0)-1.0)/(2.0*sqrt(3.0))),
             EVec((sqrt(3.0)+1.0)/(2.0*sqrt(3.0)))),
            (EVec((sqrt(5.0) - sqrt(3.0))/(2.0*sqrt(5))),
             EVec(0.5),
             EVec((sqrt(5.0) + sqrt(3.0))/(2.0*sqrt(5)))),
            (EVec(0.5+a1),EVec(0.5-a1),
             EVec(0.5+a2),EVec(0.5-a2)),
            (EVec(0.5),
             EVec(0.5*(sqrt(5.0-2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(-sqrt(5.0-2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(sqrt(5.0+2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(-sqrt(5.0+2.0*sqrt(10.0/7.0))/3.0) + 0.5))
            )
        self.weightsAll=(
            (1.0,),
            (0.5,
             0.5),
            (5.0/18.0,
             8.0/18.0,
             5.0/18.0),
            (w1,w1,w2,w2),
            (0.5*(128.0/225.0),
             0.5*(322.0+13.0*sqrt(70.0))/900.0,
             0.5*(322.0+13.0*sqrt(70.0))/900.0,
             0.5*(322.0-13.0*sqrt(70.0))/900.0,
             0.5*(322.0-13.0*sqrt(70.0))/900.0)
            )
        self.setOrder(order)

    def setOrder(self,order,domain=[0.0,1.0]):
        Q_base.setOrder(self,order)

        points = self.points
        weights = self.weights       
        self.points = []
        self.weights = []
        for i in range(order):
           self.points.append(EVec((domain[1]-domain[0])*points[i][0] +domain[0]))
           self.weights.append((domain[1]-domain[0])*weights[i])
        
class LobattoEdge(Q_base):
    """
    Gauss-Lobatto quadrature on the unit interval.
    """
    def __init__(self,order=1):
        order=1
        Q_base.__init__(self,order)
        self.pointsAll=(
            (EVec(0.0),EVec(1.0)),)
        self.weightsAll=(
            (0.5,0.5),)
        self.setOrder(order)
    def setOrder(self,k):
        self.order = 1
        self.points = self.pointsAll[0]
        self.weights = self.weightsAll[0]

class LobattoEdgeAlt(Q_base):
    """
    Gauss-Lobatto quadrature on the [-1:1] interval.
    """
    def __init__(self,order=1):
        Q_base.__init__(self,order)
        a1 = 0.2*sqrt(5.0)
        a2 = sqrt(21)/7.0
        self.pointsAll=(
            (EVec(-1.0),EVec(1.0)),
            (EVec(-1.0),EVec(0.0),EVec(1.0)),
            (EVec(-1.0),EVec(-a1),EVec(a1) ,EVec(1.0)),
            (EVec(-1.0),EVec(-a2),EVec(0.0),EVec(a2),EVec(1.0)),)
        self.weightsAll=(
            (1.0,1.0),
            (1.0/3.0, 4.0/3.0, 1.0/3.0),
            (1.0/6.0, 5.0/6.0,5.0/6.0,1.0/6.0),
            (1.0/10.0,49.0/90.0,32.0/45.0,49.0/90.0,1.0/10.0),)<commit_msg>added higher order rules for GaussLobattoEdge<commit_after>"""
A class hierarchy for numerical integration on reference domains in 1,2, and 3D.
"""
from EGeometry import *
from Profiling import logEvent
from math import *

class Q_base:
    """
    The base class for quadrature methods.
    """
    def __init__(self,order=1):
        self.points=()
        self.weights=()
        self.pointsAll=(())
        self.weightsAll=(())
        self.order=order
    def setOrder(self,k):
        self.order=k
        if self.order > len(self.pointsAll):
            logEvent("WARNING Q_base requested order=%d > max allowed=%d setting order to max" % (self.order,
                                                                                               len(self.pointsAll)))
            self.order= len(self.pointsAll)
        self.points=self.pointsAll[self.order-1]
        self.weights=self.weightsAll[self.order-1]

class GaussPoint(Q_base):
    """
    A dummy class for integrating the boundary of the unit interval (i.e. a point).
    """
    def __init__(self,order=1):
        Q_base.__init__(self,order)
        self.pointsAll=(
            (EVec(0.0),),)
        self.weightsAll=(
            (1.0,),)
        self.setOrder(order)
    def setOrder(self,k):
        self.order = k
        self.points = self.pointsAll[0]
        self.weights = self.weightsAll[0]

LobattoPoint = GaussPoint

class GaussEdge(Q_base):
    """
    Gaussian Quadrature on the unit interval.
    """
    def __init__(self,order=1):
        #mwf for convenience, see Guern and Ermond
        a1 = 0.5*sqrt((15.0+2.0*sqrt(30))/35.0)
        a2 = 0.5*sqrt((15.0-2.0*sqrt(30))/35.0)
        w1 = 0.25 - sqrt(5./6.)/12.
        w2 = 0.25 + sqrt(5./6.)/12.
        Q_base.__init__(self,order)
        self.pointsAll=(
            (EVec(0.5),),
            (EVec((sqrt(3.0)-1.0)/(2.0*sqrt(3.0))),
             EVec((sqrt(3.0)+1.0)/(2.0*sqrt(3.0)))),
            (EVec((sqrt(5.0) - sqrt(3.0))/(2.0*sqrt(5))),
             EVec(0.5),
             EVec((sqrt(5.0) + sqrt(3.0))/(2.0*sqrt(5)))),
            (EVec(0.5+a1),EVec(0.5-a1),
             EVec(0.5+a2),EVec(0.5-a2)),
            (EVec(0.5),
             EVec(0.5*(sqrt(5.0-2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(-sqrt(5.0-2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(sqrt(5.0+2.0*sqrt(10.0/7.0))/3.0) + 0.5),
             EVec(0.5*(-sqrt(5.0+2.0*sqrt(10.0/7.0))/3.0) + 0.5))
            )
        self.weightsAll=(
            (1.0,),
            (0.5,
             0.5),
            (5.0/18.0,
             8.0/18.0,
             5.0/18.0),
            (w1,w1,w2,w2),
            (0.5*(128.0/225.0),
             0.5*(322.0+13.0*sqrt(70.0))/900.0,
             0.5*(322.0+13.0*sqrt(70.0))/900.0,
             0.5*(322.0-13.0*sqrt(70.0))/900.0,
             0.5*(322.0-13.0*sqrt(70.0))/900.0)
            )
        self.setOrder(order)

    def setOrder(self,order,domain=[0.0,1.0]):
        Q_base.setOrder(self,order)

        points = self.points
        weights = self.weights       
        self.points = []
        self.weights = []
        for i in range(order):
           self.points.append(EVec((domain[1]-domain[0])*points[i][0] +domain[0]))
           self.weights.append((domain[1]-domain[0])*weights[i])
        
class LobattoEdge(Q_base):
    """
    Gauss-Lobatto quadrature on the unit interval.
    """
    def __init__(self,order=1):
        Q_base.__init__(self,order)
        a1 = 0.2*sqrt(5.0)
        a2 = sqrt(21.0)/7.0
        self.pointsAll=(
            (EVec(0.0),EVec(1.0)),
            (EVec(0.0),EVec(0.5),EVec(1.0)),
            (EVec(0.0),EVec(0.5*(1.-a1)),EVec(0.5*(1.+a1)),EVec(1.0)),
            (EVec(0.0),EVec(0.5*(1.-a2)),EVec(0.5),EVec(0.5*(1.+a2)),EVec(1.0)),
            )
        self.weightsAll=(
            (0.5,0.5),
            (1.0/6.0, 4.0/6.0, 1.0/6.0),
            (1.0/12.0, 5.0/12.0,5.0/12.0,1.0/12.0),
            (1.0/20.0,49.0/180.0,32.0/90.0,49.0/180.0,1.0/20.0)
            )
        self.setOrder(order)

class LobattoEdgeAlt(Q_base):
    """
    Gauss-Lobatto quadrature on the [-1:1] interval.
    """
    def __init__(self,order=1):
        Q_base.__init__(self,order)
        a1 = 0.2*sqrt(5.0)
        a2 = sqrt(21)/7.0
        self.pointsAll=(
            (EVec(-1.0),EVec(1.0)),
            (EVec(-1.0),EVec(0.0),EVec(1.0)),
            (EVec(-1.0),EVec(-a1),EVec(a1) ,EVec(1.0)),
            (EVec(-1.0),EVec(-a2),EVec(0.0),EVec(a2),EVec(1.0)),)
        self.weightsAll=(
            (1.0,1.0),
            (1.0/3.0, 4.0/3.0, 1.0/3.0),
            (1.0/6.0, 5.0/6.0,5.0/6.0,1.0/6.0),
            (1.0/10.0,49.0/90.0,32.0/45.0,49.0/90.0,1.0/10.0),)
{'fraction_non_alphanumeric': 0.1997787610619469, 'fraction_numerical': 0.0959070796460177, 'mean_word_length': 1.6520973892637136, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2252783', 'n_tokens_mistral': 3837, 'n_tokens_neox': 3611, 'n_words': 462}