On Wed, 27 Jul 2011, Lorenzo Alessio Botti wrote: > here is a patch to avoid inverse maps on sides reinit. > I hope that the main developers will consider it and > have some time to give it a look since it might break a lot of code...
Aha, so it might. You're adding a new method that can make DG much faster on h-conforming meshes, but also using it in FE::reinit(side) pretty much makes the entire library (well, anyone integrating boundary conditions...) dependent on it. I'll try the new version on some of our more numerically rigorous app regression tests and see if I see any breakage before I risk committing anything. You know, we could and probably *should* use a similar trick to your DG code elsewhere in the library. Instead of seeing whether element vs neighbor sides are conforming via "if (refinement_type == "p")", could we test "if (elem->level() == neighbor->level())"? Then we only have to fall back on inverse_map when the particular elem/neighbor interface is non-conforming, not when any interface anywhere in the mesh might be. I haven't tested that idea on ex21 yet; would you? We could factor the whole test-then-side-map-or-inverse-map routine into a little utility function somewhere and then it would be a drop-in replacement for inverse_map in the JumpErrorEstimator, in the current DG example, and eventually in the FEMSystem DG stuff that Truman is working on. But there is some breakage with your current patch: the library doesn't yet *compile* with --enable-everything (we need to add InfFE::side_map() or InfFE gets turned into a non-instantiable Abstract Base Class. I'm a bit conflicted about that. If we just put an InfFE::side_map that's libmesh_not_implemented(), then that won't break our InfFE example (integrating over a boundary "side" that doesn't exist would be pointless) but it would break anyone doing symmetry BCs or DG with infinite elements. Anyone doing that? If so speak up now; if not then I think I'll commit anyway and leave those features broken for now. > FEBase::last_side,last_edge are now the type of the side and not the side > number. I'm going to switch the data type to ElemType accordingly; otherwise I get compiler warnings about unsafe signed/unsigned conversions. > I think that it is safe to call > FE::init_face_shape_functions (qrule->get_points(), side.get()) > only if the type of the side is different but I might be wrong... Hmm... I think you're correct, *only* in the case that reference_side_points is generated by the same quadrature rule as was used previously. Adding a "side_p_level != _p_level" test ought to take care of that case. The part of this problem I always thought would be nasty is the relative rotation of sides, but it looks like your elem_nodes_map takes care of that neatly. > In FEBase I added a method that returns the nodes coordinates for reference > frame elements. > This method could also be used to speed-up the solution output when using > reference frame finite elements. > I did it based on FEBase::on_reference_element() > would you please check that everything is correct. This kind of coding is > error-prone and I'm not sure about the PIRAMID5. John, double-check me on this: I believe the PYRAMID5 in Lorenzo's patch is wrong - the Pyramid5 we use, going by both my own hazy memory and by the ::shape() equations in fe_lagrange_shape_3D.C, has points 0-3 in the zeta=0 plane, not zeta=-1. I checked the other elements and they looked fine, so at least if you made and missed a typo I've missed it too. > In dbg mode a warning says that the parameter tolerance of the method reinit > is not used. > The examples seems to run fine. We used to want a tolerance there to pass along to inverse_map; if we're no longer dependent on inverse_map then we no longer need that tolerance. For now I'll comment out that argument to get rid of the warning. In the long run I think we ought to do the same optimization on edge_reinit and then get rid of the tolerance parameters to both. > For dG discretizations on 3D linear meshes the time required for shape > computations on sides is reduces by 1/3. > On second order meshes higher gains are to be expected and false negative > Jacobians can be totally avoided. This is excellent. And ought to be extended further - IIRC just running KellyErrorEstimator also takes much longer than it should (i.e. it can be a sizeable minority of runtime, when it should almost always be trivial) in part because of the many inverse maps involved. I think some of the constraint and projection calculations could benefit from avoiding inverse maps on sides too. Only trouble is, if we make the library depend on this code internally, I don't want AMR to break for the less common configurations. So that would mean doing better than a libmesh_not_implemented() for InfFE: we'd need to add the InfFE entries to FEBase::get_refspace_nodes() and either move side_map up the inheritance tree or duplicate it. I don't have any time to do that myself right now, though. And otherwise the patch looks good. So (assuming it passes my app tests, and unless anyone has objections) I'm planning to commit things nearly as-is, with InfFE as libmesh_not_implemented() and with the PYRAMID5 case fixed. And if anyone wants faster jump error estimators any time soon and has more time to play with it, that should soon be low-hanging-fruit. --- Roy ------------------------------------------------------------------------------ BlackBerry® DevCon Americas, Oct. 18-20, San Francisco, CA The must-attend event for mobile developers. Connect with experts. Get tools for creating Super Apps. See the latest technologies. Sessions, hands-on labs, demos & much more. Register early & save! http://p.sf.net/sfu/rim-blackberry-1 _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
