Diffrenet results for user element and feap element

[nima_ikm]

Hello,

I wrote a simple user element for triangle element in linear quasi static environment. Unfortunately, my result is different from FEAP element.  I used the same input file; IUSER for user element and IFEAP for FEAP element, actually they are same.
However, I just changed my user element to quadrilateral element, the results are exactly same with feap element (it means my user element for quadrilaterla is working fine ). As result, I think the numerical integration, probably I used for triangle element is wrong. Could you please give me some feedbackS on it.All the files are attached.

Thank you so much, in advance.

[Prof. R.L. Taylor]

Please tell us which version of feap you are using when you post a question.  It helps to know of the problem you are experiencing is due to a known bug in the program.

It is hard to tell why your element does not work for the triangle.  You can print the result you obtain and compare to the value feap computes.  EPRI will show the last stiffness after a tangent.  You set 'p' for each quadrature point, this can be moved outside the loop -- of course there is only one point for your element.  One point should compute the exact stiffness for a plane problem.  Be careful how you form the matrices to ensure that if your input data has more than 2 dof/node specified it still works.  You tacitly assume it will always be '2', but if you used with a beam where 'ndf = 3' it would always fail to work.

You do not need to initialize arrays such as s and p under ISW = 3, feap will do it.

[nima_ikm]

Dear Prof. Taylor,
Thank you very much for your kind help and answer.
Based on EPRInt, the stiffness for last element is different for FEAP element and user element. Also, I am using FEAP 8.4 in Linux environment. Also, I attached the user element for quadrilateral which work fine and exact with FEAP element.

triangle Element                                                                Quad Element                                  Comments
  ++++++++++++++++                                   ++++++++++++++++++        ++++++++++++++++                             
 call TINT2D ( 2, LINT, GP )                                         call INT2D ( -3, LINT, GP )                     !Gauss Quadrature
 do ig=1,LINT                                                               do ig=1,LINT
         SS=GP(1:3,ig);                                                       SS=GP(1:2,ig);                                 ! Integration points         
          W=GP(4,ig);                                                           W=GP(3,ig);                                      ! Weights
          CALL TRISHP ( SS, XL, NDM, NEL, XJAC, shp )     CALL SHP2D ( SS, XL, shp, XJAC, NDM, NEL, IX, FLG )
          (….)                                                                         (….)
 Enddo                                                                         Enddo
d.o.f=6                                                                        d.o.f=8
GP(4,3)                                                                       GP(3,4)                                                    !allocation of Guass points

However, for numerical integration, I used as follows (which are same in both Quad and Triangle):
    
S=S+transpose(B)*C*B*W*det(J)

However, I am not sure this numerical integration also is correct for the Triangle (maybe I sould normalaize W by multiplying to 0.5- but still it doesn't work).

P.S. I extend it even to 3D brick element and Tet element it works fine. Then I really confused why it doesn't work.

Thank you for your help.                                                                               
                                                                               

[Prof. R.L. Taylor]

You have NEL in the argument of the call to TRISHP.  It should be the order of the element, which is "1" a linear triangle.

There are some other errors that my compiler complained about.  But I downloaded your original element so maybe you fixed them - e.g., TGP not defined integer.

[nima_ikm]

Dear Professor,

IT IS SOLVED; thank you. I forgot to mention the order of triangle (NEL), as exactly you mentioned it.

I am really appreciate.