Varying material properties based on data points

General
1

Dear Community,

I am using FEniCS to model the deformation of a functionally graded plate (actually a fiber reinforced cylinder but I’ve phrased my question for a plate to simplify things slightly) using displacement boundary conditions. EDIT: (I am using 2019.2.0.dev0 on WSL2).

My question is about how I can vary the material properties of the plate based on some data points derived from the displacement field solution. The question is almost identical to evdkmv’s question - the solution of which involved writing a customized Newton solver. I decided to make another post in case there have been any updates concerning this since the OP, sorry if this is not the correct place to ask.

I have been struggling to implement a Newton solver of my own, as my code (that has been adapted from the Hyperelasticity Demo) is very dissimilar to evdkmv’s.

My aim is to make K_0, currently a constant in my strain energy function, a spatially dependent scalar function of my first spatial coordinate x[0] (only). The difficulty is that I do not have an explicit expression for this function but instead generate data points for it’s value across the x[0] domain of the plate (a lookup table) based on my solution.

# Example of data points to define K_0 with
x_vals = [-5.0, -3.88888889, -2.77777778, -1.66666667, -0.55555556,  0.55555556,
          1.66666667,  2.77777778,  3.88888889,  5.0]
y_vals = [5.0071, 4.0214, 3.0651, 1.1977, 0.60081, 1.8251, 5.5443, 16.8423,
          51.1623, 155.4178]

My first attempt was to fit a cubic interpolation spline for K_0() using scipy.interpolate.interp1d() and then use K_0(x[0]) in my psi and sigma functions after assigning x = fe.SpatialCoordinate(mesh). But this results in “ValueError: object arrays are not supported” as in evdkmv’s question, RemDelaporteMathurin’s question and freya’s question.

I have tried fitting polynomials and rational functions and using these to define an fe.Expression() for K_0, but these do not produce a great fit without using an excessive degree (with rational functions often producing discontinuities).

I am considering fitting Chebyshev polynomials now. But I just wanted to check if there have been any updates since the posting of RemDelaporteMathurin’s question that might enable a more direct way to define my material properties from these data points (or if it may be simpler in my case)?

I have left a reduced version of my code below (that currently treats K_0 as constant).

import dolfin as fe
import numpy as np

# MESH PARAMS
LENGTH = 10
HEIGHT = 2
X_MIN, X_MAX, Y_MIN, Y_MAX = -LENGTH/2, LENGTH/2, -HEIGHT/2, HEIGHT/2
N_X, N_Y = 100, 20

# MATERIAL PARAMS
K_0 = fe.Constant(10)                  # Bulk
MU_0 = fe.Constant(10)                 # Shear

# DISPLACEMENT PARAMS
DISPLACEMENT = 5

# Cauchy stress
def sigma():
    cauchy_stress = (K_0*(fe.det(F) - 1)*I
    + MU_0*(fe.det(F)**(-5/3))*(F*F.T - (1/3)*(fe.tr(F.T*F))*I))
    return cauchy_stress

# Strain energy
def psi():    
    strain_energy = ((K_0/2)*((fe.det(F)) - 1)**2
    + (MU_0/2)*((((fe.det(F))**(-2/3))*fe.tr(F.T*F)) - 3))
    return strain_energy

# Optimization options for the form compiler
fe.parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create mesh and define function space
mesh = fe.RectangleMesh(fe.Point(X_MIN, Y_MIN), fe.Point(X_MAX, Y_MAX), 
                        N_X, N_Y, "right/left")
V = fe.VectorFunctionSpace(mesh, "Lagrange", 1)
x = fe.SpatialCoordinate(mesh)

# Mark boundary subdomians
left =  fe.CompiledSubDomain("near(x[0], side) && on_boundary", side = X_MIN)
right = fe.CompiledSubDomain("near(x[0], side) && on_boundary", side = X_MAX)

# Define Dirichlet boundary (x = X_MIN or x = X_MAX)
lhs = fe.Expression(("l_displacement", "0.0"), l_displacement = -DISPLACEMENT/2,
                     degree = 1)
rhs = fe.Expression(("r_displacement", "0.0"), r_displacement = DISPLACEMENT/2, 
                     degree = 1)

bcl = fe.DirichletBC(V, lhs, left)
bcr = fe.DirichletBC(V, rhs, right)
bcs = [bcl, bcr]

# Define functions
du = fe.TrialFunction(V)            # Incremental displacement
v  = fe.TestFunction(V)             # Test function
u  = fe.Function(V)                 # Displacement from previous iteration

# Kinematics
d = u.geometric_dimension()
I = fe.Identity(d)                  # Identity tensor
F = I + fe.grad(u)                  # Deformation gradient
C = F.T*F                           # Right Cauchy-Green tensor

# Total potential energy (ignoring body forces and traction)
Pi = psi()*fe.dx

# Compute first variation of Pi (directional derivative about u in the 
# direction of v)
LF = fe.derivative(Pi, u, v)

# Compute Jacobian of first variation
J = fe.derivative(LF, u, du)

# Solve variational problem
fe.solve(LF == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)

# Steps to generate new data points for next iteration ...

Any advice is of course greatly appreciated!

Thanks in advance,
Kit