Roy and John: Thanks for the suggestions.
>> Would it be possible to get all the values you'll need from tau_system
>> before-hand? Depending on the size of the problem, and assuming you
>> could also need multiple values per quadrature point, that might not be
>> feasible. I'm not sure how you'd store them all either.
>> The benefit would be you wouldn't have to worry about the initialization
>> states of two different systems at the same time...
In my case it would be possible to store the tau solution data in a
vector. I'm not sure of the best way to index that vector within
element_time_derivative though. Assuming the same quadrature rule is
used in both systems, I suppose one could keep a counter (as a member
variable in NavierSystem) that is incremented at each new quadrature
point, and use that to index the tau data?
>>> Hmm... even if you did it that way, you'd also have to mess around
>>> with the quadrature rule in tau_system, to ensure that you were
>>> getting values at the integration points you're using in your Navier
>>> Stokes system. Perhaps refactoring FEMSystem code is more trouble
>>> than it's worth and writing your own function really is the way to go.
>>> I'd still make it two functions - a tau_elem_init(Elem &) to handle
>>> the dof_indices(), FE::reinit(), etc. calls, so that you do batch
>>> calculations rather than trying to evaluate one point at a time.
This seems like a good approach as well. I'll let y'all know what I end
up doing :)