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
Index: src/fe/fe_boundary.C
===================================================================
--- src/fe/fe_boundary.C (revision 4878)
+++ src/fe/fe_boundary.C (working copy)
@@ -158,6 +158,8 @@
{
// initialize quadrature rule
qrule->init(side->type(), side_p_level);
+ if(qrule->shapes_need_reinit())
+ shapes_on_quadrature = false;
// FIXME - could this break if the same FE object was used
// for both volume and face integrals? - RHS
@@ -196,8 +198,13 @@
// 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, qrule->get_points(), qp);
+ this->side_map(elem, side.get(), s, *ref_qp, qp);
// compute the shape function and derivative values
// at the points qp
@@ -253,6 +260,8 @@
{
// initialize quadrature rule
qrule->init(edge->type(), elem->p_level());
+ if(qrule->shapes_need_reinit())
+ shapes_on_quadrature = false;
// We might not need to reinitialize the shape functions
if ((this->get_type() != elem->type()) ||
@@ -304,7 +313,8 @@
side_p_level = std::max(side_p_level, elem->neighbor(s)->p_level());
if (side->type() != last_side ||
- side_p_level != _p_level)
+ side_p_level != _p_level ||
+ !shapes_on_quadrature)
{
// Set the element type
elem_type = elem->type();
------------------------------------------------------------------------------
The demand for IT networking professionals continues to grow, and the
demand for specialized networking skills is growing even more rapidly.
Take a complimentary Learning@Ciosco Self-Assessment and learn
about Cisco certifications, training, and career opportunities.
http://p.sf.net/sfu/cisco-dev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
