In [1]:
import numpy as np
import sys
#sys.path.append("/home/limu/Code/Monty")
%matplotlib inline
from matplotlib import pyplot as pl
import fem_laplace as fl

conduc=2. #conductivity
influx=10 #source intensity

$M_T(0,0,0)=1$
$M_T(1,0,0)=1/4(x_a+x_b+x_c+x_d)$
$M_T(0,1,0)=1/4(y_a+y_b+y_c+y_d)$

$M_T(2,0,0)=1/10(x^2_a+x^2_b+x^2_c+x^2_d+x_a x_b+x_b x_c+x_c x_d+x_a x_c+x_a x_d+x_b x_d)$

$M_J(1,1,0)=1/20 (2(x_a y_a+x_b y_b+x_c y_c+x_d y_d)+ +x_a y_b+x_a y_c+x_a y_d+x_b y_a+x_b y_c+x_b y_d+ +x_c y_a+x_c y_b+x_c y_d+x_d y_a+x_d y_b+x_d y_c)$

In [58]:
nodes,faces=fl.load("ex3.msh","./",seltype=4)
nodes2,faces2=fl.load("ex3.msh","./",seltype=2)
faces-=1
faces2-=1
dim=3
coord1=np.r_[[nodes[f] for f in faces[2:].T]] # 
coord1.shape
Out[58]:
(1402L, 4L, 3L)
$$M_T [α, β] = \frac{1}{(d + 1)(d + 2)}\left( \sum^{d+1} x_i[α] \sum^{d+1} x_i[β] + \sum^{d+1} x_i[α] x_i[β]\right)$$
In [59]:
#momentove integraly
def M(coord,i=1,j=0,lim=10):
    d=coord.shape[2]
    if (i==0) or (j==0):
        k=i+j
        return coord[:lim,:,k-1].sum()/(d+1)
    else: 
        rep=(((coord[:lim,:,i-1]*coord[:lim,:,j-1])).sum(1))
        rep+=coord[:lim,:,i-1].sum(1)*coord[:lim,:,j-1].sum(1)
        return rep/(d+1)/(d+2)
rep2=M(coord1,1,1,None)
rep2[:10]
Out[59]:
array([ 0.03568036,  0.00551643,  0.01362373,  0.00776698,  0.01211163,
        0.00048448,  0.01007916,  0.00294401,  0.01455487,  0.00780304])
In [60]:
x_a,x_b,x_c,x_d=coord1[:,:,0].T
rep=1/10.*(x_a**2+x_b**2+x_c**2+x_d**2+x_a *x_b+x_b* x_c+x_c *x_d+x_a *x_c+x_a *x_d+x_b *x_d)
try:
    assert np.allclose(rep[:10],rep2[:10],0.01)
except:
    print(rep[:10]/rep2[:10])

overeni hranicnich ploch

In [11]:
def alfac(ve):
    ve=np.sort(ve)
    return ["%03i%03i%03i"%tuple(list(ve[np.r_[:4]!=i])) for i in range(4)]
allvod=[]
for f in faces[2:].T:
    allvod.extend(alfac(f))
len(allvod) #pocet sten vsech ctyrstenu
Out[11]:
5608
In [12]:
from collections import Counter
cc=Counter(allvod)
side=[z for z,w in cc.items() if w==1]
print(len(side))
s=set()
for a in side:
    s=s.union(set([int(a[i:i+3]) for i in range(0,9,3)]))
len(s)
pside=np.array(list(s))
xyside=nodes[pside-1] #souradnice hranicnich vrcholu

696
In [13]:
#primo z plosek zapsanych v meshi
xyside=nodes[faces2[2:].ravel()-1]
dis1=np.sqrt(xyside[:,0]**2+xyside[:,1]**2)
pl.plot(dis1,xyside[:,2],',')
Out[13]:
[<matplotlib.lines.Line2D at 0x7a1de10>]
In [21]:
edglist=list(faces2[2:].ravel()-1)
edglist=np.unique(edglist)
print("celkem %i hranice %i dolni podstava %i bodu"%(len(nodes),len(edglist),sum(nodes2[edglist,2]==0)))

celkem 418 hranice 350 dolni podstava 213 bodu

vypocet determinantu (objemu)

In [15]:
def maperm(dim=4):
    #seznam vsech permutaci a jejich parita
    import itertools as its
    orig=np.r_[:dim]
    allper=np.array([a for a in its.permutations(orig)])
    zign=np.r_[[abs(a-orig).sum() for a in allper]]
    pari=-np.ones(len(zign))
    pari[0]=0
    pval=1
    pari[zign==2]=pval
    while (np.any(pari<0)):
        pval=pval+1
        for sper in allper[pari==pval-1]:
            pari[(abs(allper-sper).sum(1)==2)*(pari<0)]=pval
    return allper,1-2*(pari%2)
vec1=coord1[:,1:,:]-coord1[:,0,:][:,np.newaxis,:]
allper,sign=maperm(vec1.shape[1])
orig=allper[0] #vec1.diagonal(0,1,2).sum(1)
vol1=vec1[:,orig,orig].prod(1)
for i in range(1,len(sign)):
    vol1+=vec1[:,orig,allper[i]].prod(1)*sign[i]
ok=pl.hist(vol1*1000,20)
In [99]:
nodind=np.zeros(len(nodes),dtype='int') #cislovani bodu zvlast pro interni a externi 
try:
    nedg=len(edglist)
    nint=len(nodes)-nedg
    nodind[edglist]=np.r_[1:nedg+1]
    nodind[nodind==0]=np.r_[:nint]
    nodind[edglist]=np.r_[:nedg]
    edgsel=np.zeros_like(nodind,dtype='bool')
    edgsel[edglist]=True
except:
    print("edgsel defined")
    nint=len(nodes)-sum(edgsel)
    nodind[edgsel==False]=np.r_[:nint]
    nodind[edgsel]=np.r_[:sum(edgsel)]
conduc=2. #conductivity
influx=10 #source intensity
In [63]:
def angles(fwk):
    orig=np.r_[:len(fwk)]
    alldet=[[np.linalg.det(fwk[:,orig!=i][orig!=j])*(2*((i+j)%2)-1) for i in range(len(fwk))] for j in range(len(fwk))]
    return np.array(alldet)
def angles2(fwk):
    orig=np.r_[:len(fwk)]
    return np.array([np.linalg.solve(fwk,orig==i) for i in orig])
fwk=np.concatenate([nodes[faces[2:,20]].T,np.ones(dim+1).reshape(1,dim+1)]).T
rep2=angles2(fwk)
rep1=angles(fwk)/vol1[20]
rep1/rep2
Out[63]:
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [-0., -0.,  1., -0.]])
In [108]:
cntedg=np.r_[[edgsel[f[2:]].sum() for f in faces.T]] #kolik vrcholu je hranicnich
facsub=np.r_[:len(cntedg)][cntedg<=1]
len(facsub)
Out[108]:
306
In [114]:
dim=3
debug=False
matI=np.zeros((nint,nint))
matX=np.zeros((nint,nedg))
stat=[]
nnint=1./(2*3*4) #integral ni nj
ipos=0
#for fx in faces[2:].T[:30]:
for fx in faces[2:,:].T:
    fwk=np.concatenate([nodes[fx].T,np.ones(dim+1).reshape(1,dim+1)]).T
    #print(fx)
    surf=abs(np.linalg.det(fwk))
    alpha,beta,gamma,delta=angles(fwk)/surf
    if debug: stat.append(list(alpha)+list(beta)+[surf,vol1[ipos]])
    #inter=[i not in edglist for i in fx] #
    inter=edgsel[fx]==False
    #print(inter)
    ipos+=1
    for i in np.r_[:dim+1][inter]:
        #if fx[i] in edglist: continue
        ni=nodind[fx[i]]
        matI[ni,ni]+=surf/6*(nnint-(alpha[i]**2+beta[i]**2+gamma[i]**2)/2.*conduc)
        for j in (np.r_[1:dim+1]+i)%(dim+1):
            val=surf/6*(nnint-(alpha[i]*alpha[j]+beta[i]*beta[j]+gamma[i]*gamma[j])/2.*conduc) #grad ni grad nj
            if inter[j]: matI[ni,nodind[fx[j]]]+=val
            else: matX[ni,nodind[fx[j]]]+=val
In [115]:
matX.sum()
Out[115]:
1.2343041670638268
In [119]:
astat=np.r_[stat]
astat.min(0)
Out[119]:
array([ -1.39583392e+02,  -1.74479982e+02,  -1.11018379e+02,
        -7.29539529e+01,  -2.37946149e+02,  -9.74767575e+01,
        -1.30777408e+02,  -7.29641805e+01,   1.45095719e-05,
         1.45095719e-05])

jemnejsi sit

In [14]:
nodes2,faces2=fl.load("ex3plus.msh","/Users/Admin/Dokumenty/Lab/",seltype=4)
coord2=np.r_[[nodes2[f-1] for f in faces2[2:].T]]
vec1=coord2[:,1:,:]-coord2[:,0,:][:,np.newaxis,:]
allper,sign=maperm(vec1.shape[1])
orig=allper[0] #vec1.diagonal(0,1,2).sum(1)
vol1b=vec1[:,orig,orig].prod(1)
for i in range(1,len(sign)):
    vol1b+=vec1[:,orig,allper[i]].prod(1)*sign[i]
ok=pl.hist(vol1b*1000,20)
In [15]:
nodes2,faces2b=fl.load("ex3plus.msh","/Users/Admin/Dokumenty/Lab/",seltype=2)
coord2b=np.r_[[nodes2[f-1] for f in faces2b[2:].T]]
In [16]:
#vypocet alpha, beta (pomoci subdeterminantu)
k=20
dim=3
fx=faces[2:,k]-1
fwk=np.array(list(nodes[fx].T)+[np.ones(4)]).T # souradnice bodu
surf=vol1[k]
slps=np.array([np.linalg.solve(fwk,np.r_[:dim+1]==i) for i in np.r_[:dim+1]])
slps.T
#delta=np.roll((np.roll(fwk[:,:dim],1,0)*fwk[:,dim-1::-1])*np.r_[1,-1],1,0).sum(1)/surf
Out[16]:
array([[  1.70104429e+01,   8.12915278e+00,  -2.51395957e+01,
         -7.05590842e-16],
       [ -1.30537406e+01,   1.52812003e+01,  -2.22745965e+00,
          4.45223890e-16],
       [ -5.75458966e+00,   4.97644521e+00,  -8.69804400e+00,
          9.47618845e+00],
       [ -3.45371253e+00,   2.16661982e+00,   2.28709271e+00,
          1.65530415e-16]])
In [17]:
surf2=np.linalg.det(fwk)
[[np.linalg.det(fwk[np.r_[:dim+1]!=i][:,np.r_[:dim+1]!=j])/surf2*(-1)**(i+j) for i in range(dim+1)] for j in range(dim+1)]
Out[17]:
[[17.010442883007972, 8.1291527753090467, -25.139595658317035, 0.0],
 [-13.053740603239614, 15.281200250908874, -2.2274596476692725, -0.0],
 [-5.7545896622231956,
  4.9764452108377935,
  -8.698043999116031,
  9.4761884505014304],
 [-3.4537125295134699, 2.1666198241642873, 2.2870927053491834, -0.0]]
In [20]:
orig=np.r_[:dim+1]
z0=orig[orig!=i]
vec1[:,z0][:,:,z0[z0!=3]].shape
Out[20]:
(23533L, 3L, 3L)
In [24]:
#efektivnejsi vypocet subdeterminantu - ve vývoji
i,j=1,2
#precislujeme puvodni sloupce na jinou kombinaci (s vynechanim sloupce i, resp. j)
#sloupec 3 jsou same 1 - v soucinu (prod(1)) se muzou vynechat
renum=np.r_[:dim+1]
orig=renum[renum!=i]
orig2=renum[renum!=j]
orig2=orig2[orig2!=3]
newper=orig[np.array(allper)]
vec1=coord1
alpx=vec1[:,orig2][:,:,orig[orig!=3]].prod(1)
for k in range(1,len(sign)):
    alpx+=vec1[:,orig2][:,:,newper[k][newper[k]!=3]].prod(1)*sign[k]