From: Adam Arbree <arbree@cs...>  20081029 02:20:09

I am trying to integrate a jump term over the interior boundary between two 3D elements. To do a sanity check on my code, I run the code snippet at the end of this mail. It fails at the last assertion and I cannot figure out why. I am sorry if I am missing something simple. Adam /************* Example ********************/ // Assume that face has a nonnull neighbor of elem. // I am using MONOMIAL finite elements, TET4 elements and Fifth order Gaussian quadrature. // Both face_fe and neighbor_fe are initialized with an identical quadrature rule. const Elem* side_neighbor = elem>neighbor(face); const unsigned int side_on_neighbor = side_neighbor>which_neighbor_am_i(elem); libmesh_assert(elem>neighbor(face) == side_neighbor); libmesh_assert(side_neighbor>neighbor(side_on_neighbor) == elem); face_fe>reinit(elem, face); neighbor_fe>reinit(side_neighbor, side_on_neighbor); const std::vector<Point>& f_point = face_fe>get_xyz(); const std::vector<Point>& n_point = neighbor_fe>get_xyz(); unsigned int f_size = f_point.size(); unsigned int n_size = n_point.size(); libmesh_assert(f_size == n_size); for (unsigned int q_pt = 0; q_pt < f_size; ++q_pt) { const Point f = f_point[q_pt]; const Point n = n_point[q_pt]; libmesh_assert(f(0) == n(0) && f(1) == n(1) && f(2) == n(2)); //FAILS! } 