Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
You can subscribe to this list here.
2003 
_{Jan}
(4) 
_{Feb}
(1) 
_{Mar}
(9) 
_{Apr}
(2) 
_{May}
(7) 
_{Jun}
(1) 
_{Jul}
(1) 
_{Aug}
(4) 
_{Sep}
(12) 
_{Oct}
(8) 
_{Nov}
(3) 
_{Dec}
(4) 

2004 
_{Jan}
(1) 
_{Feb}
(21) 
_{Mar}
(31) 
_{Apr}
(10) 
_{May}
(12) 
_{Jun}
(15) 
_{Jul}
(4) 
_{Aug}
(6) 
_{Sep}
(5) 
_{Oct}
(11) 
_{Nov}
(43) 
_{Dec}
(13) 
2005 
_{Jan}
(25) 
_{Feb}
(12) 
_{Mar}
(49) 
_{Apr}
(19) 
_{May}
(104) 
_{Jun}
(60) 
_{Jul}
(10) 
_{Aug}
(42) 
_{Sep}
(15) 
_{Oct}
(12) 
_{Nov}
(6) 
_{Dec}
(4) 
2006 
_{Jan}
(1) 
_{Feb}
(6) 
_{Mar}
(31) 
_{Apr}
(17) 
_{May}
(5) 
_{Jun}
(95) 
_{Jul}
(38) 
_{Aug}
(44) 
_{Sep}
(6) 
_{Oct}
(8) 
_{Nov}
(21) 
_{Dec}

2007 
_{Jan}
(5) 
_{Feb}
(46) 
_{Mar}
(9) 
_{Apr}
(23) 
_{May}
(17) 
_{Jun}
(51) 
_{Jul}
(41) 
_{Aug}
(4) 
_{Sep}
(28) 
_{Oct}
(71) 
_{Nov}
(193) 
_{Dec}
(20) 
2008 
_{Jan}
(46) 
_{Feb}
(46) 
_{Mar}
(18) 
_{Apr}
(38) 
_{May}
(14) 
_{Jun}
(107) 
_{Jul}
(50) 
_{Aug}
(115) 
_{Sep}
(84) 
_{Oct}
(96) 
_{Nov}
(105) 
_{Dec}
(34) 
2009 
_{Jan}
(89) 
_{Feb}
(93) 
_{Mar}
(119) 
_{Apr}
(73) 
_{May}
(39) 
_{Jun}
(51) 
_{Jul}
(27) 
_{Aug}
(8) 
_{Sep}
(91) 
_{Oct}
(90) 
_{Nov}
(77) 
_{Dec}
(67) 
2010 
_{Jan}
(25) 
_{Feb}
(36) 
_{Mar}
(98) 
_{Apr}
(45) 
_{May}
(25) 
_{Jun}
(60) 
_{Jul}
(17) 
_{Aug}
(36) 
_{Sep}
(48) 
_{Oct}
(45) 
_{Nov}
(65) 
_{Dec}
(39) 
2011 
_{Jan}
(26) 
_{Feb}
(48) 
_{Mar}
(151) 
_{Apr}
(108) 
_{May}
(61) 
_{Jun}
(108) 
_{Jul}
(27) 
_{Aug}
(50) 
_{Sep}
(43) 
_{Oct}
(43) 
_{Nov}
(27) 
_{Dec}
(37) 
2012 
_{Jan}
(56) 
_{Feb}
(120) 
_{Mar}
(72) 
_{Apr}
(57) 
_{May}
(82) 
_{Jun}
(66) 
_{Jul}
(51) 
_{Aug}
(75) 
_{Sep}
(166) 
_{Oct}
(232) 
_{Nov}
(284) 
_{Dec}
(105) 
2013 
_{Jan}
(168) 
_{Feb}
(151) 
_{Mar}
(30) 
_{Apr}
(145) 
_{May}
(26) 
_{Jun}
(53) 
_{Jul}
(76) 
_{Aug}
(33) 
_{Sep}
(23) 
_{Oct}
(72) 
_{Nov}
(125) 
_{Dec}
(38) 
2014 
_{Jan}
(47) 
_{Feb}
(62) 
_{Mar}
(27) 
_{Apr}
(8) 
_{May}
(12) 
_{Jun}
(2) 
_{Jul}
(22) 
_{Aug}
(22) 
_{Sep}

_{Oct}
(17) 
_{Nov}
(20) 
_{Dec}
(12) 
2015 
_{Jan}
(25) 
_{Feb}
(2) 
_{Mar}
(16) 
_{Apr}
(13) 
_{May}
(21) 
_{Jun}
(5) 
_{Jul}
(1) 
_{Aug}
(8) 
_{Sep}
(9) 
_{Oct}
(30) 
_{Nov}
(8) 
_{Dec}

2016 
_{Jan}
(16) 
_{Feb}
(31) 
_{Mar}
(43) 
_{Apr}
(18) 
_{May}
(21) 
_{Jun}
(11) 
_{Jul}
(17) 
_{Aug}
(26) 
_{Sep}
(4) 
_{Oct}
(16) 
_{Nov}
(5) 
_{Dec}
(6) 
2017 
_{Jan}
(1) 
_{Feb}
(2) 
_{Mar}
(5) 
_{Apr}
(4) 
_{May}
(1) 
_{Jun}
(11) 
_{Jul}
(5) 
_{Aug}

_{Sep}
(3) 
_{Oct}
(1) 
_{Nov}
(7) 
_{Dec}

2018 
_{Jan}
(8) 
_{Feb}
(8) 
_{Mar}
(1) 
_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 



1

2

3
(3) 
4
(6) 
5

6

7
(4) 
8
(12) 
9
(5) 
10

11
(1) 
12

13

14
(5) 
15

16
(2) 
17
(1) 
18

19

20

21

22
(1) 
23
(1) 
24

25
(1) 
26

27

28

29

30
(2) 
31



From: Roy Stogner <roystgnr@ic...>  20060808 15:26:45

On Tue, 8 Aug 2006, Tim Kr=F6ger wrote: >> That's exactly why the inverse_map function has the "bool secure" >> argument. If secure is false (as it is when contains_point() calls >> inverse_map() ), then Newton isn't expected to converge and >> inverse_map just exits after 10 iterations. > > Strange enough, this didn't happen last week when I had changed the cod= e such=20 > that only active elements were inserted into the tree but they were sti= ll=20 > handled only by their centroid: The code crashed because Newton didn't=20 > converge within 20 iterations (with warnings from the 11th iteration on= ). I=20 > have no explanation for that because when I look at the code, it seems = to me=20 > that you are right. Is it possible that the crash was when inverse_map() was being called from some other library function? We turn secure on in several places where we expect Newton to converge, just to doublecheck us in case one of those expectations is false. >> No good ones. All I can think to do is warn about the problem. When >> you add an interface for outofmesh MeshFunction extensions, let that >> interface turn off (well, replace) the linear search, and add an >> assert() or #ifdef DEBUG block in that code path which verifies that >> it doesn't see any quadratic mappings. > > Okay, so I think the best thing to do is: Whenever the outofmeshexte= nsion=20 > of MeshFunction is enabled, it will check Elem::has_affine_map() for al= l=20 > elements (at least in debug mode), and if any element returns true, an = error=20 > will be generated. But, however, I wonder whether Elem::has_affine_map= ()=20 > could occasionally return false due to rounding errors. What do you th= ink?=20 > Should the equality checks be replaced by some fabs()<eps stuff? As John pointed out there is an epsilon in the the point operator=3D=3D that handles that for us. That actually may cause some problems elsewhere in the library, now that I think about it. If we perturb mesh nodes by something like TOLERANCE/10, that should be enough to make a noticeable difference in a sufficiently fine finite element solution, but it won't be enough to convince the FE class that the elements are all distinct  we'd end up erroneously caching shape function & derivative calculations from element to element, as well as erroneously caching mapping derivative calculations from quadrature point to quadrature point.  Roy 
From: Roy Stogner <roystgnr@ic...>  20060808 15:17:40

On Tue, 8 Aug 2006, Karl Tomlinson wrote: > 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. Okay. If you don't need C1 geometry everywhere, then doing a pseudoHermite geometry on unstructured meshes starts to make a little more sense. I'm still not sure what the right thing to do at oddvalence nodes is, but now it's a "what's the best data structure" problem rather than a "how is this even possible" problem. Talking about a C1 geometry on an unstructured mesh was starting to confuse me. On an unstructured mesh you don't have one geometric mapping, you have an entirely different mapping for each finite element. Mapping continuity between elements usually isn't important at all; it's just that a couple of the exceptions (like combining Hermite mappings with Hermite solution spaces to handle more interesting geometries) are very interesting. >> 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? Yes  in fact if you'll look at the fe_hermite_* code, you'll see that we just give up and error() if we discover that an element isn't a rectangle parallel to the coordinate axes. In theory you could get C1 solutions on any mesh where the xi/eta/zeta directions are consistent at element interfaces, but I only needed rectilinear domains and so I left the code simple. >> 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. This is all going to take some thought. If you want to maintain a C1 map between an unstructured 2D reference grid and a deformed grid, though, you're right that using the CloughTochers should make that pretty natural. > 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. Exactly. > 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. And that makes perfect sense too  once you start perturbing the geometry it's very unlikely that you can still get away with a patchbased rather than an elementbased description. I'd like to have an implementation that can handle both sorts of situation, if possible  but as a first draft it may make sense to implement mapping information on a perelement level but try to keep the API abstract enough to allow memory optimization in the future. >> 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. It is, but only because I don't know enough about handling high polynomial degrees intelligently without getting swamped by floating point error. I've got libMesh's adaptive p refinement capped at around p=10 or p=11 for just that reason. > But I don't know what "restricted below p=3" means. Take a look at the quadratic "restricted CloughTocher" triangles for a 2D example: you can start with the cubic CT triangle, then constrain away the side flux degrees of freedom. The result has fewer degrees of freedom and p=2 accuracy like a normal quadratic element would, but it's still got some cubic terms so the quadrature rules don't get any cheaper. The tet situation is the same: you can start with a p=5 macroelement tet, then constrain away a bunch of the degrees of freedom to get something that has performance more like a cubic element. >> 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? Yes; how we do things now is to just copy a Gaussian rule over each subelement. I guess it's not a serious problem, but tripling the cost of FEM calculations isn't something to be happy about. Do you have any references to macroelementspecific quadrature rules in the literature? I'd been thinking about deriving something more efficient for the cubic CT, but I'd hate to reinvent the wheel. > If the solution variables need to C1 wrt xyz, then I think there > are some requirements for the geometry representation. It entirely depends on the finite element space  CloughTochers will build a C1 solution space on top of any geometry; Hermites need consistent local coordinate directions at nodes.  Roy 
From: <tim@ce...>  20060808 14:07:59

Dear John and Roy, On Tue, 8 Aug 2006, John Peterson wrote: > Tim Kr=F6ger writes: > > Okay, so I think the best thing to do is: Whenever the > > outofmeshextension of MeshFunction is enabled, it will check > > Elem::has_affine_map() for all elements (at least in debug mode), and > > if any element returns true, an error will be generated. But, > > however, I wonder whether Elem::has_affine_map() could occasionally > > return false due to rounding errors. What do you think? Should the > > equality checks be replaced by some fabs()<eps stuff? > > I think the point=3D=3D tests do have fabs() tests in them. You are right, thank you (except that the doxygen docu tells something=20 wrong about Point::operator=3D=3D and Point::operator!=3D). Hence,=20 everything is clear now unless someone objects to the aforementioned=20 suggestion. Best Regards, Tim 
From: John Peterson <peterson@cf...>  20060808 13:45:02

Tim Kr=F6ger writes: > Dear Roy, >=20 > On Mon, 7 Aug 2006, Roy Stogner wrote: >=20 > > On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > > > >> On Fri, 4 Aug 2006, Roy Stogner wrote: > >>=20 > >>> The inverse map iteration definitely won't converge for some ele= ments > >>> if the physical point is outside the element, but that shouldn't= cause > >>> a problem for Elem::constains_point() anymore; if its call to > >>> inverse_map() doesn't converge, the result is a distant point th= at > >>> will make sure on_reference_element() is false. > >>=20 > >> No: If the nonlinear solver doesn't converge, it will call error(= ) and thus=20 > >> crash. At least, this call to error() should be removed  perha= ps by=20 > >> giving the solver class an optional argument. > > > > That's exactly why the inverse_map function has the "bool secure" > > argument. If secure is false (as it is when contains_point() call= s > > inverse_map() ), then Newton isn't expected to converge and > > inverse_map just exits after 10 iterations. >=20 > Strange enough, this didn't happen last week when I had changed the=20= > code such that only active elements were inserted into the tree but=20= > they were still handled only by their centroid: The code crashed=20 > because Newton didn't converge within 20 iterations (with warnings=20= > from the 11th iteration on). I have no explanation for that because= =20 > when I look at the code, it seems to me that you are right. >=20 > >> Any suggestions? > > > > No good ones. All I can think to do is warn about the problem. W= hen > > you add an interface for outofmesh MeshFunction extensions, let = that > > interface turn off (well, replace) the linear search, and add an > > assert() or #ifdef DEBUG block in that code path which verifies th= at > > it doesn't see any quadratic mappings. >=20 > Okay, so I think the best thing to do is: Whenever the=20 > outofmeshextension of MeshFunction is enabled, it will check=20 > Elem::has_affine_map() for all elements (at least in debug mode), an= d=20 > if any element returns true, an error will be generated. But,=20 > however, I wonder whether Elem::has_affine_map() could occasionally=20= > return false due to rounding errors. What do you think? Should the= =20 > equality checks be replaced by some fabs()<eps stuff? I think the point=3D=3D tests do have fabs() tests in them. J 
From: <tim@ce...>  20060808 06:45:10

Dear Roy, On Mon, 7 Aug 2006, Roy Stogner wrote: > On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > >> On Fri, 4 Aug 2006, Roy Stogner wrote: >>=20 >>> The inverse map iteration definitely won't converge for some elements >>> if the physical point is outside the element, but that shouldn't cause >>> a problem for Elem::constains_point() anymore; if its call to >>> inverse_map() doesn't converge, the result is a distant point that >>> will make sure on_reference_element() is false. >>=20 >> No: If the nonlinear solver doesn't converge, it will call error() and t= hus=20 >> crash. At least, this call to error() should be removed  perhaps by= =20 >> giving the solver class an optional argument. > > That's exactly why the inverse_map function has the "bool secure" > argument. If secure is false (as it is when contains_point() calls > inverse_map() ), then Newton isn't expected to converge and > inverse_map just exits after 10 iterations. Strange enough, this didn't happen last week when I had changed the=20 code such that only active elements were inserted into the tree but=20 they were still handled only by their centroid: The code crashed=20 because Newton didn't converge within 20 iterations (with warnings=20 from the 11th iteration on). I have no explanation for that because=20 when I look at the code, it seems to me that you are right. >> Any suggestions? > > No good ones. All I can think to do is warn about the problem. When > you add an interface for outofmesh MeshFunction extensions, let that > interface turn off (well, replace) the linear search, and add an > assert() or #ifdef DEBUG block in that code path which verifies that > it doesn't see any quadratic mappings. Okay, so I think the best thing to do is: Whenever the=20 outofmeshextension of MeshFunction is enabled, it will check=20 Elem::has_affine_map() for all elements (at least in debug mode), and=20 if any element returns true, an error will be generated. But,=20 however, I wonder whether Elem::has_affine_map() could occasionally=20 return false due to rounding errors. What do you think? Should the=20 equality checks be replaced by some fabs()<eps stuff? Best Regards, Tim 
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. 
From: Roy Stogner <roystgnr@ic...>  20060807 20:18:39

On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > On Fri, 4 Aug 2006, Roy Stogner wrote: > >> The inverse map iteration definitely won't converge for some elements >> if the physical point is outside the element, but that shouldn't cause >> a problem for Elem::constains_point() anymore; if its call to >> inverse_map() doesn't converge, the result is a distant point that >> will make sure on_reference_element() is false. > > No: If the nonlinear solver doesn't converge, it will call error() and = thus=20 > crash. At least, this call to error() should be removed  perhaps by=20 > giving the solver class an optional argument. That's exactly why the inverse_map function has the "bool secure" argument. If secure is false (as it is when contains_point() calls inverse_map() ), then Newton isn't expected to converge and inverse_map just exits after 10 iterations. Of course, inverse_map() is *supposed* to exit by returning an outofbounds point, not just whatever its last iterate was. Derek, didn't we fix this? Or did you send me a patch I forgot to commit to CVS? > The other problem that I see when leaving the linear search in is that = when=20 > I will later implement the possibiliy to evaluate the MeshFunction outs= ide=20 > the domain covered by the grid, the linear search will slow this down=20 > considerabely. That is undesirable. > What about adding some heuristic safety margin around the bounding box = in=20 > the case of quadratic mappings? Well, not really satisfying ... As long as the margin is a proveable upper bound, I'll be happy with it. Trying to figure out what would be an upper bound on a quadratic Hex makes my head spin, though. > Perhaps better: Enable the linear search only if quadratic mappings are= =20 > involved. That sounds good too, but deciding whether or not quadratic mappings are involved is easier said than done. I'd go so far as to say that secondorder geometric elements with linear mappings is one of the most common usage patterns in libMesh, since we need the extra nodes in the secondorder elements to store degree of freedom numbers for any p>1 finite element. > Any suggestions? No good ones. All I can think to do is warn about the problem. When you add an interface for outofmesh MeshFunction extensions, let that interface turn off (well, replace) the linear search, and add an assert() or #ifdef DEBUG block in that code path which verifies that it doesn't see any quadratic mappings.  Roy 
From: Roy Stogner <roystgnr@ic...>  20060807 20:05:20

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. > 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. What I still can't wrap my head around is what you do for really 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. > > 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. 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. > 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. >> 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. >>> 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'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. > (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. >> 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 perMesh mapping FEType. > > For the most part I'm happy with one perMesh mapping FEFamily. > The possible exception is CloughTocher and Hermite FEFamilies. > > 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. 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. 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. > 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. > Are CloughTocher and Hermite FEFamilys similar enough to be > different ElemTypes of the same FEFamily? No; even if the code worked, the resulting function space would only be C0 conforming. 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.  Roy 
From: <tim@ce...>  20060807 09:43:38

Dear Roy, On Fri, 4 Aug 2006, Roy Stogner wrote: > On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > >> I did these changes now; see the attached patch. Please commit it to th= e=20 >> repository (provided that you are satisfied with it). > > One problem I see: in the case of quadratic mappings (or any more > complicated mappings, if the University of Auckland folks decide to > code them!) it's possible for an element's nodes to all be inside a > bounding box even though part of the element interior extends outside > the box. You are right, I did not think about this. > I don't know of any fast way to get the real bounding box for even a > quadratic element. This will be a very uncommon problem, though  > perhaps an adequate solution for now is just to leave in the fallback > linear search. > The inverse map iteration definitely won't converge for some elements > if the physical point is outside the element, but that shouldn't cause > a problem for Elem::constains_point() anymore; if its call to > inverse_map() doesn't converge, the result is a distant point that > will make sure on_reference_element() is false. No: If the nonlinear solver doesn't converge, it will call error() and=20 thus crash. At least, this call to error() should be removed =20 perhaps by giving the solver class an optional argument. The other problem that I see when leaving the linear search in is that=20 when I will later implement the possibiliy to evaluate the=20 MeshFunction outside the domain covered by the grid, the linear search=20 will slow this down considerabely. What about adding some heuristic safety margin around the bounding box=20 in the case of quadratic mappings? Well, not really satisfying ... Perhaps better: Enable the linear search only if quadratic mappings=20 are involved. Any suggestions? Best Regards, Tim 
From: Karl Tomlinson <k.tomlinson@au...>  20060807 06:36:54

Thanks, Roy, for your positive comments and thoughts. On Fri, 4 Aug 2006 16:09:16 0500 (CDT), Roy Stogner wrote: > On Sat, 5 Aug 2006, Karl Tomlinson wrote: > >> However, we have many existing meshes that use Hermite >> interpolation for the geometric variables (xyz). [snip] >> 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? I now see that libMesh is storing xyzbased derivatives for the dofs in Hermite meshes, whereas I'm used to thinking of xibased derivatives. 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. > 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.) A similar issue exists where an element is adjacent on one side to 2 elements of higher hrefinement level (with a hanging node). I'll have to look more at the libMesh (xyz) implementation. > 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. 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? >> 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. > > Even if we do something more memory efficient (separate perFEType > 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. (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.) > 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 perMesh mapping FEType. For the most part I'm happy with one perMesh mapping FEFamily. The possible exception is CloughTocher and Hermite FEFamilies. 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)? Maybe using only CloughTocher elements might be solution here. Are CloughTocher and Hermite FEFamilys similar enough to be different ElemTypes of the same FEFamily? Karl. 
From: Roy Stogner <roystgnr@ic...>  20060804 21:27:21

On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > I did these changes now; see the attached patch. Please commit it to t= he=20 > repository (provided that you are satisfied with it). One problem I see: in the case of quadratic mappings (or any more complicated mappings, if the University of Auckland folks decide to code them!) it's possible for an element's nodes to all be inside a bounding box even though part of the element interior extends outside the box. I don't know of any fast way to get the real bounding box for even a quadratic element. This will be a very uncommon problem, though  perhaps an adequate solution for now is just to leave in the fallback linear search. I'll add that search back in (but changed to loop over only active elements) and commit the patch to CVS this afternoon. > I also removed the fallback linear search since in my opinion it is no= t=20 > necessary any more. Anyway, if the linear search includes nonactive=20 > elements (as it did), it would again give wrong results, and if not, it= is=20 > very likely to crash, because the involved inverse map seems to be very= =20 > instable in the case that the point is very far away from the element. = At=20 > least, this is true (and not surprising) for distortet HEX8 elements. The inverse map iteration definitely won't converge for some elements if the physical point is outside the element, but that shouldn't cause a problem for Elem::constains_point() anymore; if its call to inverse_map() doesn't converge, the result is a distant point that will make sure on_reference_element() is false. > Additionally, I added a function get_point_locator() to MeshFunction, b= ecause=20 > in my application, I want to reuse the MeshFunction's PointLocator for=20 > additional purposes. Looks good.  Roy 
From: Roy Stogner <roystgnr@ic...>  20060804 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 easeoflearning. 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 CloughTocher 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 nonLagrange 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. Perelement 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 perElem pointer to a mapping FEType object would be a bit much. Even if we do something more memory efficient (separate perFEType 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 perMesh mapping FEType.  Roy Stogner 
From: <tim@ce...>  20060804 20:50:46

Dear Roy, On Fri, 4 Aug 2006, Roy Stogner wrote: > On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > >> This could be fixed by letting Tree::Tree() loop over active elements on= ly. > > I don't think that's a complete fix. The problem with mesh refinement > / coarsening isn't just that the mesh is more complicated, it's that > the mesh keeps *changing*. Okay, I did not think about this. I don't use MeshFunction with=20 changing grids  i.e., I rebuild the MeshFunction anyway each time=20 the mesh might change. > You'd be welcome (and thanked!) for the fixes. I'd like it if we > could figure out how to make fix 1 a little more robust, but even just > looping over active elements would be better than nothing. I'd rather > place a warning in the MeshFunction documentation than let this class > hold up a 0.6.0 release, and it would be better for that warning to > say "rebuild this object after every mesh change" and not "don't use > this object on a refined mesh". I did these changes now; see the attached patch. Please commit it to=20 the repository (provided that you are satisfied with it). I did set up a test scenario in which the old implementation gave=20 wrong results and the new one got it right. I also removed the fallback linear search since in my opinion it is=20 not necessary any more. Anyway, if the linear search includes=20 nonactive elements (as it did), it would again give wrong results,=20 and if not, it is very likely to crash, because the involved inverse=20 map seems to be very instable in the case that the point is very far=20 away from the element. At least, this is true (and not surprising)=20 for distortet HEX8 elements. Additionally, I added a function get_point_locator() to MeshFunction,=20 because in my application, I want to reuse the MeshFunction's=20 PointLocator for additional purposes. >> Another point is that I would like MeshFunction to be evaluable outside = the=20 >> mesh domain and return a useddefined constant there (instead of crashin= g).=20 >> I think I would be able to do these changes as well (using optional=20 >> parameters so that existing code will not change its behaviour). Any=20 >> comments on this? > Adding an optional parameter to the constructor would be fine. Adding > a method that sets that parameter might be better  although that's > kind of a matter of taste. I did not implement this yet, but I will probably do this later=20 (unless someone else does) because I am likely to need it. Best Regards, Tim 
From: Karl Tomlinson <k.tomlinson@au...>  20060804 19:16:06

Hi. The Bioengineering Institute at the University of Auckland is looking to replace its current FEM modelling software. The main reasons are probably: * We want to develop open source code. * Although our existing software has OpenMP multithreading for shared memory systems, we want software to run on clusters also, and we were intimidated by the idea of introducing message passing to legacy Fortran code full of common blocks and data structures that are accessed directly from everywhere. (See http://www.cmiss.org/ if you want more info.) We've been impressed by what we've seen of libMesh and it already implements many of the features that we are interested in. 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...). Hopefully we will also be able to contribute additional functionality to libMesh. 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). 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? Would it be reasonable to implement this? 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). 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? Thanks, Karl. 
From: Roy Stogner <roystgnr@ic...>  20060804 11:46:01

On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > On Thu, 3 Aug 2006, Roy Stogner wrote: > > 1.) Unless I missed something, not only active elements are stored in t= he=20 > Tree object, so that, as you mentioned, the PointLocator class might re= turn a=20 > coarse element instead of a descendant. > > This could be fixed by letting Tree::Tree() loop over active elements o= nly. I don't think that's a complete fix. The problem with mesh refinement / coarsening isn't just that the mesh is more complicated, it's that the mesh keeps *changing*. We need to make sure that new elements are inserted into the Tree after every refinement step and removed after every coarsening  and just as importantly for the modularity of the code, we ought to have that happen without messing too badly with the Mesh and Elem classes, and it ought to work even if the user never calls anything but MeshFunction::operator() after the AMR/C. > 2.) Elements are arranged by the location of their centroid only, altho= ugh=20 > they in fact might happen to have nonempty intersection with the bound= ing=20 > boxes of more than one of the children of a TreeNode object. After fix= ing=20 > Problem 1, this will cause Tree::find_element() to occasionally return = NULL=20 >  i.e. if the point to look for lies in a different TreeNode's boundin= g box=20 > than the centroid of the corresponding element. In this case,=20 > PointLocatorTree::operator() will fall back to a linear search and thus= still=20 > find the correct element but be slow. > > I think that this problem also can be solved quite easily: One would ju= st=20 > have to insert such elements in all the TreeNode's children that they h= ave=20 > nonempty intersection with. That sounds like the best way to do things. > I would not refuse to make these changes to the code myself, but before= , I=20 > would like to wait for your comments. You'd be welcome (and thanked!) for the fixes. I'd like it if we could figure out how to make fix 1 a little more robust, but even just looping over active elements would be better than nothing. I'd rather place a warning in the MeshFunction documentation than let this class hold up a 0.6.0 release, and it would be better for that warning to say "rebuild this object after every mesh change" and not "don't use this object on a refined mesh". > Another point is that I would like MeshFunction to be evaluable outside= the=20 > mesh domain and return a useddefined constant there (instead of crashi= ng).=20 > I think I would be able to do these changes as well (using optional=20 > parameters so that existing code will not change its behaviour). Any=20 > comments on this? Adding an optional parameter to the constructor would be fine. Adding a method that sets that parameter might be better  although that's kind of a matter of taste.  Roy 
From: <tim@ce...>  20060804 07:16:41

Dear Roy and all, On Thu, 3 Aug 2006, Roy Stogner wrote: > On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > >>> You might want look over the previous recent discussion on the mailing >>> lists before you come to rely on MeshFunction too heavily. There's a >>> lot of suspicious code in there that makes me believe it won't work >>> on an adaptively hrefined mesh. >>=20 >> Well, it seems to work for my meshes (which are hrefined but not=20 >> prefined) with the mentioned modifications. > > You might want to set up some validation problems and make sure that > "seems to work" is really "works exactly". I don't use MeshFunction > myself, so it's possible that I'm just being paranoid, but it looks to > me like the PointLocator classes might return a coarse element when they > should be returning one of it's descendant elements. If that happens, > then as long as you're using Lagrange elements MeshFunction may still > give you a solution without crashing, but the solution will be > slightly inaccurate. Thank you for your hint; I just looked at the code more closely and=20 now agree that there might be two problems in the code: 1.) Unless I missed something, not only active elements are stored in=20 the Tree object, so that, as you mentioned, the PointLocator class=20 might return a coarse element instead of a descendant. This could be fixed by letting Tree::Tree() loop over active elements=20 only. 2.) Elements are arranged by the location of their centroid only,=20 although they in fact might happen to have nonempty intersection with=20 the bounding boxes of more than one of the children of a TreeNode=20 object. After fixing Problem 1, this will cause Tree::find_element()=20 to occasionally return NULL  i.e. if the point to look for lies in a=20 different TreeNode's bounding box than the centroid of the=20 corresponding element. In this case, PointLocatorTree::operator()=20 will fall back to a linear search and thus still find the correct=20 element but be slow. I think that this problem also can be solved quite easily: One would=20 just have to insert such elements in all the TreeNode's children that=20 they have nonempty intersection with. I would not refuse to make these changes to the code myself, but=20 before, I would like to wait for your comments. Another point is that I would like MeshFunction to be evaluable=20 outside the mesh domain and return a useddefined constant there=20 (instead of crashing). I think I would be able to do these changes as=20 well (using optional parameters so that existing code will not change=20 its behaviour). Any comments on this? Best Regards, Tim 
From: Roy Stogner <roystgnr@ic...>  20060803 15:31:39

On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: >> You might want look over the previous recent discussion on the mailing >> lists before you come to rely on MeshFunction too heavily. There's a >> lot of suspicious code in there that makes me believe it won't work >> on an adaptively hrefined mesh. > > Well, it seems to work for my meshes (which are hrefined but not pref= ined)=20 > with the mentioned modifications. You might want to set up some validation problems and make sure that "seems to work" is really "works exactly". I don't use MeshFunction myself, so it's possible that I'm just being paranoid, but it looks to me like the PointLocator classes might return a coarse element when they should be returning one of it's descendant elements. If that happens, then as long as you're using Lagrange elements MeshFunction may still give you a solution without crashing, but the solution will be slightly inaccurate.  Roy 
From: <tim@ce...>  20060803 15:19:48

Dear Roy, On Thu, 3 Aug 2006, Roy Stogner wrote: > On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > You might want look over the previous recent discussion on the mailing > lists before you come to rely on MeshFunction too heavily. There's a > lot of suspicious code in there that makes me believe it won't work > on an adaptively hrefined mesh. Well, it seems to work for my meshes (which are hrefined but not=20 prefined) with the mentioned modifications. Without these=20 modifications, however, it was impossible to use MeshFunction at all=20 (unless using some workarounds). Best Regards, Tim 
From: Roy Stogner <roystgnr@ic...>  20060803 14:54:27

On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > Dear libMesh developers, You ought to send stuff like this to libmeshdevel, unless it's really time critical and Sourceforge's mailing lists are acting up again. > Could you please commit the attached patch to the libMesh cvs repositor= y? It=20 > consists of the following changes: > > 1.) Three changes in the MeshFunction class (see discussion in libmesh= users=20 > dated July 31, 2006); You might want look over the previous recent discussion on the mailing lists before you come to rely on MeshFunction too heavily. There's a lot of suspicious code in there that makes me believe it won't work on an adaptively hrefined mesh. > 2.) EquationSystems::read() now also recognizes the type ExplicitSystem. It looked good to me, so I'll commit it just as soon as I'm sure gcc agrees. I don't use MeshFunction myself, however  could Derek or someone who does doublecheck that part of the patch? Bearing in mind that if you've already had to work around MeshFunction's bugs yourself, Tim's fix is likely to break your workarounds.  Roy 