Re: [Libmesh-devel] Loop subdivision surface elements

[Roman V. <vet...@et...>]

>   And if you're not comfortable generating the pull request I could help - if you send along the patch I can try and break it up into number of smaller changes on a new branch, and we can start discussion.

I'm very new to git, and it's all still a bit confusing, so I think I'll 
gladly accept the help. :)

Here's the patch. You'll need to run ./bootstrap as I didn't include the 
autoconf/automake/libtool output.

A few more explanations on the code:

1. As already mentioned (also by Norbert), the solution on a given 
element depends on the 1-ring of that element (i.e., the direct neighbor 
elements). Obviously, boundary elements lack such neighbors. Therefore, 
subdivision meshes require special preparation, where an additional 
layer of elements is added along the mesh boundaries. These additional 
elements do not belong to the physical domain.

2. The added boundary elements are named "ghost" elements in the 
literature. I guess this name somewhat clashes with the "ghost" elements 
that lie on processor domain boundaries.

3. Any program using subdivision surface elements will need to operate 
only on non-ghost elements. To this end, a new element type TRI3SD is 
added, providing the is_ghost() member.

4. Nodes that have valence 6 are called "regular", otherwise they are 
"irregular". Elements are regular iff all their nodes are regular. 
Regular elements have exactly 12 shape functions (quartic box splines). 
Irregular elements, however, need special treatment.

5. Even though the subdivision surface shape functions are quartic, one 
Gauss point per element is perfectly sufficient. The traditional 
theoretical lower bound (predicting a minimum of 6 points for triangular 
elements with quartic shape fcts) doesn't apply due to the overlapping 
subdiv shape functions.

6. For efficiency, we internally reorder the nodes of TRI3SD elements in 
such a way that if the elements is irregular, the irregular vertex is 
always the first (see also 
http://sourceforge.net/mailarchive/message.php?msg_id=20808509) and only 
the first. Clearly, this means that only meshes with elements with at 
most one irregular vertex are supported. (The mesh preparator will 
properly complain otherwise.)

7. The custom DofMap originally mentioned by Norbert 
(http://sourceforge.net/mailarchive/message.php?msg_id=20778890) was 
avoided in the present implementation by generalizing the way DofMap 
computes the dof_indices.

8. Due to 1., boundary conditions are non-trivial to impose. Have a look 
at miscellaneous_ex11 or the original paper (Int. J. Numer. Meth. Engng. 
47, 2039-2072 (2000)) for one way to impose them. See Int. J. Numer. 
Meth. Eng. 61, 380–405 (2004) for an alternative constraint approach.

9. I don't know what happens if a user attempts to use AMR, libmesh's 
boundary functionality, or constraints on the subdivision surface 
elements. I simply haven't tried.

Feel free to ask more questions.

-Roman