From: John Peterson <jwpeterson@gm...>  20081029 02:51:12

On Tue, Oct 28, 2008 at 9:20 PM, Adam Arbree <arbree@...> wrote: > 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! This could fail because of floating point comparisons. Try a relative_fuzzy_equals comparison instead...  John 