Copy of ufl.algebra.product

1

Hi,

I am currently trying to establish a spectral decomposition of 2D and 3D tensors in fenics for an Ogden-type material model. The code is based on analytic evaluations of the characteristic polynomial of the tensor and has proven to work well. Since I am assuming that the eigenvalues are of algebraic multplicity 1, I have to take care of this fact after evaluating them: I do this by adding some numerical pertubations to the individual eigenvalues, if they are equal within a certain tolerance.

My problem is that, to check if the eigenvalues have changed after a conditional expression, I need a copy of the eigenvalue before it might have been changed. If I simply make a copy of this object, it seems to be only a shallow copy. Does anybody have an idea on how to make a deep copy, i.e. a new object that simpy stores all properties of my eigenvalue before I change it using a conditional expression?

Here is a snippet of the code, I am working with:

###########################
# load required libraries #
###########################
import dolfin as df                                                             # dolfin module
import ufl                                                                      # ufl module

####################################################
# Method for analytical calculation of eigenvalues #
####################################################
def eigenValues(T, tol=df.DOLFIN_EPS):
    '''
    Method to calculate the (real) eigenvalues of a second order tensor T using
    its principal invariants.
    '''
    dim=T.geometric_dimension()                                                 # get geometric dimension of tensor T to decide which subroutine to use for eigenvalue computations

    # choose projection tensor routine based on dimensionality of problem
    if dim == 2:                                                                # 2D-problem
        return _eigVal2D(T,tol)
    else:                                                                       # 3D-problem
        return _eigVal3D(T,tol)

########################################################################
# Method for the analytical calculation of eigenvalues for 2D-Problems #
########################################################################
def _eigVal2D(T, tol):
    '''
    Analytically calculate eigenvalues for a two-dimensional tensor with a
    characteristic polynomial equation of the form

                    lambda^2 - I1*lambda + I3 = 0   .

    Since the characteristic polynomial is a reduced quadratic equation, the
    eigenvalues can be determined using the p-q-formula.

    NOTE:
    The method implemented here, implicitly assumes that the polynomial has
    only real roots, since imaginary ones should not occur in this use case.
    '''

    # determine perturbation from tolerance
    pert = 2*tol

    # get required invariants
    I1 = df.tr(T)                                                               # trace of tensor
    I3 = df.det(T)                                                              # determinant of tensor

    # determine discriminant
    D = I1-4*I3                                                                 # preliminary value for discriminant
    D=ufl.conditional(ufl.lt(D,tol),D+pert,D)                                   # add numerical perturbation to discriminant, if close to zero

    # calculate polynomial roots
    # -> Here, the polynomial solution formula is only used to calculate the
    # -> "bigger" root, since - if I1 and sqrt(D) are of similar magnitude - a
    # -> loss of significance could lead to erroneous results. For the second
    # -> root, Vietas formula is used.
    lam1 = 0.5*(I1+ufl.sqrt(D))
    lam2 = I3/lambda1

    # use Miehe's method to ensure algebraic multiplicity of 1
    # ("Computation of isotropic tensor functions" by C. Miehe (1993))
    maxLambda = ufl.max_value(lam1, lam2)
    lambda1=ufl.conditional(ufl.lt(abs(lam2-lam1)/maxLambda,tol),lam1*(1+pert),lam1)
    lambda2=ufl.conditional(ufl.ne(lam1,lambda1),lam2*(1-pert),lam2)

    # return tuple of roots (eigenvalues)
    return lambda1, lambda2

Especially, in the 3D case (not shown here, since the 2D case already shows the problem), the comparison of the old and new values lam1 and lambda1 is where I would need the copy of the ufl object, because multiple changes could be necessary if pairs of eigenvalues are nearly equal.

The object type of the individual eigenvalues is ufl.algebra.product and I have absolutely no idea, how to make a copy of them…

1 Like