Commit afe6ac2a by Romain Hild

### start interpolation with picture

parents
interpol.py 0 → 100644
 # -*- coding: utf-8 -*- """ Spyder Editor This is a temporary script file. """ import tensorflow as tf import numpy as np import random as ran import time import matplotlib.pyplot as plt model=tf.constant(3.0,dtype=tf.float32) print(model) def test_function(x,y): z = x + y return z def random_mesh(x,h): theta=ran.uniform(-1.0,1.0) x[0] = x[0] +theta*h[0]/5 x[1] = x[1] +theta*h[1]/5 return x def id_mesh(x,h): return x class mesh: def __init__(self, Lx, Ly, Nx, Ny): self.Nx = Nx self.Ny = Ny self.Lx = Lx self.Ly = Ly self.Nc = Nx*Ny self.Nv = (Nx+1)*(Ny+1) self.hx = Lx/Nx self.hy = Ly/Ny self.h = min(self.hx,self.hy) self.nodes = np.zeros((self.Nv,2),float) self.centers = np.zeros((self.Nv,2),float) self.cells = np.zeros((self.Nc,4),int) self.stencil = [] def create(self,newmesh): k=0 for j in range(0,self.Ny+1): for i in range(0,self.Nx+1): self.nodes[k,:]=[i*self.hx,j*self.hy] if i !=0 and j !=0 and i!= self.Nx and j!= self.Ny: newmesh(self.nodes[k,:],(self.hx,self.hy)) k = k+1 k=0 for j in range(0,self.Ny): for i in range(0,self.Nx): self.cells[k,0:4] = [(self.Nx+1)*j+i,(self.Nx+1)*j+i+1,(self.Nx+1)*(j+1)+i+1,(self.Nx+1)*(j+1)+i] k = k+1 self.create_stencil() def local_nodes(self, ncell): x = M.nodes[M.cells[ncell,:],0] y = M.nodes[M.cells[ncell,:],1] return x,y def create_stencil(self): self.stencil=[1,-1,self.Nx,-self.Nx,self.Nx-1,self.Nx+1,-self.Nx-1,-self.Nx+1] class projection: def __init__(self, deg): self.deg = deg self.point = np.zeros(deg+1) self.weight = np.zeros(deg+1) def create(self): if self.deg==0: self.point[0]=0.0 self.weight[0]=2.0 if self.deg==1: self.point[0]=-0.5773502691896257 self.point[1]=0.5773502691896257 self.weight[0]=1.0 self.weight[1]=1.0 if self.deg==2: self.point[0]=-0.7745966692414834 self.point[1]=0.0 self.point[2]=0.7745966692414834 self.weight[0]=0.5555555555555556 self.weight[1]=0.8888888888888888 self.weight[2]=0.5555555555555556 def local_projection(self, M, init_f,Pi,x1,x2,x3,x4): res = 0.0 v = 0.0 for k1 in range(0,self.deg+1): for k2 in range(0,self.deg+1): xl=self.point[k1] yl=self.point[k2] s1=0.25*(1.0-xl)*(1.0-yl) s2=0.25*(1.0+xl)*(1.0-yl) s3=0.25*(1.0+xl)*(1.0+yl) s4=0.25*(1.0-xl)*(1.0+yl) x = x1[0]*s1+x2[0]*s2+x3[0]*s3+x4[0]*s4 y = x1[1]*s1+x2[1]*s2+x3[1]*s3+x4[1]*s4 j11=0.25*((1+yl)*(x3[0]-x4[0]) +(1-yl)*(x2[0]-x1[0])) j12=0.25*((1+yl)*(x3[1]-x4[1]) +(1-yl)*(x2[1]-x1[1])) j21=0.25*((1-xl)*(x4[0]-x1[0]) +(1+xl)*(x3[0]-x2[0])) j22=0.25*((1-xl)*(x4[1]-x1[1]) +(1+xl)*(x3[1]-x2[1])) j = j11*j22-j21*j12 v = v +j res= res + self.weight[k1]*self.weight[k2]*abs(j)*init_f(x,y) Pi=res/v return Pi def projection(self, M, init_f,Pi): for i in range(0,M.Nc): x1=M.nodes[M.cells[i,0],:] x2=M.nodes[M.cells[i,1],:] x3=M.nodes[M.cells[i,2],:] x4=M.nodes[M.cells[i,3],:] res = 0.0 v = 0.5*((x1[0]*x2[1]-x1[1]*x2[0])+(x2[0]*x3[1]-x2[1]*x3[0])+(x3[0]*x4[1]-x3[1]*x4[0])+(x4[0]*x1[1]-x4[1]*x1[0])) for k1 in range(0,self.deg+1): for k2 in range(0,self.deg+1): xl=self.point[k1] yl=self.point[k2] s1=0.25*(1.0-xl)*(1.0-yl) s2=0.25*(1.0+xl)*(1.0-yl) s3=0.25*(1.0+xl)*(1.0+yl) s4=0.25*(1.0-xl)*(1.0+yl) x = x1[0]*s1+x2[0]*s2+x3[0]*s3+x4[0]*s4 y = x1[1]*s1+x2[1]*s2+x3[1]*s3+x4[1]*s4 j11=0.25*((1+yl)*(x3[0]-x4[0]) +(1-yl)*(x2[0]-x1[0])) j12=0.25*((1+yl)*(x3[1]-x4[1]) +(1-yl)*(x2[1]-x1[1])) j21=0.25*((1-xl)*(x4[0]-x1[0]) +(1+xl)*(x3[0]-x2[0])) j22=0.25*((1-xl)*(x4[1]-x1[1]) +(1+xl)*(x3[1]-x2[1])) j = j11*j22-j21*j12 res= res + self.weight[k1]*self.weight[k2]*abs(j)*init_f(x,y) Pi[i]=res/v return Pi class Picture: def __init__(self, raf, M): self.raf = raf self.Nx = M.Nx self.Ny = M.Ny self.Lx = M.Lx self.Ly = M.Ly self.picture = np.zeros((raf*M.Nx,raf*M.Ny)) def test_point(self, qx, qy, x, y): j = 3 c = False for i in range(0,4): if (((qy[i]>y) != (qy[j]>y)) and (x < ((qx[j]-qx[i]) * (y - qy[i]) / (qy[j]-qy[i]) + qx[i]) )): c = not c j = i return c def point_ref_cell(self, ipoint): icell_ref = int((ipoint/self.raf) % self.Nx) + int(ipoint/(self.raf**2*self.Nx))*self.Nx return icell_ref def find_cell(self, M, ipoint, x, y ): res = False icell_ref=self.point_ref_cell(ipoint) ic = icell_ref [qx,qy]=M.local_nodes(icell_ref) res= self.test_point(qx,qy,x,y) k = 0 while not res: ic = icell_ref + M.stencil[k] if ic >= 0 and ic < M.Nc: [qx,qy]=M.local_nodes(ic) res = self.test_point(qx,qy,x,y) k =k +1 return ic def construct_picture(self, M, Pi): for i in range(0,self.raf*M.Nx): for j in range(0,self.raf*M.Ny): pixel = i*self.raf*M.Ny+j ix = (j+1)* M.Lx/(self.raf*M.Nx+1) iy = (i+1)* M.Ly/(self.raf*M.Ny+1) ic=self.find_cell(M,pixel,ix,iy) self.picture[i,j]= Pi[ic] t0=time.clock() M = mesh(1.0,1.0,200,200) M.create(random_mesh) t1=time.clock() P = projection(2) P.create() Pi = np.zeros(M.Nc) P.projection(M,test_function,Pi) t2=time.clock() Pic = Picture(4,M) Pic.construct_picture(M,Pi) t3=time.clock() print " mesh " + str(t1-t0) + " projection " + str(t2-t1) + " image " + str(t3-t2) x = np.linspace(M.Lx/(Pic.raf*M.Nx+1), Pic.raf*M.Nx* M.Lx/(Pic.raf*M.Nx+1), Pic.raf*M.Nx) y = np.linspace(M.Ly/(Pic.raf*M.Ny+1), Pic.raf*M.Ny* M.Ly/(Pic.raf*M.Ny+1), Pic.raf*M.Ny)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!