From: Roy S. <roy...@ic...> - 2006-08-04 21:09:21
|
On Sat, 5 Aug 2006, Karl Tomlinson wrote: > So, rather than rewrite everything again, we would like to use > libMesh as the core finite element library for our projects > (although there may still be some internal debates over whether we > should be using Fortran...). Aren't there always internal debates about one language vs. another? The only reasons to prefer Fortran over C++ are performance and ease-of-learning. For a code like libMesh which has been developed as a learning and research tool rather than for production engineering, C++ was clearly better: any performance problems are better solved by throwing cheap hardware at a C++ code than by throwing expensive graduate student time at a Fortran code, and the original developers were already past the hard part of the C++ learning curve. > Hopefully we will also be able to contribute additional > functionality to libMesh. That's the open source idea! These things sort of snowball. I got involved when I realized it would be much easier to add Clough-Tocher elements and W^{2,p} calculations to libMesh than it would be to add AMR and parallelism to my own Matlab code. Neither Ben nor John had any need for Hermite elements, but now that I've been lured into writing them it might attract more developers? Good stuff. > However, we have many existing meshes that use Hermite > interpolation for the geometric variables (xyz). > > I see that the Hermite FEFamily has been added since 0.5.0, but > AFAICT from FE< Dim, T >::init_shape_functions, libmesh always > uses the Lagrange FEFamily for geometric variables (the map). This is true. > FE< Dim, T >::reinit has the comment: > > "In the future we can specify different types of maps" > > Does this indicate that more general interpolations for the > geometric variables were intended? There are all sorts of features that have always been intended, but some will be easier to implement than others... > Would it be reasonable to implement this? I suspect implementing other mappings to the point where you can run all the standard examples may be reasonable; no harder than adding AMR to non-Lagrange variables was, at least. The big problem will be getting all the odd combinations of features right. Right now adaptive p refinement works for most calculations, for example, but you can't write p refinement code capable of restarts because none of the Mesh I/O classes are aware of element p_level. A couple other users are now discovering that nobody ever updated the MeshFunction support classes to work with adaptive h refinement. It's hard to keep everything in sync. > Any thoughts on how this might be done? > > 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? 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. 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. Per-element NURBS information could work but would waste RAM. I'll think about it some more, but no bright ideas are coming to mind at the moment. > 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 per-Elem pointer to a mapping FEType object would be a bit much. Even if we do something more memory efficient (separate per-FEType std::vector<Elem *>s in Mesh? A tree storing the first and last element ID number in each contiguous range with the same FEType?), mixing different mapping FETypes in the same mesh is a little weird. I guess more flexibility is better than less, but to simplify code elsewhere (the FE caching I just added to CVS, for instance) it would be nice if we could assume there was just one per-Mesh mapping FEType. --- Roy Stogner |