Commit d90c004d by Nicolas Fressengeas

### First Commit after I discovered Git

Début du projet pour les étudiants Supélec
parents
This diff is collapsed.
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### #This file implements the finite_difference class #The FiniteDifference class is ## based on the SymbolicArithmetic class ## useful to find neighborhood or dependency #An instance of this class is the description of the equation te be solved #I do not believe it can be either compiled or parallelized # #File FiniteDifference.sage #Author : Nicolas Fressengeas # # Each FiniteDifference object is one equation def neighborhood_two(expr,v): """neighorood_two(expr,v) parses espr to compute its neigborhood in terms of v. The neighborhhod in terms of v is a matrix of all the arguments of v appearing in self. Does the same as neigborhood_one though with sage rather that maxima. """ if expr.operator()==None: return [] elif expr.operator()==v: return vector(expr.operands()) else: return list(neighborhood_two(i,v) for i in expr.operands()) class FiniteDifference(Expression): #__init__ instances the class as SymbolicArithmetic, sustracting rhs to lhs of equation if necessary def __init__(self,eq): eq=SR(eq) if eq.is_relational(): Expression.__init__(self,SR,eq.left()-eq.right()) else: Expression.__init__(self,SR,eq) #start maxima and define a function used for the neighborhood #To my opinion, the following should work. Howeverer recursivity does not work #self.get_neighbors=maxima.function('f,v','( local_neighbors(x):=get_neighbors(x,v), if atom(f) then [] elseif op(f)=v then args(f) else map(local_neighbors,args(f)) )') #However, we can do it in Maxima directly #maxima('get_neighbors_f(f,v):=(local_neighbors(x):=get_neighbors_f(x,v), if atom(f) then [] elseif op(f)=v then (f) else map(local_neighbors,args(f)))') #maxima('get_neighbors(f,v):=map(args,flatten(get_neighbors_f(f,v)))') #Another maxima function is need to translate the arguments of a given list of functions #The commented ont is written recursively #The other one is done using pattern matching #It is hopefully faster #maxima('translate(expr,flist,v):=(translation(x):=translate(x,flist,v), if atom(expr) then expr elseif member(op(expr),flist) then apply(op(expr),args(expr)+v) else apply(op(expr),map(translation,args(expr))) )') #maxima('translate(expr,listf,v):=(matchdeclare(ff,lambda([x],if atom(x) then false else member(op(x),listf)),gg,lambda([x],if atom(x) then false else op(x)=booz)),defrule(r1,ff,booz[op(ff),args(ff)+v]),defrule(r2,gg,apply(args(gg)[1],args(gg)[2])),apply1(expr,r1,r2))') # a cform function to apply a set of rules to comply with the requirements of booz and C #maxima('cform(expr,maxlength):=(matchdeclare(ffc,true,eec,lambda([x],(x#1 and x#0)),fftmp,lambda([x],if atom(x) then false else (member(op(x),["+","*"]) and length(args(x))>maxlength))),defrule(rc,ffc^eec,pow(ffc,eec)),defrule(rtmp,fftmp,simpftmpiter(op(fftmp),args(fftmp))),apply1(expr,rc,rtmp))') maxima('cform(expr,maxargs):=(matchdeclare(ffc,true,eec,lambda([x],(x#1 and x#0)),fcplx,lambda([f],if atom(f) then false else (length(args(f))>maxargs))),defrule(rc,ffc^eec,pow(ffc,eec)),defrule(rcplx,fcplx,simpfarobas(fcplx)),applyb1(expr,rc,rcplx))') #maxima('localize(expr,flist,args):=(matchdeclare(ff,lambda([x],(member(x,flist)))),defrule(r,ff,apply(ff,args)),apply1(expr,r))') # def __repr__(self): # return SymbolicArithmetic.__repr__(self) def neighborhood_one(self,v): """Neighorood_one(v) parses self to compute its neigborhood in terms of v. The neighborhhod in terms of v is a matrix of all the arguments of v appearing in self. This is done using maxima since sage cannot (yet?) do it. """ #n=Pattern(sageobj(maxima('get_neighbors('+self.__repr__()+','+v.__repr__()+')'))) #This should work but maxima hangs #The follwong works but is supposedly slower and maxima.eval is supposedly less stable n=Pattern(sage_eval(maxima.eval('get_neighbors('+self.__repr__()+','+v.__repr__()+')'))) #This matrix can contain several identical lines #This makes ne sense in terms of patterns and has to be removed #return(n.set_ify()) return (n) def neighborhood(self,vlist): """Neighborhhod(v) parses self to compute its neighborhood in terms of the function list argument. This neighborhood is the union of all neighborhoods """ var("simpfvariable") n=Pattern([]) for v in vlist: n.extend(neighborhood_two(self,v)) # n.extend(self.neighborhood_one(v)) n=Pattern(flatten(n)) return(n.set_ify().tuple_ify()) def dependency(self,v): """The dependency is the convolution of a neighborhood by its opposite. This is purely numerical array manipulation and could be compiled using spyx files. We will not do it because it belongs to the FiniteDifference class which is not compiled. Unless there is a solution to do it nonetheless. """ dep=Pattern([]) n=self.neighborhood(v) for i in n: for j in n: dep.append(vecteurs.diffv(i,j)) return (dep.set_ify()) def translate_one(self,expr,v,fifi): """The translate_one(v,f) method translates by the amount of v the arguments of f in self. """ def trfunc(*extra): return apply(fifi,vecteurs.sumv(v,extra)) # dico={} # dico[fifi]=trfunc return expr.substitute_function(fifi,trfunc) def translation(self,v,flist): """The translation(v,flist) method translates by the amount of v the arguments of all the function appearing both in self and in flist in self. """ res=self for f in flist: res=self.translate_one(res,v,f) return res #return FiniteDifference(sageobj(maxima('translate('+self.__repr__()+','+flist.__repr__()+','+list(v).__repr__()+')'))) #Could also be (though does not work if unknown functions are used in self) : #t=sage_eval(("map(function,flist)","maxima.eval('translate('+self.__repr__()+','+flist.__repr__()+','+list(v).__repr__()+')')"),locals={'flist':flist,'self':self,'v':v}) #return (FiniteDifference(sage_eval(t))) #But using maxima is way too slow.... is sage better ? def cform(self,maxargs): """Applies a set of rules to comply with booz and C. """ return FiniteDifference(sageobj(maxima('cform('+self.__repr__()+','+maxargs.__repr__()+')'))) # def localize(self,args,flist): # """In Self, applies all functions of flist to args. # """ # return FiniteDifference(sageobj(maxima('localize('+self.__repr__()+','+flist.__repr__()+','+list(args).__repr__()+')'))) def save(self,filename,compress=True): f=open(filename,"r") f.write(self.__repr__()) f.close() Expression.save(self,filename,compress) def simplify(self): self=self.full_simplify() return self
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### # This file implements the FiniteDifferenceSystem class # It basically is a list of FiniteDifference # It thus inherits from the list class # Each FiniteDifferenceSystem object is one differential system to be solved on one point or one area load(os.path.join(os.path.expanduser(EscapadeInstallDir),'FiniteDifference.sage')) class FiniteDifferenceSystem(list): #It constructs as a list #Each item will be constructed itself with the FiniteDifference constructor def __init__(self,l): list.__init__(self,map(FiniteDifference,l)) #The representation __repr__ of class list will do the printing job def neighborhood(self,vlist): """Neighborhhod(v) parses self to compute its neighborhood in terms of the function list argument. This neighborhood is the union of all neighborhoods """ n=Pattern([]) for p in self: n.extend(p.neighborhood(vlist)) return(n.set_ify()) def dependency(self,v): """The dependency is the convolution of a neighborhood by its opposite. This is purely numerical array manipulation and could be compiled using spyx files. We will not do it because it belongs to the FiniteDifference class which is not compiled. Unless there is a solution to do it nonetheless. """ dep=Pattern([]) n=self.neighborhood(v) for i in n: for j in n: dep.append(vecteurs.diffv(i,j)) return (dep.set_ify()) def squared_norm(self): """Returns the squared norm of the system: i.e. the sum of all squared components """ norm=0 for equ in self: norm+=equ**2 return norm def translation(self,v,flist): """The translation(v,flist) method translates by the amount of v the arguments of all the function appearing both in the elements of self and in flist. """ system=[] for equ in self: system.append(equ.translation(v,flist)) return FiniteDifferenceSystem(system) def diff(self,x): """Differentiation of the elements of the list. """ diff=[] for eq in self: diff.append(eq.diff(x)) return diff def cform(self): """Applies a set of rules to comply with booz and C. """ return FiniteDifferenceSystem(sageobj(maxima('cform('+self.__repr__()+')'))) def localize(self,args,flist): """In Self, applies all functions of flist to args. """ return FiniteDifferenceSystem(sageobj(maxima('localize('+self.__repr__()+','+flist.__repr__()+','+list(args).__repr__()+')'))) def simplify(self): system=[] for equ in self: system.append(equ.simplify()) return(system)
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### #This file implements the FiniteDifferenceSystemList class # It basically is a list of FiniteDifferenceSystem # It thus inherits from the list class # The FiniteDifferenceSystemList object is the list of all differential systems to be solved over the different areas load(os.path.join(os.path.expanduser(EscapadeInstallDir),'FiniteDifferenceSystem.sage')) class FiniteDifferenceSystemList(list): #It constructs as a list #Each item will be constructed itself with the FiniteDifferenceSystem constructor def __init__(self,l): list.__init__(self,map(FiniteDifferenceSystem,l)) #The representation __repr__ of class list will do the printing job def neighborhood(self,vlist): """Neighborhhod(v) parses self to compute its neighborhood in terms of the function list argument. This neighborhood is the union of all neighborhoods """ n=Pattern([]) for p in self: n.extend(p.neighborhood(vlist)) return(n.set_ify()) def dependency(self,v): """The dependency is the convolution of a neighborhood by its opposite. This is purely numerical array manipulation and could be compiled using spyx files. We will not do it because it belongs to the FiniteDifference class which is not compiled. Unless there is a solution to do it nonetheless. """ dep=Pattern([]) n=self.neighborhood(v) for i in n: for j in n: dep.append(vecteurs.diffv(i,j)) return (dep.set_ify()) def simplify(self): systemlist=[] for system in self: systemlist.append(system.simplify()) return(systemlist)
Makefile 0 → 100644
Mesh.spyx 0 → 100644
Pattern.spyx 0 → 100644
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Copyright (C) 2009 Hubert Frauensohn # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### # This file implements the Pattern class, based on mist # It is supposed to be a liste of vectors # # As this class mostly does a numerical job, it can probably be compiled either is spyx or C++ import vecteurs class Pattern(list): #Construct as a list def __init__(self,l): list.__init__(self,l) def set_ify(self): """Remove identical elements as they do not make sense in a pattern. """ for i in self: while self.count(i)>1: self.remove(i) return(self) def tuple_ify(self): """Make each element a tuple. """ v=Pattern([]) for i in self: v.append(tuple(i)) return(v) def points(self,pos): """Returns a list of all points addressed by self when centered on pos. It basically is a translation. """ pointlist=Pattern([]) for point in self: pointlist.append(vecteurs.sumv(pos,point)) return pointlist def member(self,v): """Returns True if the intersection of self and v is not empty. """ for i in self: for j in v: if (i==j): return(True) return(False) def save(self,filename): """Saves the string representation of the pattern into the file named filename. """ f=open(filename,"w") f.write(self.__repr__()) f.close()
dd.sage 0 → 100644
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### # Functions to compute gradient and Hessian using maxima... or not so that it does not crash (!) # All explicit references to maxima have been removed at the expense of some substitutions # It seems that maxima works in the background. This may be removed in future sage # I wonder if this is faster def superdiff(f,x,n=1): # return sageobj(maxima.diff(f,x,n)) var("simpfderivee") return diff(f.subs_expr(x==simpfderivee),simpfderivee,n).subs_expr(simpfderivee==x) def grad(f,v): return vector(list(superdiff(f,x) for x in v)) def hessian(f,v): # return sageobj(maxima.hessian(f,list(v))) l=len(v) hess=matrix(SR,l) for i in range(l): hess[i,i]=superdiff(f,v[i],2) for j in range(i): hess[i,j]=superdiff(superdiff(f,v[i]),v[j]) hess[j,i]=hess[i,j] return hess def diaghessian(f,v): return (diagonal_matrix(list(1/superdiff(f,x,2) for x in v))) # A few function to help discretize a continuous differential problem # Arguments : # f is a function # x is the variable with respect to wchich the differenciation is done # dx is the discretization step #Centered derivative def ndc(f,x,dx): return((f.subs_expr(x==x+1)-f.subs_expr(x==x-1))/(2*dx)) #Left derviative def ndl(f,x,dx): return((f-f.subs_expr(x==x-1))/(dx)) #Right derivative def ndr(f,x,dx): return((f.subs_expr(x==x+1)-f)/(2*dx)) #Centered second derivative def nd2(f,x,dx): return((f.subs_expr(x==x+1)+f.subs_expr(x==x-1)-2*f)/(dx**2))
disloc.sage 0 → 100644
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### #Vector differential operators in 3D #curls def NROTC(v,r,dr): return[ ndc(v[2],r[1],dr[1])-ndc(v[1],r[2],dr[2]), ndc(v[0],r[2],dr[2])-ndc(v[2],r[0],dr[0]), ndc(v[1],r[0],dr[0])-ndc(v[0],r[1],dr[1]) ] def NROTD2C(m,r,dr): return matrix(list(NROTC(l,r,dr) for l in m.rows())) def NROTG2C(m,r,dr): return matrix(list(NROTC(l,r,dr) for l in m.columns())) #matrix cross products def MatrixCrossD(m,v): return(matrix(list(l.cross_product(v) for l in m.rows()))) def MatrixCrossG(m,v): return(matrix(list(l.cross_product(v) for l in m.columns()))) #Divergence def NDIVC(v,r,dr): return sum(list(ndc(v[i],r[i],dr[i]) for i in range(3))) def NDIVD2C(m,r,dr): return matrix(list(NDIVC(l,r,dr) for l in m.rows())) def NDIVG2C(m,r,dr): return matrix(list(NDIVC(l,r,dr) for l in m.columns())) #Les dislocs a=function("a") var("v x y z t dx dy dz dt") r=[x,y,z] dr=[dx,dy,dz] A=matrix([ [0,a(x,t),0], [0,0,0], [0,0,0]]) V=vector([v,0,0]) equ=((A-A.subs(t=t-1))/dt+NROTD2C(MatrixCrossD(A,V),r,dr)) equreal=equ[0,1] equlx=equreal.subs_expr(x-1==x,dx==dx/2) equhx=equreal.subs_expr(x+1==x,dx==dx/2) Sys=[[equreal(x=0,t=0)],[equlx(x=0,t=0)],[equhx(x=0,t=0)],[0]] N=100 mesh=numpy.zeros([N,N],dtype=int) mesh[0,:]=1 mesh[N-1,:]=2 mesh[:,0]=3 disloc=Escapade(Sys,mesh,[a]) disloc.method=1 disloc.params=[[dx,n(1/N)],[dt,1/N],[v,1]] #disloc.write_simpffiles() #Now the Init zerofunction=lambda x:0 reducedcos=lambda x:cos(x*pi/2) dislocfunction=piecewise([[(-Infinity,-1),zerofunction],[(-1,1),reducedcos],[(1,Infinity),zerofunction]]) a_init=numpy.zeros([N,N],dtype=float) relativewidth=1/4 a_init[:,0]=list(dislocfunction(((i/N)-(1/2))*2/relativewidth) for i in range(N)) disloc.simpfinit([a_init])
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Hubert Frauensohn # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### #include(os.path.join(os.path.expanduser(EscapadeInstallDir),'Pattern.spyx')) #include(os.path.join(os.path.expanduser(EscapadeInstallDir),'Mesh.spyx')) #include(os.path.join(os.path.expanduser(EscapadeInstallDir),'tosimpf.spyx')) #include(os.path.join(os.path.expanduser(EscapadeInstallDir),'tobooz.spyx')) include "Pattern.spyx" include "Mesh.spyx" include "tosimpf.spyx" include "tobooz.spyx"
setup.py 0 → 100644
 ############################################################################### # ESCAPADE # Ergonomic Solver using Cellular Automata for PArtial Differential Equation # Copyright (C) 2009 Nicolas Fressengeas # Copyright (C) 2009 Hubert Frauensohn # Distributed under the terms of the GNU General Public License (GPL), # version 2 or any later version. The full text of the GPL is available at: # http://www.gnu.org/licenses/ ############################################################################### #!/usr/local/sage-4.0.1/local/bin/python2.5 # -*- coding: cp1252 -*- # compilation: /usr/local/sage-4.0.1/local/bin/python2.5 setup.py build_ext --inplace # or ./makemodule from distutils.core import setup from distutils.core import Extension setup(name = 'vecteurs', #version = '1.0', ext_modules = [Extension('vecteurs', ['vecteurs.c']) ], #url='localhost', #author='Hub',