Dear Roy,

Thanks a lot for your reply. I will try to write a basic solver first and test those approaches. Then I can discuss with you guys what I find there. 

Sincerely Yours,

Lei Shi

On Wed, Aug 29, 2012 at 4:45 PM, Roy Stogner <> wrote:

On Wed, 29 Aug 2012, Lei Shi wrote:

I will do some experiments to compare the results based on computing
every cheap thing when needed and store everything needed.

Let us know the results!  Speculating about optimization ideas is a
poor substitute for profiling optimization ideas.

The other question is the element iterator. When we write ++iter,
what happened behind this operator? Does the current element need to
climb along the quad forest and go down to the lowest child? Can we
just construct a map with all active elementIndex -> element *
mappings? and directly use this map or vector to loop the elements.

In a SerialMesh, ++iter just increments through an STL vector,
potentially testing predicates (e.g. whether an element is active)
along the way.  ParallelMesh is similar but with a STL map under the

Searches through a quadtree/octree only happen when using a
PointLocator object, and most library code doesn't use those.  The big
exception is with PeriodicBoundary support; Derek or the other INL
people could tell you whether or not that turns out to be a
performance hit.

I know there is no anisotropic refinement in  the lib. How much
effort I need to put to add anisotropic? I mean how many percent the
current code based on those isotropic assumption?

For starters, "egrep '(->|\.)level\(*/*/*.{C,h} | wc" shows about a
hundred lines of code, the majority of which are implicitly assuming
that the idea of an isotropic "refinement level" makes sense.  Fixing
that to enable non-conforming anisotropic refinement would be a
serious effort.

Conforming anisotropic refinement might not be too bad.  Add the right
edge splitting algorithms for triangles and tets, don't bother
creating a tree structure, and you could get solutions without
changing anything else.  You'd want to create a new project_vector()
codepath to get efficient solutions.