I finally got my FEM EM code working. I profiled it and this function uses up a big hunk of time. It performs gaussian integration over a triangle. I am trying to figure out how to slice the arrays so as to push it down into the C level. Does anyone have any ideas? Thanks, Rob. ps. it looks to be intractible to me. Maybe I need to look at writing a C extension. I've never done that before. ##*************************************************************************** ##Prototype: void ComputeGaussQuadPoint(int QuadPoint, int *a.TrngleNode, ## double *SrcPointCol) ##Description: To compute the coordinates of 7-point Gauss nodes of ## a triangular patch. ##Input value: ## int QuadPoint --- node index, it can be from 0 to 6. ## int *a.TrngleNode --- the three nodes of a tringular patch. ## double *SrcPointCol --- where to store the results ##Return value: none ##Global value used: a.NodeCord, a.Qpnt ##Global value modified: none ##Subroutines called: none ##Note: Not very sure. **************************************************************************** def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a): SrcPointCol=zeros((3),Float) SrcPointCol[0] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0] SrcPointCol[1] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1] SrcPointCol[2] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2] return SrcPointCol -- The Numeric Python EM Project www.members.home.net/europax
On Sun, 26 Aug 2001, Rob wrote:
I finally got my FEM EM code working. I profiled it and this function uses up a big hunk of time. It performs gaussian integration over a triangle. I am trying to figure out how to slice the arrays so as to push it down into the C level. Does anyone have any ideas? Thanks, Rob.
ps. it looks to be intractible to me. Maybe I need to look at writing a C extension. I've never done that before. [...]
I'm not a good enough ufunc / array function hacker to come up with anything, but I doubt it, and I wonder how much time it would save anyway, given the size of the arrays involved. From the look of your comments before the function, it looks like a) you're a C programmer, not quite comfortable with Python yet, or b) you wrote it in C first, then moved it into Python. If the latter, you might want to try wrapping that C function with SWIG, though to be honest I'm not sure the overhead of learning to use SWIG is any less than for writing a C extension manually (but less error-prone I expect, and less work if you have a lot of C to wrap). I think the Scipy project has some typemaps for passing Numeric arrays; if not, I've seen some collections of SWIG typemaps for Numeric somewhere... BTW, the best place to put those comments is in a docstring. Here is a slightly more Pythonically-formatted version (I must be bored today): def ComputeGaussQuadPoint(QuadPoint, a): """Return coordinates of 7-point Gauss nodes of a triangular patch. QuadPoint: node index, from 0 to 6 a: triangular patch? """ SrcPointCol=zeros((3),Float) tn = a.TrngleNode # the three nodes of the triangular patch SrcPointCol[0] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],0] + a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],0] + a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],0]) SrcPointCol[1] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],1] + a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],1] + a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],1]) SrcPointCol[2] = (a.Qpnt[QuadPoint,0]*a.NodeCord[tn[0],2] + a.Qpnt[QuadPoint,1]*a.NodeCord[tn[1],2] + a.Qpnt[QuadPoint,2]*a.NodeCord[tn[2],2]) return SrcPointCol John
def cgqp(QuadPoint, TrngleNode, a): s = a.Qpnt[Quadpoint,:] c = Numeric.take(a.NodeCord, TrngleNode) return Numeric.add.reduce(s * c, axis=1) This may or may not be right. The data structures I would have to set up to test it are too much for Sunday morning. -----Original Message----- From: numpy-discussion-admin@lists.sourceforge.net [mailto:numpy-discussion-admin@lists.sourceforge.net] On Behalf Of Rob Sent: Sunday, August 26, 2001 8:39 AM To: numpy-discussion@lists.sourceforge.net. Subject: [Numpy-discussion] Can this function by Numpy-ized? I finally got my FEM EM code working. I profiled it and this function uses up a big hunk of time. It performs gaussian integration over a triangle. I am trying to figure out how to slice the arrays so as to push it down into the C level. Does anyone have any ideas? Thanks, Rob. ps. it looks to be intractible to me. Maybe I need to look at writing a C extension. I've never done that before. ##********************************************************************** ***** ##Prototype: void ComputeGaussQuadPoint(int QuadPoint, int *a.TrngleNode, ## double *SrcPointCol) ##Description: To compute the coordinates of 7-point Gauss nodes of ## a triangular patch. ##Input value: ## int QuadPoint --- node index, it can be from 0 to 6. ## int *a.TrngleNode --- the three nodes of a tringular patch. ## double *SrcPointCol --- where to store the results ##Return value: none ##Global value used: a.NodeCord, a.Qpnt ##Global value modified: none ##Subroutines called: none ##Note: Not very sure. ************************************************************************ **** def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a): SrcPointCol=zeros((3),Float) SrcPointCol[0] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0] SrcPointCol[1] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1] SrcPointCol[2] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2] return SrcPointCol -- The Numeric Python EM Project www.members.home.net/europax _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net http://lists.sourceforge.net/lists/listinfo/numpy-discussion
Thanks Paul, I will try it out. That at least gives me some direction. I've been pouring over the Numpy manual, but there are so many different functions! I would like to try to write an extension that includes some of these FEM routines, but thats for later. BTW, at work I have to use Windows NT, but I've grown to love the calldll/DynWin extension. You can even interact with the windows kernel with that one. Rob. "Paul F. Dubois" wrote:
def cgqp(QuadPoint, TrngleNode, a): s = a.Qpnt[Quadpoint,:] c = Numeric.take(a.NodeCord, TrngleNode) return Numeric.add.reduce(s * c, axis=1)
This may or may not be right. The data structures I would have to set up to test it are too much for Sunday morning.
-----Original Message----- From: numpy-discussion-admin@lists.sourceforge.net [mailto:numpy-discussion-admin@lists.sourceforge.net] On Behalf Of Rob Sent: Sunday, August 26, 2001 8:39 AM To: numpy-discussion@lists.sourceforge.net. Subject: [Numpy-discussion] Can this function by Numpy-ized?
I finally got my FEM EM code working. I profiled it and this function uses up a big hunk of time. It performs gaussian integration over a triangle. I am trying to figure out how to slice the arrays so as to push it down into the C level. Does anyone have any ideas? Thanks, Rob.
ps. it looks to be intractible to me. Maybe I need to look at writing a C extension. I've never done that before.
##********************************************************************** ***** ##Prototype: void ComputeGaussQuadPoint(int QuadPoint, int *a.TrngleNode, ## double *SrcPointCol) ##Description: To compute the coordinates of 7-point Gauss nodes of ## a triangular patch. ##Input value: ## int QuadPoint --- node index, it can be from 0 to 6. ## int *a.TrngleNode --- the three nodes of a tringular patch. ## double *SrcPointCol --- where to store the results ##Return value: none ##Global value used: a.NodeCord, a.Qpnt ##Global value modified: none ##Subroutines called: none ##Note: Not very sure. ************************************************************************ **** def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
SrcPointCol=zeros((3),Float)
SrcPointCol[0] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]
SrcPointCol[1] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]
SrcPointCol[2] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
return SrcPointCol
-- The Numeric Python EM Project
www.members.home.net/europax
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net http://lists.sourceforge.net/lists/listinfo/numpy-discussion
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net http://lists.sourceforge.net/lists/listinfo/numpy-discussion
-- The Numeric Python EM Project www.members.home.net/europax
The good news is that your routine reduces execution time by 30%. The bad news is that the wrong numbers are coming out of it. I'm trying some otner permutations of "take" and "add.reduce" in the function to see if I can get it. One method that worked was splitting up "c" into c1,c2,c3, such that: c1=Numeric.take(a.NodeCord[:,0],TrngleNode) etc and then using Numeric.add.reduce(s*c1,1) etc This gives the right results, but is slower than plain Python. I'll keep at it. Thanks again. "Paul F. Dubois" wrote:
def cgqp(QuadPoint, TrngleNode, a): s = a.Qpnt[Quadpoint,:] c = Numeric.take(a.NodeCord, TrngleNode) return Numeric.add.reduce(s * c, axis=1)
This may or may not be right. The data structures I would have to set up to test it are too much for Sunday morning.
-----Original Message----- From: numpy-discussion-admin@lists.sourceforge.net [mailto:numpy-discussion-admin@lists.sourceforge.net] On Behalf Of Rob Sent: Sunday, August 26, 2001 8:39 AM To: numpy-discussion@lists.sourceforge.net. Subject: [Numpy-discussion] Can this function by Numpy-ized?
I finally got my FEM EM code working. I profiled it and this function uses up a big hunk of time. It performs gaussian integration over a triangle. I am trying to figure out how to slice the arrays so as to push it down into the C level. Does anyone have any ideas? Thanks, Rob.
ps. it looks to be intractible to me. Maybe I need to look at writing a C extension. I've never done that before.
##********************************************************************** ***** ##Prototype: void ComputeGaussQuadPoint(int QuadPoint, int *a.TrngleNode, ## double *SrcPointCol) ##Description: To compute the coordinates of 7-point Gauss nodes of ## a triangular patch. ##Input value: ## int QuadPoint --- node index, it can be from 0 to 6. ## int *a.TrngleNode --- the three nodes of a tringular patch. ## double *SrcPointCol --- where to store the results ##Return value: none ##Global value used: a.NodeCord, a.Qpnt ##Global value modified: none ##Subroutines called: none ##Note: Not very sure. ************************************************************************ **** def ComputeGaussQuadPoint(QuadPoint,TrngleNode,a):
SrcPointCol=zeros((3),Float)
SrcPointCol[0] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],0]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],0]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],0]
SrcPointCol[1] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],1]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],1]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],1]
SrcPointCol[2] = a.Qpnt[QuadPoint,0]*a.NodeCord[TrngleNode[0],2]\ + a.Qpnt[QuadPoint,1]*a.NodeCord[TrngleNode[1],2]\ + a.Qpnt[QuadPoint,2]*a.NodeCord[TrngleNode[2],2]
return SrcPointCol
-- The Numeric Python EM Project
www.members.home.net/europax
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net http://lists.sourceforge.net/lists/listinfo/numpy-discussion
-- The Numeric Python EM Project www.members.home.net/europax
participants (3)
-
John J. Lee -
Paul F. Dubois -
Rob