From: Steffen Petersen <steffen.petersen@tu...>  20060426 20:59:33

Roy, I just tried to clean up some shape function files and code the Bernstein shapes in a more generic way (including arbitrary orders for quadrilaterals). Although I have adapted the corresponding fe_ElemType.C (i.e. the n_dofs functions), I got this error for p > 6 ERROR: unsigned char too small to hold Elem::_p_level! Recompile with Elem:_p_level set to something bigger! [0] /home/mubsp/tmp/libmesh/include/geom/elem.h, line 1325, compiled Apr 21 2006 at 12:21:08 For different element types, the order seems to be fixed at different levels. Can you perhaps give me a quick hint? Steffen 
From: Steffen Petersen <steffen.petersen@tu...>  20060426 20:59:33

Roy, I just tried to clean up some shape function files and code the Bernstein shapes in a more generic way (including arbitrary orders for quadrilaterals). Although I have adapted the corresponding fe_ElemType.C (i.e. the n_dofs functions), I got this error for p > 6 ERROR: unsigned char too small to hold Elem::_p_level! Recompile with Elem:_p_level set to something bigger! [0] /home/mubsp/tmp/libmesh/include/geom/elem.h, line 1325, compiled Apr 21 2006 at 12:21:08 For different element types, the order seems to be fixed at different levels. Can you perhaps give me a quick hint? Steffen 
From: Roy Stogner <roystgnr@ic...>  20060426 21:08:53

On Wed, 26 Apr 2006, Steffen Petersen wrote: > I just tried to clean up some shape function files > and code the Bernstein shapes in a more generic way > (including arbitrary orders for quadrilaterals). > Although I have adapted the corresponding fe_ElemType.C > (i.e. the n_dofs functions), > I got this error for p > 6 > > ERROR: unsigned char too small to hold Elem::_p_level! > Recompile with Elem:_p_level set to something bigger! > [0] /home/mubsp/tmp/libmesh/include/geom/elem.h, line 1325, compiled Apr 21 > 2006 at 12:21:08 > > For different element types, the order seems to be > fixed at different levels. Can you perhaps give > me a quick hint? Because no element type can support arbitrarily high p (in fact, floating point and matrix conditioning seem to kill us around p = 11), I've put a max_order function in the FEInterface code to tell the DoFMap just how high each element can go. If the base polynomial degree + the p_level elevation exceeds this max order, the DoFMap bumps the p_level down again. This obviously isn't a complete solution (if p refinement is impossible, 99% of the time you'd prefer h refinement to nothing) but it's a start. The problem comes in when the base polynomial degree exceeds that max_order even with p_level() == 0  then the DofMap tries to drop the p_level to 1 (which casts turn into 2^32  1 or so) and the p_level() function freaks. I put an assert in dof_map.C this morning to prevent that misleading error message from cropping up; should I add a more informative error message to dof_map.C to take its place? Oh, and to summarize all the rambling: In src/fe/fe_interface.C, edit the FEInterface::max_order() function at the bottom of the file to return a higher number for the BERNSTEIN/QUAD* case.  Roy 
From: Steffen Petersen <steffen.petersen@tu...>  20060426 21:28:17

> > On Wed, 26 Apr 2006, Steffen Petersen wrote: > > > I just tried to clean up some shape function files > > and code the Bernstein shapes in a more generic way > > (including arbitrary orders for quadrilaterals). > > Although I have adapted the corresponding fe_ElemType.C > > (i.e. the n_dofs functions), > > I got this error for p > 6 > > > > ERROR: unsigned char too small to hold Elem::_p_level! > > Recompile with Elem:_p_level set to something bigger! > > [0] /home/mubsp/tmp/libmesh/include/geom/elem.h, line 1325, > compiled Apr 21 > > 2006 at 12:21:08 > > > > For different element types, the order seems to be > > fixed at different levels. Can you perhaps give > > me a quick hint? > > Because no element type can support arbitrarily high p (in fact, > floating point and matrix conditioning seem to kill us around p = 11), Jep, and (depending on the problem you solve) the poor convergence of the iterative solvers may also somewhat blow the benifit of higher orders. > I've put a max_order function in the FEInterface code to tell the > DoFMap just how high each element can go. If the base polynomial > degree + the p_level elevation exceeds this max order, the DoFMap > bumps the p_level down again. This obviously isn't a complete > solution (if p refinement is impossible, 99% of the time you'd prefer > h refinement to nothing) but it's a start. > > The problem comes in when the base polynomial degree exceeds that > max_order even with p_level() == 0  then the DofMap tries to drop the > p_level to 1 (which casts turn into 2^32  1 or so) and the p_level() > function freaks. > > I put an assert in dof_map.C this morning to prevent that misleading > error message from cropping up; should I add a more informative error > message to dof_map.C to take its place? Perhaps a message with a hint to the max_order() function. I was expecting something like this, but obviously needed a more direct hint. Thanks, Steffen > > Oh, and to summarize all the rambling: > In src/fe/fe_interface.C, edit the FEInterface::max_order() function > at the bottom of the file to return a higher number for the > BERNSTEIN/QUAD* case. >  > Roy > 
From: Roy Stogner <roystgnr@ic...>  20060426 21:23:26

On Wed, 26 Apr 2006, Steffen Petersen wrote: > I just tried to clean up some shape function files > and code the Bernstein shapes in a more generic way > (including arbitrary orders for quadrilaterals). Oh, and out of curiosity  can you make the high order Bernstein bases hierarchic? (i.e. such that the basis functions at each node for order p+1 begin with the basis functions at that node for order p) Currently we only support padaptive constraints on elements with hierarchic bases on side nodes, System::project_vector will only handle p refinement correctly on elements with hierarchic bases on every node, and I have no plans to do the work necessary to handle more general cases anytime soon. ;) Of course, nonhierarchic highp elements are still useful, and just getting the correct function numbering on a hierarchic basis can be tricky (check out my shameful hack at src/utils/number_lookups.h sometime), but I'd hate for you to start getting bad adaptive p results without warning you about why. In fact, should I add a "hierarchic_bases" FEInterface function to let the library know which elements do and don't have that characteristic? It would allow the code to switch to the more expensive constraint and projection calculations that nonhierarchic elements require in the future, and would allow us to barf with a warning when someone tries to do p adaptivity and p projections on elements which don't support them today.  Roy 
From: Steffen Petersen <steffen.petersen@tu...>  20060427 12:16:20

> On Wed, 26 Apr 2006, Steffen Petersen wrote: > > > I just tried to clean up some shape function files > > and code the Bernstein shapes in a more generic way > > (including arbitrary orders for quadrilaterals). > > Oh, and out of curiosity  can you make the high order Bernstein bases > hierarchic? (i.e. such that the basis functions at each node for > order p+1 begin with the basis functions at that node for order p) Unfortunately, that is not possible :( The reason why I still prefer the Bernstein shapes is that they give much better performance for acoustic computations at high wave numbers. However, there is probably no meaningfull reason to use nonhierarchic shapes with local p refinement. > Currently we only support padaptive constraints on elements with > hierarchic bases on side nodes, System::project_vector will only > handle p refinement correctly on elements with hierarchic bases on > every node, and I have no plans to do the work necessary to handle > more general cases anytime soon. ;) > > Of course, nonhierarchic highp elements are still useful, and just > getting the correct function numbering on a hierarchic basis can be > tricky (check out my shameful hack at src/utils/number_lookups.h > sometime), but I'd hate for you to start getting bad adaptive p > results without warning you about why. > > In fact, should I add a "hierarchic_bases" FEInterface function to let > the library know which elements do and don't have that characteristic? Seems reasonable to me. Steffen > It would allow the code to switch to the more expensive constraint and > projection calculations that nonhierarchic elements require in the > future, and would allow us to barf with a warning when someone tries > to do p adaptivity and p projections on elements which don't support > them today. >  > Roy > 