Sliding with gaps and friction parameters, INF and NAN issues
issueid=563 07-10-2012 06:13 AM
Sliding with gaps and friction parameters, INF and NAN issues
Sliding with gaps and friction parameters, INF and NAN issues

I'm trying to implement friction in a model of indentation of soft tissue. I used the following contact parameters:

<contact type="sliding_with_gaps">
<laugon>0</laugon>
<tolerance>1</tolerance>
<penalty>1</penalty>
<two_pass>0</two_pass>
<auto_penalty>1</auto_penalty>
<fric_coeff>0.57</fric_coeff>
<fric_penalty>1</fric_penalty>
<search_tol>0.01</search_tol>
<seg_up>0</seg_up>
<ktmult>1</ktmult>

I also tried:
<contact type="sliding_with_gaps">
<laugon>0</laugon>
<tolerance>1</tolerance>
<penalty>10</penalty>
<two_pass>0</two_pass>
<auto_penalty>1</auto_penalty>
<fric_coeff>0.57</fric_coeff>
<fric_penalty>10</fric_penalty>
<search_tol>0.01</search_tol>

For the attached model FEBio produces INF and NAN in the convergence norms and in PostView only the faces that contain rigidly supported nodes are visible.

What is the cause for these errors? Am I using the wrong type of parameters?

I am actually interested in creating a glued interfaces such that no sliding is allowed to occur between the indentor and the skin surface once contact is occurs. How do I go about implementing that? Is this the correct approach?


Kevin
Issue Details
Issue Number 563
Project FEBio
Category Unknown
Status Fixed
Priority Unknown
Affected Version Unknown
Fixed Version (none)
Users able to reproduce bug 0
Users unable to reproduce bug 0
Assigned Users (none)
Tags (none)




07-17-2012 02:56 AM
I am really stuck with this.

Does anybody know the cause for these errors? Am I using the wrong contact parameters?



Kevin

07-20-2012 02:36 PM
Lead Code Developer
Hi Kevin,

Have you tried turning of the friction? I wonder how the model runs whithout friction. Also, try turning of the autopenalty. Your indentor is a rigid body and the autopenalty doesn't always work well with rigid bodies. I will try to take another look at the model early next week and let you know if have any insights.

Cheers,

Steve.

07-24-2012 04:21 AM
Hey Steve,

Yeah the model runs fine without friction. Unfortunately I really need to implement friction (currently there is a lot of unrealistic sliding occurring under the indentor).

I based the auto penalty option on this forum discussion:
http://mrlforums.sci.utah.edu/forums...ontact-problem

However that was without friction. I tried to run the model without the auto penalty option but that does not change anything and I still have the NAN and INF issues (see attached temp_contact4.feb model).

Thanks, for your help. I hope we can fix this soon,



Kevin

08-02-2012 04:54 AM
Hey Steve,

I am still stuck with this so my FEBio work is currently on hold.

Did you have a look at this model yet? Any thoughts on fixing the friction issues?

Thanks,


Kevin

08-20-2012 03:26 AM
Does anybody have any feedback on this?

Kevin

08-22-2012 06:35 AM
Moderator
Hi Kevin,

I was able to run your temp_contact4.zip file without encountering this NAN or INF issue. The analysis ran until time 0.1813 then terminated due to a negative jacobian. I am using the latest development version of FEBio (which I believe you have), running it on a Mac. Have you tried running this model on the latest version?

Best,
Gerard

08-22-2012 08:47 AM
Hi Gerard,

Thanks for looking into this. I tried that model (temp_contact4) again this time with the latest source code (1.5.1.3711) which I compiled in Windows 7 64-bit and it now fails at 0.137768 as it detects the NAN's.



Kevin

08-22-2012 01:30 PM
Moderator
Hi Kevin,

I am not sure why we are getting different outcomes, but I looked more carefully at the temp_contact4 model and here are my suggestions, which produce a much better outcome:

1) For the Mooney-Rivlin material you are using c1=0.197 and k=0.1. Generally speaking, to enforce near-incompressibility, k should be much larger than c1 or c2. I suggest k=200.
2) For the penalty parameters in the contact interface, you are using penalty=1 and fric_penalty=1. Your shell thickness is set to 0.91 and the elements under the shell layer have a thickness of ~2. So it makes more sense to use a penalty value of ~c1/thickness. I suggest penalty=fric_penalty=0.1
3) For some reason some of your shell elements have zero thickness. I suggest specifying a non-zero thickness (I searched and replaced 0 with 0.91).

With these changes, I get the analysis to converge much more rapidly, up to time t=0.4766. The analysis terminates on a different error now:

Assertion failed: (TMT >= 0), function ContactNodalForce, file FESlidingInterface.cpp, line 818.

This error means that the magnitude of the frictional traction is not positive. It would take me a little while to debug this error, but in the meantime why don't you try these settings on your end and see what you get.

Best,

Gerard

08-27-2012 09:16 AM
Hey Gerard,

Thanks for this! I addressed the points you raised but still get the errors, also at t=0.4766. I had not set the thickness for the rigid indentor part. In reality that was a solid part but I model it as a shell for simplicity. What thickness would I set for the rigid indentor? Currently you've set it the same as the skin layer, I could do the same but I hope the results are independent of the choice of the thickness.

If I run your temp_contact4b.feb model then it also runs until t=0.4766. However I obtain a NaN error (I take it you were using debug mode) and the log file shows this:

convergence norms : INITIAL CURRENT REQUIRED
residual 9.281968e+000 -1.#IND00e+000 0.000000e+000
energy 9.780373e+000 1.#QNAN0e+000 9.780373e-002
displacement 5.242316e+002 2.010453e+005 1.890978e-001
6
*************************************************************************
* FATAL ERROR *
* *
* NAN Detected. Run aborted. *
*************************************************************************

Let me know if you figure anything out or if you have more tips for me.


Kevin

11-14-2012 09:39 AM
This is still an issue for me. Are there any updates on the NAN errors?

I get the impression it is due to either the bulk modulus choices or the ktmult option.

Attached are some files which fail or run successfully depending on the choice of these parameters.


Kevin

11-19-2012 07:03 AM
I've attached another model which worked initially but then I increased the friction coefficient (to 1) and now I get NaN errors suddenly at around 50%.

These are the contact settings:

<contact type="sliding_with_gaps">
<penalty>2.500000e-02</penalty>
<auto_penalty>0.000000e+00</auto_penalty>
<two_pass>0.000000e+00</two_pass>
<laugon>0.000000e+00</laugon>
<tolerance>1.000000e+00</tolerance>
<gaptol>0.000000e+00</gaptol>
<minaug>0.000000e+00</minaug>
<maxaug>1.000000e+01</maxaug>
<fric_coeff>1.000000e+00</fric_coeff>
<fric_penalty>2.500000e-02</fric_penalty>
<seg_up>0.000000e+00</seg_up>
<search_tol>1.000000e-02</search_tol>





Kevin

11-19-2012 07:19 AM
NaN's no longer occur on the NAN_at_50.zip model if I reduce the penalty parameters. However now I get a lot of penetration.

Hope that somebody can help me learn what the source of the NaN errors is.



Kevin

11-24-2012 07:15 AM
Moderator
Hi Kevin,

The source of the NaN error is that an internal calculation of the metric tensor for one of the patches on the contact surface produces a determinant of zero. A zero determinant would imply that the surface collapsed to a curve (an unexpected outcome), causing the NaN in subsequent calculations. There is no clear physical reason for that to occur, and it is surprising that the determinant is exactly zero, so we are a bit mystified by this outcome.

As far as we can tell, there is no bug in the code that would cause this specific outcome. This is an unusual circumstance that, so far, we only seem to encounter in this specific problem that you are trying to solve. If you could reproduce this error in a simpler contact geometry with fewer elements, that would be extremely helpful to help us determine whether the problem is a bug in the code or a problem with the model (maybe the mesh?). This analysis takes so long to run in the debugger that we are not able to play with multiple scenarios to find the bug, if there is one.

Best,

Gerard

11-26-2012 03:27 AM
Hi Gerards,

Thanks for your reply.

Do you have a patch or element ID? Do you know if it is a patch on the master or the slave surface?

I think you might be right that this is specific to this case. I create this model and its .feb file directly from MRI data using Matlab. So far simple models have not produced this error so distilling this to a smaller model might be difficult. But perhaps the bug is actually in one of my Matlab functions. However what puzzles me is that for quite a number of parameter choices the model runs fine, then suddenly when a certain parameter is altered it produces the NaN error. So if it is the model layout or the mesh you'd expect it to always produce the error independent of the parameters right? (I need flexibility in the parameters because this is for inverse FEA optimisation towards MRI derived deformations&force)

Also I've tried setting the contact definitions by hand instead using PreView but then it still produces NaN errors. Therefore it does not seem to be my patch/contact definitions in Matlab (but it could still be a problem of the mesh which at present PreView does not detect).

I am going to check the mesh generator I created and also the export function to .FEB files (the mesh seems fine but perhaps there is a bug in my .feb file export function, although the jacobians seem fine in postview for the model...). It would help to know what element/patch or part we think is responsible.

Could you explain how you found out about this error information? If I want to study the problem in the same way, how do I go about doing that? Do I simply install the debug source code and run it in debug mode and study the error message/the line where it gets stuck?

Thanks a lot,


Kevin

11-28-2012 07:18 AM
Hi Gerard,

I've attached a faster model (with a terribly coarse mesh) called NAN_quick, which produces the NaN error at t=0.3195 and takes about 1 min 46 sec to get there.

I've studied the meshing algorithms I created in Matlab and studied the areas of the elements and it seems fine to me. Do you have any information on what face is responsible? I suppose it would have to be on the slave surface since the master is rigid and can therefore not collapse to a curve as you say. I think the mesh starts out fine (FEBio probably would produce a warning if it didn't) so not sure what to look for at this point.

Hope the faster model helps, thanks,


Kevin

11-28-2012 10:39 PM
Moderator
Hi Kevin,

Thanks for sending this reduced-size model. It was very helpful.

When I ran it with the parameters you selected, I did eventually encounter a NAN exception. However, I noticed that it occurred after very many iterations, when it seemed that the solution would not converge anyway (whether or not NAN was encountered). This made me realize that my last response to you did not address the real problem.

I reran this NAN_quick problem with a friction coefficient of zero and it converged very nicely until t=1. However, I did notice that the indenter had a large penetration into the tissue, indicating that the penalty factor you used (0.01, with auto_penalty turned off) might have been too small. I restored the friction coefficient to 0.9, switched auto_penalty on, and set the scale factor for the penalty to 1 (for both <penalty> and <fric_penalty>) and the analysis also converged nicely to t=1, but with a lot of indenter penetration. Increasing the two penalty factors to 10 also produced very nice convergence to t=1, while reducing the penetration, though it still seemed a little too much penetration. I further increased the two penalty factors to 100, but now the analysis no longer converged beyond t=0.4. So I switched from BFGS to full Newton (<max_ups> set to zero) and convergence was now achieved very nicely to t=1 (50 time steps, 3.3 average iterations per time step), with very little penetration, 3 1/2 minutes elapsed time, see the attached input file.

I draw two conclusions from this experience: (1) Contact with friction converges much better with the full Newton scheme than BFGS (I realize this seems a sweeping generalization, but I also find it true in other contact problems). (2) I don't think there is a bug in the contact algorithm, nor is there a problem with your mesh. The only tricky issue with your model is that you have a shell surrounding hex elements and you define your contact interface on the shell; but FEBio associates the contact surface with the hex elements and therefore uses the properties of those hex elements (not the shells) to evaluate the penalty factor when auto_penalty is on. If the hex elements are significantly softer than the shell elements (as happens in your model), using a <penalty> scale factor of 1 is too small to prevent significant penetration. That's why I needed to increase the <penalty> scale factor to a much larger value to prevent penetration. For consistency, I used the same value for <fric_penalty>.

In any case, I hope that this approach also works well with your finer mesh, let me know how it comes out.

Best,

Gerard

p.s. In response to your earlier question about debugging, this is done by compiling the code using settings suitable for debugging, then running the code within the debugger to automatically trace the errors and examine the values of all the variables at the location where the error occurs. I use XCode as the compiler on my Mac and it has a very nice interface for compiling in either debug or release mode, and for displaying the content of variables in any function.

11-29-2012 06:37 AM
Excellent thanks a lot for this Gerard! I'll try to do the same for the denser mesh. I saw that this simulation slowly gets into trouble before it produces the error. Sometimes though the simulation does not even start properly and immediately produces a NaN error which I hope does not mean that this quick simulation produces the NaN for different reasons. But I am going to implement your suggestions and give it a go for my normal meshes.

Thanks again,


Kevin

07-13-2015 10:15 AM
Lead Code Developer
I'm cleaning up some old posts and I was wondering if this issue has been resolved? I'll consider it fixed unless I hear otherwise.

For the record, I also want to point that as of FEBio 2.3 (or was it 2.3.1, I forget), FEBio treats a NAN as a negative Jacobian: instead of aborting the entire run, it retries with a smaller time step size. The source of occasional NANs in contact is still unclear, but this trick has often allowed us to a point where a NAN was detected.

Cheers,

Steve

07-14-2015 06:34 AM
I haven't encountered this problem anymore. Thanks,

Kevin

+ Reply