On Tue, Oct 28, 2008 at 9:20 PM, Adam Arbree <[EMAIL PROTECTED]> 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 non-null 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
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
