In [3]:
import numpy as np
import sys
#sys.path.append("/home/limu/Code/Monty")
#import fem_laplace as fl
%matplotlib inline
from matplotlib import pyplot as pl
from numpy import r_
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
In [4]:
nodes,faces=load("ex2fin.msh","./",seltype=2)
dist=np.sqrt(nodes[:,0]**2+nodes[:,1]**2)
amax=dist.max()
edgsel=dist>amax-0.01
#amax=nodes[:,0].max()
#edgsel=nodes[:,0]>amax-0.01*scal
edgval=5 #okrajova teplota
conduc=2. #faktor vodivosti
influx=10 #
dist=np.sqrt(nodes[:,0]**2+(nodes[:,1]+0.2)**2)
heasel=dist<0.03 #heat source
#testovani "nepripojenych" uzlu
from collections import Counter
cc=Counter(list(faces[-3:].ravel()-1))
ocnt=np.zeros(len(nodes))
ocnt[cc.keys()]=cc.values()
print("nepripojeno %i bodu"%sum(ocnt==0))
faces[2:]-=1 #opravime na indexovani od 0, nikoli od 1
edgsel[ocnt==0]=False #nepripojene vyradime
nint=sum(ocnt>0)-sum(edgsel)
nodind=np.zeros(len(nodes),dtype='int') #cislovani bodu zvlast pro interni a externi (okrajove) vrcholy
nodind[(edgsel==False)*(ocnt>0)]=np.r_[:nint]
nodind[edgsel]=np.r_[:sum(edgsel)]
#----------------
matI=np.zeros((nint,nint))
matX=np.zeros((nint,sum(edgsel)))
stat=[]
nnint=1./(2*3*4) #integral ni nj
for fx in faces[2:].T:
fwk=nodes[fx].copy()
fwk[:,2]=1
surf=abs(np.linalg.det(fwk))
alpha,beta=-np.roll((np.roll(fwk[:,:2],1,0)-fwk[:,:2]),1,0).T/surf
stat.append(list(alpha)+list(beta)+[surf])
inter=edgsel[fx]==False #
for i in np.r_[:3][inter]:
ni=nodind[fx[i]]
matI[ni,ni]+=surf/6*(nnint-(alpha[i]**2+beta[i]**2)/2.*conduc)
for j in [(i+1)%3,(i+2)%3]:
val=surf/6*(nnint-(alpha[i]*alpha[j]+beta[i]*beta[j])/2.*conduc) #grad ni grad nj
if inter[j]: matI[ni,nodind[fx[j]]]+=val
else: matX[ni,nodind[fx[j]]]+=val
Rside=-matX.dot(np.ones(sum(edgsel))*edgval)
Rside[nodind[heasel]]-=influx #influx
from scipy import sparse as sp
from scipy.sparse import linalg
matS=sp.csc_matrix(matI)
try:
Csol=linalg.splu(matS)
Cred=Csol.solve(rhs=Rside)
except:
dia=matI.diagonal()
print("diagonala:",dia.min(),dia.max(),len(dia))
ix,iy=matS.nonzero()
iix=ix[ix!=iy]
iiy=iy[ix!=iy]
qsel=matS[iix,iiy]
print("mimo ni:",qsel.max(),qsel.min(),len(iix))
else:
#--------------------------------------
dx,dy=nodes[(ocnt>0),:2].T
tmpall=np.zeros(dx.shape)
tmpall[edgsel[ocnt>0]==False]=Cred
tmpall[edgsel[ocnt>0]]=edgval
pl.scatter(dx,dy,c=tmpall)
pl.colorbar()
#innd=nodes[edgsel==False,:2]
In [5]:
size=[100,100]
xran=[dx.min(),dx.max()]
yran=[dy.min(),dy.max()]
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]]
tmap=griddata(nodes[:,:2], tmpall, (grid_x,grid_y), method='linear').T
pl.imshow(tmap,extent=np.r_[xran,yran],origin="lower")
pl.colorbar()
Out[5]:
In [ ]: