negative jacobians in DEMO_febio_0044_mammography_01

[gibboncode]

negative jacobians in DEMO_febio_0044_mammography_01

[Kevin-Mattheus-Moerman]

ANDRE! Nice to "see" you here.

Did you alter the demo in any way or is this the demo "out of the box"? If you made changes, what did you change? Where does the simulation stop, i.e. at what point?

The "standard" demo runs fine on my machine. I think this is not a MacOS/Ubuntu or MATLAB 2017b/MATLAB 2019b issue, I think this is likely an "issue" with the differences in your hardware and FEBio. I think perhaps my demo is a bit sensitive to some numerical differences in the way the FEBio model runs on your hardware and on mine. Here are some hopefully useful comments. Let me know if you learn anything based on them. If you are still stuck we can have a Jitsi/TeamViewer/Skype/.. meeting some time to trouble shoot this issue on your machine.

Perhaps the initial mesh is bad

Negative Jacobian warning/errors imply that FEBio cannot find a solution to the next part of the simulation without distorting an element so much that it creates a negative volume. If the element is initially of a bad quality (its so skewed its easy to "bump it" into a negative volume) we need to look at improving the mesh quality, more elements could do the trick. So in this case we can change the point spacing or the volume factor for instance. FEBio might list the violating elements too so you can plot/check out where they are and what they look like initially. I think it might actually be a contact or time stepping issue though (see next points).

Take it slow, i.e. smaller time steps may be required

Negative Jacobians may also occur if we are forcing FEBio to take a big leap to the next time increment. For instance if the whole motion is done at once it is such a huge step it may only come up with "solutions" where elements are inverted or distorted excessively, leading to negative Jacobians. Conversely we could ask it to take a nanosecond step (e.g. if the full simulation is 1s) and it would likely instantly solve as FEBio is like "oh thats easy the solution is nearly the same as where we are now". You have to find a happy medium between jumping to the final solution (e.g. at 1s) and taking tiny steps. I usually use 10 steps initially (resulting in a 1/10 time step size), I then try 20, 50 etc., so you play with increasingly smaller time steps to see if it does solve.

numTimeSteps=10; %Number of time steps desired

I would run FEBio in internal mode so you can see the text of the log file in the command prompt. Besides the time step size we also want the "optimal number of iterations" to make sense. After choosing time stepping settings you monitor the log file and see how many steps it takes to solve one time increment. If you have the following:

opt_iter=10; %Optimum number of iterations

and for instance 1/10 as a time step size then after 10 iterations we should solve for t=0.1. If it goes into 12, 15 or 20 iterations then the problem appears more difficult than we think (as we are telling FEBio that 10 iterations sounds like the ideal case). BTW if more are used FEBio will reduce the time step size automatically (if you let it), and if less are used FEBio will increase the time step size automatically (if you let it). So getting this number right helps. You can get it right by looking at the look file and seeing how many it is using to solve the steps you set. Play around with this number so that it refines the time steps when needed and restores to being fast when needed too. It is a dark art to get these right 😄.

Contact is difficult

Contact problems are notoriously difficult to solve. My guess is the most sensitive parameter is:

contactPenalty=50; %Contact penalty

It determines how much force to used to push nodes away when penetrations occur during contact. I like to image this in the following way. If the penalty is too low the forces are too low and excessive penetration will occur. A bit of penetration is okay but at some point I feel the situation becomes unstable as the geometry starts to make less and less sense. The model might then fail because of this. Conversely if the penalty parameter is too high you can think of the excessive force being applied which is kinda "violent" and might distort the elements when contact occurs. Like nodes being hit by a baseball bat (rather than being pushed) at the interface. This might cause the nodes to be pulled away so far that they end up inverting an element (negative Jacobian).
I tend to start with a contact penalty parameter that is too low. It will likely solve quite far and you can see the amount of contact penetration (one shape will "spill/bleed through" the other). The magnitude of this number depends on the stiffness and node spacing of the surfaces involved so it can be tricky to get right. I think for most of my problems (e.g. units in mm, stiffness in kPa) starting at 0.1, then 1, then 10 etc sounds like a good way to start off. For this model I ended up with 50 which is perhaps too high. It might be that there is little difference between 10, 20, 30 (e.g. all of these enforce no or little penetration) however if you keep increasing it then eventually you start getting the unstable contact stuff I describe above. If I change to 20 for this model it still runs fine but I have a bit of contact penetration near the chest. (BTW some times it may look like a shape is penetrating but the problem is more geometric, e.g. if all the slave nodes are outside/on the master surface but an edge might geometrically have to pass through the master surface).

I do think this model has some ugly contact associated with it. I wonder if perhaps I am loading too close to the chest as well. If we were to position the compressing plates further away there would be less of the complex tissue "pulling" close to the chest. Perhaps the compression magnitude is also somewhat excessive.

Let me know if this helps.

[Kevin-Mattheus-Moerman]

I'll close this for now but do let me know if you have any remaining questions on this.

[sergioblancoib]

Dear Kevin,
I am writting you because I have the same issue with the mammography problem that had André. If I run the code as it is in gibboncode, I am able only to run four steps before it crashes (around time 0.4). If I see the postprocess (i am sending the picture) I see a total penetration inside the breast. I've tried with previous versions of febio but I got the same results. I've tried also to change the contact penalty with the similar results.

Could you please help me? Is there samething that I am missing?

Thank you very much,

Sergio