dodaltuin / soft-tissue-pignn

Physics-informed graph neural network (GNN) emulation of soft-tissue mechanics
MIT License
21 stars 4 forks source link

Material Parameters Problem #2

Open A474852023 opened 2 months ago

A474852023 commented 2 months ago

Dear David, In the TwistingCube part, the E is [1,25] in your code(as well as your article) , what's the unit of E ? We are tring to perform the steel matertial with E = 210e9. The problem is the train_loss is very huge image How to solve this problem?

Another question is the loading of each part seems the same. If I want to change the loading, what can I do? Should I retrain the model, which is timeconsuming.

Thank you very much!

A474852023 commented 2 months ago

By the way, since there are only three files in the simulation folder (displacement, theta, pe-values), whether the size of the model was changed when testing the model?

dodaltuin commented 2 months ago

Hi there,

I'm not actually sure what units / scale are - I took this example straight from the FEniCS documentation here, where the units/scale aren't specified.

Regarding the loading value - the best thing to do is to give this value as an input to the emulator along with the material parameters, so that the GNN can generalise to different values without retraining. One thing to note on that - in the paper, I included material parameters after the message passing stage. You might want to try including this information (and the loading value) with the node features, which should give better accuracy.

It was however an unfortunate oversight on my behalf that I wrote the code assuming the loading value would be fixed (see for instance the fixed value external_forces.surface_force on line 166 of this file). You will need to make updates to allow this value to vary during training.

Finally, yes, I did change the size of the model, but I only put one dataset online. Simulations for different mesh sizes can be run by a slight adjustment of the FEniCS code at the same link from above. I have also made available post-processing code to get raw FEniCS simulation data into a format suitable for the GNN in this notebook here.

Thanks, David

A474852023 commented 2 months ago

Thank you very much for your reply, but I still have a couple questions that I haven't figured out yet:

  1. when I need to change the model size, do I need to retrain the model?

  2. on page 11, section 3.1 part of your paper, you set the Young's modulus to [5-10] kPa. But in the constuctive_law.py of your code (TwiceClamepedBeam), the range of E is [5e5, 10e5]. I would like to ask what should be the units here?

  3. I tried the following simulation: a square film made of a linear elastic material with four fixed surfaces, front, back, left and right, with a uniformly distributed pressure applied to the upper surface. I tried to modify the constitutive_law.py file to use the linear-elastic material, and modeled the data processing after your FenicsMeshProcessingLV.ipynbfile, which generates pressure-surface-facet.npy. But the error tends to drop from a few tens of millions as I train, the and jumps a lot throughout, and also the final prediction is completely different from the real value. Have you encountered similar problems?

image

Here's part of the code I modified in the constitutional_law.py file:


def LinearElastic(params, F, J, fibres=None):
    """Linear Elastic strain energy density function

    Inputs:
    -----------
    params: jnp.array
       Material parameter vector (E, nu)
    F     : jnp.array
       Deformation gradient
    J     : jnp.array
       Determinant of F
    fibres: None
       Not used here as we assume isotropic material

    Returns:
    ----------
    sed: float
       Strain energy density
    """

    # Extract material parameters
    E, nu = params

    # Convert to Lame parameters
    mu = E / (2 * (1 + nu))
    lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))

    # Right Cauchy-Green tensor
    C = jnp.matmul(F.T, F)

    # Green-Lagrange strain tensor
    E_tensor = 0.5 * (C - jnp.eye(3))

    # First invariant of the strain tensor
    I1 = jnp.trace(E_tensor)

    # Second invariant of the strain tensor (needed for the quadratic term)
    I2 = 0.5 * (jnp.trace(jnp.matmul(E_tensor, E_tensor)) - I1 ** 2)

    # Strain energy density function for linear elasticity
    sed = 0.5 * lambda_ * I1 ** 2 + mu * I2

    return sed

constitutive_law = LinearElastic
A474852023 commented 2 months ago

By the way, could you please give me a copy of the FEniCS code for the TwiceClampedBeam model? I would like to try to reproduce your model, but it is missing for example the density parameter.