Contact in igafeap

[resammc]

Dear Prof. Taylor / Prof. Govindjee / FEAP Admin,

I'm beginning to work with the contact module in igafeap. I was wondering what is the reference (book, article, ...) used for the implementation.

Any information in this regard would be greatly appreciated.

Regards,

Resam

[Prof. R.L. Taylor]

The primary reference is:
@article{dimitri14,
author    = {R. Dimitri and L. De Lorenzis and M.A. Scott and P. Wriggers and R.L. Taylor and G. Zavarise},
title     = {Isogeometric large deformation frictionless contact using {T-splines}},
journal   = cmame,
volume    = {269},
year      = {2014},
pages     = {394--414} }

At some time this needs to be improved to use a "mortar" type surface treatment. 

[resammc]

Thank you very much for the reference and also for the hint.

[Axuni]

Hi Resam

We have been using contact in the last months. Actually I modified a bit your output subroutines in order to visualize the detection points.
The main thing was to include 'c_pair.h', take a look at the driver number ndrv, as well as putting whatever information we wanted (gap, status) into a module and calling those arrays in uparaview_nurb.f
Probably you have done it already. It would be a great addition to a further revision of your NPVI.

Cheers,

[resammc]

Hi Axuni,

This looks very interesting and helpful. Thanks a lot for sharing the information. Unfortunately, I have not done this yet. So, if you like, you can directly post the necessary changes in the NPVI post, or simply send them to me, and I will add it to the module. It will definitely save a significant time from us and the other users!

Regards,

Resam

[Axuni]

Hi Resam,
sorry for my late response.
Here you (and anybody else interested) have my proposal for adding contact output to NPVI.

I started creating a module (1) where I will store all the information. I am sure you can do this using palloc and ualloc. Remember that my implemenattion was only a prototype.
Then allocating the arrays in celmt06.f (2), where I found that csw=14 was doing the trick. Again, this can be placed in a more convenient place. In csw=3,6 I save the information I require, while I am still inside the loops.
Finally, in uparaview_nurbs.f (3), I added some variables in order to write down the output only when there is a valid contact drive.
At this point, I am pretty sure that you may come up with a more elegant solution, as you have a better overview of NPVI.
Sorry that I am only adding schemes of the changes, as my contact 3D elements have a lot of other stuff that I would have to clean if I wanted to send it to you.
However, for the uparaview file it would be far easier. So please let me know if you need it.

Remarks:
a. As the number of gauss points in a contact element can be defined arbitrarily, I implemented the output based on POINTS only, (VTU Type=2). I also tried outputing the elements as shells (type=23). The former option is robust and shows all output I require so far. The latter looks good, but only works for 9 points, and even then it is not using the middle point to interpolate results.
b. I tried creating a different PIECE in paraview where I could have only contact output. It did not work as expected, because of the way Paraview uses those pieces internally. Pitty. Therefore, I had to add dummy contact output to all other points, and dummy "solid" output to all contact points.
c. it would be great if you can give me a hand with pressure output

1. The module

Code: [Select]

      module contact_output
      implicit none
      save

      integer iCont_data      ! Flag for initialization of contact data
      integer c_pairs
      double precision,allocatable :: c_gpcds(:,:,:) ! coordinates of gauss points in contact elements (ndf,n_gp,c_pair)
      double precision,allocatable :: c_gpdsp(:,:,:) ! displacements of gauss points in contact elements (ndf,n_gp,c_pair)
      integer,allocatable :: c_status(:,:) !contact status at gauss point c_status(n_gp,c_pair)
      double precision,allocatable :: c_gn(:,:) !penetration at gauss point c_gn(n_gp,c_pair)
c     for IGA contact
      integer iCont_points   ! Flag for output as points (arbitrary number of points)
      integer c_nodes         ! number of nodes per element
      integer,allocatable :: c_elpair(:) !number of elements per contact pair
     
      data iCont_data/0/
      data iCont_points/1/
      data c_pairs/0/

      end module

2. A scheme of the changes in Celmt06.f

Code: [Select]

subroutine celmt06 (...)
use contact_output
[...]
case(3,6)
        cel_tmp  = 0
        c_elpair(npair) = neps_s(1)*neps_s(2)
        do i= 1,npair-1
          cel_tmp = cel_tmp + c_elpair(i)
        end do

        do ifac2 = 1,neps_s(2)

          n1(2) = k_order1(2) + 1

          do ifac1 = 1,neps_s(1)
            n1(1) = k_order1(1) + 1
            nng = 1
            do ng2 = 1,ngauss
              do ng1 = 1,ngauss

[...]

c               save data for output purposes
                if (iCont_points.eq.0) then
                  p = p_id(ng1,ng2)
                  if (ng1.eq.2 . AND. ng2.eq.2) goto 999 !Do nothing for the middle point
                else if (iCont_points.eq.1) then
                  p = nng
                end if !iCont_points

                cel_id = iifac + cel_tmp

                if (iCont_data.eq.1) then
                  do i=1,3
                    c_gpcds(i,p,cel_id) = xS(i)
                  end do
                  if (nng.eq.nngaus .AND. cel_id.eq.c_pairs) then
                    iCont_data = 2
                  end if
                end if
                do i=1,3
                  c_gpdsp(i,p,cel_id) = xS(i)-c_gpcds(i,p,cel_id)
                end do
                c_status(p,cel_id) = istgn !+istgt(p)
                c_gn(p,cel_id)     = ch2(p1(4),npt)!gn_int
                if (c_gn(p,cel_id).gt.0) c_gn(p,cel_id) = 0.d0  !write zero when contact

[...]

      case(14)

[...]

c       allocate variables for user output
        if (iCont_data.eq.0) then
          c_pairs = c_pairs + neps_s(1)*neps_s(2)
          if (npair.eq.numcp) then
            if (iCont_points.eq.1) then   ! for output as points only (default)
                c_nodes = ngauss*ngauss
            end if
            allocate ( c_elpair(numcp) )
            allocate ( c_gpcds(3,c_nodes,c_pairs) )
            allocate ( c_gpdsp(3,c_nodes,c_pairs) )
            allocate ( c_status(c_nodes,c_pairs) )
            allocate ( c_gn(c_nodes,c_pairs) )
            c_elpair(:)  = 0
            c_gpcds(:,:,:) = 0.0d0
            c_gpdsp(:,:,:) = 0.0d0
            c_status(:,:)  = 0
            iCont_data=1
          end if !numcp
        end if !iCont


3. A scheme of changes in uparaview_nurbs.f

Code: [Select]

      subroutine uparaview_nurb(...)
      use contact_output
[...]
      if (ndrv.eq.2 .OR. ndrv.eq.6) then
        write(plu,1020) numnp+c_nodes*c_pairs,numel+c_pairs !add contact elements
      else
        c_pairs = 0
        write(plu,1020) numnp,numel               ! Start Mesh/Piece Section
      end if
[...]
      write(plu,1010) '<Points>'                ! Start Point/Node data
      write(plu,1030) 'Float64','nodes',3
      do i = 1,numnp
        write(plu,2010) (x(ii,i),ii = 1,ndm) ,(0.0d0,ii = ndm+1,3)
        write(plu,*)
      end do ! i
      if (ndrv.eq.2 .OR. ndrv.eq.06) then
        do j=1,c_pairs
          do i = 1,c_nodes
            write(plu,2010) (c_gpcds(ii,i,j),ii = 1,ndm)
          end do ! i
        end do !j
      end if
      write(plu,1010) '</DataArray>'            ! Close Node data
      write(plu,1010) '</Points>'               ! Close Points section
[...]
        write(plu,1030) 'Float64','Principal Stress',7 ! Start Prin Stresses
        do i = 1,numnp
          do ii =1,7
            write(plu,2000) psig(i,ii)
          end do ! ii
          write(plu,*)
        end do ! i
        if (ndrv.eq.2 .OR. ndrv.eq.6) then !dummy
          do j=1,c_pairs
            do i = 1,c_nodes
              do ii = 1,7
              write(plu,2000) 0.d0
              end do ! ii
            end do ! i
          end do
        end if
        write(plu,1010) '</DataArray>'               ! Close Prin Stresses
[...]
c     Output contact variables
c ======================================================================
      if (ndrv.eq.2 .OR. ndrv.eq.6) then
        write(plu,1030) 'Int32','C_Status',1       ! Start contact status
        do i = 1,numnp
          write(plu,2020) 0 !dummy
        end do ! i
        do ii = 1,c_pairs
          do i=1,c_nodes
            write(plu,2020) c_status(i,ii)
          end do
        end do
        write(plu,1010) '</DataArray>'            ! Close contact status
       
        write(plu,1030) 'Float64','C_gap',1       ! Start contact gap
        do i = 1,numnp
          write(plu,2000) 0.d0
        end do ! i
        do ii = 1,c_pairs
          do i=1,c_nodes
            write(plu,2000) c_gn(i,ii)
          end do
        end do
        write(plu,1010) '</DataArray>'            ! Close contact gap
      end if


[resammc]

Dear Axuni,

Thanks a lot for sharing your work. I will go through it (asap). Unfortunately, I still have not got time to look deep into the contact module routines.

Best,

Resam

[1user]

Dear all,

Here is another question regarding the contact formulation in IGA. I need to plot the contact pressure at the detection point. So, I assume I should print the value of the normal force (lambda_hint) divided by the Jacobian (ch3(p3(1),npt)). However, the values I get are a bit lower than in FEM. I would be glad to get any information here.

Best regards

[Prof. R.L. Taylor]

It appears you are doing the right thing.  Differences in solutions could be due to penetration differences. 
What problem are  you solving?

[1user]

Dear Prof. Taylor,

I am really grateful for your response.
My aim was to observe the contact pressure distribution over the contact region and on the border to compare it with the oscillating results in FEM with 2nd order Lagrangian elements (Gauss point to the rigid surface formulation). 
I attach my input file for IGA.
I checked my gap values. They are matching for FEM and IGA. I also get the same reaction force.
I am wondering if I should use some additional multiplier for the area in IGA. A lot of the are calculated in the routine cont_normint06.f. Meanwhile I use only value of the Jacobian (ch3(p3(1),npt)).

Thank you in advance!
Best regards

[Prof. R.L. Taylor]

I am not sure what you want us to do with the files.  They contain commands you have developed and also refer to a user material model for the contact.

[1user]

Dear Prof. Taylor,

I am sorry for the confusion. I just wanted to demonstrate the geometry of the task.
My question to the developers of the program was of course about the area of the integration point.
I got confused with the nature of the coefficient fac3. All in all, I just need to be sure that I divide my normal force obtained in the integration point (penalty coefficient * gap) by the correct area.
If there is any articles/reports I could get more information from it would be also useful.

Thank you for your help!
Best regards

[1user]

Dear all,

I know understand my mistake. Jacobian is the overall area of the element here and fac3 is the area of the particular detection point.

Thank you for your time!
Best regards

[resammc]

I think, in c_surfnurb2d.f, line 130:

        if(fac.eq.2) then
          jj = len1*len2 - len2 + 1
        else
          jj = 1
        endif

should be changed to:

        if(fac.eq.2) then
          jj = len1*len2 - len1+ 1
        else
          jj = 1
        endif

or simply

        if(fac.eq.2) then
          jj = len1*(len2 - 1) + 1
        else
          jj = 1
        endif

as the difference between face 5 and face 2 is in the 2nd direction (knot vector 2) and therefore one has to pass all the sides from 1 to len1*(len2-1) to reach the first side of the face 2.

[Prof. R.L. Taylor]

Thank you.