Commit a2b69c9a authored by Romain Hild's avatar Romain Hild

vectorialize projection

parent afe6ac2a
......@@ -5,15 +5,15 @@ Spyder Editor
This is a temporary script file.
"""
import tensorflow as tf
# 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)
# model=tf.constant(3.0,dtype=tf.float32)
print(model)
# print(model)
def test_function(x,y):
z = x + y
......@@ -125,41 +125,36 @@ class projection:
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]))
def projection(self, M, init_f):
xx1 = M.nodes[M.cells[:,0],:]
xx2 = M.nodes[M.cells[:,1],:]
xx3 = M.nodes[M.cells[:,2],:]
xx4 = M.nodes[M.cells[:,3],:]
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
res = np.zeros(M.Nc)
v = 0.5*((xx1[:,0]*xx2[:,1]-xx1[:,1]*xx2[:,0])+(xx2[:,0]*xx3[:,1]-xx2[:,1]*xx3[:,0])+(xx3[:,0]*xx4[:,1]-xx3[:,1]*xx4[:,0])+(xx4[:,0]*xx1[:,1]-xx4[:,1]*xx1[:,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 = xx1[:,0]*s1+xx2[:,0]*s2+xx3[:,0]*s3+xx4[:,0]*s4
y = xx1[:,1]*s1+xx2[:,1]*s2+xx3[:,1]*s3+xx4[:,1]*s4
j11=0.25*((1+yl)*(xx3[:,0]-xx4[:,0]) +(1-yl)*(xx2[:,0]-xx1[:,0]))
j12=0.25*((1+yl)*(xx3[:,1]-xx4[:,1]) +(1-yl)*(xx2[:,1]-xx1[:,1]))
j21=0.25*((1-xl)*(xx4[:,0]-xx1[:,0]) +(1+xl)*(xx3[:,0]-xx2[:,0]))
j22=0.25*((1-xl)*(xx4[:,1]-xx1[:,1]) +(1+xl)*(xx3[:,1]-xx2[:,1]))
j = j11*j22-j21*j12
res = res + self.weight[k1]*self.weight[k2]*abs(j)*init_f(x,y)
Pi = res/v
return Pi
......@@ -219,14 +214,15 @@ M.create(random_mesh)
t1=time.clock()
P = projection(2)
P.create()
Pi = np.zeros(M.Nc)
P.projection(M,test_function,Pi)
Pi = P.projection(M,test_function)
# print(Pi)
t2=time.clock()
Pic = Picture(4,M)
Pic.construct_picture(M,Pi)
print "coucou " + str(Pic.picture)
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)
# 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)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment