Element eigenvalues

[Colin McAuliffe]

Dear all,

I am writing a umacr to compute element eigenpairs for our unsymmetric user elements. The umacr follows the same procedure as the code for the command EIGE, for symmetric elements:

         do el = 1,numel
             call formfe(np(40),np(26),np(26),np(26),.false.,.false.
     &            ,.false.,.false.,3,el,el,1)
             call elmteig(hr(np(36)),storel,el)
         enddo

Where the subroutine elmteig computes the eigenvalues of the local tangent using the lapack routine dgeev. For now I am testing the routine on a mesh with 1 element using either a symmetric feap element or one of our unsymmetric user elements and I have noticed the following strange behavior:

1. If I use either the feap command EIGE or my element eigenvalue  routine without first going through the normal solution commands, TANG/UTAN,,1, the eigenvalues turn out wrong compared to exporting the tangent and computing the eigenvalues in matlab.

2. Applying boundary conditions changes the eigenvalues of the local tangent if parfeap is used, but not in serial feap.

Any help in explaining either of these would be appreciated.

All the best
Colin

[FEAP_Admin]

For Q1:  Is your problem non-linear?  If so solution steps will alter your tangent.  To debug further examine what is contained in 's'
after the call to elmlib( ) in program/pform.f [serial feap] and compare this to what is being sent to peige( ) in pmacr3.f [again just testing serial FEAP with the EIGE command.]

For Q2: the issue is that in parFEAP we assemble a modified version of the stiffness matrix where prescribed displacement dofs get their row and columns zeroed out and a '1' is placed on the diagonal -- likewise the right-hand side entry in 'p' is set to the prescribed values; i.e. we assemble the entire matrix 'nneq' in size with all dofs (unknown as well as known).  The routines that are doing this are parfeap/smodify.f and parfeap/pmodify.f.  I believe this only occurs when you use blocked equations.  For unblocked equations smodify( ) gets skipped in parfeap/pform.f.

[Colin McAuliffe]

1. The same thing happens for a linear problem, I attached a sample input file. The s that sld2d1 gives is different than that the one that peige gets if I do not have the TANG command. If the TANG command is used s is the same in both routines. I attached the tangents given by these routines as well. This is not really crucial to figure out since we will nearly always be using the eigenvalue command after a UTAN when using our nonlinear elements, I am just trying to understand what feap is doing here.

2. Ok thanks, it should be strait forward for me to make a switch to skip smodify when using the umacr

[Prof. R.L. Taylor]

Re: Q1:

Your input file is requesting a transient solution.  There could be a 'bug' in the element routine as it may not know what to do about the transient term.  Do you expect to get the eigenvalues of the stress linearization?  In a usual solution mode the program computes: S = c1*K + c2*C + c3*M where K,C,M are stiffness, damping mass.  We need to check now what is done when you ask for 'S' using 'eige' -- probably there are some missing flags that mean you are getting the term augmented by the mass.  We need to check more what is happening before we can give a complete answer.  You may also be able to follow the path from 'pmacr3.f' where the [eige] is executed to see what is happening by the time you get to where the element stiffness is called.  The call path should be 'formfe' then 'pform' -- which may explain why parfeap is different too.

Thanks for the question --

rlt

[Prof. R.L. Taylor]

Can you try the following:

In pmacr3.f (/program directory) before the call to 'formfe' in the [eige] command add the following statements.

        hflgu = .false.
        h3flg = .false.
        ctan(1) = 1.0d0
        ctan(2) = 0.0d0
        ctan(3) = 0.0d0

You will need to also add an include statement:

        include   'eltran.h'

Let us know if this corrects your program.  We still need to do some additional checks but it should fix the stiffness eigenvalues (they will not include any terms from 'C' or 'M' type terms and will not disturb the history variables).

rlt

[Prof. R.L. Taylor]

Sorry it is 'h3flgu'  not 'h3flg'

[Colin McAuliffe]

Ok interesting. When I comment out the TANG command, in pmacr3.f at eige, ctan(1:3) is 1.0,0.0,0.0 so it is computing eigenvalues of the stiffness only. On the other hand if the TANG is included, ctan(2) and (3) are set and so that hr(np(36)) contains M+C+K at the time of eigenvalue computation. So I guess values for ctan are set at some point during the call to TANG/UTAN?

Also, regarding question 2, I tried an example with unblocked format and there is still a difference between the stiffness matrix given by feap vs parfeap for elements containing restrained dofs, so I don't think smodify was the issue after all. I suppose the difference will be somewhere in program/pform.f vs parfeap/pform.f ?

Thanks