I do not have a lot of first hand experience with this, as I've never took the time to implement a proper solution for this problem. Also, I don't know a lot about FEM methods. So what I am going to say is mainly based on observing how various Mathematica functions work.
Don't use Sow/Reap for this
I don't think
Reap are designed with this application in mind.
It is recommended to always use tags (seconds argument of
Reap. Otherwise things like this may happen:
f[x_] := Reap[x; Sow]
Reap[f[Sow]] (* users doing this may not know about the internals of f *)
Reap catches everything (incorrectly) and leaves nothing for the outer one.
Having functions which do not
Sowed values can be similarly unsafe, even if using tags. Also, having to use tags already decreases the user-friendliness of the approach.
Bundling in lists is simple
The simplest way is to return a list, but of course it's possible to do better when writing a polished package.
In version 10, you can also consider bundling the information in an
Association. The result will be very similar to using "properties", as
FittedModel does, except that there can't be a default behaviour of the output (similar to how
FittedModel can be used as a numerical function).
Objects with properties
I recommend using solution (2): try to imitate basic object oriented programming and return an expression which:
Doesn't show its (complex and irrelevant) internal structure by default. Use
Format to achieve this.
Can be directly used for the main intended purpose, e.g. an
InterpolatingFunction or a
FittedModel can be directly used as numerical functions.
Has a number of properties or methods, which can be accessed as
This is very common with built-in Mathematica functions, but usually it is not (well) documented. You can use
MeshRegion, etc. this way. The
NDSolve plugin framework also uses this style.
Usually, properties or methods can be queried as
object["Methods"] (depending the type of the object). If you stick to this convention, then the use of your objects will already be familiar to many users.
To make these objects even more user friendly, you may choose to define functions for retrieving these properties. For example,
StiffnessMatrix[x_FEMSolution] := x["StiffnessMatrix"]
This is what the
DifferentialEquations`InterpolatingFunctionAnatomy` package does (take a look at its source).
As I see it, the main disadvantage of having top level functions like this to retrieve properties is that is requires introducing a large number of top-level functions (and choosing a large number of names, each of which has a potential to conflict with other packages). Using the
object["StiffnessMatrix"] syntax has the advantage of tying the name
"StiffnessMatrix" to this one special object type. The name cannot conflict with anything except other properties or methods of the same type of object.
This technique is really emulating classes from object oriented programming (without things like inheritance), and solves the same type of problem classes do: organization. If a certain function is only going to be used with a certain type of data structre, then it's best to formally define the data structure and tie the name of the function to it.
So there are many advantages to using solution (2), but what are the disadvantages?
One possible disadvantage is that if you bundle a lot of rarely used ancillary information into one object, it is going to take up a lot of space. Depending on the particular application, it may not be worth permanently storing all this information when most of the time the user only wants to use the main result only, which might be considerably smaller than the ancillary information.
For this situation I recommend the following:
A task like solving PDEs using a FEM method can usually be broken down into several steps. You may want to implement these steps as separate fuctions. These functions will work with a data structure that contains full information. You can make this (large) data structure accessible through the methods/properties calling style used by
You can also define a higher level function which calls these lower level functions in the correct order and computes the result in one go. This high level function can then return a data structure which contains less ancillary information.
Advanced users of your package can (if they wish to) use the lower level functions and work with the data structure that contains full information (but is more difficult to use).
Others can use the easy-to-use higher level functions and understand that the price to pay for this convenience is that they only have access to some ancillary information.
Mathematica 10 contains a FEM framework. You may want to read up on how it works. I did not yet study it in detail, but it seems to work in a similar way to what I described above. See here and related tutorials. You can also study the advanced NDSolve documentation and the NDSolve method plugin framework to get ideas.