Design of output of computationally expensive parameter function

[Dominik Steinberger]

Dear all,

I am currently implementing an Eigenstrain solver. For this purpose I have currently defined an FESystem via

FESystem<3> fe(FE_Q<3>(1), 3);


with the FEValuesExtractors

constFEValuesExtractors::Vector displacements(0);


The DoFs are subsequently distributed using the regular

DoFHandler<3> dof_handler(triangulation);
dof_handler.distribute_dofs(fe);


At the moment I have a function which computes the eigenstrain at a given point

class EigenstrainFunction : public TensorFunction<2, 3>
{
public:
    virtual Tensor<2, 3> value(const Point<3> &p) const {
        Tensor<2, 3> eigenstrain;
       
        // some rather expensive computation

        return plastic_eigenstrain;
    }
};

and during the assembly of my right hand side this function gets called at the quadrature points.
So far so good.

But now I would like to have the eigenstrain in my VTK output as well, exactly as the program "sees" it, i.e. in a way that the values at the quadrature points are preserved.
To my understanding the correct way to do this is to use one of the VectorTools::project routines to get respective values at the degrees of freedom.
Unfortunately the eigenstrain does not have any corresponding degrees of freedom.
And projecting the values would again lead to an evaluation of this rather expensive function at the exact same quadrature points that I have already done an evaluation on.

What is the most elegant solution to this problem under the constraint, that the function should only be called once per quadrature point?

Options that come to mind, which I am absolutely not sure about whether they are even feasible:

• As my eigenstrain matrix is symmetric I could pack it into two solution vectors and reuse the DoFHandler from the 3 component displacement to get the projected components, then use them in the assembly of the right hand side to interpolate their values at the quadrature points using the shape functions.
• I could set up another FESystem
  
  FESystem<3> elasticity_element(FE_Q<3>(1), 3);
  FESystem<3> eigenstrain_element(FE_Q<3>(1), 6);
  FESystem<3> fe(elasticity_element, 1, eigenstrain_element, 1);
  
  but could I then use FEValuesExtractors::SymmetricTensor on the eigenstrain element? And could I only project my function onto the eigenstrain element?

Best regards,

Dominik

[Jean-Paul Pelteret]

Dear Dominik,

I believe that your question is directly along the lines of this recent discussion:

https://groups.google.com/forum/#!topic/dealii/FW9Jy5_CUGY

Does the information contained in that post help you? 

Regards,

Jean-Paul