From: Karl Tomlinson <k.tomlinson@au...>  20060808 03:24:01

Roy Stogner writes: > On Mon, 7 Aug 2006, Karl Tomlinson wrote: > >> On Fri, 4 Aug 2006 16:09:16 0500 (CDT), Roy Stogner wrote: >> >>> On Sat, 5 Aug 2006, Karl Tomlinson wrote: >>> >>>> I'm used to thinking of x,y,z as variables in themselves, like >>>> dependent (system) variables in libMesh, each with their own >>>> FEType. But I think it is reasonable that x,y,z should all have >>>> the same FEType (but not necessarily the same as that of each of >>>> the dependent variables). >>> >>> This makes sense when the domain is transformed from a cartesian grid, >>> but what do you do on more general unstructured meshes? >> >> I now see that libMesh is storing xyzbased derivatives for the >> dofs in Hermite meshes, whereas I'm used to thinking of xibased >> derivatives. > > That's right. Basically I wanted to make it easier to handle > selectively refined (whether by libMesh's hierarchical AMR or by a > priori mesh grading) meshes. Once you do that, d/dxi is no longer > consistent between neighboring elements, but d/dx still is. I can see advantages in this, but it might provide some challenges for Hermite geometries as the dx/dxi scale factors are not yet available. > >> When thinking in terms of xyzderivatives, there needs to be some >> way of calculating the scale factors (dx/dxi) for the geometry. >> This can be done with edgelengths or something similar but it may be >> iterative and may not be simple. >> >> When the dofs are xiderivatives this step is not an issue. > > Well, if your mesh is topologically equivalent to a subset of a > cartesian grid, there's no problem: your "master elements" can all be > the same, and then making a distinction between d/dx and d/dxi (for > mapping functions) doesn't make sense. Once you do a little > hierarchic refinement, then there's a scaling factor that comes into > play, but it's still simple. Are you talking about prefinement here (I haven't thought about these issues) or hrefinement (I can see use of scaling factors here)? > What I still can't wrap my head around > is what you do for really unstructured meshes: Our meshes do try to be as topologically similar as possible to a subset of a cartesian grid. There are places where this is not the case but we don't have C1 continuity there  we don't require C1 for our second order problems but it makes geometries look nice and provides a dofcheap highorder interpolation. >>> Imagine you've got three Hermite quads meeting at one point, for >>> example  what are the mapping variables at that node? If xi and eta >>> are different between elements, derivatives with respect to xi and eta >>> won't be uniquely defined on element interfaces. >> >> This situation does make things more complicated for the >> xiderivative implementation. Some linear constraints would be needed >> to require C1 continuity. (The scale factor matrix provides this >> in the xyzderivative implementation.) > > I'm not convinced it's possible to get C1 continuous mapping functions > with Hermite elements in this situation, regardless of how you code > them. Yes I don't think it is possible in general. It can be C1 at the node but only almost C1 away from the node. But the same issue applies to solution variables with Hermite elements, doesn't it? > The Hermite elements get away with fewer degrees of freedom > than you would expect a C1 quad or hex to need, because they take > advantage of the fact that you're using meshes where the mixed > derivatives (e.g. d^2/dxideta in 2D) are in the same "direction" on > all neighboring elements. I don't see how that can happen unless you > have four quads (or 8 hexes) meeting at every node. The Hermite elements normally share their derivatives with adjacent elements but this need not be the case. The zeroth order derivatives are always going to be consistent (and shared) between elements, but, in the 3 quads at a point case, the firstorder derivatives could be expressed as linear combinations of one pair of derivatives to ensure C1 continuity at the node (but not necessarily on the edges). I don't know whether their is an appropriate way to tie mixed derivatives to each other or whether this would help with continuity at all (probably not). This can be difficult to implement. The way it has been done here is by having a separate node for each of the three elements and then tie the appropriate dofs together, but we have only done this to provide the necessary C0 continuity. (It could also be done by having the derivatives as element parameters so only the zeroth derivatives are automatically shared, but this would mean that constraints would normally be involved even for the standard 4 quads at a node case.) These are the issues that I was/am hoping ClouchTocher elements might solve. > >> A similar issue exists where an element is adjacent on one side to >> 2 elements of higher hrefinement level (with a hanging node). > > This is much easier to handle by comparison; sure there's some scaling > involved, but as long as the mapping is consistent (all the mapping > values on the hanging node correspond to the correct linear > combinations of the values on the neighboring vertices), the result > should still be C1 and should still be straightforward to compute. Yes. The same could also apply for the derivative values at nodes that are shared (not hanging) if using d/dxi, but there may be issues in mapping nodal dofs to Hermite (element) dofs. > >>> Also, it's easy to conceive of situations where the mapping variables >>> take up a different amount of data than any solution variable. It >>> might be nice to handle NURBS geometries via a two step map: a first >>> order Lagrange mapping from master elements to a few unit boxes, then >>> a few NURBS points to define a map from those into physical space. >>> Perelement NURBS information could work but would waste RAM. >> >> I'm not familiar enough with NURBS to comment here. > > Imagine a quadratic mapping, then. Say you want to discretize the > domain from 1<x<1 and x^21<y<x^2+1. You can describe a mapping into > that domain by using triangles which each have 6 nodes and a quadratic > map, but you could also do it by mapping a bunch of 3 node master > triangles onto a unit square, then using a single 9 node quadratic map > to transform the unit square onto the physical domain. In the former > case you end up using approximately 1 mapping node per 2 triangles, in > the latter you use approximately 3. > > This isn't a huge memory savings until you get up to high p mappings > or complicated NURBS, I know, but I don't want to preclude it unless > we have to. > >> I was thinking that geometric variables should be as similar as >> possible to solution (or any other  e.g. material parameter) >> variables, but are you saying that geometric variables may be more >> complicated than solution variables? > > It's not so much that they're more or less complicated, it's that > they're "differently" complicated. I think that for geometric > variables you may often find that the most memoryefficient way to > store the mapping is as a composition of many simple (perhaps linear > Lagrange) mappings from master elements into "geometric patches" and a > few complicated (perhaps NURBS) mappings from those patches into the > physical domain. There's no analog to such a process among solution > variables. I think I'm seeing the issues here: It is inefficient to refine a geometry that doesn't change. By providing a 2 step mapping, only the simplest mapping need be refined. And there is no need for a high order representation of the internal element boundaries in a homogeneous domain, when usually representing the boundary is the only reason for the high order representation. Much of our legacy code was designed for Lagrangian finitedeformation (nonlinear) solid mechanics. In this problem, solution variables are the geometry variables. Similarly for problems involving fitting meshes to boundary data points. > >>>> Should the shape functions for the map be tied more closely to an >>>> FE for the Mesh (or an Elem) rather than to the FEs for each of >>>> the dependent variables in each System on the Mesh? >>> >>> I'd rather not have a mapping FE associated with each Elem unless we >>> can think of a cunning data structure; I felt bad about just adding a >>> byte for p_level until I realized most compilers would have turned >>> that byte into alignment padding anyway. Adding a perElem pointer to >>> a mapping FEType object would be a bit much. > > I'm willing to reconsider this. Looking over elem.h, it's clear that > we're going to have something like a dozen pointers per Elem on even > 2D meshes; adding another one won't kill us. I was wondering if this was the case, but if we use this argument too many times we'll soon be using lots of memory. > I'm still not sure what good it would be, though. Generally when we > have finite elements that are compatible with each other, we put them > in the same FEFamily, even though one is a quad and the other a > triangle. I can't imagine situations where I'd want to mix and match > mapping FEFamilies; it wouldn't be safe. I'm happy to leave this for now at least. Our meshes usually only have one FEFamily, and solution variables in libmesh only have one FEFamily. What I'm primarily interested in is providing similar flexibility to independent variables as what is already provided for dependent variables. > >> (As I understand it, prefinement essentially changes the FEType >> through the order and so it is just having different FEFamilys >> that is the issue. Correct me if I'm wrong.) > > You're correct. Also we currently only support adaptive p refinement > of hierarchictype basis functions, which makes mixing different p > levels in the same mesh a little easier. Sounds sensible. > >> >> I'm not familiar with CloughTocher but they look pretty clever. >> Hermites don't provide TRI or TET ElemTypes so having some >> CloughToucher to provide this may be useful. Is this a feasible >> combination (with hanging node constraints perhaps)? > > I had hoped it would be when I started coding them, but no, that > doesn't work. To make the interface C1, you'd need to constrain the > flux on the Hermite side to be quadratic instead of cubic  but > there's no way to do that without the constraint "spilling over" into > neighboring elements, and that can ruin the locality or even the > approximation accuracy of your basis functions. OK then. Thanks for your investigation. > > There is a quad macroelement which is compatible with the > CloughTochers and which is likely a little more efficient in a hybrid > mesh than dissecting the quads and using CT alone, but I don't think > the slight gain is worth the coding time. Sounds reasonable. > > Also, keep in mind that there's no Clough Tet elements in libMesh, > either. C1 tets are ugly; you need something like p=11 to build them > without macroelements, and the best macroelements suitable for any Tet > mesh require p=5 polynomials and can't be restricted below p=3. My > main motivation for writing the Hermite class in the first place was > that I wanted to get some C1 3D results without spending a year or two > bug hunting through a hundred macroelement terms. Sounds scary. But I don't know what "restricted below p=3" means. Is this something to do with restricting a fine mesh to a coarse mesh? Or is is something to do with the second order derivative parameters involved? > >> Maybe using only CloughTocher elements might be solution here. > > That would work, but I'd be worried about what CloughTocher mappings > might do to your quadrature rules. Granted, any kind of nonaffine > map can mess up your nice exact Gaussian quadrature, but the > CloughTocher basis functions have internal subelement boundaries > which might be even worse. The quadrature rules should be applied over each of the subelements individually (or through a macroquadrature rule). Isn't this also an issue when integrating solution variables or their shape functions? > ... I guess having a space which was C1 almost > everywhere might be good for some applications, but my fourth order > problems would start producing bad results. If the solution variables need to C1 wrt xyz, then I think there are some requirements for the geometry representation. du/dx = du/dxi * dxi/dx Scaling factors are chosen so that this quantity is continuous at the nodes, but, if dxi/dx is different in adjacent elements and varies within an element, then there may not be C1 continuity along the whole edge (or face). (IIRC a necessary condition for C1 continuity is that the ratios of scale factors in adjacent elements must be equal, which is not necessarily true of a Lagrange mesh.) The requirements are satisfied if a C1 representation is used for the geometry. 