node = hou.pwd()
geo = node.geometry()

from numpy import *
from compiler.ast import flatten

E=5
miu=2.0/3.0;
D=E/(1-miu**2)*mat([[1,miu,0],[miu,1,0],[0,0,(1-miu)/2]])


alist=[]
BxyList=[]
KofElementList=[]
PointInPrim=[]
Kassemble=[[0]*3 for i in range(3)]
KassembleList=[]


BdeMatrix=mat(zeros((3,6)))
AdeMatrix=mat(zeros((3,3)))
dim=2*len(geo.points())

K=mat(zeros((dim,dim)))

attrib=geo.addAttrib(hou.attribType.Prim, "AdeMatrix",flatten([[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]))

for prim in geo.prims():
    i=0
    ReleventVTX=[prim.vertex(0).point().number(),prim.vertex(1).point().number(),prim.vertex(2).point().number()]
    PointInPrim.append(ReleventVTX)
    for vtx in prim.vertices():
        AdeMatrix[i,0]=1.0
        AdeMatrix[i,1]=vtx.point().position().x()
        AdeMatrix[i,2]=vtx.point().position().y()        
        i=i+1
    i=0
      
    AdeMatrix=AdeMatrix.I
    prim.setAttribValue(attrib,flatten(AdeMatrix.tolist()))
    alist.append(AdeMatrix)
      
    BdeMatrix[0,0]=AdeMatrix[1,0]
    BdeMatrix[0,2]=AdeMatrix[1,1]
    BdeMatrix[0,4]=AdeMatrix[1,2]
    BdeMatrix[1,1]=AdeMatrix[2,0]
    BdeMatrix[1,3]=AdeMatrix[2,1]
    BdeMatrix[1,5]=AdeMatrix[2,2]
    BdeMatrix[2,0]=AdeMatrix[2,0]
    BdeMatrix[2,2]=AdeMatrix[2,1]
    BdeMatrix[2,4]=AdeMatrix[2,2]
    BdeMatrix[2,1]=AdeMatrix[1,0]
    BdeMatrix[2,3]=AdeMatrix[1,1]
    BdeMatrix[2,5]=AdeMatrix[1,2]
    
    BxyList.append(BdeMatrix)
    
    KdeMatrix=BdeMatrix.T*D*BdeMatrix
    KofElementList.append(KdeMatrix)
    
    Kassemble[0][0]=KdeMatrix[0:2,0:2]
    Kassemble[0][1]=KdeMatrix[0:2,2:4]
    Kassemble[0][2]=KdeMatrix[0:2,4:6]
    Kassemble[1][0]=KdeMatrix[2:4,0:2]
    Kassemble[1][1]=KdeMatrix[2:4,2:4]
    Kassemble[1][2]=KdeMatrix[2:4,4:6]
    Kassemble[2][0]=KdeMatrix[4:6,0:2]
    Kassemble[2][1]=KdeMatrix[4:6,2:4]
    Kassemble[2][2]=KdeMatrix[4:6,4:6]
    
    KassembleList.append(Kassemble)
    
    for i in range(3):
        for j in range(3):
            K[ReleventVTX[i]*2:(ReleventVTX[i]*2+2),ReleventVTX[j]*2:(ReleventVTX[j]*2+2)]+=Kassemble[i][j]

fixedpoints=[]
fixedrows=[]
fixednumbers=len(geo.findPointGroup("fixed").points())
pointgroups=geo.findPointGroup("fixed").points()
forcelist=list(geo.floatListAttribValue("forcelist"))

for i in range(fixednumbers):
    fixedpoints.append(pointgroups[i].number())
    fixedrows.append(2*fixedpoints[i])
    fixedrows.append(2*fixedpoints[i]+1)

notfixed=range(len(geo.points()))

notfixed=list(set(notfixed)-set(fixedpoints))



KKernal=delete(K,fixedrows,axis=0)
KKernal=delete(KKernal,fixedrows,axis=1)



forcelist=array(list(geo.floatListAttribValue("forcelist")))
forcelist=delete(forcelist,fixedrows)



displacement=KKernal.I*mat(forcelist).T


j=0
for i in notfixed:
    point=geo.iterPoints()[i]
    vector=hou.Vector3(displacement[2*j,0],displacement[2*j+1,0],0)
    point.setPosition(point.position()+vector)
    j=j+1

int ptnum=i@totalnum;
float forcelist[];
for(int i=0;i<ptnum;i++)
{
    float force[]=point(0,"force",i);
    forcelist[2*i]=force[0];
    forcelist[2*i+1]=force[1];
}

f[]@forcelist=forcelist;