Return type of ScalarReferenceFiniteElement::EvalReferenceShapeFunctions and Gradients

[craffael]

I know we have discussed this a while ago, but given the new code in the branch fe_tools I think it's worth to reconsider the return type of the following two methods (both belong to ScalarReferenceFiniteElement):

virtual std::vector<Eigen::RowVectorXd> EvalReferenceShapeFunctions(const Eigen::MatrixXd& refcoords) const = 0;
virtual std::vector<Eigen::MatrixXd> GradientsReferenceShapeFunctions(const Eigen::MatrixXd& refcoords) const = 0;

I see three problems with the current return types:

1. Computing the local element matrices is unnecessary complex and inefficient: At the moment, both methods are exclusively called by LocCompLagrFEPreprocessor which transforms the returned data into a more suitable form.
2. Given a FE solution vector, it is rather complicated and inefficient to evaluate the solution function at a given point on the mesh. (for loop over local shape functions is necessary).
3. We solved a similar problem, namely the return type of Geometry::Jacobians() differently. This is not very consistent.

I suggest that the return types of both methods become Èigen::MatrixXd matrices of size NumRefShapeFunctions() x NumQuadPoints(), respectively NumRefShapeFunctions() x (DimLocal*NumQuadPoints). This solves all three problems mentioned above:

1. Computing e.g. the stiffness matrix is rather simple: GradientsScalarReferenceShapeFunctions()*W*GradientsScalarReferenceShapeFunctions().transpose() where w is a diagonal (assuming diffusion tensor is the identity) square matrix of size (DimLocal*NumQuadPoints) that has every quadrature weight DimLocal times on the diagonal. The mass matrix or boundary integrals can be computed very similarly.
2. If we have a local solution vector v of length NumRefShapeFunctions(), we can compute the value of the fe-solution at quadrature points qp by evaluating v.transpose()*EvalScalarReferenceShapeFunctions(qp).
3. Is clear.

I know that the my proposed change makes the code at first a little bit less readable, but as soon as someone has to deal with assembly of local element matrices, it's much easier this way... And I think in 90% of the use cases, these methods are called to compute local element matrices...

[hiptmair]

Since, a similar scheme is in use already for returning Jacobians, adopting the same convention for ScalarReferenceFiniteElement should not be too confusing. I endorse the proposed changes.

[craffael]

done.