I'm going to forward this to the libmesh-users mailing list; I hope
you don't mind. I'll explain why at the end of my reply.
On Tue, 31 Mar 2009, Brent Kraczek wrote:
> Sorry, please use this copy of poly_main.C
Here's (hopefully the only) two of your problems:
> // Element iterators are a nice way to iterate through
> // all the elements, or all the elements that have some property.
> // There are many types of element iterators, but here we will
> // use the most basic type, the const_elem_iterator. The iterator
> // el will iterate from the first to the last element. The
> // iterator end_el tells us when to stop. It is smart to make
> // this one const so that we don't accidentally mess it up!
> // const_elem_iterator el (mesh.elements_begin());
> // const const_elem_iterator end_el (mesh.elements_end());
>
> MeshBase::const_element_iterator el = mesh.elements_begin();
> const MeshBase::const_element_iterator end_el = mesh.elements_end();
During system assembly, looping over elements with "the most basic
type" of iterator was appropriate in example 3, where every single
element was to be assembled. It's not appropriate if you're using
adaptivity, in which case only the refined elements (what we call the
"active elements" should contribute to the matrix, and the coarser
("parent" or "ancestor") elements are just there to aid constraint
calculations and to enable future coarsening. It's also not
appropriate if you're running in parallel, in which case each
processor should only be responsible for assembling the contribution
from its local elements.
See the other examples; in its most general form the assembly loop is
from active_local_elements_begin() to end().
> // The element matrix and right-hand-side are now built
> // for this element. Add them to the global matrix and
> // right-hand-side vector. The SparseMatrix::add_matrix()
> // and NumericVector::add_vector() members do this for us.
> system.matrix->add_matrix (Ke, dof_indices);
> system.rhs->add_vector (Fe, dof_indices);
The second problem is not what's here, but what *isn't* here. Because
example 3 isn't adaptive, we don't bother constraining element
matrices before adding them to the global system. If you've got an
adaptive geometrically nonconforming mesh (i.e. nearly any adapted
mesh in libMesh) you need to handle the constrained degrees of freedom
on hanging nodes; the way we do this is with a constraint matrix pre-
and post-multiplication at the element level. Search through some of
the other examples and you'll find:
dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
which does this.
Now, why I sent this to the mailing list:
1. Because other users should see it. I guessed before looking what
bugs you were going to have, because you're not the first person I've
seen with these exact same bugs. Maybe having the SourceForge mailing
list archives come up with more search engine hits for: "libMesh
adaptivity bug problem error" or some such will help.
2. Because the other developers should see it. I thought we had
agreed *years* ago (after a faculty member visiting UT had the same
bugs) that we'd change all the library's example programs to use the
infintesimally-less-efficient-but-vastly-more-robust APIs. I guess we
never got around to actually making the change? That's too bad. All
the accumulated nanoseconds of CPU time saved by not using active
local iterators when they're not necessary aren't worth even one
user forgetting to switch to them when they are necessary.
So, I'm going to add lines for the active iterators and the constraint
operations to all our examples in SVN. The old iterator lines will be
commented out with explanation so that users who really want to avoid
every last JNE operation know how to perform a tiny optimization on
uniform meshes in serial... but then the consequences are on their
heads, not ours.
I'm also going to turn on local iterators in every example. This is
an even clearer choice: code with no adaptivity in it can't be run
adaptively, but parallel execution can be a execution-time decision,
so there's no reason not to have every example prepared to use it.
---
Roy
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
