     cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), values.size());
     A.row_indices = row_indices; A.column_indices = column_indices; A.values = values;
 
-    if (globalbcz == DIR_NEU ) std::cerr << "DIR_NEU NOT IMPLEMENTED "<<std::endl;
-    if (globalbcz == NEU_DIR ) std::cerr << "NEU_DIR NOT IMPLEMENTED "<<std::endl;
-    if (globalbcz == dg::PER ) std::cerr << "PER NOT IMPLEMENTED "<<std::endl;
+    if (boundary == DIR_NEU ) std::cerr << "DIR_NEU NOT IMPLEMENTED "<<std::endl;
+    if (boundary == NEU_DIR ) std::cerr << "NEU_DIR NOT IMPLEMENTED "<<std::endl;
+    if (boundary == dg::PER ) std::cerr << "PER NOT IMPLEMENTED "<<std::endl;
     
     return A;
 }

    cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), values.size());
    A.row_indices = row_indices; A.column_indices = column_indices; A.values = values;

    if (globalbcz == DIR_NEU ) std::cerr << "DIR_NEU NOT IMPLEMENTED "<<std::endl;
    if (globalbcz == NEU_DIR ) std::cerr << "NEU_DIR NOT IMPLEMENTED "<<std::endl;
    if (globalbcz == dg::PER ) std::cerr << "PER NOT IMPLEMENTED "<<std::endl;
    
    return A;
}

            for(unsigned k=0; k<pyF.size(); k++)
                for( unsigned l=0; l<pxF.size(); l++)
                    pxyz[k*g.n()+l]= 1.*pyF[k]*pxF[l];
            if (boundary == dg::DIR)
            {
                if ( x[i]==g.x0() || x[i]==g.x1()  ||y[i]==g.y0()  || y[i]==g.y1())
                {
