From: Boyce G. <gri...@ci...> - 2011-10-05 23:31:14
|
Is FE::reinit() in fe_boundary.C intended to work when qrule is NULL? Line 203 of fe_boundary.C calls qrule->get_points() whether or not qrule has been set to anything: // Find where the integration points are located on the // full element. std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); this->side_map(elem, side.get(), s, qrule->get_points(), qp); Can this revert to using inverse_map if there are user-supplied points? Or is this not the right solution? Thanks, -- Boyce |
From: Lorenzo A. B. <bot...@gm...> - 2011-10-19 09:03:03
|
> > Looking at the implementation of FE<Dim,T>::side_map, it seems like it should be possible to avoid calling inverse_map() in the reinit(elem,side,...), even when there are user-supplied points, by doing something like: > > // Find where the integration points are located on the > // full element. > const std::vector<Point>* ref_qp; > if (pts != NULL) > ref_qp = pts; > else > ref_qp = &qrule->get_points(); > std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); > this->side_map(elem, side.get(), s, *ref_qp, qp); > > One other change is that side_map() now needs to reinit values if shapes_on_quadrature is false. > > Does this seem correct? It appears to work in my tests. > > -- Boyce Yes, it is correct. I think that it makes sense to add if(qrule->shapes_need_reinit()) shapes_on_quadrature = false; after line 251 and 156 of fe_boundary.C (that is after qrule->init()) so that the reinit on side behaves like the reinit on elements. Thus you should be able to fix side_map simply by adding a !shape_on_quadrature condition. This of coarse means that the pts vector should be a vector of points on the reference side. Lorenzo |
From: Boyce G. <gri...@ci...> - 2011-10-19 23:09:03
Attachments:
fe_boundary.C.diff
|
On 10/19/11 5:02 AM, Lorenzo Alessio Botti wrote: > >> >> Looking at the implementation of FE<Dim,T>::side_map, it seems like it should be possible to avoid calling inverse_map() in the reinit(elem,side,...), even when there are user-supplied points, by doing something like: >> >> // Find where the integration points are located on the >> // full element. >> const std::vector<Point>* ref_qp; >> if (pts != NULL) >> ref_qp = pts; >> else >> ref_qp =&qrule->get_points(); >> std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); >> this->side_map(elem, side.get(), s, *ref_qp, qp); >> >> One other change is that side_map() now needs to reinit values if shapes_on_quadrature is false. >> >> Does this seem correct? It appears to work in my tests. >> >> -- Boyce > > Yes, it is correct. > I think that it makes sense to add > > if(qrule->shapes_need_reinit()) > shapes_on_quadrature = false; > > after line 251 and 156 of fe_boundary.C (that is after qrule->init()) so that > the reinit on side behaves like the reinit on elements. > > Thus you should be able to fix side_map > simply by adding a !shape_on_quadrature condition. > > This of coarse means that the pts vector should be > a vector of points on the reference side. OK -- attached is a diff with these changes. Can someone with commit privileges please take a look and, if it seems OK, commit to the repo? Thanks! -- Boyce |
From: Lorenzo A. B. <bot...@gm...> - 2011-10-06 07:59:38
|
Hi Boyce, > // Find where the integration points are located on the > // full element. > std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); > this->side_map(elem, side.get(), s, qrule->get_points(), qp); > > Can this revert to using inverse_map if there are user-supplied points? > Or is this not the right solution? > > I introduced the above change in order avoid inverse maps and speed up quadrature evaluation on sides. It should be reverted to inverse map if the quadrature points you are providing are defined in the physical frame (xyz), if the quadrature points are defined on the reference side the code should work. However, if your points are in the physical frame the best option (you don't need to modify the reinit behavior) is probably to find the reference points in your code with inverse map and then use reinit providing a vector of points FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, your_xyz_points, your_reference_frame_points); fe_neighbor_face->reinit(neighbor, &your_reference_frame_points); This assuming that you are using basis functions defined in the reference frame. Lorenzo > ------------------------------------------------------------------------------ > All the data continuously generated in your IT infrastructure contains a > definitive record of customers, application performance, security > threats, fraudulent activity and more. Splunk takes this data and makes > sense of it. Business sense. IT sense. Common sense. > http://p.sf.net/sfu/splunk-d2dcopy1 > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: Boyce G. <gri...@ci...> - 2011-10-18 18:20:50
|
On 10/6/11 3:59 AM, Lorenzo Alessio Botti wrote: > Hi Boyce, > >> // Find where the integration points are located on the >> // full element. >> std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); >> this->side_map(elem, side.get(), s, qrule->get_points(), qp); >> >> Can this revert to using inverse_map if there are user-supplied points? >> Or is this not the right solution? >> >> > > I introduced the above change in order avoid inverse maps and speed up > quadrature evaluation on sides. > It should be reverted to inverse map if the quadrature points you are providing > are defined in the physical frame (xyz), if the quadrature points are defined > on the reference side the code should work. OK --- two issues with the current implementation of the side version of reinit in fe_boundary.C are (1) it breaks if qrule is NULL, and (2) it seems to ignore the argument "pts" (even if pts is specified, the routine will try to use the initialization points specified by the quadrature rule). Looking at the code, I'm not sure what the correct modification should be. I originally thought that lines 199--200 should be changed to something like: // Find where the integration points are located on the // full element. std::vector<Point> qp; if (pts != NULL) this->inverse_map(elem, xyz, qp, tolerance); else this->side_map(elem, side.get(), s, qrule->get_points(), qp); or possibly: // Find where the integration points are located on the // full element. std::vector<Point> qp; if (!shapes_on_quadrature) this->inverse_map(elem, xyz, qp, tolerance); else this->side_map(elem, side.get(), s, qrule->get_points(), qp); However, now I'm not sure if that is correct. I think that the call to inverse_map() (which is currently commented out) is actually redundant, because xyz should be initialized by the call to this->compute_face_map(). Right? Should the code simply be: // Find where the integration points are located on the // full element. std::vector<Point> qp; if (pts == NULL) this->side_map(elem, side.get(), s, qrule->get_points(), qp); ??? > However, if your points are in the physical frame the best option > (you don't need to modify the reinit behavior) is probably to find the reference points > in your code with inverse map and then use reinit providing a vector of points > > FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, your_xyz_points, your_reference_frame_points); > fe_neighbor_face->reinit(neighbor,&your_reference_frame_points); > > This assuming that you are using basis functions defined in the reference frame. OK --- I will give this a try, thanks. -- Boyce |
From: Boyce G. <gri...@ci...> - 2011-10-18 18:59:40
|
On 10/6/11 3:59 AM, Lorenzo Alessio Botti wrote: > However, if your points are in the physical frame the best option > (you don't need to modify the reinit behavior) is probably to find the reference points > in your code with inverse map and then use reinit providing a vector of points > > FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, your_xyz_points, your_reference_frame_points); > fe_neighbor_face->reinit(neighbor,&your_reference_frame_points); > > This assuming that you are using basis functions defined in the reference frame. The points that I'm passing to reinit are reference coordinates on the side element. To avoid calling reinit(elem,side,tolerance,pts), I'm now trying: AutoPtr<Elem> side_elem = elem->build_side(side); fe_boundary->reinit(side_elem.get(), &reference_coords); but I get an error: ERROR: Unsupported 2D element type!: 0 If I print out the element type for the side proxy element, it seems like I always get 0, regardless of the parent element type. Are proxy elements not supposed to know their type? -- Boyce |
From: Boyce G. <gri...@ci...> - 2011-10-18 19:12:46
|
On 10/18/11 3:02 PM, Boyce Griffith wrote: > > > On 10/6/11 3:59 AM, Lorenzo Alessio Botti wrote: >> However, if your points are in the physical frame the best option >> (you don't need to modify the reinit behavior) is probably to find the >> reference points >> in your code with inverse map and then use reinit providing a vector >> of points >> >> FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, >> your_xyz_points, your_reference_frame_points); >> fe_neighbor_face->reinit(neighbor,&your_reference_frame_points); >> >> This assuming that you are using basis functions defined in the >> reference frame. > > The points that I'm passing to reinit are reference coordinates on the > side element. > > To avoid calling reinit(elem,side,tolerance,pts), I'm now trying: > > AutoPtr<Elem> side_elem = elem->build_side(side); > fe_boundary->reinit(side_elem.get(), &reference_coords); > > but I get an error: > > ERROR: Unsupported 2D element type!: 0 Grrr, I forgot to switch dimension of the FE object... -- Boyce |
From: Boyce G. <gri...@ci...> - 2011-10-18 19:30:11
|
On 10/18/11 1:48 PM, Boyce Griffith wrote: >> However, if your points are in the physical frame the best option >> (you don't need to modify the reinit behavior) is probably to find the >> reference points >> in your code with inverse map and then use reinit providing a vector >> of points >> >> FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, >> your_xyz_points, your_reference_frame_points); >> fe_neighbor_face->reinit(neighbor,&your_reference_frame_points); >> >> This assuming that you are using basis functions defined in the >> reference frame. Aha, I remember now why I need FE::reinit(elem,side,tolerance,pts) for this code --- I need to be able to compute surface normals at the points, and it seems like that must be done via this version of reinit(). --- Boyce |
From: Boyce G. <gri...@ci...> - 2011-10-18 19:59:03
|
On 10/18/11 3:32 PM, Boyce Griffith wrote: > > > On 10/18/11 1:48 PM, Boyce Griffith wrote: >>> However, if your points are in the physical frame the best option >>> (you don't need to modify the reinit behavior) is probably to find the >>> reference points >>> in your code with inverse map and then use reinit providing a vector >>> of points >>> >>> FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), elem, >>> your_xyz_points, your_reference_frame_points); >>> fe_neighbor_face->reinit(neighbor,&your_reference_frame_points); >>> >>> This assuming that you are using basis functions defined in the >>> reference frame. > > Aha, I remember now why I need FE::reinit(elem,side,tolerance,pts) for > this code --- I need to be able to compute surface normals at the > points, and it seems like that must be done via this version of reinit(). Looking at the implementation of FE<Dim,T>::side_map, it seems like it should be possible to avoid calling inverse_map() in the reinit(elem,side,...), even when there are user-supplied points, by doing something like: // Find where the integration points are located on the // full element. const std::vector<Point>* ref_qp; if (pts != NULL) ref_qp = pts; else ref_qp = &qrule->get_points(); std::vector<Point> qp; //this->inverse_map (elem, xyz, qp, tolerance); this->side_map(elem, side.get(), s, *ref_qp, qp); One other change is that side_map() now needs to reinit values if shapes_on_quadrature is false. Does this seem correct? It appears to work in my tests. -- Boyce |