import urllib,os from numpy import r_ srpt=0 xran=[-1.5,1.5] yran=[-0.5,0] def load(ifile="sample2D_x.msh",indir="work/",seltype=2): if indir[:4]=='http': #nahrajeme soubor z webove adresy lines=urllib.request.urlopen(indir+ifile).readlines() else: lines=open(indir+ifile).readlines() #nalezeni pocatku seznamu vrcholu ibeg=[i for i in range(len(lines)) if str(lines[i]).strip()[-5:]=='Nodes'] # pole nodes obsahuje souradnice vrcholu site if len(ibeg)<2: print("only %i Nodes sign"%len(ibeg)) return None,None nodes=r_[[[float(b) for b in a.split()[1:]] for a in lines[ibeg[0]+2:ibeg[1]]]] inex=[i for i in range(len(lines)) if str(lines[i]).strip()[-8:]=='Elements'] #elementy if len(inex)<2: print("only %i Elements sign"%len(inex)) return None,None cones=[[int(b) for b in a.split()] for a in lines[inex[0]+2:inex[1]]] # vybereme z nich trojuhelniky # pole faces obsahuje indexy vrcholu a tagy faces=r_[[c for c in cones if c[1]==seltype]].transpose()[3:] return nodes,faces from numpy import where,array,zeros from numpy.linalg import det #vypocet plochy 3uhelnika def area(a,b,c): b1,c1=b-a,c-a return det([b1,c1])/2. def make_mat(nodes,faces,side,srfc=10,exval=1.): '''element from integrating grad Ni grad Nj side: is border vertex? ''' from scipy import sparse as sp from scipy import ndimage zexter=[sum([(f==b+1) for b in where(side)[0]]) for f in faces[2:5]] #elementy s hranicnimy body (urcenymi [side]) print("pocet elementu s alesp. 1 externim uzlem:",sum(zexter).sum()) coord1=[nodes[f-1][:,:2] for f in faces[2:5]] # jen 2 souradnice jsou nenulove ve 2D relco1=array([coord1[(i+1)%3]-coord1[(i+2)%3] for i in range(3)]) dist1=(relco1**2).sum(2) #ctverce vzdalenosti bodu v 3-uhelniku diag1=zeros(len(nodes)) #zfac=r_[:faces.shape[1]] for i in range(3): sel=zexter[i]==0 ids=faces[i+2][sel]-1 # index protilehleho vrcholu diag1+=ndimage.sum(dist1[i][sel]/4./surf[sel],ids,r_[:len(nodes)]) #lze take pridavat primo do coo_matrix formatu - viz nize A=sp.dia_matrix((diag1,[0]),shape=(len(nodes),len(nodes))).tocsr() B=zeros(len(nodes)) B[srpt]=srfc #(zdroj) for i in range(3): sel=zexter[i]==0#+zexter[j]==0 for j in [(i+1)%3,(i+2)%3]: vals=(relco1[i][sel]*relco1[j][sel]).sum(1)/surf[sel]/4. selj=zexter[j][sel]==0 A += sp.coo_matrix((vals[selj],(faces[i+2][sel][selj]-1,faces[j+2][sel][selj]-1)),shape=A.shape).tocsr() B[faces[i+2][sel][selj==False]-1] += vals[selj==False]*exval return A[side==False].transpose()[side==False], B[side==False] def make_mat_linear(nodes,faces,side,srfc=10,srpt=None,exval=1.,omega=0): '''element from integrating Ni grad Nj ''' surf=r_[([area(*nodes[f[2:]-1,:2]) for f in faces.transpose()])] coord1=r_[[nodes[f-1][:,:2] for f in faces[2:5]]] # jen 2 souradnice jsou nenulove I10,I01=coord1.sum(0).transpose()/3. S2=(coord1**2).sum(0) S2+=(coord1[r_[1,2,0]]*coord1).sum(0) #S2+=(coord1[r_[2,0,1]]*coord1).sum(0) I20,I02=S2.transpose()/6. I11=coord1.prod(2).sum(0)/6. I11+=(coord1[r_[1,2,0],:,0]*coord1[:,:,1]).sum(0)/12. I11+=(coord1[r_[1,2,0],:,1]*coord1[:,:,0]).sum(0)/12. relco1=coord1[r_[1,2,0]]-coord1[r_[2,0,1]] #vektor protilehle strany 3uhelnika relco2=coord1[r_[1,2,0]]*coord1[r_[2,0,1],:,::-1] relmix=relco2[:,:,0]-relco2[:,:,1] #smisene souciny #(abs(relco2-relco1).max(1)) if omega: Ldiag=relco1[:,:,1]**2*I20+relco1[:,:,0]**2*I02 Ldiag-=2*relco1.prod(2)*I11 Ldiag+=2*(relco1[:,:,1]*I10-relco1[:,:,0]*I01)*relmix Ldiag+=relmix**2 else: Ldiag=0 Ldiag=omega*Ldiag-(relco1**2).sum(2) from numpy import zeros,where,ones from scipy import sparse as sp D=sp.dia_matrix((zeros(len(nodes)),[0]),shape=(len(nodes),len(nodes))).tocsr() #side=nodes[:,2]==yran[0] zexter=[sum([(f==b+1) for b in where(side)[0]]) for f in faces[2:5]] B=zeros(len(nodes)) if srpt==None: srpt=where((nodes[:,0]==0)*(nodes[:,1]==0))[0][0] B[srpt]=srfc for i in range(3): sel=zexter[i]==0 ids=faces[i+2][sel]-1 D+=sp.coo_matrix((-Ldiag[i][sel]/surf[sel]/4,(ids,ids)),shape=D.shape).tocsr() for j in [(i+1)%3,(i+2)%3]: if omega>0: Ndiag=(relco1[i][sel]*relco1[j][sel]*r_[[I02[sel],I20[sel]]].transpose()).sum(1) Ndiag+=(relco1[i][sel,0]*relco1[j][sel,1]+relco1[i][sel,1]*relco1[j][sel,0])*I11[sel] Ndiag-=(relmix[i][sel]*relco1[j][sel,1]+relmix[j][sel]*relco1[i][sel,1])*I10[sel] Ndiag-=(relmix[i][sel]*relco1[j][sel,0]+relmix[j][sel]*relco1[i][sel,0])*I01[sel] Ndiag+=(relmix[i][sel]*relmix[j][sel]) #I00 else: Ndiag=0 Ndiag=Ndiag*omega-(relco1[i][sel]*relco1[j][sel]).sum(1) selj=zexter[j][sel]==0 #Ndiag[selj]/surf[sel][selj]/4. D += sp.coo_matrix((-Ndiag[selj]/surf[sel][selj]/4.,(faces[i+2][sel][selj]-1,faces[j+2][sel][selj]-1)),shape=D.shape).tocsr() B[faces[i+2][sel][selj==False]-1] += -Ndiag[selj==False]/surf[sel][selj==False]/4.*exval return D[side==False].transpose()[side==False], B[side==False] def solve(A,B,nodes,side,size=[200,100],intermethod='linear',exval=-1): '''solving sparse matrix interpolating on 2D grid with "size" dimensions xran,yran are global variables ''' from scipy.sparse import linalg as la Csol=la.splu(A) Cred=Csol.solve(rhs=B) from numpy import ones values=r_[Cred.transpose(),exval*ones(sum(side))] xg,yg=[r_[nodes[side==False,1],nodes[side,0]],r_[nodes[side==False,2],nodes[side,1]]] from scipy.interpolate import griddata from numpy import mgrid,random grid_x,grid_y = mgrid[xran[0]:xran[1]:1j*size[0], yran[0]:yran[1]:1j*size[1]] #random.random_integers(1,10,(10,2)).shape return griddata(r_[[xg,yg]].transpose(), values, (grid_x,grid_y), method=intermethod) def main(): global srpt nodes,faces=load() side=nodes[:,1]==yran[0] surf=array([area(*nodes[f[2:]-1,:2]) for f in faces.transpose()]) srpt=where((nodes[:,0]==0)*(nodes[:,1]==0))[0][0] A,B=make_mat(nodes,faces,side,10.,srpt,1.)