neumannEquationBC.h

Go to the documentation of this file.

c File generated by ogmg/neumann.maple
c Fourth-order BC for mixed boundary condition using the normal derivative of the PDE 
c BC:    an1*ur+an2*us+an0*u = g 
c PDE:   c11*urr + c12*urs + c22*uss + c1*ur + c2*us + c0 *u = f 
c PDE.r: urrr = b0*u + b1*us + b2*uss + b3*usss + bf  

c real alpha1,alpha2,a1,a2,a0
c real rxi,ryi,sxi,syi,rxr,rxs,sxr,sxs,ryr,rys,syr,sys
c real rxx,ryy,sxx,syy
c real rxrr,rxrs,rxss,ryrr,ryrs,ryss
c real sxrr,sxrs,sxss,syrr,syrs,syss
c real rxx,ryy,sxx,syy
c real rxxr,ryyr,rxxs,ryys, sxxr,syyr,sxxs,syys
c real rxNormI,rxNormIs,rxNormIss,rxNormIr,rxNormIrr
c real n1,n1s,n1ss,n2,n2s,n2ss
c real an1,an1s,an1ss,an2,an2s,an2ss 
c real an1r,an1rr,an2r,an2rs 
c real ff,ffs,ffr,g,gs,gss 
c real c11,c11r,c11s,c12,c12r,c12s,c22,c22r,c22s,c1,c1r,c1s,c2,c2r,c2s,c0,c0r,c0s
c real b0,b1,b2,b3,bf

#beginMacro fourthOrderNeumannEquationBC(DIR,DIM)
 #If #DIM == "3" 
   stop 45
 #End
 rxi=rx(i1,i2,i3)
 ryi=ry(i1,i2,i3)
 sxi=sx(i1,i2,i3)
 syi=sy(i1,i2,i3)

 rxr=rxr2(i1,i2,i3)
 rxs=rxs2(i1,i2,i3)
 ryr=ryr2(i1,i2,i3)
 rys=rys2(i1,i2,i3)

 sxr=sxr2(i1,i2,i3)
 sxs=sxs2(i1,i2,i3)
 syr=syr2(i1,i2,i3)
 sys=sys2(i1,i2,i3)

 rxx=rxx22(i1,i2,i3)
 ryy=ryy22(i1,i2,i3)
 rxrr=rxrr2(i1,i2,i3)
 rxrs=rxrs2(i1,i2,i3)
 rxss=rxss2(i1,i2,i3)
 ryrr=ryrr2(i1,i2,i3)
 ryrs=ryrs2(i1,i2,i3)
 ryss=ryss2(i1,i2,i3)

 sxx=sxx22(i1,i2,i3)
 syy=syy22(i1,i2,i3)
 sxrr=sxrr2(i1,i2,i3)
 sxrs=sxrs2(i1,i2,i3)
 sxss=sxss2(i1,i2,i3)
 syrr=syrr2(i1,i2,i3)
 syrs=syrs2(i1,i2,i3)
 syss=syss2(i1,i2,i3)

 rxxr=rxi*rxrr+rxr*rxr + sxi*rxrs + sxr*rxs
 ryyr=ryi*ryrr+ryr*ryr + syi*ryrs + syr*rys

 rxxs=rxi*rxrs+rxs*rxr + sxi*rxss + sxs*rxs
 ryys=ryi*ryrs+rys*ryr + syi*ryss + sys*rys

 sxxr=rxi*sxrr+rxr*sxr + sxi*sxrs + sxr*sxs
 syyr=ryi*syrr+ryr*syr + syi*syrs + syr*sys

 sxxs=rxi*sxrs+rxs*sxr + sxi*sxss + sxs*sxs
 syys=ryi*syrs+rys*syr + syi*syss + sys*sys

 alpha1=a1*nsign  ! nsign=2*side-1
 alpha2=a1*nsign

 #If #DIR == "R"
  rxNormI=1./sqrt(rxi**2+ryi**2)
  rxNormIs=-(rxi*rxs+ryi*rys)*rxNormI**3
  rxNormIss=-(rxi*rxss+ryi*ryss+rxs*rxs+rys*rys)*rxNormI**3 -3.*(rxi*rxs+ryi*rys)*rxNormI**2*rxNormIs

  n1=rxi*rxNormI
  n1s=rxs*rxNormI + rxi*rxNormIs
  n1ss=rxss*rxNormI + 2.*rxs*rxNormIs + rxi*rxNormIss
  n2=ryi*rxNormI
  n2s=rys*rxNormI + ryi*rxNormIs
  n2ss=ryss*rxNormI + 2.*rys*rxNormIs + ryi*rxNormIss
  an1=alpha1*(n1*rxi+n2*ryi)
  an2=alpha2*(n1*sxi+n2*syi)
  an1s=alpha1*(n1*rxs+n2*rys+n1s*rxi+n2s*ryi)
  an1ss=alpha1*(n1*rxss+n2*ryss+2.*(n1s*rxs+n2s*rys)+n1ss*rxi+n2ss*ryi)
  an2s=alpha2*(n1*sxs+n2*sys + n1s*sxi+n2s*syi)
  an2ss=alpha2*(n1*sxss+n2*syss + 2.*(n1s*sxs+n2s*sys) + n1ss*sxi+n2ss*syi)
 #Else
  sxNormI=1./sqrt(sxi**2+syi**2)
  sxNormIr=-(sxi*sxr+syi*syr)*sxNormI**3
  sxNormIrr=-(sxi*sxrr+syi*syrr+sxr*sxr+syr*syr)*sxNormI**3 -3.*(sxi*sxr+syi*syr)*sxNormI**2*sxNormIr

  n1=sxi*sxNormI
  n1r=sxr*sxNormI + sxi*sxNormIr
  n1rr=sxrr*sxNormI + 2.*sxr*sxNormIr + sxi*sxNormIrr
  n2=syi*sxNormI
  n2r=syr*sxNormI + syi*sxNormIr
  n2rr=syrr*sxNormI + 2.*syr*sxNormIr + syi*sxNormIrr
  an1=alpha1*(n1*rxi+n2*ryi)
  an2=alpha2*(n1*sxi+n2*syi)
  an1r=alpha1*(n1*rxr+n1r*rxi + n2*ryr+n2r*ryi)
  an1rr=alpha1*(n1*rxrr+n2*ryrr+ 2.*(n1r*rxr+n2r*ryr) + n1rr*rxi+n2rr*ryi)
  an2r=alpha2*(n1*sxr+n1r*sxi + n2*syr+n2r*syi)
  an2rr=alpha2*(n1*sxrr+n2*syrr+ 2.*(n1r*sxr+n2r*syr) + n1rr*sxi+n2rr*syi)
 #End

 c11=rxi**2+ryi**2
 c11r=2.*(rxi*rxr+ryi*ryr)
 c11s=2.*(rxi*rxs+ryi*rys)
 c12=2.*(rxi*sxi+ryi*syi)
 c12r=2.*(rxr*sxi+rxi*sxr + ryr*syi+ ryi*syr)
 c12s=2.*(rxs*sxi+rxi*sxs + rys*syi+ ryi*sys)
 c22=sxi**2+syi**2
 c22r=2.*(sxi*sxr+syi*syr)
 c22s=2.*(sxi*sxs+syi*sys)
 c1=rxx+ryy
 c1r=rxxr+ryyr
 c1s=rxxs+ryys
 c2=sxx+syy
 c2r=sxxr+syyr
 c2s=sxxs+syys
 c0=0.
 c0r=0.
 c0s=0.

c Coefficients of u(i1-2,i2,i3) u(i1-1,i2,i3) u(i1,i2,i3) u(i1+1,i2,i3) u(i1+2,i2,i3)
c um2=-1/2/dr**3
c um1=1/dr**3
c u0 =-1/ds**2*(b0*ds**2-2*b2)
c up1=-1/dr**3
c up2=1/2/dr**3
 #If #DIR == "R"
 b0=(-c12*c1s*an1**2*c11*a0-c12*c11s*an1**3*c0-c11s*an1*c12**2*an1s*a0+c12**2*c11*an1ss*an1*a0+c1*an1**3*c11*c0-c0r*an1**3*c11**2+c12*c11s*an1**2*c1*a0-c22*c11**2*an1ss*an1*a0+c12*c12s*an1*c11*an1s*a0+2*c22*c11**2*an1s**2*a0-c2*an1*c11**2*an1s*a0+c12*c0s*an1**3*c11+c0*an1**2*c11**2*a0+c11r*an1*c11*c12*an1s*a0-2*c12**2*c11*an1s**2*a0-c1**2*an1**2*c11*a0+2*c12*c1*an1*c11*an1s*a0-c11r*an1**2*c11*c1*a0-c12r*an1*c11**2*an1s*a0+c11r*an1**3*c11*c0+c1r*an1**2*c11**2*a0)/an1**3/c11**3

 b1=(c22*c11**2*an2ss*an1**2-c12*c11s*an1**3*c2+c12*c2s*an1**3*c11+c12*c0*an1**3*c11+c11s*an1**2*c12**2*a0+c11s*an1**2*c12**2*an2s-2*c12**2*c11*an1s**2*an2-c12**2*c11*an2ss*an1**2+c11r*an1**3*c11*c2-c12*c12s*an1**2*c11*a0-c12*c12s*an1**2*c11*an2s-2*c12*c1*an1**2*c11*a0-2*c12*c1*an1**2*c11*an2s+c11r*an1*c11*c12*an1s*an2+c12**2*c11*an1ss*an1*an2+2*c12**2*c11*an1s*a0*an1+2*c12**2*c11*an1s*an2s*an1-c11r*an1**2*c11*c1*an2-c11r*an1**2*c11*c12*an2s-c1**2*an1**2*c11*an2+2*c22*c11**2*an1s**2*an2+c12*c12s*an1*c11*an1s*an2-c12r*an1*c11**2*an1s*an2+2*c12*c1*an1*c11*an1s*an2-c12*c1s*an1**2*c11*an2-c2*an1*c11**2*an1s*an2-c22*c11**2*an1ss*an1*an2-2*c22*c11**2*an1s*a0*an1-2*c22*c11**2*an1s*an2s*an1+c12*c11s*an1**2*c1*an2-c11s*an1*c12**2*an1s*an2+c0*an1**2*c11**2*an2+c2*an1**2*c11**2*a0+c2*an1**2*c11**2*an2s+c1*an1**3*c11*c2+c1r*an1**2*c11**2*an2+c12r*an1**2*c11**2*a0+c12r*an1**2*c11**2*an2s-c11r*an1**2*c11*c12*a0-c2r*an1**3*c11**2)/an1**3/c11**3

 b2=(-c11r*an1**2*c11*c12*an2-2*c12**2*c11*an2s*an1**2-c12*c11s*an1**3*c22-2*c12*c1*an1**2*c11*an2+c12*c2*an1**3*c11-c12*c12s*an1**2*c11*an2+c11s*an1**2*c12**2*an2-c22r*an1**3*c11**2+2*c22*c11**2*an2s*an1**2+c22*c11**2*a0*an1**2+c12*c22s*an1**3*c11-2*c22*c11**2*an1s*an2*an1+c11r*an1**3*c11*c22+c1*an1**3*c11*c22+c2*an1**2*c11**2*an2-c12**2*c11*a0*an1**2+c12r*an1**2*c11**2*an2+2*c12**2*c11*an1s*an2*an1)/an1**3/c11**3

 b3=(c12*c22*an1**3*c11-c12**2*c11*an2*an1**2+c22*c11**2*an2*an1**2)/an1**3/c11**3

 bf=(-c11r*an1*c11*c12*an1s+c11r*an1**2*c11*c1+2*c12**2*c11*an1s**2-c12**2*c11*an1ss*an1+c11s*an1*c12**2*an1s-2*c12*c1*an1*c11*an1s-c12*c12s*an1*c11*an1s+c12*c1s*an1**2*c11-c12*c11s*an1**2*c1-2*c22*c11**2*an1s**2+c22*c11**2*an1ss*an1+c1**2*an1**2*c11+c2*an1*c11**2*an1s-c0*an1**2*c11**2+c12r*an1*c11**2*an1s-c1r*an1**2*c11**2)/an1**3/c11**3*g+(-c11s*an1**2*c12**2-c2*an1**2*c11**2+2*c12*c1*an1**2*c11+2*c22*c11**2*an1s*an1+c11r*an1**2*c11*c12-c12r*an1**2*c11**2-2*c12**2*c11*an1s*an1+c12*c12s*an1**2*c11)/an1**3/c11**3*gs+(-c22*c11**2*an1**2+c12**2*c11*an1**2)/an1**3/c11**3*gss-c12/c11**2*ffs+1/c11*ffr+(-c11r*an1**3*c11*ff-c1*an1**3*c11*ff+c12*c11s*an1**3*ff)/an1**3/c11**3

 #Else
 b0=(-c12*c2r*an2**2*c22*a0-c12*c22r*an2**3*c0-c22r*an2*c12**2*an2r*a0+c12**2*c22*an2rr*an2*a0+c2*an2**3*c22*c0-c0s*an2**3*c22**2+c12*c22r*an2**2*c2*a0-c11*c22**2*an2rr*an2*a0+c12*c12r*an2*c22*an2r*a0+2*c11*c22**2*an2r**2*a0-c1*an2*c22**2*an2r*a0+c12*c0r*an2**3*c22+c0*an2**2*c22**2*a0+c22s*an2*c22*c12*an2r*a0-2*c12**2*c22*an2r**2*a0-c2**2*an2**2*c22*a0+2*c12*c2*an2*c22*an2r*a0-c22s*an2**2*c22*c2*a0-c12s*an2*c22**2*an2r*a0+c22s*an2**3*c22*c0+c2s*an2**2*c22**2*a0)/an2**3/c22**3

 b1=(c11*c22**2*an1rr*an2**2-c12*c22r*an2**3*c1+c12*c1r*an2**3*c22+c12*c0*an2**3*c22+c22r*an2**2*c12**2*a0+c22r*an2**2*c12**2*an1r-2*c12**2*c22*an2r**2*an1-c12**2*c22*an1rr*an2**2+c22s*an2**3*c22*c1-c12*c12r*an2**2*c22*a0-c12*c12r*an2**2*c22*an1r-2*c12*c2*an2**2*c22*a0-2*c12*c2*an2**2*c22*an1r+c22s*an2*c22*c12*an2r*an1+c12**2*c22*an2rr*an2*an1+2*c12**2*c22*an2r*a0*an2+2*c12**2*c22*an2r*an1r*an2-c22s*an2**2*c22*c2*an1-c22s*an2**2*c22*c12*an1r-c2**2*an2**2*c22*an1+2*c11*c22**2*an2r**2*an1+c12*c12r*an2*c22*an2r*an1-c12s*an2*c22**2*an2r*an1+2*c12*c2*an2*c22*an2r*an1-c12*c2r*an2**2*c22*an1-c1*an2*c22**2*an2r*an1-c11*c22**2*an2rr*an2*an1-2*c11*c22**2*an2r*a0*an2-2*c11*c22**2*an2r*an1r*an2+c12*c22r*an2**2*c2*an1-c22r*an2*c12**2*an2r*an1+c0*an2**2*c22**2*an1+c1*an2**2*c22**2*a0+c1*an2**2*c22**2*an1r+c2*an2**3*c22*c1+c2s*an2**2*c22**2*an1+c12s*an2**2*c22**2*a0+c12s*an2**2*c22**2*an1r-c22s*an2**2*c22*c12*a0-c1s*an2**3*c22**2)/an2**3/c22**3

 b2=(-c22s*an2**2*c22*c12*an1-2*c12**2*c22*an1r*an2**2-c12*c22r*an2**3*c11-2*c12*c2*an2**2*c22*an1+c12*c1*an2**3*c22-c12*c12r*an2**2*c22*an1+c22r*an2**2*c12**2*an1-c11s*an2**3*c22**2+2*c11*c22**2*an1r*an2**2+c11*c22**2*a0*an2**2+c12*c11r*an2**3*c22-2*c11*c22**2*an2r*an1*an2+c22s*an2**3*c22*c11+c2*an2**3*c22*c11+c1*an2**2*c22**2*an1-c12**2*c22*a0*an2**2+c12s*an2**2*c22**2*an1+2*c12**2*c22*an2r*an1*an2)/an2**3/c22**3

 b3=(c12*c11*an2**3*c22-c12**2*c22*an1*an2**2+c11*c22**2*an1*an2**2)/an2**3/c22**3

 bf=(-c22s*an2*c22*c12*an2r+c22s*an2**2*c22*c2+2*c12**2*c22*an2r**2-c12**2*c22*an2rr*an2+c22r*an2*c12**2*an2r-2*c12*c2*an2*c22*an2r-c12*c12r*an2*c22*an2r+c12*c2r*an2**2*c22-c12*c22r*an2**2*c2-2*c11*c22**2*an2r**2+c11*c22**2*an2rr*an2+c2**2*an2**2*c22+c1*an2*c22**2*an2r-c0*an2**2*c22**2+c12s*an2*c22**2*an2r-c2s*an2**2*c22**2)/an2**3/c22**3*g+(-c22r*an2**2*c12**2-c1*an2**2*c22**2+2*c12*c2*an2**2*c22+2*c11*c22**2*an2r*an2+c22s*an2**2*c22*c12-c12s*an2**2*c22**2-2*c12**2*c22*an2r*an2+c12*c12r*an2**2*c22)/an2**3/c22**3*gr+(-c11*c22**2*an2**2+c12**2*c22*an2**2)/an2**3/c22**3*grr-c12/c22**2*ffr+1/c22*ffs+(-c22s*an2**3*c22*ff-c2*an2**3*c22*ff+c12*c22r*an2**3*ff)/an2**3/c22**3

 #End

#endMacro