Re: [Libmesh-users] Initialize DenseSubMatrix with m x n (non-const) size

[John P. <jwp...@gm...>]

On Tue, Aug 10, 2021 at 7:34 PM edgar <edg...@cr...> wrote:

> Hello!
>
> I am trying to do something like this (systems of equations, example 6):
>
>      DenseSubMatrix<Number> Ke_var[3][3] =
>        {
>          {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
> DenseSubMatrix<Number>(Ke)},
>          {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
> DenseSubMatrix<Number>(Ke)},
>          {DenseSubMatrix<Number>(Ke), DenseSubMatrix<Number>(Ke),
> DenseSubMatrix<Number>(Ke)}
>        };
>
> Except that I would like to do it like this (does not work):
>
>      const unsigned int mesh_dim = mesh.mesh_dimension ();
>      DenseSubMatrix<Number> Ke_var[mesh_dim][mesh_dim];
>      for (unsigned int i=0; i<mesh_dim; i++)
>          for (unsigned int j=0; j<mesh_dim; i++)
>              Ke_var[i][j] = DenseSubMatrix<Number>(Ke);
>
> - One of the reasons for which this does not work is that there is no
> empty constructor defined for =DenseSubMatrix<T>=. In other words, the
> second line in the above code requires parametrized initialization.
> - I considered using double pointers (**) and =new= (dynamic
> allocation), but that would be a bomb for memory.
> - Since =mesh_dim= is =const=, I think that even the ISO C++ would be
> fine with =Ke_var[mesh_dim][mesh_dim]=.
> - I also thought of playing with =new_m= and =new_n= from the definition
> of the constructor (something like =DenseSubMatrix<Number> Ke_var (Ke,
> 0, 0, mesh_dim, mesh_dim)=), but I don't know how I would do that.
>
> The point is to be able to generate the matrix with specific dimensions.
> Can you help me, please? Thanks!
>

I don't know of a way to achieve this using 2D arrays, but if you are
willing to use a std::vector of DenseSubMatrices, you can do something like:

std::vector<DenseSubMatrix<Number>> Ke_vec(mesh_dim*mesh_dim,
DenseSubMatrix<Number>(Ke));

You could then get the (i,j) entry at index k = mesh_dim*j + i. And of
course you could wrap all this up in a helper class if you wanted the
syntax to be nicer.

--
John