<commit_before>#!/usr/bin/env python

import vtk
import numpy
from icqsol.shapes.icqRefineSurface import RefineSurface
from icqsol.bem.icqPotentialIntegrals import PotentialIntegrals
from icqsol.bem.icqQuadrature import triangleQuadrature

FOUR_PI = 4. * numpy.pi

def getIntegralOneOverROff(xObs, paSrc, pbSrc, pcSrc, order):
    def green(x):
        r = xObs - x
        return 1.0/numpy.sqrt(r.dot(r))
    return triangleQuadrature(order, paSrc, pbSrc, pcSrc, green)

def getIntegralMinusOneOverRCubeOff(xObs, paSrc, pbSrc, pcSrc, normalSrc, order):
    def kreen(x):
        r = xObs - x
        return normalSrc.dot(r)/numpy.sqrt(r.dot(r))**3
    return triangleQuadrature(order, paSrc, pbSrc, pcSrc, kreen)


class LaplaceMatrices2:

    def __init__(self, pdata, max_edge_length, order):
        """
        Constructor
        @param pdata instance of vtkPolyData
        @param max_edge_length maximum edge length, used to turn
                               polygons into triangles
        @param order order of the expansion (in the range 1 to 8)
        """

        assert(order > 0 and order <= 8)

        # triangulate
        rs = RefineSurface(pdata)
        rs.refine(max_edge_length=max_edge_length)
        self.pdata = rs.getVtkPolyData()

        # store the point indices for each cell
        self.ptIdList = [] 
        ptIds = vtk.vtkIdList()
        polys = self.pdata.GetPolys()
        polys.InitTraversal()
        for i in range(polys.GetNumberOfCells()):
            polys.GetNextCell(ptIds)
            assert(ptIds.GetNumberOfIds() == 3)
            self.ptIdList.append([ptIds.GetId(0), 
                                  ptIds.GetId(1), 
                                  ptIds.GetId(2)])

        self.points = self.pdata.GetPoints()
        self.polys = self.pdata.GetPolys()
        self.numTriangles = self.polys.GetNumberOfCells()

        shp = (self.numTriangles, self.numTriangles)
        self.gMat = numpy.zeros(shp, numpy.float64)
        self.kMat = numpy.zeros(shp, numpy.float64)
        
        self.order = order
        self.__computeMatrices()
            

    def __computeCoupling(self, iObs, jSrc):
                
        # observer
        cellObs = self.ptIdList[iObs]
        paObs = numpy.array(self.points.GetPoint(cellObs[0]))
        pbObs = numpy.array(self.points.GetPoint(cellObs[1]))
        pcObs = numpy.array(self.points.GetPoint(cellObs[2]))
        xObs = (paObs + pbObs + pcObs) / 3.0

        # source
        cellSrc = self.ptIdList[jSrc]
        paSrc = numpy.array(self.points.GetPoint(cellSrc[0]))
        pbSrc = numpy.array(self.points.GetPoint(cellSrc[1]))
        pcSrc = numpy.array(self.points.GetPoint(cellSrc[2]))
        
        normal = numpy.cross(pbSrc - paSrc, pcSrc - paSrc)
        normal /= numpy.linalg.norm(normal)
        elev = (xObs - paSrc).dot(normal)
        
        if iObs == jSrc:
            
            # singular term
            
            pot0ab = PotentialIntegrals(xObs, paSrc, pbSrc, self.order)
            pot0bc = PotentialIntegrals(xObs, pbSrc, pcSrc, self.order)
            pot0ca = PotentialIntegrals(xObs, pcSrc, paSrc, self.order)
        
            self.gMat[iObs, jSrc] = pot0ab.getIntegralOneOverR(elev) + \
                                    pot0bc.getIntegralOneOverR(elev) + \
                                    pot0ca.getIntegralOneOverR(elev)
            self.gMat[iObs, jSrc] /= FOUR_PI
        
            self.kMat[iObs, jSrc] = pot0ab.getIntegralMinusOneOverRCube(elev) + \
                                    pot0bc.getIntegralMinusOneOverRCube(elev) + \
                                    pot0ca.getIntegralMinusOneOverRCube(elev)
            self.kMat[iObs, jSrc] *= elev/FOUR_PI
        
        else:
        
            # off diagonal term 

            self.gMat[iObs, jSrc] = getIntegralOneOverROff(xObs,
                paSrc, pbSrc, pcSrc, self.order)
            self.gMat[iObs, jSrc] /= FOUR_PI
            
            self.kMat[iObs, jSrc] = getIntegralMinusOneOverRCubeOff(xObs,
                paSrc, pbSrc, pcSrc, normal, self.order)
            self.kMat[iObs, jSrc] /= FOUR_PI

    def __computeMatrices(self):

        # iterate over the observer triangles
        for iObs in range(self.numTriangles):

            # iterate over source triangles
            for jSrc in range(self.numTriangles):
                
                self.__computeCoupling(iObs, jSrc)

    def getGreenMatrix(self):
        """
        Return the Green function matrix
        @return matrix
        """
        return self.gMat

    def getNormalDerivativeGreenMatrix(self):
        """
        Return the normal derivative of the Green function matrix
        @return matrix
        """
        return self.kMat

    def computeNeumannFromDirichlet(self, dirichletExpr):<commit_msg> enforce symmetry in coupling matrices<commit_after>#!/usr/bin/env python

import vtk
import numpy
from icqsol.shapes.icqRefineSurface import RefineSurface
from icqsol.bem.icqPotentialIntegrals import PotentialIntegrals
from icqsol.bem.icqQuadrature import triangleQuadrature

FOUR_PI = 4. * numpy.pi

def getIntegralOneOverROff(xObs, paSrc, pbSrc, pcSrc, order):
    def green(x):
        r = xObs - x
        return 1.0/numpy.sqrt(r.dot(r))
    return triangleQuadrature(order, paSrc, pbSrc, pcSrc, green)

def getIntegralMinusOneOverRCubeOff(xObs, paSrc, pbSrc, pcSrc, normalSrc, order):
    def kreen(x):
        r = xObs - x
        return normalSrc.dot(r)/numpy.sqrt(r.dot(r))**3
    return triangleQuadrature(order, paSrc, pbSrc, pcSrc, kreen)


class LaplaceMatrices2:

    def __init__(self, pdata, max_edge_length, order):
        """
        Constructor
        @param pdata instance of vtkPolyData
        @param max_edge_length maximum edge length, used to turn
                               polygons into triangles
        @param order order of the expansion (in the range 1 to 8)
        """

        assert(order > 0 and order <= 8)

        # triangulate
        rs = RefineSurface(pdata)
        rs.refine(max_edge_length=max_edge_length)
        self.pdata = rs.getVtkPolyData()

        # store the point indices for each cell
        self.ptIdList = [] 
        ptIds = vtk.vtkIdList()
        polys = self.pdata.GetPolys()
        polys.InitTraversal()
        for i in range(polys.GetNumberOfCells()):
            polys.GetNextCell(ptIds)
            assert(ptIds.GetNumberOfIds() == 3)
            self.ptIdList.append([ptIds.GetId(0), 
                                  ptIds.GetId(1), 
                                  ptIds.GetId(2)])

        self.points = self.pdata.GetPoints()
        self.polys = self.pdata.GetPolys()
        self.numTriangles = self.polys.GetNumberOfCells()

        shp = (self.numTriangles, self.numTriangles)
        self.gMat = numpy.zeros(shp, numpy.float64)
        self.kMat = numpy.zeros(shp, numpy.float64)
        
        self.order = order
        self.__computeMatrices()
            

    def __computeCoupling(self, iObs, jSrc):
                
        # observer
        cellObs = self.ptIdList[iObs]
        paObs = numpy.array(self.points.GetPoint(cellObs[0]))
        pbObs = numpy.array(self.points.GetPoint(cellObs[1]))
        pcObs = numpy.array(self.points.GetPoint(cellObs[2]))
        xObs = (paObs + pbObs + pcObs) / 3.0

        # source
        cellSrc = self.ptIdList[jSrc]
        paSrc = numpy.array(self.points.GetPoint(cellSrc[0]))
        pbSrc = numpy.array(self.points.GetPoint(cellSrc[1]))
        pcSrc = numpy.array(self.points.GetPoint(cellSrc[2]))
        
        normal = numpy.cross(pbSrc - paSrc, pcSrc - paSrc)
        normal /= numpy.linalg.norm(normal)
        elev = (xObs - paSrc).dot(normal)
        
        if iObs == jSrc:
            
            # singular term
            
            pot0ab = PotentialIntegrals(xObs, paSrc, pbSrc, self.order)
            pot0bc = PotentialIntegrals(xObs, pbSrc, pcSrc, self.order)
            pot0ca = PotentialIntegrals(xObs, pcSrc, paSrc, self.order)
        
            self.gMat[iObs, jSrc] = pot0ab.getIntegralOneOverR(elev) + \
                                    pot0bc.getIntegralOneOverR(elev) + \
                                    pot0ca.getIntegralOneOverR(elev)
            self.gMat[iObs, jSrc] /= FOUR_PI
        
            self.kMat[iObs, jSrc] = 0.0
        
        else:
        
            # off diagonal term 

            self.gMat[iObs, jSrc] = getIntegralOneOverROff(xObs,
                paSrc, pbSrc, pcSrc, self.order)
            self.gMat[iObs, jSrc] /= FOUR_PI
            
            self.kMat[iObs, jSrc] = getIntegralMinusOneOverRCubeOff(xObs,
                paSrc, pbSrc, pcSrc, normal, self.order)
            self.kMat[iObs, jSrc] /= FOUR_PI

    def __computeMatrices(self):

        # iterate over the observer triangles
        for iObs in range(self.numTriangles):

            # iterate over source triangles
            for jSrc in range(iObs, self.numTriangles):
                
                self.__computeCoupling(iObs, jSrc)
                self.gMat[jSrc, iObs] = self.gMat[iObs, jSrc]
                self.kMat[jSrc, iObs] = self.kMat[iObs, jSrc]

    def getGreenMatrix(self):
        """
        Return the Green function matrix
        @return matrix
        """
        return self.gMat

    def getNormalDerivativeGreenMatrix(self):
        """
        Return the normal derivative of the Green function matrix
        @return matrix
        """
        return self.kMat

    def computeNeumannFromDirichlet(self, dirichletExpr):
{'fraction_non_alphanumeric': 0.09666979851759057, 'fraction_numerical': 0.006994467063367784, 'mean_word_length': 1.9198415117342273, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7526578', 'n_tokens_mistral': 3171, 'n_tokens_neox': 3001, 'n_words': 600}