Option for parameters as Quadrature values

[srosenbu]

We currently have the option for different material laws that are defined on submeshes. This allows for using material models with differently shaped history together. However, for cases where the same model but with different parameters is used, it might make more sense to have the option of defining the parameters as values on the quadrature points instead of globally for the whole (sub) mesh.

This can be useful for:

1. Layers of different age in additiv emanufacturing
2. Statistically perturbated material properties in a domain

In order for this to work, we would need:

• A unified way to define parameters for the constitutive models
• In the IncrSmallstrainModel in the main loop, one would need to make the distinction between looping through an array of parameters or using the same parameters for each quadrature point.

[joergfunger]

If this is really required, couldn't you just define another scaling factor on the quadrature space for each parameter and then set these values (either a random field, or from the history when it was printed, ..)?

[srosenbu]

@joergfunger , I am not sure what you mean. Maybe I can clarify how I would want this to work. When I write a constitutive model, I usually have a loop that calls an evaluate_ip for each quadrature point. It might look similar to this

def evaluate(
        self,
        del_t: float,
        grad_del_u: np.ndarray,
        mandel_stress: np.ndarray,
        tangent: np.ndarray,
        history: np.ndarray | dict[str, np.ndarray] | None,
) -> None:
    n_gauss = mandel_stress.size//6
    for ip in range(n_gauss):
        evaluate_ip(ip, del_t, grad_del_u, ..., self.parameters)

Here self.parameters could be an array containing the parameters that are the same for the whole domain.

We could now add an optional parameter quadrature_parameters

def evaluate(
        self,
        del_t: float,
        grad_del_u: np.ndarray,
        mandel_stress: np.ndarray,
        tangent: np.ndarray,
        history: np.ndarray | dict[str, np.ndarray] | None,
        quadrature_parameters: np.ndarray | None,
) -> None:
    n_gauss = mandel_stress.size//6
    if quadrature_parameters is None:
        for ip in range(n_gauss):
            evaluate_ip(ip, del_t, grad_del_u, ..., self.parameters)
    else:
        for ip in range(n_gauss):
            parameters = quadrature_parameters[ip]
            evaluate_ip(ip, del_t, grad_del_u, ..., parameters)

The changes in the code would be minimal for each constitutive model. I am however not sure, if this option should be enforced in the interface. I could also imagine that we write an example model for this and if someone wanted to use this option they would need to write their own solver.

[joergfunger]

What I mean is that you could just use another entry in the history variables dict that stores these parameters.
So when you call evaluate_ip(ip, del_t, grad_del_u, ..., parameters), you pass these parameters automatically with the history variables (and your material model has to know that there is another value in the history variables to be used e.g. for scaling the parameters or directly setting them). So in the evaluate_ip method you could say

if history['E'] is None:
        E = parameters['E']
else:
        E = history['E']

Alternatively, you could store a general array in the history variables, but I guess this is very similar to your approach (just avoiding another array, and "hide" it in the history dict.

def evaluate(
        self,
        del_t: float,
        grad_del_u: np.ndarray,
        mandel_stress: np.ndarray,
        tangent: np.ndarray,
        history: np.ndarray | dict[str, np.ndarray] | None,
        quadrature_parameters: np.ndarray | None,
) -> None:
    n_gauss = mandel_stress.size//6
    if history['quadrature_parameters'] is None:
        for ip in range(n_gauss):
            evaluate_ip(ip, del_t, grad_del_u, ..., self.parameters)
    else:
        for ip in range(n_gauss):
            parameters = history['quadrature_parameters'][ip]
            evaluate_ip(ip, del_t, grad_del_u, ..., parameters)

[srosenbu]

Ok, that's possible as well. I did a similar approach with optional output/history variables in my implementation of the JH2 model: optional_history.

From an implementation point of view, I would prefer to have the if-condition outside of the evaluate_ip, because I usually implement the loop for all models inside a package and only rewrite the evaluate_ip. Then I would not need to think about this extra functionality when I add a new model.

For me it comes down to how flexible we want this approach to be: Do we want this to work with almost every model that we implement without too much additional code, then I would vote for my proposal. This approach would have the disadvantage that we would need to save all parameters on each quadrature point even if we only want one of them to be different.

Your suggestion would allow for fine-grained control over which parameters are varied throughout the domain, however it would however require more additional code for each implemented material model.

[joergfunger]

What would be the examples that we currently envision where this would be needed?

[srosenbu]

@aradermacher wants to use constitutive models where the youngs modulus is different for the layers in an additive manufacturing simulation.
I also remember a conversation that I once had with @danielandresarcones where he mentioned that he may be interested to use constitutive models with varying parameters in his work. But I don't remember his exact application.

But maybe they can comment on that themselves

[aradermacher]

To comment ;): the question is more, if we would like to force everyone to write the material laws in a way that field parameters are possible or not. That we can adapt the models for that need is clear.

[danielandresarcones]

I do use sometimes stochastic material properties such as Young's modulus as Gaussian Random Fields, but I directly take an adaptation of the implementation at FenicsXConcrete. There, the material parameters are defined as a Function over the subspace, then the vector values are calculated independently and assigned to the dolfinx Function object. I have not found any limitations up to now with this procedure, but my current cases are quite simple material-wise (static linear elasticity with only one subdomain). Having control over the quadrature points will probably be advantageous at some point when I try to approach the bias in the constitutive law, but I have not really looked into it and don't know if I can directly solve it with FenicsX (or with the new 8.0 version).