Kratos: [StructuralMechanicsApplication] SpectralDecomposition and GaussSeidelEigenSystem Unexpected Results

[philipkalkbrenner]

Could it be that in constitutive_law_utilities -> SpectralDecomposition the eigenvectors and eigenvalues called from GaussSeidelEigenSystem in the math_utils.h come switched? So that the first entry of the eigen_values_matrix does not refer to the first entry of the eigen_vector_matrix.
eigen_values_matrix(0,0) has not eigen_vectors_matrix(0) as its corresponding eigenvector.
I am having some issues with that sice I am using SpectralDecomposition for my d+d- damage law for masonry.
@AlejandroCornejo told me to tag @marandra
Thanks,
Philip

[AlejandroCornejo]

@marandra you told me some days ago that there was some issue regarding this EigenSystems right? :)

[philipkalkbrenner]

I found something that could be related since SpectralDecomposition was changed:

4383

[AlejandroCornejo]

maybe @loumalouomega knows something!

[loumalouomega]

The previous version works exactly the same, the new version has altered the order (tranpose) to adapt to the criteria of @marandra , which is the one usually used in literature

[philipkalkbrenner]

ok,
well I have doubts using it, since I am calculating the decomposition of an easy 2d case stress state by hand and I compare it with the results obtained by SpectralDecomposition and the first entry of the eigenvectors is not the eigenvector which belongs to the first eigenvalue. Even if the transpose is used in SpectralDecomposition.

[loumalouomega]

Can you be more detailed?

[philipkalkbrenner]

I realized a problem while I was using the SpectralDecomposition of constitutive_law_utilities in the masonry 2d damage law. The returned stress_tensor_tension and stress_tensor_compression were not the same I obtained by an easy computation by hand (for a 2d case). So, I started to check the SpectralDecomposition and figured out that the eigen_vectors and eigen_values coming from the math utilities are not belonging together, e.g. if we take the Eigenvalue from position 1 of the eigen_values_matrix the corresponding Eigenvector is not contained at position 1 of the eigen_vectors_matrix.
And as soon as we compute the sigma_tension_tensor we use the function to compute the tensile part of the stress by considering that the eigenvalues and corresponding eigenvectors are at the same position (here it's i):
sigma_tension_tensor = eigen_values_matrix(i, i) * outer_prod(eigen_vectors_container[i], eigen_vectors_container[i]);

I hope the problem I have becomes clearer now...

[loumalouomega]

The old utility has exactly the same order, and is deprecated, the new one is transposed. But the methods are fine, are checked by several tests, the only problem can be the confusion of using it. The methods are documented in the code and explain the input/output

[philipkalkbrenner]

Ok but then I am still wondering why SpectralDecomposition is returning the wrong tensile and compressive parts of the stress tensor, at least for my 2D case, since I am just calling the function.
I put some KRATOS_WATCH in SpectralDecomposition and I enter with this stress tensor:
stress_tensor = (0.0 , 0.0 , -1549.17) ( in voigt)
the
eigen_values_matrix = [2,2] ( ( 1549.17,0 ) , ( 0 , -1549.17 ) )
and the
eigen_vectors_container = [ [2] (0.707107 , 0.707107) , [2] (-0.707107,0.707107) ]
all inside SpectralDecomposition.

so you can now check by hand, that the first eigen_vector in the container does not belong to the first eigenvalue of the matrix.
But the function pasted in my comment above assumes it.
The container should be
eigen_vectors_container = [ [2] (-0.707107 , 0.707107) , [2] (0.707107,0.707107) ] ,
just both the vectors switched.

[loumalouomega]

Can you check that if you multiply: A = LDL', if not then there is an error. I don't know how the tension/compression works, this was done by @AlejandroCornejo . I updated when I updated the function, but it was supposed to be consistent. In https://github.com/KratosMultiphysics/Kratos/pull/4383 we added the BDBtProductOperation operation, can you check that the methods you are using were changed consistently?

https://github.com/KratosMultiphysics/Kratos/blob/02f667e378bd346c2844b29ffde973cc244b5fec/kratos/utilities/math_utils.h#L1668 as it is on the documentation of the utility

[philipkalkbrenner]

I dont know to which tensors you refer by mentioning A, L, D

[loumalouomega]

I mean the matrix after the decomposition : A is your original matrix, D is the eigenvalues matrix, L is eigenvector matrix. I think that the documentation I attached explains it fine

[philipkalkbrenner]

I checked it, the matrices coming from GaussSeidelEigenSystem fullfill A=LDL'.
But later in SpectralDecomposition the assembly of the auxiliar_vector should be constructed like this to read the correct eigenvector:
line 782 in constitutive_law_utilities should be changed to this:
auxiliar_vector[j] = eigen_vectors_matrix(j,i);
before it was:
auxiliar_vector[j] = eigen_vectors_matrix(i, j);

[loumalouomega]

OK, but if the order has change the operation should change?

El vie., 29 mar. 2019 13:51, philipkalkbrenner [e-mail]
escribió:

I checked it, the matrices coming from GaussSeidelEigenSystem fullfill
A=LDL'.
But later in SpectralDecomposition the assembly of the auxiliar_vector
should be constructed like this to read the correct eigenvector:
line 782 in constitutive_law_utilities should be changed to this:
auxiliar_vector[j] = eigen_vectors_matrix(j,i);
before it was:
auxiliar_vector[j] = eigen_vectors_matrix(i, j);

—
You are receiving this because you were assigned.
Reply to this email directly, view it on GitHub
https://github.com/KratosMultiphysics/Kratos/issues/4528#issuecomment-477909522,
or mute the thread
https://github.com/notifications/unsubscribe-auth/ATEMgD1IdaVSMkf48j7KbEdfgwzKLdwKks5vbcw1gaJpZM4cG-XO
.

[loumalouomega]

I commented in the original PR, the change was done consistently, could be that the original order was wrong?

[philipkalkbrenner]

Could be that the original one was wrong. I just checked the current one and did not write the old one.

[loumalouomega]

That's what I suggested, in that case ping @AlejandroCornejo

[philipkalkbrenner]

I have opened a branch with the changes I did in SpectralDecomposition.
https://github.com/KratosMultiphysics/Kratos/tree/discussing_spectral_decomposition
@AlejandroCornejo could you check it?
I think before merging anything, the tests which use the SpectralDecomposition function should be checked. I commented out the tests for the masonry CLs (2d and 3d)
But as far as I know, your (@AlejandroCornejo) CLs are also using the SpectralDecomposition. Please add others which use this function.

[AlejandroCornejo]

I'll compile that branch and run some cases in order to check everything. THX

[AlejandroCornejo]

It seems that everything is okkey in my CL @philipkalkbrenner

[philipkalkbrenner]

Ok, then I will start the PR