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
------------------------------------------------------------------------------
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-d2d-oct
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
