<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 gaussPtsAndWeights

FOUR_PI = 4. * numpy.pi

class LaplaceMatrices:

    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 Gauss quadrature scheme
        """

        assert(order > 0 and order <= 5)
        
        self.normalDerivativeName = 'normal_derivative'
        self.potentialName = 'potential'

        # 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 setPotentialName(self, name):
        self.potentialName = name

    def setNormalDerivativeName(self, name):
        self.normalDerivativeName = name

    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 
            pb2 = pbSrc - paSrc
            pc2 = pcSrc - paSrc
            areaSrc = numpy.linalg.norm(numpy.cross(pb2, pc2))
            normDistance = numpy.linalg.norm((paSrc + pbSrc + pcSrc)/3. - xObs) \
                         / numpy.sqrt(areaSrc)
            offDiagonalOrder = min(8, max(1, int(8 * 2 / normDistance)))

            # Gauss quadrature weights
            gpws = gaussPtsAndWeights[offDiagonalOrder]

            # number of Gauss points
            npts = gpws.shape[1]

            # triangle positions and weights
            xsis, etas, weights = gpws[0, :], gpws[1, :], gpws[2, :]
            
            self.gMat[iObs, jSrc] = 0.0
            self.kMat[iObs, jSrc] = 0.0
            for k in range(npts):
                dr = xObs - paSrc - pb2*xsis[k] - pc2*etas[k]
                drNorm = numpy.sqrt(dr.dot(dr))
                self.gMat[iObs, jSrc] += weights[k] / drNorm
                self.kMat[iObs, jSrc] += weights[k] * normal.dot(dr) / drNorm**3

            self.gMat[iObs, jSrc] *= 0.5 * areaSrc / FOUR_PI
            self.kMat[iObs, jSrc] *= 0.5 * areaSrc / FOUR_PI
            
    def getVtkPolyData(self):
        """
        Get the (modified) vtkPolyData object
        @return object
        """
        return self.pdata

    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):
        """
        Get the Neumann boundary values from the Dirichlet boundary conditions
        @param dirichletExpr expression for the potential values
        @return response
        """
        from math import pi, sin, cos, log, exp, sqrt
        
        kMat = self.getNormalDerivativeGreenMatrix()
        n = kMat.shape[0]
        
        v = numpy.zeros((n,), numpy.float64)
        
        # subtract residue
        for i in range(n):
            kMat[i, i] -= 0.5
        
        # evaluate potential on each triangle center
        ptIds = vtk.vtkIdList()
        polys = self.pdata.GetPolys()
        points = self.pdata.GetPoints()
        polys.InitTraversal()
        for i in range(n):
            
            polys.GetNextCell(ptIds)
            ia, ib, ic = ptIds.GetId(0), ptIds.GetId(1), ptIds.GetId(2)
            pa = numpy.array(points.GetPoint(ia))
            pb = numpy.array(points.GetPoint(ib))
            pc = numpy.array(points.GetPoint(ic))
            pMid = (pa + pb + pc)/3.
            x, y, z = pMid
            
            # set the Dirichlet value
            v[i] = eval(dirichletExpr)
        
        gMat = self.getGreenMatrix()
        
        # solve
        gM1 = numpy.linalg.inv(gMat)
        
        normalDerivative = gM1.dot(kMat).dot(v)
        
        # add field
        cellData = self.pdata.GetCellData()
        
        potentialData = vtk.vtkDoubleArray()
        potentialData.SetNumberOfComponents(1)
        potentialData.SetNumberOfTuples(n)
        potentialData.SetName(self.potentialName)

        normalDerivData = vtk.vtkDoubleArray()
        normalDerivData.SetNumberOfComponents(1)
        normalDerivData.SetNumberOfTuples(n)
        normalDerivData.SetName(self.normalDerivativeName)

        for i in range(n):
            potentialData.SetTuple(i, [v[i],])
            normalDerivData.SetTuple(i, [normalDerivative[i],])

        cellData.AddArray(potentialData)
        cellData.AddArray(normalDerivData)
        
        return normalDerivative 


###############################################################################


def testSingTriangle():

    "Single triangle"

    h = 0.1
    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(3)
    points.SetPoint(0, [1., -1.*h/3., -1.*h/3.])
    points.SetPoint(1, [1., 2.*h/3., -1.*h/3.])
    points.SetPoint(2, [1., -1.*h/3., 2.*h/3.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.Allocate(1, 1)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

def testTwoTrianglesCoplanar():

    "Two triangles"

    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(4)
    points.SetPoint(0, [0., 0., 0.])
    points.SetPoint(1, [1., 0., 0.])
    points.SetPoint(2, [0., 1., 0.])
    points.SetPoint(3, [1., 1., 0.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    
    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 1)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

def testTwoTriangles():

    "Two triangles"

    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(4)
    points.SetPoint(0, [0., 0., 0.])
    points.SetPoint(1, [1., 0., 0.])
    points.SetPoint(2, [0., 1., 0.])
    points.SetPoint(3, [0., 0., 1.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    
    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 3)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

if __name__ == '__main__':
    #testSingleTriangle()
    #testTwoTrianglesCoplanar()
    testTwoTriangles()
<commit_msg> down to 21 secs from 46s initially<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 gaussPtsAndWeights

FOUR_PI = 4. * numpy.pi

class LaplaceMatrices:

    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 Gauss quadrature scheme
        """

        assert(order > 0 and order <= 5)
        
        self.normalDerivativeName = 'normal_derivative'
        self.potentialName = 'potential'

        # 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 setPotentialName(self, name):
        self.potentialName = name

    def setNormalDerivativeName(self, name):
        self.normalDerivativeName = name                
            
    def getVtkPolyData(self):
        """
        Get the (modified) vtkPolyData object
        @return object
        """
        return self.pdata

    def __computeMatrices(self):

        # iterate over the source triangles
        for jSrc in range(self.numTriangles):

            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]))
        
            normalSrc = numpy.cross(pbSrc - paSrc, pcSrc - paSrc)
            normalSrc /= numpy.linalg.norm(normalSrc)

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

                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]))

                # observer is at mid point
                xObs = (paObs + pbObs + pcObs) / 3.0

                elev = (xObs - paSrc).dot(normalSrc)

                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 
                    pb2Src = pbSrc - paSrc
                    pc2Src = pcSrc - paSrc
                    areaSrc = numpy.linalg.norm(numpy.cross(pb2Src, pc2Src))
            
                    normDistance = numpy.linalg.norm((paSrc + pbSrc + pcSrc)/3. - xObs) \
                         / numpy.sqrt(areaSrc)
                    offDiagonalOrder = min(8, max(1, int(8 * 2 / normDistance)))

                    # Gauss quadrature weights
                    gpws = gaussPtsAndWeights[offDiagonalOrder]

                    # number of Gauss points
                    npts = gpws.shape[1]

                    # triangle positions and weights
                    xsis, etas, weights = gpws[0, :], gpws[1, :], gpws[2, :]
            
                    self.gMat[iObs, jSrc] = 0.0
                    self.kMat[iObs, jSrc] = 0.0
                    for k in range(npts):
                        dr = xObs - paSrc - pb2Src*xsis[k] - pc2Src*etas[k]
                        drNorm = numpy.sqrt(dr.dot(dr))
                        self.gMat[iObs, jSrc] += weights[k] / drNorm
                        self.kMat[iObs, jSrc] += weights[k] * normalSrc.dot(dr) / drNorm**3

                    self.gMat[iObs, jSrc] *= 0.5 * areaSrc / FOUR_PI
                    self.kMat[iObs, jSrc] *= 0.5 * areaSrc / FOUR_PI


    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):
        """
        Get the Neumann boundary values from the Dirichlet boundary conditions
        @param dirichletExpr expression for the potential values
        @return response
        """
        from math import pi, sin, cos, log, exp, sqrt
        
        kMat = self.getNormalDerivativeGreenMatrix()
        n = kMat.shape[0]
        
        v = numpy.zeros((n,), numpy.float64)
        
        # subtract residue
        for i in range(n):
            kMat[i, i] -= 0.5
        
        # evaluate potential on each triangle center
        ptIds = vtk.vtkIdList()
        polys = self.pdata.GetPolys()
        points = self.pdata.GetPoints()
        polys.InitTraversal()
        for i in range(n):
            
            polys.GetNextCell(ptIds)
            ia, ib, ic = ptIds.GetId(0), ptIds.GetId(1), ptIds.GetId(2)
            pa = numpy.array(points.GetPoint(ia))
            pb = numpy.array(points.GetPoint(ib))
            pc = numpy.array(points.GetPoint(ic))
            pMid = (pa + pb + pc)/3.
            x, y, z = pMid
            
            # set the Dirichlet value
            v[i] = eval(dirichletExpr)
        
        gMat = self.getGreenMatrix()
        
        # solve
        gM1 = numpy.linalg.inv(gMat)
        
        normalDerivative = gM1.dot(kMat).dot(v)
        
        # add field
        cellData = self.pdata.GetCellData()
        
        potentialData = vtk.vtkDoubleArray()
        potentialData.SetNumberOfComponents(1)
        potentialData.SetNumberOfTuples(n)
        potentialData.SetName(self.potentialName)

        normalDerivData = vtk.vtkDoubleArray()
        normalDerivData.SetNumberOfComponents(1)
        normalDerivData.SetNumberOfTuples(n)
        normalDerivData.SetName(self.normalDerivativeName)

        for i in range(n):
            potentialData.SetTuple(i, [v[i],])
            normalDerivData.SetTuple(i, [normalDerivative[i],])

        cellData.AddArray(potentialData)
        cellData.AddArray(normalDerivData)
        
        return normalDerivative 


###############################################################################


def testSingTriangle():

    "Single triangle"

    h = 0.1
    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(3)
    points.SetPoint(0, [1., -1.*h/3., -1.*h/3.])
    points.SetPoint(1, [1., 2.*h/3., -1.*h/3.])
    points.SetPoint(2, [1., -1.*h/3., 2.*h/3.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.Allocate(1, 1)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

def testTwoTrianglesCoplanar():

    "Two triangles"

    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(4)
    points.SetPoint(0, [0., 0., 0.])
    points.SetPoint(1, [1., 0., 0.])
    points.SetPoint(2, [0., 1., 0.])
    points.SetPoint(3, [1., 1., 0.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    
    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 1)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

def testTwoTriangles():

    "Two triangles"

    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(4)
    points.SetPoint(0, [0., 0., 0.])
    points.SetPoint(1, [1., 0., 0.])
    points.SetPoint(2, [0., 1., 0.])
    points.SetPoint(3, [0., 0., 1.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    
    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 3)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = LaplaceMatrices(pdata,
                               max_edge_length=1000.,
                               order=order)
        print 'order = ', order
        print 'g matrix: ', lslm.getGreenMatrix()
        print 'k matrix: ', lslm.getNormalDerivativeGreenMatrix()

if __name__ == '__main__':
    #testSingleTriangle()
    #testTwoTrianglesCoplanar()
    testTwoTriangles()

{'fraction_non_alphanumeric': 0.11521194650536364, 'fraction_numerical': 0.016823401540569917, 'mean_word_length': 1.8264992653933485, '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': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2817925', 'n_tokens_mistral': 7191, 'n_tokens_neox': 6775, 'n_words': 1342}