SimVascular / svFSI

A multi-physics finite element solver for patient-specific blood flow simulation including fluid-structure interaction and cardiac electrophysiology
Other
31 stars 49 forks source link

Problems with only minor symmetry stiffness matrix? #39

Open mrp089 opened 3 years ago

mrp089 commented 3 years ago

Part of #32. Our G&R model creates a tangential stiffness matrix that only has minor and no major symmetries. It thus has 36 unique entries (instead of 21 with major symmetries). I extended CCTOVOIGT to handle that in https://github.com/mrp089/svFSI/commit/5796f6e9b95ba0d597876ea514e048bb960f258d.

@vvedula22 Are you aware of any other places where svFSI might assume a stiffness matrix with major symmetries?

vvedula22 commented 3 years ago

@mrp089 I don't recall all instances. Do a simple `grep -i cctovoigt' in the folder to search for all function calls and add the extra logical argument to them.

mrp089 commented 3 years ago

It looks like there are no other uses of CCTOVOIGT. I was wondering if we make any assumptions after the stiffness tensor has been converted to a 6x6 matrix.

vvedula22 commented 3 years ago

Yes, the matrix is a 6x6 assuming major and minor symmetries. Lack of major symmetry is not a major problem but lack of minor symmetry can pose issues. I am attaching a couple of snapshots of how the material stiffness tensor is assembled to the virtual work equation

Screenshot from 2021-07-15 20-28-09

Screenshot from 2021-07-15 20-27-02

As you see, when you don't have major symmetry, you need to assemble the full **C** matrix (or Dm in the code). Make sure that the entries **C**_{ABCD} in the matrix are correctly located. Eventually, you should make sure that the component-wise product of the variation in virtual E, C matrix, and the variation in E, should make sense and evaluate to a scalar.

mrp089 commented 3 years ago

Thanks for the overview!

Our Dm still has minor symmetries, i.e. 36 unique entries. I wrote code in https://github.com/mrp089/svFSI/commit/5796f6e9b95ba0d597876ea514e048bb960f258d to assemble the lower triangle into Dm (instead of copying it from the upper triangle as with major symmetry). So as far as Dm is concerned, it always supported minor-symmetry-only problems due to its 6x6 dimension, it just assumed major symmetry during assembly.

My concern is that other parts of the code "downstream" still assume that Dm has major symmetries (e.g. multiplications, linear solver, ...).

vvedula22 commented 3 years ago

Dm is used to assemble the local stiffness matrix in STRUCT.f and USTRUCT.. We also use it in IB.f and POST.f but those shouldn't be a concern for you now. There is no other place where you make any assumption on the structure of the matrix.

From my previous comment, once you assemble the full 6x6 Dm matrix, you need to check the expressions for Bt D B in the following subourines:

subroutines STRUCT3D and STRUCT2D from STRUCT.f

subroutines USTRUCT3D_M and USTRUCT2D_M from USTRUCT.f`

and make sure that there are correctly applied. For e.g., check if the correct row-column component of Bm is acting on D and vice-versa.

From my understanding, simply creating a full Dm matrix should be sufficient. But you may verify the expressions just to make sure.

vvedula22 commented 3 years ago

I looked at your commit mrp089@5796f6e and it appears that the 3D assembly for Dm is correct. However, the 2D one looks incorrect to me. It should be,

            Dm(2,1) = CC(2,2,1,1)
            Dm(3,1) = CC(1,2,1,1)
            Dm(3,2) = CC(1,2,2,2)

Please verify.

mrp089 commented 3 years ago

Thanks for checking @vvedula22, you're absolutely correct! And thank you for pointing out the lines for the multiplications. I'll check if they generalize for a full Dm.