text: static inline void coarsen_quadratic (Point point, scalar s)
  {
    double sum = 0.;
  #if dimension == 1
    foreach_child()
      sum += 9*s[]-s[child.x,0,0];
    s[] = sum/(16);
    
  #elif dimension==2
    foreach_child()
      sum += 10*s[]-s[child.x,0,0]-s[0,child.y,0];
    s[] = sum/(32);
    
  #else //dimension ==3 
    foreach_child()
      sum += 11*s[]-s[child.x,0,0]-s[0,child.y,0]-s[0,0,child.z];
    s[] = sum/(64);
  #endif
  }
  
  static inline void refine_quadratic (Point point, scalar s)
  {
    #if dimension == 1
    foreach_child()
      s[]=(30*coarse(s,0)+5*coarse(s,child.x)-3*coarse(s,-child.x))/32;
    #elif dimension == 2
      foreach_child()
      s[]=(296.*coarse(s,0,0) +5.*coarse(s,child.x,child.y) +
  	 146.*(coarse(s,child.x,0) + coarse(s,0,child.y)) +
  	 98.*(coarse(s,-child.x, 0) + coarse(s,0,-child.y)) -
  	 61.*(coarse(s,child.x,-child.y) + coarse(s,-child.x,child.y))-
  	 91.*(coarse(s,-child.x,-child.y)))/576.;
    #else //dimension == 3
      foreach_child()
        s[]=(412*coarse(s,0,0,0)+
  	   262.*(coarse(s,child.x,0,0)+coarse(s,0,child.y,0)+coarse(s,0,0,child.z))+
  	   121.*(coarse(s,child.x,child.y,0)+coarse(s,child.x,0,child.z)+
  		 coarse(s,0,child.y,child.z))+
  	   214.*(coarse(s,-child.x,0,0)+coarse(s,0,-child.y,0)+coarse(s,0,0,-child.z))-
  	   11.*(coarse(s,child.x,child.y,child.z))+
  	   55.*(coarse(s,child.x,-child.y,0)+coarse(s,-child.x,child.y,0)+
  		coarse(s,child.x,0,-child.z)+coarse(s,-child.x,0,child.z)+
  		coarse(s,0,child.y,-child.z)+coarse(s,0,-child.y,child.z))-
  	   95.*(coarse(s,child.x,child.y,-child.z)+coarse(s,child.x,-child.y,child.z)+
  		coarse(s,-child.x,child.y,child.z))+
  	   25.*(coarse(s,-child.x,0,-child.z)+coarse(s,0,-child.y,-child.z)+
  		coarse(s,-child.x,-child.y,0))-
  	   143.*(coarse(s,-child.x,-child.y,child.z)+coarse(s,-child.x,child.y,-child.z)+
  		 coarse(s,child.x,-child.y,-child.z))-
  	   155.*coarse(s,-child.x,-child.y,-child.z))/1728;
      #endif
  }
  
  static inline void refine_3x3_linear(Point point, scalar s){
    double a = 8.*(s[1,0]+s[1,-1]-s[-1,-1]-s[-1,0]-s[-1,1]+s[1,1]);
    double b = 8.*(-s[1,-1]-s[0,-1]-s[-1,-1]+s[-1,1]+s[0,1]+s[1,1]);
    double c = 3.*(-s[1,-1]+s[-1,-1]-s[-1,1]+s[1,1]);
    double d = s[];
    foreach_child()
      s[]=d+(a*child.x+b*child.y+c*(child.x*child.y))/192.;
  }
  
  
  # if dimension == 3
  static inline void refine_3x3x3_linear(Point point, scalar s){
    double a = 64.*(s[1,0,0]+s[1,-1,0]-s[-1,-1,0]-s[-1,0,0]-s[-1,1,0]+s[1,1,0]
  		  +s[1,0,-1]+s[1,-1,-1]-s[-1,-1,-1]-s[-1,0,-1]-s[-1,1,-1]+s[1,1,-1]
  		  +s[1,0,1]+s[1,-1,1]-s[-1,-1,1]-s[-1,0,1]-s[-1,1,1]+s[1,1,1]);
    double b = 64.*(-s[1,-1,0]-s[0,-1,0]-s[-1,-1,0]+s[-1,1,0]+s[0,1,0]+s[1,1,0]
  		  -s[1,-1,-1]-s[0,-1,-1]-s[-1,-1,-1]+s[-1,1,-1]+s[0,1,-1]+s[1,1,-1]
  		  -s[1,-1,1]-s[0,-1,1]-s[-1,-1,1]+s[-1,1,1]+s[0,1,1]+s[1,1,1]);
    double c = 64.*(-s[0,0,-1]-s[1,1,-1]-s[-1,-1,-1]-s[1,0,-1]-s[-1,0,-1]-s[0,1,-1]
  		  -s[0,-1,-1]-s[1,-1,-1]-s[-1,1,-1]+s[0,0,1]+s[1,1,1]+s[-1,-1,1]
  		  +s[1,0,1]+s[-1,0,1]+s[0,1,1]+s[0,-1,1]+s[1,-1,1]+s[-1,1,1]);
    double d = 24.*(-s[1,-1,0]+s[-1,-1,0]-s[-1,1,0]+s[1,1,0]-s[1,-1,-1]+s[-1,-1,-1]
  		  -s[-1,1,-1]+s[1,1,-1]-s[1,-1,1]+s[-1,-1,1]-s[-1,1,1]+s[1,1,1]) ; 
    double e = 24.*(s[1,-1,-1]+s[0,-1,-1]+s[-1,-1,-1]-s[-1,1,-1]-s[0,1,-1]-s[1,1,-1]
  		  -s[1,-1,1]-s[0,-1,1]-s[-1,-1,1]+s[-1,1,1]+s[0,1,1]+s[1,1,1]) ; 
    double f = 24.*(-s[1,0,-1]-s[1,-1,-1]+s[-1,-1,-1]+s[-1,0,-1]+s[-1,1,-1]-s[1,1,-1]
  		  +s[1,0,1]+s[1,-1,1]-s[-1,-1,1]-s[-1,0,1]-s[-1,1,1]+s[1,1,1]) ; 
    double g = 9.*(s[1,-1,-1]-s[-1,-1,-1]+s[-1,1,-1]-s[1,1,-1]-s[1,-1,1]+s[-1,-1,1]-s[-1,1,1]+s[1,1,1]); 
    double h = s[];
    foreach_child()
      s[]=h+(a*child.x+b*child.y+c*child.z+d*(child.x*child.y)+e*(child.y*child.z)+f*(child.x*child.z)+g*(child.x*child.y*child.z))/4608.;
  }
  static inline void refine_quad(Point point, scalar s){
    double a = 3.*448.*(s[1,0,0]+s[1,-1,0]-s[-1,-1,0]-s[-1,0,0]-s[-1,1,0]+s[1,1,0]
  		  +s[1,0,-1]+s[1,-1,-1]-s[-1,-1,-1]-s[-1,0,-1]-s[-1,1,-1]+s[1,1,-1]
  		  +s[1,0,1]+s[1,-1,1]-s[-1,-1,1]-s[-1,0,1]-s[-1,1,1]+s[1,1,1]);
    double b = 3.*448.*(-s[1,-1,0]-s[0,-1,0]-s[-1,-1,0]+s[-1,1,0]+s[0,1,0]+s[1,1,0]
  		  -s[1,-1,-1]-s[0,-1,-1]-s[-1,-1,-1]+s[-1,1,-1]+s[0,1,-1]+s[1,1,-1]
  		  -s[1,-1,1]-s[0,-1,1]-s[-1,-1,1]+s[-1,1,1]+s[0,1,1]+s[1,1,1]);
    double c = 3.*448.*(-s[0,0,-1]-s[1,1,-1]-s[-1,-1,-1]-s[1,0,-1]-s[-1,0,-1]-s[0,1,-1]
  		  -s[0,-1,-1]-s[1,-1,-1]-s[-1,1,-1]+s[0,0,1]+s[1,1,1]+s[-1,-1,1]
  		  +s[1,0,1]+s[-1,0,1]+s[0,1,1]+s[0,-1,1]+s[1,-1,1]+s[-1,1,1]);
    double d = 3.*168.*(-s[1,-1,0]+s[-1,-1,0]-s[-1,1,0]+s[1,1,0]-s[1,-1,-1]+s[-1,-1,-1]
  		  -s[-1,1,-1]+s[1,1,-1]-s[1,-1,1]+s[-1,-1,1]-s[-1,1,1]+s[1,1,1]) ; 
    double e = 3.*168.*(s[1,-1,-1]+s[0,-1,-1]+s[-1,-1,-1]-s[-1,1,-1]-s[0,1,-1]-s[1,1,-1]
  		  -s[1,-1,1]-s[0,-1,1]-s[-1,-1,1]+s[-1,1,1]+s[0,1,1]+s[1,1,1]) ; 
    double f = 3.*168.*(-s[1,0,-1]-s[1,-1,-1]+s[-1,-1,-1]+s[-1,0,-1]+s[-1,1,-1]-s[1,1,-1]
  		   +s[1,0,1]+s[1,-1,1]-s[-1,-1,1]-s[-1,0,1]-s[-1,1,1]+s[1,1,1]) ; 
    double g = 3.*63.*(s[1,-1,-1]-s[-1,-1,-1]+s[-1,1,-1]-s[1,1,-1]-s[1,-1,1]+s[-1,-1,1]-s[-1,1,1]+s[1,1,1]);
    double h = 3.*(240.*(s[1,0,0]+s[-1,0,0])
  	      +144.*(s[1,-1,0]+s[-1,-1,0]+s[1,1,0]+s[1,0,-1]+s[-1,0,-1]+s[1,0,1]+s[-1,0,1])
  	      -96.*(s[0,-1,0]+s[0,1,0]+s[0,0,-1]+s[0,0,1])
  	      +48.*(s[1,-1,-1]+s[-1,-1,-1]+s[-1,1,-1]+s[1,1,-1]+s[1,-1,1]+s[-1,-1,1]+s[-1,1,1]+s[1,1,1])
  	      -192.*(s[0,-1,-1]+s[0,1,-1]+s[0,-1,1]+s[0,1,1]));
    double i = 0.;
    double j = 0. ;
    double k = 3584.*(7.*s[]
  		    +4.*(s[1,0,0]+s[-1,0,0]+s[0,1,0]+s[0,-1,0]+s[0,0,1]+s[0,0,-1])
  		    +(s[1,1,0]+s[-1,1,0]+s[1,-1,0]+s[-1,-1,0]
  		      +s[1,0,1]+s[-1,0,1]+s[1,0,-1]+s[-1,0,-1]
  		      +s[0,1,1]+s[0,-1,1]+s[0,1,-1]+s[0,-1,-1])
  		    -2.*(s[1,1,1]+s[-1,1,1]+s[1,-1,1]+s[1,1,-1]
  			 +s[-1,-1,1]+s[-1,1,-1]+s[1,-1,-1]+s[-1,-1,-1]));
  						    
     
    foreach_child()
      s[]=(k+a*child.x+b*child.y+c*child.z+d*(child.x*child.y)+e*(child.y*child.z)
  	 +f*(child.x*child.z)+g*(child.x*child.y*child.z)+h*(child.x*child.x)+
  	 i*(child.y*child.y)+j*(child.z*child.z))/(32256.*3);
  }
  
  # endif
  static inline void coarsen1D_twoleaf_quad (Point point, scalar s){
    double sum = 0.;
    foreach_child()
      sum += 16*s[];
    sum+=-(s[-1]+s[1]);
    s[] = sum/(30);
  }
  
  static inline void coarsen1D_leftleaf_quad (Point point, scalar s){
    double sum = 0.;
    foreach_child(){
      if (child.x==-1)
        sum += 591.2*s[];
      else
        sum+=(((490.5*s[])+(-32.3*s[child.x])));
    }
    sum+=-49.4*(s[-1]);
    s[] = sum/(1000);
  }
  
  static inline void coarsen1D_rightleaf_quad (Point point, scalar s){
    double sum = 0.;
    foreach_child(){
      if (child.x==1)
        sum += 591.2*s[];
      else
        sum+=(((490.5*s[])+(-32.3*s[child.x])));
    }
    sum+=-49.4*(s[1]);
    s[] = sum/(1000);
    }
  
  static inline void coarsen2D_oneleaf_quad(Point point,scalar s,int n)
  {
    double w[6]= {0.321645313553608, 0.282670262980445, -0.010788941335132,
  		-0.052461227242077, -0.013486176668914,-0.055158462575860};
    if (n == 1){
      s[]=(w[0]*(fine(s,0,0)+fine(s,0,1))+
  	 w[1]*(fine(s,1,0)+fine(s,1,1))+
  	 w[2]*(fine(s,2,1)+fine(s,2,0))+
  	 w[3]*(fine(s,1,2)+fine(s,1,-1))+
  	 w[4]*(fine(s,0,2)+fine(s,0,-1))+
  	 w[5]*s[-1]);
    }
    else if (n == 2){
      s[]=(w[0]*(fine(s,0,0)+fine(s,1,0))+
  	 w[1]*(fine(s,0,1)+fine(s,1,1))+
  	 w[2]*(fine(s,1,2)+fine(s,0,2))+
  	 w[3]*(fine(s,2,1)+fine(s,-1,1))+
  	 w[4]*(fine(s,2,0)+fine(s,-1,0))+
  	 w[5]*s[0,-1]);
    }
    else if (n == 4)
      s[]=(w[0]*(fine(s,1,1)+fine(s,1,0))+
  	 w[1]*(fine(s,0,0)+fine(s,0,1))+
  	 w[2]*(fine(s,-1,0)+fine(s,-1,1))+
  	 w[3]*(fine(s,0,-1)+fine(s,0,2))+
  	 w[4]*(fine(s,1,-1)+fine(s,1,2))+
  	 w[5]*s[1]);
    else // (n == 8)
      s[]=(w[0]*(fine(s,1,1)+fine(s,0,1))+
  	 w[1]*(fine(s,0,0)+fine(s,1,0))+
  	 w[2]*(fine(s,0,-1)+fine(s,1,-1))+
  	 w[3]*(fine(s,-1,0)+fine(s,2,0))+
  	 w[4]*(fine(s,-1,1)+fine(s,2,1))+
  	 w[5]*s[0,1]);
  }
  
  static inline void coarsen2D_twoleaf_quad(Point point,scalar s,int n)
  {
    if (n==5 || n==10)
      {
        if (n==5)
  	s[]=(9*(fine(s,0,0)+fine(s,0,1)+
  		fine(s,1,1)+fine(s,1,0))-
  	     (s[-1]+s[1]+
  	      fine(s,1,2)+fine(s,0,2)+
  	      fine(s,1,-1)+fine(s,0,-1)))/30;
        
        else if (n == 10)
  	s[]=(9*(fine(s,0,0)+fine(s,0,1)+
  		fine(s,1,1)+fine(s,1,0))-
  	     (s[0,-1]+s[0,1]+
  	      fine(s,2,1)+fine(s,2,0)+
  	      fine(s,-1,1)+fine(s,-1,0)))/30;
      }
    else
      {
        double w[6]={0.343774069319640, 0.287034659820282, 0.256482670089859,
  		   -0.016174582798460, -0.020539152759949, -0.050449293966624};
        if (n == 3)
  	s[]=(w[0]*fine(s,0,0)+
  	     w[1]*(fine(s,1,0)+fine(s,0,1))+
  	     w[2]*fine(s,1,1)+
  	     w[3]*(fine(s,0,2)+fine(s,2,0))+
  	     w[4]*(fine(s,1,2)+fine(s,2,1))+
  	     w[5]*(s[-1]+s[0,-1]));
        else if (n == 6)
  	s[]=(w[0]*fine(s,1,0)+
  	     w[1]*(fine(s,1,1)+fine(s,0,0))+
  	     w[2]*fine(s,0,1)+
  	     w[3]*(fine(s,1,2)+fine(s,-1,0))+
  	     w[4]*(fine(s,0,2)+fine(s,-1,1))+
  	     w[5]*(s[1]+s[0,-1]));
        else if (n== 9)
  	s[]=(w[0]*fine(s,0,1)+
  	     w[1]*(fine(s,1,1)+fine(s,0,0))+
  	     w[2]*fine(s,1,0)+
  	     w[3]*(fine(s,0,-1)+fine(s,2,1))+
  	     w[4]*(fine(s,1,-1)+fine(s,2,0))+
  	     w[5]*(s[-1]+s[0,1]));
        else // (n==12)
  	s[]=(w[0]*fine(s,1,1)+
  	     w[1]*(fine(s,0,1)+fine(s,1,0))+
  	     w[2]*fine(s,0,0)+
  	     w[3]*(fine(s,1,-1)+fine(s,-1,1))+
  	     w[4]*(fine(s,0,-1)+fine(s,-1,0))+
  	     w[5]*(s[1]+s[0,1]));
  	}
  }
  
  static inline void coarsen2D_threeleaf_quad(Point point,scalar s,int n)
  {
    double w[5]= {0.312119914346895, 0.2674946466809427,
  		-0.050963597430407, -0.035032119914347, -0.01910064239828};
    if (n == 7)
      s[]=(w[0]*(fine(s,0,0)+fine(s,1,0))+
  	 w[1]*(fine(s,0,1)+fine(s,1,1))+
  	 w[3]*(s[-1]+s[1])+w[2]*s[0,-1]+
  	 w[4]*(fine(s,0,2)+fine(s,1,2)));
    else if (n == 11)
      s[]=(w[0]*(fine(s,0,0)+fine(s,0,1))+
  	 w[1]*(fine(s,1,1)+fine(s,1,0))+
  	 w[3]*(s[0,-1]+s[0,1])+w[2]*s[-1]+
  	 w[4]*(fine(s,2,0)+fine(s,2,1)));
    else if (n == 13)
      s[]=(w[0]*(fine(s,0,1)+fine(s,1,1))+
  	 w[1]*(fine(s,1,0)+fine(s,0,0))+
  	 w[3]*(s[-1]+s[1])+w[2]*s[0,1]+
  	 w[4]*(fine(s,1,-1)+fine(s,0,-1)));
    else // (n == 14)
      s[]=(w[0]*(fine(s,1,1)+fine(s,1,0))+
  	 w[1]*(fine(s,0,0)+fine(s,0,1))+
  	 w[3]*(s[0,-1]+s[0,1])+w[2]*s[1]+
  	 w[4]*(fine(s,-1,1)+fine(s,-1,0)));
  }
  
  
  static inline void coarsen2D_fourleaf_quad(Point point, scalar s){
    double sum = 0;
    foreach_child()
      sum+=8*s[];
    sum-=(s[-1]+s[1]+s[0,1]+s[0,-1]);
    s[] = sum/28; 
  }
  
  static inline void coarsen3D_oneleaf_quad(Point point, scalar s, int n){
  }
  static inline void coarsen3D_twoleaf_quad(Point point, scalar s, int n){
  }
  static inline void coarsen3D_threeleaf_quad(Point point, scalar s, int n){
  }
  static inline void coarsen3D_fourleaf_quad(Point point, scalar s, int n){
  }
  
  static inline void coarsen3D_fiveleaf_quad(Point point, scalar s,int n){
    double w[5]={0.181093394077449, 0.119589977220957, 0.018223234624146,
  	       -0.039863325740319, -0.116173120728929};
    fprintf(ferr,"Lets use this variable %g\n",w[2]);
  }
  
  static inline void coarsen3D_sixleaf_quad(Point point, scalar s){
    double sum = 0;
    foreach_child()
      sum+=4*s[];
    sum-=(s[-1]+s[1]+s[0,1]+s[0,-1]+s[0,0,1]+s[0,0,-1]);
    s[] = sum/28; 
  }
  
  static inline int caser(Point point)
  {
    int n=0;
    int j=0;
    for (int i = -1;i<=1;i+=2){
      if (is_leaf(neighbor(i,0,0)))
        n+=pow(2,j);
      j++;
  #if dimension > 1 
      if (is_leaf(neighbor(0,i,0)))
        n+=pow(2,j);
      j++;
  #endif
  #if dimension > 2 
      if (is_leaf(neighbor(0,0,i)))
        n+=pow(2,j);
      j++;
  #endif
    }
    return n;
  }
  
  
  static inline void quadratic_coarsening(Point point, scalar s)
  {
    int n = caser(point);
    if (n==0)
      coarsen_quadratic(point, s); 
  #if dimension == 1
    else if (n==1) 
      coarsen1D_leftleaf_quad(point,s);
    else if (n==2)
      coarsen1D_rightleaf_quad(point,s);
    else if (n==3)
      coarsen1D_twoleaf_quad(point,s);
    
  #elif dimension == 2
    else if (n==1||n==2||n==4||n==8)
      coarsen2D_oneleaf_quad(point,s,n);
    else if (n==3||n==5||n==6||n==9||n==10||n==12)  
      coarsen2D_twoleaf_quad(point,s,n); 
    else if(n==7||n==11||n==13||n==14) 
      coarsen2D_threeleaf_quad(point,s,n);
    else if (n==15)
      coarsen2D_fourleaf_quad(point,s);
    
  #else  // dimension == 2
    else if (n==1||n==2||n==4||n==8||n==16||n==32)
      coarsen3D_oneleaf_quad(point,s,n);
    else if (n==3||n==5||n==9||n==17||n==33||n==6||n==10||n==18||
  	   n==34||n==12||n==20||n==36||n==24||n==40||n==48)
      coarsen3D_twoleaf_quad(point,s,n);
    else if(n==7||n==11||n==19||n==35||n==13||n==21||n==37||
  	  n==25||n==41||n==49||n==14||n==22||n==38||n==26||
  	  n==42||n==50||n==28||n==44||n==56)
      coarsen3D_threeleaf_quad(point,s,n);
    else if(n==15||n==23||n==27||n==29||n==30||n==39||n==43||
  	  n==45||n==46||n==51||n==53||n==54||n==57||n==58||n==60) 
      coarsen3D_fourleaf_quad(point,s,n);
    else if (n==31||n==47||n==55||n==59||n==61||n==62)
      coarsen3D_fiveleaf_quad(point,s,n);
    else if (n==63)
      coarsen3D_sixleaf_quad(point,s);
  #endif
  }
  
  static inline void quadratic_boundary_restriction(Point point, scalar s)
  {
    int n = caser(point);
    if (n==0)
      restriction_average(point, s); 
  #if dimension == 1
    else if (n==1) 
      coarsen1D_leftleaf_quad(point,s);
    else if (n==2)
      coarsen1D_rightleaf_quad(point,s);
    else if (n==3)
      coarsen1D_twoleaf_quad(point,s);
    
  #elif dimension == 2
    else if (n==1||n==2||n==4||n==8)
      coarsen2D_oneleaf_quad(point,s,n);
    else if (n==3||n==5||n==6||n==9||n==10||n==12)  
      coarsen2D_twoleaf_quad(point,s,n); 
    else if(n==7||n==11||n==13||n==14) 
      coarsen2D_threeleaf_quad(point,s,n);
    else if (n==15)
      coarsen2D_fourleaf_quad(point,s);
  #else 
    else if (n==1||n==2||n==4||n==8||n==16||n==32)
      coarsen3D_oneleaf_quad(point,s,n);
    else if (n==3||n==5||n==9||n==17||n==33||n==6||n==10||n==18||
  	   n==34||n==12||n==20||n==36||n==24||n==40||n==48)
      coarsen3D_twoleaf_quad(point,s,n);
    else if(n==7||n==11||n==19||n==35||n==13||n==21||n==37||
  	  n==25||n==41||n==49||n==14||n==22||n==38||n==26||
  	  n==42||n==50||n==28||n==44||n==56)
      coarsen3D_threeleaf_quad(point,s,n);
    else if(n==15||n==23||n==27||n==29||n==30||n==39||n==43||n==45||
  	  n==46||n==51||n==53||n==54||n==57||n==58||n==60) 
      coarsen3D_fourleaf_quad(point,s,n);
    else if (n==31||n==47||n==55||n==59||n==61||n==62)
      coarsen3D_fiveleaf_quad(point,s,n);
    else if (n==63)
      coarsen3D_sixleaf_quad(point,s);
  #endif
  }
  