anchor.diff:      ergrad%pos=0
       allocate(ergrad%f(1:mesh%nv))
       do i=1, mesh%nv
  -       ergrad%f(i)=0._r8
  +      ergrad%f(i)=0._r8
       end do
   
       !Laplacian error
  
positive.code:   end subroutine laplacian_truncation
  
  
  
    subroutine relaxac(mesh)
      !--------------------------------------------------------
  
negative.code:     lap%name="lap_est"
      allocate(lap%f(1:lap%n))
      do i=1, mesh%nv
        lap%f(i)=0._r8
        ergrad%f(i)=0._r8
        do j=1, mesh%v(i)%nnb
          !Edge index
          k=mesh%v(i)%ed(j)
          !Hexagonal edge midpoint
          p=mesh%edhx(k)%c%p
          !ExactGrad=ExactGradVector*NormalVectorEdge*CorrectionOutHx
          gradex=dot_product(df(p),mesh%edhx(k)%nr)*mesh%hx(i)%nr(j)
          !Estimated Gradient Normal component
          gradest=func%f(mesh%v(i)%nb(j))-func%f(i)
          len=mesh%ed(k)%leng
          !arclen(p, mesh%v(i)%p) + &
          !arclen(p, mesh%v(mesh%v(i)%nb(j))%p)
          gradest=gradest/ len !mesh%ed(k)%leng
          !Maximum gradient error for cell
          ergrad%f(i)=max(abs(gradest-gradex), ergrad%f(i))
          !Updade Laplacian
          lap%f(i)=lap%f(i)+gradest*mesh%edhx(k)%leng
        end do
        lap%f(i)=lap%f(i)/mesh%hx(i)%areag
        !Error in laplacian
        error%f(i)=lap_ex%f(i)-lap%f(i)
         !print"(i4, 3f16.8)", i, lap%f(i), lap_ex%f(i), error%f(i)
      end do
  
  