On Mon, 14 May 2012, Derek Gaston wrote:
Ok - one more thing on this.... a question actually:
Why do we need to "contract()" for SerialMesh? Shouldn't it be ok
to have holes in the element / node numbering?
Hmmm... originally it was absolutely not okay. For starters, we had
a ton of code that assumed (mesh.max_elem_id() == mesh.n_elem()), and
in fact our code had no choice, since "MeshBase::max_elem_id()" didn't
even exist yet.
On the other hand, max_elem_id() *does* exist now, we're already using
it in places that want arrays indexed over all elements like the
ErrorEstimator subclasses, and in the long run we don't really want
arrays indexed over all elements anywhere anyway.
On the gripping hand: if we have a numbering full of holes, what is
SerialMesh::add_elem() supposed to do? It would currently ignore them
and append to the end, leaking memory until contract() is called. If
instead we tried to fill holes in in one at a time, that would be more
expensive in CPU than just letting a contract() bring everything
together at once. Even in ParallelMesh, where we expect the local
numberings to have holes, we renumber local elements contiguously to
keep the hole count down, letting us use a partially contiguous data
structure to keep asymptotic complexity down.
Wait a minute (we're four digressions in - anyone still following my
derailed train of thought?) - I never actually updated the mapvector
implementation to make use of mostly-contiguous storage. It would
still be nice to keep the option open, though.
Right now we kind of have weird behavior where nodes will be
changing numbers all the time as we adapt the mesh (even if a node
isn't destroyed and recreated it can be renumbered through
contract() ). That seems like weird behavior that is kind of
anti-user. For instance, what if you want to track the value of a
variable at a node that is created through refinement? You probably
can't do that with our current system because that node number will
be changing all the danged time.... right?
Right, but how would this be improved if we did try to avoid
renumbering? Unless we effectively don't allow coarsening, any node
created through refinement could still be destroyed in a subsequent
coarsening and then recreated with a different number later.
Even if we did come up with a way to track vanishing and reappearing
nodes, nodes are an artifact of the discretization - caring about a
coefficient at node N is like caring about the float at memory address
0xC0DE. Track the value by Point location rather than by Node index
and you get a sane result even if the node gets coarsened away, even
if you switch from a Lagrange to a non-nodal basis set, etc.
...
I write the above paragraph hesitantly, since I suspect that "track
the value at a node" is a hypothetical that's been oversimplified for
benefit of discussion and that your actual node-tracking needs aren't
so heedlessly dismissable. IIRC you've often been the first person to
chime in with "add auxilliary nodal or element-wise data to a separate
ExplicitSystem" when people have asked about node-indexed data in the
past, so trying to imagine what problem has you falling back on fixed
indices is frightening. Feel free to surprise me with "of course,
tracking by Point fixes everything, thanks!" rather than "here's a
more complicated case where your ideas totally fail", though. ;-)
Also, on a related note, it looks like in equation_systems.C
around line 194 sets "mesh_changed = true"... no matter what. So
calling eq->reinit() will cause a mesh.contract() _no matter
what_... even if you haven't modified the mesh.
It looks like this is in there because of backwards compatibility...
Correct. And it's been in there for quite a while, and we haven't
even fixed the example codes which contain examples of the API usage
requiring the backwards compatibility.
The idea (from Ben IIRC) was that we could often cut down on max
memory requirements in an AMR run by doing a refine-and-coarsen as a
two-step process: coarsen first, free up some memory after the
projections are done, then refine second. To make this work, we added
coarsen and refinement steps within eq->reinit().
The trouble is that all our user code was already manually calling
refine_and_coarsen_elements(), disabling any memory savings and also
making it impossible for eq->reinit() to examine the MeshRefinement
method return values to tell whether or not the mesh had been changed.
We ought to fix this situation. Maybe add a bool
"mesh_already_changed" argument to reinit(), defaulting to true for
backwards compatibility, but allowing users to pass in false if they
have an unchanged mesh and/or if they want to do the two-pass AMR?
Just my first idea, anyone have a better one?
Meaning that just calling eq->reinit() you will get your nodes and
elements renumbered!
True. We never really expect anyone to call reinit() except after
changing the mesh, though, do we? I suppose maybe if your mesh is
staying fixed but your stencil is changing you might be calling it to
reinitialize ImplicitSystem matrices?
---
Roy
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
