0

I am using the Eigen library to solve a FEM problem in which I use several similar sparse matrices for the different kinds of derivatives I calculate. To build the sparse matrix to solve the system I would like to use the comma initializer, but that is not supported for sparse matrices (https://eigen.tuxfamily.org/dox/group__TutorialSparse.html and https://stackoverflow.com/a/43532618/15811117). In that answer, Henri Menke does suggest the Eigen unsupported Kronecker product, though I don't think that will work here.

With an NxM mesh, I calculate several derivative matrices (which are sparse), Dx2, Dxz, and Dz2 of dimension N*M x N*M. My solver matrix is as follows, (K's are just constants, and 0 is just filler--a zero matrix, N*M x N*M)

K3*Dx2 + K1*Dz2   K2*Dxz            0              0                 0
K2*Dxz            K3*Dz2 + K1*Dx2   0              0                 0
0                 0                 K1*(Dx2+Dz2)   0                 0
0                 0                 0              K3*Dx2 + K1*Dz2   K2*Dxz
0                 0                 0              K2*Dxz            K3*Dz2 + K1*Dx2

As of right now, I build it by converting the sparse into dense, using the comma initializer, and then converting back to sparse by .sparseView(). The only problem with this is that if the mesh is much larger than 16x16, I get a std::bad_alloc error (which is reasonably expected, since it creates a huge matrix: 5*N*M x 5*N*M), and I need a mesh at least 100x100 (and there will not always be as many 0 matrices as currently shown, though it will remain a sparse matrix).

What is the best way to build this matrix? I have thought of trying to use Triplets to do so, similar to this (Convert an Eigen matrix to Triplet form C++).

Edit: I don't think that using Triplets will be that easy (though I am currently attempting it) because I edit several of the values in each D_ _ matrix for boundary conditions...unless there is a way to find a Triplet in an std::vector that has given row and column number?

Update: Both the Triplet and the Kronecker product methods are feasible, but I have run into several problems:

With the Kronecker product, I get the following error (about a hundred times or so) at compilation:

.../eigen/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h:225:97: error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<int, std::complex<double>, Eigen::internal::scalar_product_op<int, std::complex<double> > >’

With the Triplet method, I defined some operators and a function to 'shift' the matrix to the right place (basically what the Kronecker product is doing). I just have to figure out an index out of bounds issue.

Izek H
  • 39
  • 5
  • Are the derivative matrices sparse? – CJR Jul 13 '21 at 12:16
  • @CJR yes, they are sparse; they were originally made using the setFromTriplets function. I'm now trying to see if I can keep them as vectors of triplets until I build this matrix. – Izek H Jul 13 '21 at 16:12
  • 2
    Why do you think the Kronecker product is not a viable solution? With the Kronecker product you can create the matrices obtained repeating Dx2, Dxz and Dz2 and then sum them to obtain the final matrix. – LGrementieri Jul 13 '21 at 20:49
  • @LGrementieri I guess I misunderstood how to use the Kronecker product--I agree, it should work. Sadly, I get hundreds of errors when I compile :( See update. – Izek H Jul 13 '21 at 21:59

1 Answers1

0

The problem comes from the fact that the matrices I used in the Kronecker product were of type SparseMatrix<int> and SparseMatrix<complex<double>>. There are no operator overloads for int and complex<double>, so by making both matrices of compatible types, the code compiles and runs as expected.

Izek H
  • 39
  • 5