## [Libmesh-users] Question about integration on sides

 [Libmesh-users] Question about integration on sides From: Adam Arbree - 2008-10-29 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 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& f_point = face_fe->get_xyz(); const std::vector& 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! } ```

### Thread view

 [Libmesh-users] Question about integration on sides From: Adam Arbree - 2008-10-29 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 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& f_point = face_fe->get_xyz(); const std::vector& 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! } ```
 Re: [Libmesh-users] Question about integration on sides From: John Peterson - 2008-10-29 02:51:12 ```On Tue, Oct 28, 2008 at 9:20 PM, Adam 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 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& f_point = face_fe->get_xyz(); > const std::vector& 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 ```
 Re: [Libmesh-users] Question about integration on sides From: Roy Stogner - 2008-10-29 02:54:03 ```On Tue, 28 Oct 2008, Adam 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. Just because the two neighboring elements are using an identical quadrature rule doesn't mean they're going to get matching points. For a simple example: Imagine two cubes touching each other. Now, rotate one of them 90 degrees around their common axis. The side quadrature points rotate with it, so even if they matched the other cube's quadrature points before, they won't after. And yet both meshes are perfectly valid. > I am sorry if I am missing something simple. Simple but disappointing: we can make no guarantees, even on conforming meshes, that quadrature points will match the way you want them to. The robust thing to do is the inverse maps that I suggested in a recent mailing list entry. Do it the way JumpErrorEstimator does it and your code will work even for non-conforming interfaces. If you're determined to avoid the inverse maps you could try to get the quadrature points on one element in terms of face xi-eta coordinates, figure out which coordinate (xi, eta, 1-xi, whatever) each of those corresponds to on the face of the adjoining element, and then use the result as a face quadrature rule there. But I wouldn't recommend it; the code I'm imagining wouldn't have to do any Newton solves but it would be a huge morass of switch statements and linear transformations. --- Roy ```
 Re: [Libmesh-users] Question about integration on sides From: Adam Arbree - 2008-10-29 13:47:40 ```Thanks again. I will go with the inverse transformation. Adam ________________________________ From: Roy Stogner [mailto:roystgnr@...] Sent: Tue 10/28/2008 10:53 PM To: Adam Arbree Cc: libmesh-users@... Subject: Re: [Libmesh-users] Question about integration on sides On Tue, 28 Oct 2008, Adam 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. Just because the two neighboring elements are using an identical quadrature rule doesn't mean they're going to get matching points. For a simple example: Imagine two cubes touching each other. Now, rotate one of them 90 degrees around their common axis. The side quadrature points rotate with it, so even if they matched the other cube's quadrature points before, they won't after. And yet both meshes are perfectly valid. > I am sorry if I am missing something simple. Simple but disappointing: we can make no guarantees, even on conforming meshes, that quadrature points will match the way you want them to. The robust thing to do is the inverse maps that I suggested in a recent mailing list entry. Do it the way JumpErrorEstimator does it and your code will work even for non-conforming interfaces. If you're determined to avoid the inverse maps you could try to get the quadrature points on one element in terms of face xi-eta coordinates, figure out which coordinate (xi, eta, 1-xi, whatever) each of those corresponds to on the face of the adjoining element, and then use the result as a face quadrature rule there. But I wouldn't recommend it; the code I'm imagining wouldn't have to do any Newton solves but it would be a huge morass of switch statements and linear transformations. --- Roy ```