# # SageMath code accompaniment to the manuscript # Computational Serendipity and Tensor Product Finite Element Differential Forms'' # by Andrew Gillette, Tyler Kloefkorn Victoria Sanders # Version date: May 29, 2018 # x,y,z=var('x,y,z') U = CoordinatePatch((x,y,z)) F = DifferentialForms(U) ZZ=IntegerRing() from sage.tensor.differential_form_element import d import itertools import numpy #################### # # standard spanning set construction tools # #################### def k(form): form._cleanup() kappa_form = DifferentialForm(F, form._degree - 1) for comp in form._components: fun = form._components[comp].factor() i = 1 for n, coord in enumerate(F.base_space().coordinates()): if n in comp: newcomp = tuple(a for a in comp if a != n) kappa_form[newcomp] += fun*coord*(-1)^(i+1) i += 1 kappa_form._cleanup() return kappa_form def HrLkRn(r,k,n): U = CoordinatePatch(var(n)) F = DifferentialForms(U) monomials = [] for mon in itertools.product(var(n),repeat=r): monomials.append(prod(mon)) allforms = [] for dxs in itertools.combinations(range(0,F.ngens()),k): for mon in monomials: form = DifferentialForm(F,k) form[dxs] = mon if form not in allforms: allforms.append(form) return allforms def PrLkRn(r,k,n): full = [] for i in range(0,r+1): full.extend(HrLkRn(i,k,n)) return full def JrLkRn(r,j,n): Hs = [] for l in range(1,len(var(n))-j): eH = HrLkRn(r+l-1,j+1,n) if j+1 != 0: for form in eH: for comp in form._components: ldeg = 0 fun = form[comp] for w,coord in enumerate(F.base_space().coordinates()): if w not in comp and fun.degree(coord) == 1: ldeg += 1 if ldeg >= l: Hs.append(form) final = [] for form1 in Hs: new = k(form1) if new not in final: final.append(new) return final def dJrLkRn(r,j,n): allforms = [] for element in JrLkRn(r,j,n): new = d(element) allforms.append(new) return allforms def SrLkRn(a,b,c): full=PrLkRn(a,b,c) full.extend(JrLkRn(a,b,c)) for element in JrLkRn(a+1,b-1,c): new=d(element) new._cleanup() full.append(new) return full def trimmedPrLkRn(a,b,c): full = PrLkRn(a-1,b,c) for form in HrLkRn(a-1,b+1,c): new = k(form) new._cleanup() if new not in full: full.append(new) return full def trimmedSrLkRn(a,b,c): full = trimmedPrLkRn(a,b,c) full.extend(JrLkRn(a,b,c)) for form in JrLkRn(a,b-1,c): new = d(form) new._cleanup() full.append(new) return full def Prst(r,s,t): full=[] for i in range(0,r+1): for j in range(0,s+1): for k in range(0,t+1): polynomial=x^i*y^j*z^k if polynomial not in full: full.append(polynomial) return full def QrLkRn(r,k,n): full=[] U = CoordinatePatch(var(n)) F = DifferentialForms(U) for alternator in itertools.combinations(range(0,F.ngens()),k): list1=[0,0,0] for i in range(0,F.ngens()): list1[i]=r for i in alternator: list1[i]-=1 for polynomial in Prst(list1[0],list1[1],list1[2]): form=DifferentialForm(F,k) form[alternator]=polynomial full.append(form) return full #################### # # computational basis subspace builders # #################### def vertices_k0(): full=[] var = F.base_space().coordinates() for sign1 in [1,-1]: for sign2 in [1,-1]: sign3=numpy.array([+1,-1]) form1,form2=(var[0]+sign1)*(var[1]+sign2)*(var[2]+sign3) full.extend((form1,form2)) return full def edges_k0(i): full=[] var = F.base_space().coordinates() for j in range(0,3): k=ZZ((j+1)%3) l=ZZ((j+2)%3) for sign1 in [1,-1]: sign2=numpy.array([+1,-1]) form1,form2=var[j]^i*(var[k]+sign1)*(var[l]+sign2)*(var[j]^2-1) full.extend((form1,form2)) return full def faces_tensor_k0(i): full=[] sign=numpy.array([+1,-1]) var=F.base_space().coordinates() for a in range(0,3): b=ZZ((a+1)%3) c=ZZ((a+2)%3) for j in range(0,i-1): k=i-1 form1,form2=var[a]^j*var[b]^k*(var[c]+sign)*(var[a]^2-1)*(var[b]^2-1) form3,form4=var[a]^k*var[b]^j*(var[c]+sign)*(var[a]^2-1)*(var[b]^2-1) full.extend((form1,form2,form3,form4)) form5,form6=var[a]^(i-1)*var[b]^(i-1)*(var[c]+sign)*(var[a]^2-1)*(var[b]^2-1) full.extend((form5,form6)) return full def faces_k0(i): full=[] var=F.base_space().coordinates() for v1 in range(0,3): v2=ZZ((v1+1)%3) v3=ZZ((v1+2)%3) for j in range(0,i-3): k=i-4-j sign=numpy.array([+1,-1]) form1,form2=var[v1]^j*var[v2]^k*(var[v3]+sign)*(var[v1]^2-1)*(var[v2]^2-1) full.extend((form1,form2)) return full def interior_tensor_k0(i): full=[] for j in range(0,i): for k in range(0,i): l=i-1 form1=x^j*y^k*z^l*(x^2-1)*(y^2-1)*(z^2-1) form2=x^l*y^j*z^k*(x^2-1)*(y^2-1)*(z^2-1) form3=x^k*y^l*z^j*(x^2-1)*(y^2-1)*(z^2-1) for form in [form1,form2,form3]: if form not in full: full.append(form) return full def interior_k0(i): full=[] for j in range(0,i-5): for k in range(0,i-5): l=i-6-j-k if l>=0: form=x^j*y^k*z^l*(x^2-1)*(y^2-1)*(z^2-1) full.append(form) return full def edges_k1(i): full=[] var = F.base_space().coordinates() sign1=numpy.array([+1,-1]) sign2=numpy.array([+1,-1]) for j in range(0,3): k=ZZ((j+1)%3) l=ZZ((j+2)%3) for sign in sign1: form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form1[j],form2[j]=var[j]^i*(var[k]+sign)*(var[l]+sign2) full.extend((form1,form2)) return full def edges_tilde_k1(i): full=[] var=F.base_space().coordinates() for j in range(0,3): k=ZZ((j+1)%3) l=ZZ((j+2)%3) for sign1 in [1,-1]: for sign2 in [1,-1]: form=DifferentialForm(F,1) form[j]=(i+1)*var[j]^i*(var[k]+sign1)*(var[l]+sign2) form[k]=var[j]^(i-1)*(var[l]+sign2)*(var[j]^2-1) form[l]=var[j]^(i-1)*(var[k]+sign1)*(var[j]^2-1) full.append(form) return full def faces_k1(i): full=[] sign=numpy.array([+1,-1]) var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) for j in range(0,i-1): k=i-2-j form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form3=DifferentialForm(F,1) form4=DifferentialForm(F,1) form1[u],form2[u]=var[u]^j*var[v]^k*(var[w]+sign)*(var[v]^2-1) form3[u],form4[u]=var[u]^j*var[w]^k*(var[v]+sign)*(var[w]^2-1) full.extend((form1,form2,form3,form4)) return full def faces_tilde_k1(i): full=[] if i<2: return full sign=numpy.array([+1,-1]) var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form3=DifferentialForm(F,1) form4=DifferentialForm(F,1) form1[u],form2[u]=var[v]^(i-2)*(var[w]+sign)*(var[v]^2-1) form3[u],form4[u]=var[w]^(i-2)*(var[v]+sign)*(var[w]^2-1) full.extend((form1,form2,form3,form4)) for j in range(1,i-1): k=i-2-j form5=DifferentialForm(F,1) form6=DifferentialForm(F,1) form5[0],form6[0]=(var[0]^j)*(var[1]^k)*(var[2]+sign)*(var[1]^2-1) form5[1],form6[1]=-(var[0]^(j-1))*(var[1]^(k+1))*(var[2]+sign)*(var[0]^2-1) full.extend((form5,form6)) for j in range(1,i-1): k=i-2-j form5=DifferentialForm(F,1) form6=DifferentialForm(F,1) form5[0],form6[0]=(var[0]^j)*(var[2]^k)*(var[1]+sign)*(var[2]^2-1) form5[2],form6[2]=-(var[0]^(j-1))*(var[2]^(k+1))*(var[1]+sign)*(var[0]^2-1) full.extend((form5,form6)) for j in range(1,i-1): k=i-2-j form5=DifferentialForm(F,1) form6=DifferentialForm(F,1) form5[1],form6[1]=(var[1]^j)*(var[2]^k)*(var[0]+sign)*(var[2]^2-1) form5[2],form6[2]=-(var[1]^(j-1))*(var[2]^(k+1))*(var[0]+sign)*(var[1]^2-1) full.extend((form5,form6)) return full def faces_hat_k1(i): full=[] if i<2: return full sign1=numpy.array([+1,-1]) sign2=numpy.array([-1,-1]) var = F.base_space().coordinates() for k in range(0,3): other=[0,1,2] other.pop(other.index(k)) for l in range(len(other)): m=other[l] n=other[ZZ((l+1)%2)] form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form1[k],form2[k]=var[m]^(i-2)*(var[n]+sign1)*(var[m]^2-1) for j in range(1,i-1): form3=DifferentialForm(F,1) form4=DifferentialForm(F,1) form3[k],form4[k]=(i+1)*var[k]^j*var[m]^(i-j-2)*(var[n]+sign1)*(var[m]^2-1) form3[n],form4[n]=var[k]^(j-1)*var[m]^(i-j-2)*(var[k]^2-1)*(var[m]^2+sign2) full.extend((form3,form4)) full.extend((form1,form2)) return full def interior_k1(i): full=[] var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) for j in range(0,i-3): for k in range(0,i-j-3): l=i-4-j-k poly1=var[u]^j*var[v]^k*var[w]^l poly2=var[u]^l*var[v]^j*var[w]^k poly3=var[u]^k*var[v]^l*var[w]^j polys=numpy.array([poly1,poly2,poly3]) form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form3=DifferentialForm(F,1) form1[u],form2[u],form3[u]=polys*(var[v]^2-1)*(var[w]^2-1) # full.extend((form1,form2,form3)) full.append(form1) return full def interior_tilde_k1(i): full=[] if i<4: return full var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) form1=DifferentialForm(F,1) form2=DifferentialForm(F,1) form1[u]=var[v]^(i-4)*(var[v]^2-1)*(var[w]^2-1) form2[u]=var[w]^(i-4)*(var[v]^2-1)*(var[w]^2-1) full.extend((form1,form2)) if i==4: full.pop() for j in range(1,i-3): k=i-4-j form3=DifferentialForm(F,1) form3[0]=(var[0]^j)*(var[1]^k)*(var[1]^2-1)*(var[2]^2-1) form3[1]=-(var[0]^(j-1))*(var[1]^(k+1))*(var[0]^2-1)*(var[2]^2-1) full.append(form3) for j in range(1,i-3): k=i-4-j form3=DifferentialForm(F,1) form3[0]=(var[0]^j)*(var[2]^k)*(var[1]^2-1)*(var[2]^2-1) form3[2]=-(var[0]^(j-1))*(var[2]^(k+1))*(var[0]^2-1)*(var[1]^2-1) full.append(form3) for j in range(1,i-3): k=i-4-j form3=DifferentialForm(F,1) form3[1]=(var[1]^j)*(var[2]^k)*(var[0]^2-1)*(var[2]^2-1) form3[2]=-(var[1]^(j-1))*(var[2]^(k+1))*(var[0]^2-1)*(var[1]^2-1) full.append(form3) if i==5: full.pop() return full def faces_k2(i): full=[] sign=numpy.array([+1,-1]) var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) for j in range(0,i+1): k=i-j form1=DifferentialForm(F,2) form2=DifferentialForm(F,2) if v>w: v,w = w,v form1[v,w],form2[v,w]=var[v]^j*var[w]^k*(var[u]+sign) form1._cleanup() form2._cleanup() full.extend((form1,form2)) return full def faces_tilde_k2(i): full=[] if i<1: return full sign=numpy.array([+1,-1]) var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) form1=DifferentialForm(F,2) form2=DifferentialForm(F,2) form3=DifferentialForm(F,2) form4=DifferentialForm(F,2) if v>w: v,w = w,v form1[v,w],form2[v,w]=(i+1)*var[v]^i*(var[u]+sign) form1[u,w]=var[v]^(i-1)*(var[v]^2-1) form2[u,w]=var[v]^(i-1)*(var[v]^2-1) form3[v,w],form4[v,w]=(i+1)*var[w]^i*(var[u]+sign) form3[u,v]=-var[w]^(i-1)*(var[w]^2-1) form4[u,v]=-var[w]^(i-1)*(var[w]^2-1) full.extend((form1,form2,form3,form4)) for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) for j in range(1,i): k=i-j form5=DifferentialForm(F,2) form6=DifferentialForm(F,2) form5[v,w],form6[v,w]=(i+2)*var[v]^j*var[w]^k*(var[u]+sign) form5[u,w]=var[v]^(j-1)*var[w]^k*(var[v]^2-1) form6[u,w]=var[v]^(j-1)*var[w]^k*(var[v]^2-1) form5[u,v]=-var[v]^j*var[w]^(k-1)*(var[w]^2-1) form6[u,v]=-var[v]^j*var[w]^(k-1)*(var[w]^2-1) if u==1: form5=-form5 form6=-form6 full.extend((form5,form6)) return full def interior_k2(i): full=[] var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) for j in range(0,i-1): for k in range(0,i-j-1): l=i-2-j-k poly1=var[u]^j*var[v]^k*var[w]^l form1=DifferentialForm(F,2) form1[v,w]=poly1*(var[u]^2-1) if u==1: form1=-form1 full.append(form1) return full def interior_tilde_k2(i): full=[] if i<2: return full var=F.base_space().coordinates() for u in range(0,3): v=ZZ((u+1)%3) w=ZZ((u+2)%3) form1=DifferentialForm(F,2) form1[v,w]=var[u]^(i-2)*(var[u]^2-1) if u==1: form1=-form1 full.append(form1) for j in range(1,i-1): k=i-j-2 form2=DifferentialForm(F,2) form3=DifferentialForm(F,2) form4=DifferentialForm(F,2) form2[1,2]=var[0]^k*var[1]^j*(var[0]^2-1) form2[0,2]=-var[0]^(k+1)*var[1]^(j-1)*(var[1]^2-1) form3[1,2]=var[0]^k*var[2]^j*(var[0]^2-1) form3[0,1]=var[0]^(k+1)*var[2]^(j-1)*(var[2]^2-1) form4[0,2]=var[2]^j*var[1]^k*(var[1]^2-1) form4[0,1]=-var[2]^(j-1)*var[1]^(k+1)*(var[2]^2-1) full.extend((form2,form3,form4)) for j in range(0,i-1): for k in range(1,i-1-j): l=i-2-j-k form5=DifferentialForm(F,2) form5[1,2]=var[0]^j*var[1]^k*var[2]^l*(var[0]^2-1) form5[0,2]=-var[0]^(j+1)*var[1]^(k-1)*var[2]^l*(var[1]^2-1) form5[0,1]=var[0]^(j+1)*var[1]^k*var[2]^(l-1)*(var[2]^2-1) full.append(form5) if l==0: full.pop() return full #################### # # commands to build computational bases # Sprk3(r) = serendipity k-forms in 3D, order r given as input # Smrk3(r) = trimmed serendipity k-forms in #d, order r given as input # #################### def Spr03(r): full=vertices_k0() for i in range(0,r-1): full.extend(edges_k0(i)) for i in range(4,r+1): full.extend(faces_k0(i)) for i in range(6,r+1): full.extend(interior_k0(i)) return full def Spr13(r): full=[] for i in range(0,r): full.extend(edges_k1(i)) full.extend(edges_tilde_k1(r)) for i in range(2,r): full.extend(faces_k1(i)) full.extend(faces_hat_k1(r)) for i in range(4,r+1): full.extend(interior_k1(i)) return full def Smr13(r): full=[] for i in range(0,r): full.extend(edges_k1(i)) for i in range(2,r): full.extend(faces_k1(i)) full.extend(faces_tilde_k1(r)) for i in range(4,r): full.extend(interior_k1(i)) full.extend(interior_tilde_k1(r)) return full def Spr23(r): full=[] for i in range(0,r): full.extend(faces_k2(i)) full.extend(faces_tilde_k2(r)) for i in range(2,r+1): full.extend(interior_k2(i)) return full def Smr23(r): full=[] for i in range(0,r): full.extend(faces_k2(i)) for i in range(2,r): full.extend(interior_k2(i)) full.extend(interior_tilde_k2(r)) return full MR3.=QQbar[];MR3 def find_max_k0(basis): deg=0 basis_list=sage_eval(str(basis),locals={'x':x, 'y':y, 'z':z}) for element in basis_list: if element != 1: degrees=element.degrees() if max(degrees)>=deg: deg=max(degrees) return deg def find_max_k1(basis): deg=0 for form in basis: for i in range(0,3): polynomial=MR3(form[i]) degrees=polynomial.degrees() if max(degrees)>=deg: deg=max(degrees) return deg def find_max_k2(basis): deg=0 for form in basis: for i in [(1,0),(1,2),(0,2)]: polynomial=MR3(form[i]) degrees=polynomial.degrees() if max(degrees)>=deg: deg=max(degrees) return deg def find_coeffs_k0(basis_list,deg): basis_list=sage_eval(str(basis_list),locals={'x':x, 'y':y, 'z':z}) final_list=[] for form in basis_list: coeff_list=[] for xdeg in range(0,deg+1): for ydeg in range(0,deg+1): for zdeg in range(0,deg+1): poly=MR3(form) c=poly.coefficient({x:xdeg,y:ydeg,z:zdeg}) coeff_list.append(c) final_list.append(coeff_list) return final_list def find_coeffs_k1(basis_list,deg): final_list=[] for form in basis_list: coeff_list=[] for i in range(0,3): for xdeg in range(0,deg+1): for ydeg in range(0,deg+1): for zdeg in range(0,deg+1): poly=MR3(form[i]) c=poly.coefficient({x:xdeg,y:ydeg,z:zdeg}) coeff_list.append(c) final_list.append(coeff_list) return final_list def find_coeffs_k2(basis_list,deg): final_list=[] for form in basis_list: coeff_list=[] for i in [(1,0),(1,2),(0,2)]: for xdeg in range(0,deg+1): for ydeg in range(0,deg+1): for zdeg in range(0,deg+1): poly=MR3(form[i]) c=poly.coefficient({x:xdeg,y:ydeg,z:zdeg}) coeff_list.append(c) final_list.append(coeff_list) return final_list def coeff_matrix_k3(basis_list): deg=0 final_list=[] for form in basis_list: polynomial=MR3(form[0,1,2]) degrees=polynomial.degrees() if max(degrees)>=deg: deg=max(degrees) for form in basis_list: coeff_list=[] for xdeg in range(0,deg+1): for ydeg in range(0,deg+1): for zdeg in range(0,deg+1): poly=MR3(form[0,1,2]) c=poly.coefficient({x:xdeg,y:ydeg,z:zdeg}) coeff_list.append(c) final_list.append(coeff_list) A=matrix(QQ,final_list) return A #################### # # sample basis check command # #################### def basis_check(): print("Checking regular serendipity, r=1, k=1, n=3") rr=1 print("Computing standard spanning set") standard=SrLkRn(rr,1,'x,y,z') print("Standard spanning set length = "+str(len(standard))) print("Computing computational basis") comp=Spr13(rr) print("Computational basis length = "+str(len(comp))) degS=find_max_k1(standard) degC=find_max_k1(comp) deg=max(degS,degC) print("Setting deg to be ",deg) print("Coeffs for standard spanning set") std_coeffs=find_coeffs_k1(standard,deg) print("Coeffs for computational basis") comp_coeffs=find_coeffs_k1(comp,deg) A=matrix(QQ,std_coeffs) B=matrix(QQ,comp_coeffs) both=std_coeffs+comp_coeffs C=matrix(QQ,both) print(comp) print("Rank of standard spanning set: "+str(A.rank())) print("Rank of computational basis: "+str(B.rank())) print("Rank of both: "+str(C.rank())) #################### # # debugging tools # #################### # these commands take in a set, print its length, builds the coefficient matrix, then prints and computes its rank. They can be used to check linear independence of a set of 1-forms or 2-forms def lin_ind_set_check_k1(spanset): print("Input set length: "+str(len(spanset))) deg=find_max_k1(spanset) coeffs=find_coeffs_k1(spanset,deg) A=matrix(QQ,coeffs) print("Rank of input set: "+str(A.rank())) def lin_ind_set_check_k2(spanset): print("Input set length: "+str(len(spanset))) deg=find_max_k2(spanset) coeffs=find_coeffs_k2(spanset,deg) A=matrix(QQ,coeffs) print("Rank of input set: "+str(A.rank()))