I want to use Libmesh in a multiphysics context where different solvers, including solvers not written in Libmesh, share information.

The interface between solvers is very simple: they send and receive constant or writeable references to vectors of doubles.  I want to send both actual solutions and gradients.


For example, in the case of a first order electrostatic (Poisson) solver , I want to generate :

a vector <double> V ( containing the containing the nodal values of the potential and numbered the same indices as the nodes of the mesh )

three vector<double> Ei ( containing the containing the values of the electric field components in the centroids and numbered the same indices as the elements of the mesh )


Now here are my questions:


1)       I can get the values ordered per DOF via:

            std::vector<double>      v;

             equation_systems.build_discontinuous_solution_vector (v);


The indexing of the values in the vector v is not exactly the way I want it to be , but I figured out that by looping through the elements I can recover the nodal values for each element. This works, but it is not very elegant ( or efficient) as each node is being treated several times.


Question : Do you have any suggestions how to do this better ? Ideally, I want to loop through the nodes and find the corresponding value in v.


2)       How do I recover the gradients ? I can go through the whole exercise of building again an FE object, and recuperate the derivatives of the shape functions in the qpoints, which should allow to calculate the gradients ( along the lines of the construction of grad_u_old in ex.9). That seems a bit overkill to me, but I cannot think of a simpler way.


Question: What is the simplest way to do this? If going through building FE is the only way, can you please write a simple example ? There is plenty of opportunity to make errors while converting between the different indices, and I ‘d feel more comfortable if I get an example.



All the best,


Jacques Kools