Using data from external matrix as source for calculations done with deal.II - best approach?

[Maxi Miller]

As source term for calculations I am using data which is written in an external matrix (because it has been calculated by an external function), together with the x- and y-values for each matrix entry. My current approach for using this data is to use an interpolation function from GSL, and call it for each cell while looping over them. This gives me the value at that position, interpolated.

Now, though, I would like to improve that approach. I was thinking about using the project-function and interpolating the matrix onto a component in my FESystem, but as far as I can see the project-function does not allow me to only change one component in my solution while leaving the other components untouched. Of course, I also could loop once over all cells before starting the main loop and solve the equation u = f, with f the value at each position, but I am not sure how efficient that is, especially when I can use the matrix-free vectorized approach for solving for the other components.

Therefore, are there other (maybe more efficient) approaches which I can use for transferring and interpolating the values from the original matrix to something which I then can use in my deal.II-functions?

Thanks!

[Wolfgang Bangerth]

On 1/22/20 11:17 AM, 'Maxi Miller' via deal.II User Group wrote:
> As source term for calculations I am using data which is written in an
> external matrix (because it has been calculated by an external
> function), together with the x- and y-values for each matrix entry. My
> current approach for using this data is to use an interpolation function
> from GSL, and call it for each cell while looping over them. This gives
> me the value at that position, interpolated.

You might also want to look at the InterpolatedTensorProductGridData and
InterpolatedUniformGridData classes. Since these are derived from
Function<dim>, you can send them into VectorTools::project if you want.

> Now, though, I would like to improve that approach. I was thinking about
> using the project-function and interpolating the matrix onto a component
> in my FESystem, but as far as I can see the project-function does not
> allow me to only change one component in my solution while leaving the
> other components untouched. Of course, I also could loop once over all
> cells before starting the main loop and solve the equation u = f, with f
> the value at each position, but I am not sure how efficient that is,
> especially when I can use the matrix-free vectorized approach for
> solving for the other components.

Why would you want to project when you can interpolate? That seems
unnecessarily expensive :-)

But regardless, if you can't find a way to deal with just component,
embed your data into a vector-valued function where the other components
are all zero, and then project/interpolate the whole vector-valued
function onto your FE space. You might be interested in the
VectorFunctionFromScalarFunctionObject class for this purpose.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

Am Mittwoch, 22. Januar 2020 19:23:52 UTC+1 schrieb Wolfgang Bangerth:

On 1/22/20 11:17 AM, 'Maxi Miller' via deal.II User Group wrote:
> As source term for calculations I am using data which is written in an
> external matrix (because it has been calculated by an external
> function), together with the x- and y-values for each matrix entry. My
> current approach for using this data is to use an interpolation function
> from GSL, and call it for each cell while looping over them. This gives
> me the value at that position, interpolated.

You might also want to look at the InterpolatedTensorProductGridData and
InterpolatedUniformGridData classes. Since these are derived from
Function<dim>, you can send them into VectorTools::project if you want.


> Now, though, I would like to improve that approach. I was thinking about
> using the project-function and interpolating the matrix onto a component
> in my FESystem, but as far as I can see the project-function does not
> allow me to only change one component in my solution while leaving the
> other components untouched. Of course, I also could loop once over all
> cells before starting the main loop and solve the equation u = f, with f
> the value at each position, but I am not sure how efficient that is,
> especially when I can use the matrix-free vectorized approach for
> solving for the other components.

Why would you want to project when you can interpolate? That seems
unnecessarily expensive :-)

I would like to have the input data (which is interpolated) as an additional component, which then is printed together with the result of the other calculations. Therefore I am not sure if there is a simpler approach (after interpolation) to get the data into the component. I either can solve the equation u = f for that, or use project(). Or is there another one?

[Wolfgang Bangerth]

On 1/23/20 7:29 AM, 'Maxi Miller' via deal.II User Group wrote:
>
> Why would you want to project when you can interpolate? That seems
> unnecessarily expensive :-)
>
> I would like to have the input data (which is interpolated) as an additional
> component, which then is printed together with the result of the other
> calculations. Therefore I am not sure if there is a simpler approach (after
> interpolation) to get the data into the component. I either can solve the
> equation u = f for that, or use project(). Or is there another one?

What I'm saying is that you can solve for u_h=f in many different ways,
depending on how exactly you interpret the equality. In general, u_h is in
some finite element space, but f is not, so the two will not be equal
everywhere. Projection finds a u_h so that
(u_h,v_h) = (f,v_h) for all v_h
which is one way of interpreting the equality. Interpolation finds a u_h so that
u_h(x_j) = f(x_j) at all node locations x_j
which is another way of interpreting the equality. Neither is better than the
other, but the second approach is *far* cheaper to compute.

[Maxi Miller]

I implemented it the following way:

    Functions::InterpolatedTensorProductGridData<dim> interpolator(this->coordinate_values,
                                                                         this->data_storage);

    VectorFunctionFromScalarFunctionObject<dim> interpolator_object(&interpolator.value,
                                                                    target_component,
                                                                    n_components_to_use);

    VectorTools::internal::project_matrix_free<n_components_to_use, fe_degree_to_use>(MappingQGeneric<dim>(1),
                                                                    dof_handler,
                                                                    hanging_node_constraints,
                                                                    INTERPOLATION_EQUATION<dim>(fe.degree + gauss_size),
                                                                    interpolator_object,
                                                                    present_solution,
                                                                    false,
                                                                    (dim > 1 ?
                                                                         INTERPOLATION_EQUATION<dim - 1>(fe_degree_to_use)
                                                                       : Quadrature<dim - 1>(0)),
                                                                    false);

but the compilation fails with

error: ISO C++ forbids taking the address of a bound member function to form a pointer to member function.  Say ‘&dealii::Functions::InterpolatedTensorProductGridData<2>::value’. Therefore, how can I approach that problem in a better way?

Thanks!

[Wolfgang Bangerth]

> |
> Functions::InterpolatedTensorProductGridData<dim> interpolator(this->coordinate_values,
> this->data_storage);
>
> VectorFunctionFromScalarFunctionObject<dim> interpolator_object(&interpolator.value,

The error...

> but the compilation fails with
> error: ISO C++ forbids taking the address of a bound member function to
> form a pointer to member function.  Say
> ‘&dealii::Functions::InterpolatedTensorProductGridData<2>::value’.

...comes from the line above. The problem is that
&object.functionname
is simply not something that C++ understands. What you need is a
function object that takes a Point and returns a double. You can use
std::bind(&Functions::InterpolatedTensorProductGridData<dim>::value,
&interpolator)
as this argument, or maybe more modern, write
[&interpolator] (const Point<dim> &p) { return interpolator.value(p); }
in place of the first constructor argument to interpolator_object.

[Maxi Miller]

Yes, the lambda worked, thanks!

[Maxi Miller]

After additional tests I noticed that I'm not only overwriting the target_component-component with the new values, but also the other components, but with 0 instead. But instead I do not want to touch them. How can I do that (if possible)?

Thanks!

Am Donnerstag, 27. Februar 2020 17:00:12 UTC+1 schrieb Wolfgang Bangerth:

[Wolfgang Bangerth]

On 2/28/20 12:49 AM, 'Maxi Miller' via deal.II User Group wrote:
> After additional tests I noticed that I'm not only overwriting the
> target_component-component with the new values, but also the other
> components, but with 0 instead. But instead I do not want to touch them.
> How can I do that (if possible)?

Let interpolate_b_v put its results into a temporary vector and then
only copy over those elements you are interested in. If you are using
block vectors, then that means just copying the block that matters.
Otherwise you'll have to identify for each unknown which component it
corresponds to (I think that there are functions for this in DoFTools)
and then copy only those elements.