BAMresearch / fenics-constitutive

Complex constitutive models beyond the FEniCS UFL.
https://bamresearch.github.io/fenics-constitutive
MIT License
13 stars 2 forks source link

WIP Gradient damage #6

Closed joergfunger closed 5 months ago

joergfunger commented 3 years ago

I know that this is work in progress, but I had a look at it already and think that it seems already a very good approach (couldn't test all of it because the mesh file for the test test_gradient_damage was not uploaded yet and the bending example had some errors - WIP). Some remarks

TTitscher commented 3 years ago

Should we foresee to allow the users to decide whether they want the hessian as well - for the explicit simulations it seems overkill to always compute the hessian - and similar - when you perform a line search in a Newton type of procedure, you might not necessarily have to recompute the complete hessian. (and not even allocate the memory)

KRATOS, for example, passes flags like "COMPUTE_STRESS" or "COMPUTE_CONSTITUTIVE_TENSOR" to the 'evaluate' method. Currently, this is mainly relevant for you, @srosenbu, so maybe you could think about it. In any ways, this would be trivially possible in the inferface

void Evaluate(InputList, OutputList)

but does not fit the

std::pair<OutputStress, OutputDStress> Evaluate(Input) 

How is this approach working with MPI-parallelization - I guess that FEniCS is providing only an ip-vector for each MPI-domain, thus we don't have to worry about this, right?

In the example, there was a tiny bug related to dolfin log levels. I fixed that and it runs smoothly with mpi.

How do we handle non-convergence in a constitutive formulation? I think it would be nice to distinguish between an error that is really fatal and non-convergence. (This might happen e.g. in a multisurface-plasticity model). I guess so far this is not an issue and also not necessary to implement, just wanted to mention it

By default, every c++ exception is translated to the python RuntimeError by pybind, but you could also bind your own ReturnMappingError. I am, however, not sure how the dolfin.NewtonSolver reacts to exceptions during evaluating the residuals/hessians. If all exceptions (from the solver and the ip loop) cause a reduction in time step, we can treat them equally:

try:
   converged, number_of_iterations = newton_solver.solve(...)
except ReturnMappingError:
   return False, 42 # "false" will trigger a restore of the previous DOF values as x0 for the next NR step with smaller time step.

In this FixIp why could it be that there is a material allocated and the ip-vector is empty?

If there is exactly one law, the library can trivially figure out the IPs. Just convenience for the user.

I think there is a copy paste error here. The degree of the interpolation for the nonlocal strains is set with the variable related to the displacements

Thanks!

Maybe constraints is not a self explaining term to distinguish the dimension (and plane strain vs plane stress), what about dimension

I think that "dimension" rather unambiguously refers to integers 1,2,3 and is also not optimal. I checked wikipedia and other FEM libs and they use terms like "plane strain case" or "plane stress assumption". To make it self explaining, we probably have to make it longer like "off_dimension_assumption". Could just defining/documenting the meaning of "constraint" (or any other rather short word) be sufficient?