FEAP User Forum
FEAP => nurbFEAP => Topic started by: resammc on September 26, 2018, 01:56:55 AM
-
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
-
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.
-
Thank you very much for the reference and also for the hint.
-
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,
-
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
-
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
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
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
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
-
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
-
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
-
It appears you are doing the right thing. Differences in solutions could be due to penetration differences.
What problem are you solving?
-
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
-
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.
-
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
-
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
-
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.
-
Thank you.
-
Resammc, thank you for sharing this information.
I have one more question on the contact formulation within Feap IGA.
I was testing simple models using basis functions of different orders. I noticed that starting with the 4th order contact (in all directions) for 3d case (celmt06) stops working.
Turned out that the subroutine celmt06, in this case, is not called. Only for csw=0. Is there any restriction for the basis functions in contact formulation?
Thank you in advance!
Best regards
-
It may not be the contact. In general the elements in feap are only coded to third order polynomials. It may be that the displacement elements can go higher but not mixed ones. So first try just the element you want to use and be sure it works up tot he order you want to test contact.
I will also check the contact to see if it has limitations -- since it was not originally coded by us.
-
Dear Prof. Taylor,
Thank you a lot for the hint. Indeed, in Feap 8.4 solid element is not working for the 4th order NURBS. However, in Feap 8.5 it works.
Unfortunately, I can't manage to use contact in 8.5. I tried the input file you provided in another topic http://feap.berkeley.edu/forum/index.php?topic=1504.0 and it worked for me.
I am trying a very simple task - a cube and a rectangular surface. I prescribe the displacement to one of the bodies as a proportional load. From time step 0 to 1 - it moves towards the second body. If I add a second time-step, without changing the load ( time step 1 to 2, the load is equal to 1), the second body finally deforms but in a weird way.
In Feap 8.4 it works. Maybe, I need to define the load or the BCs in a different way in the new version. I attach the input file.
Thank you in advance!
Best regards
-
I have located two bugs which I hope will fix the problem.
1. In cpair06.f add the line:
cp0(13,-1) = 0.0d0 ! force generation of surface connections
2. In c_shp2d_nurb.f correct the pointers to read
! Compute Bezier using extraction operator
point = np(289) + mr(np(273)+nknot(1)-1) + (p+1)**2*(ks(1)-1)
call derbezier1d(sg(1),hr(point), p, Nshp)
! call outc_e(hr(point),p+1, ks(1), nknot(1))
point = np(289) + mr(np(273)+nknot(2)-1) + (q+1)**2*(ks(2)-1)
I hope this will fix the problem
-
Dear Prof. Taylor,
It did solve the problem. Contact is working now. Thank you so much for your help!
Best regards
-
Dear all,
I've been trying to perform the contact analysis for a relatively complex geometry for a few weeks and my efforts were not successful. Therefore, I would like to ask for some hints/help in finding the problem.
Model description:
Two wires (screen_01.png) are put very close to each other. They do not touch each other and there is a gap between them (I have tried to show this in the attached figure, screen_zoomed.png), therefore I assume no contact in the beginning. Both of the wires are fixed using Dirichlet boundary conditions at their ends. A displacement shall be applied to one end of one of the wires, in order to get the wires in contact.
The cross section is circular and I model each individual wire with two patches, one core and one skin (therefore, in the input files, you see two MATErial tags).
What is the problem?
I cannot get the model to converge even when I have applied no displacements/loads!
My observations:
- Zero displacement applied (parameter di = 0): IGAFEAP does not converge when CONTACT is ON (non-zero residual). CONTACT OFF produces the expected results (zero residuals...).
- Displacement applied (parameter di != 0): IGAFEAP does not converge when CONTACT is ON even with very low Penalty values (Penalty=0 works, of course). NURBFEAP converges with CONTACT ON, but seems like the wires do not see each other and they pass through each other.
- The problem seems to be related to the somewhat complex geometry. When I try two straight wires (built with the same strategy as above), IGAFEAP detects the contact even with very low values of the Penalty (e.g. =1.0). I can provide the input files for those cases, if needed. Even with (very) coarse meshes, IGAFEAP gives an acceptable result as long as I keep the number of quadrature points high (e.g. 10 in each direction).
In general, this is not clear to me why IGAFEAP detects contact even when the wires are not close to each other. Next problem would be the convergence problem of the contact (i.e. even if it detects contact since the wires are within a defined range/gap, why the routine cannot converge?).
Any help in this regard or suggestions to look deeper into some parts of the code would be greatly appreciated.
Attached the input files (parameter DI controls the applied displacement value).
Many thanks and kind regards,
RM
-
Hi Resam,
you have an interesting question there.
I could run your input file and reproduced the error you obtained.
Using di=0 shows no convergence, despite no contact should be detected.
I can see many points showing erroneously penetration. The see figures attached show the penetration (gap) at every detection point.
(It is tricky to see the 3D fig, but one is a "normal to contact view" and the other one is the same results but from a "lateral" view)
No detection/contact is given by gap=0 (red color). As you can see, there is penetration defined for some points, especially for points away from the master surface.
If I have the time I will take a look at the arrays... but printing some og the gaps could be a good way to start debuging.
Good luck!
Axuni
-
Hi Axuni,
Many thanks for providing these informative figures. This definitely helps in finding the source of this problem. As you mentioned, I have to take a look at these gaps and the routines corresponding to them.
-
Here is my guess on what needs to be done:
I think the surface nodes for the contact have not been adjusted for the tie. I looked to see what happens for nurbs when you have a tie and I do not see any adjustment in the routine cksurf -- for any of the many contact types -- so this may not be the issue.
I suggest you prepare a small problem that has patches and check to see if the numbering is correct after the tie.
I will also work on making the outm work for nurbs -- it seems to be a challenge as the surface facets are storing the extraction numbers too. It may take some time to fix.
-
Dear Prof. Taylor,
Thank you very much for your valuable insights and also for looking into the problem. I will check to see if the numbering is correct or not after the tie command.
-
I just took a look that, if there is no TIE command, it still shows the same behavior.
It seems that TIE is not the (main?) problem here
-
What happens if you prepare a small problem that models the two wires and see if there is enough contact zone to pick up the interactions. The discretization needs to be similar to what the wires have. If the tieing is not the problem then it has to relate to the contact treatment. The current form is the Gauss-point form of De Lorenzis -- which may be sensitive to very small contact areas.
-
I tried it with a smaller problem and straight lines (the same wire discretization (the same cross-section radius also), but a straight line for the trajectory) .
It seems like the detection routine doesn't have any problems in this case. Attached the results for a 7x7 quadrature points distribution and a 10x10 one for the same nurbs model. (the nonsymmetric gaps are because of the loading).
-
Here are some possibilities:
1.There is some tolerance in the contact that is not working.
2. The problem is ill conditioned by the long, thin mesh of wires.
3. There is a programming bug still somewhere.
Do you have any other things that may cause the problem before I look more?
-
Many thanks for your help. At the moment I have nothing to add. Attached is the input file for the straight wires, in case it is needed. I will also keep investigating to find the source of this problem.
-
Can you describe how the wires are loaded? Is it from the end only? I note when running without the contact and loading by end displacement the wire bends only for a short distance from the end. I am concerned on conditioning. This may not be the contact problem but will make convergence hard.
-
Many thanks for your remark. This is actually a very small part of a more complex model and I have tried to isolate small parts of these wires. The boundary conditions (and material model used) are therefore artificial. At least one end of each wire is completely fixed and as the aspect ratio of the individual elements are also in an acceptable range, I was not expecting any problems with the conditioning. I think I tested this model (geometry and boundary conditions) with neo-hookean and finite deformations, and the deformed shape (also stress distribution) was better, but I didn’t want to add to the non-linearity of the problem (yet).