## Re: [Libmesh-users] Subdivision surface based FEM

 Re: [Libmesh-users] Subdivision surface based FEM From: Norbert Stoop - 2008-11-11 16:56:28 ```Roy Stogner wrote: >>> Interesting. I've heard of subdivision elements being used for >>> surface mesh refinement, but in a context where the subdivision >>> surface was only C1 in the limiting sense; each discrete mesh was >>> still faceted. We could do something like that relatively easily, but >>> of course it wouldn't be as accurate unless your formulation is >>> insensitive to domain boundary discontinuities. >> >> Hm, I'm not sure if I understand your last comment, > > That's because I didn't understand your context - I was picturing a > volume domain for which you wanted a C1 boundary (which can be > critical for some fourth order problems); you're talking about a 2D > manifold domain in 3D. Yeah, sorry... Let me specify the physical problem: It is all about the simulation of thin shells, where the shell is just a 2D manifold embedded in 3D. For the physical description the curvature must be square integrable, requiring the manifold to be C1. >> but yes, the surface is C1 only in the limit case. > > In that case you don't need an exotic mapping at all; all you need is > a function to correctly "snap" points to their subdivision-defined > positions when you refine the surface. That's something we've wanted > to do (for boundaries, if not manifolds) for a while, and it may be > hard to come up with a satisfactory user API for it, but once the > API's defined the implementation would be relatively simple. > >> These mappings are expressed in the phi_map, dphidxi_map etc., >> right? So, as a dirty hack, *in principle* one could overwrite the >> FE::init_shape_functions and others to do the mapping right for this >> particular subdivision element. Is my understanding correct? > > If you want a C1 surface even in the discretized case, then fixing the > _map values would be necessary. Ok, that's what I assumed. If you're happy with a C0 discretized > surface that converges (in some norms) to a C1 surface, then all you > need to do is snap midedge points (and midface points on quads) to > their proper places after all the elements touching them have been > refined. That should properly be done by giving Elem::refine() access > to some abstract base class that can return enough information about > the manifold or domain boundary shape, but you could do it (as a dirty > hack) just by looping over nodes on active elements after each > refinement and adjusting them for consistency with the subdivision > surface defined by their parent elements. Ah, wait - the whole idea is to have an initial mesh which is *never* refined. It merely acts as a "control mesh" for the limiting surface. Each node of the control mesh corresponds to a (analytically) known point on the limit surface. This does not have to be the same spacial location, as some subdivision schemes change the position of existing nodes. As shown below, calculations will be made in the limit surface. But since this nodal correspondence is one-to-one, nodal properties can as well be treated as if they belonged to the control points. In fact, this map can also be constructed for any point in a triangle on the control mesh, needed for the quadrature points. It is given by x(xi,eta) = sum_nb N_nb(xi,eta) * x_nb where the sum is over the neighbor control nodes of the triangle in question and x_nb is the position of these control points. N_nb are the shape functions of the limit surface, which are given analytically depending on the chosen subdivision scheme. E.g., x(0,0) gives the limit surface position of the triangle's lower left control point, x(0.5,0.5) of the barycentric middle etc. In other words, we have a parametrization of the limit surface in terms of only the position of neighboring control point positions x_nb. This parametrization is like any isoparametric map from the standard triangle to the manifold. C1 continuity is given since the parametrization is "non-local" in the sense that neighboring control points are included in the sum above. To summarize, we abuse the traditional mesh as a control mesh. We define a parametrization of the limit surface which naturally gives us the needed mathematical objects such as derivatives, surface patches etc. Since we *know* where each control point converges to in the limit, we can assign back calculated nodal values to the control points and assemble the system as usual. Hope this helps to clarify... -Norbert ```