Re: [Libmesh-users] p refinement and FE::reinit(elem,side)

[Roy S. <roy...@ic...>]

On Thu, 28 Feb 2008, Lorenzo Botti wrote:

> I've found the problem affecting the p refinement in my elliptic_dg
> code...
> In fe_boundary.C the method void FE<Dim,T>::reinit(const Elem* elem,
> const unsigned int s, const Real tolerance)
> will initialize the shape functions if ((this->get_type() != elem-
> >type()) || (s != last_side) || this->shapes_need_reinit())...
> I think we need to add a condition  like if (this->get_p_level() !=
> elem->p_level())... I've tried it creating a method
> FeBase::get_p_level() and it works
> but maybe you have a better solution.

Good catch, and I think your fix is the right way to go.  If you've
got a patch, send it; otherwise I'll get to it myself this weekend.

> Another option needed by the method FE<Dim,T>::reinit(const Elem*
> elem, const unsigned int s, const Real tolerance) would be a parameter
> extraorder to
> increase the quadrature order if the element neighbor is at an higher
> p level.
> Something like
> if (neighbor->p_level() > elem->p_level())
> {
>      unsigned int extraorder =  neighbor->p_level() - elem->p_level();
> }
> fe->reinit(elem, side, extraorder)
>
> and
> FE<Dim,T>::reinit(const Elem* elem, const unsigned int s,unsigned int
> extraorder, const Real tolerance)
> {
>
> ....
> qrule->init(side->type(), elem->p_level() + extraorder);
> ...
> }
>
> Again, maybe you have a better idea...

Hmm... I don't like making the user worry about it.  It's already too
easy to write a libMesh code that assumes you're not using feature X
and then wonder why it breaks when you turn X on.  I'd rather not make
"adaptive p refinement" more likely to be such a feature.

The side version of reinit() already has an Elem* and a side; why not
have it use the max p_level between the elem and it's neighbor by
default?  Then if the user really need to do something else, changing
the quadrature rule on the fly is still an option.  We could also make
it easier to change just quadrature rule order on the fly, but I think
that should be a method in QBase, not FEBase.
---
Roy