Using Numpy and Scipy into loops in an iterative algorithm


249
views
1
6 months ago by

Good evening.

I am trying to implement the Glowinsky algorithm to solve problem for an incompressible finite elastic material.
To do so I have to compute several mathematical operations for a discontinuous field (one problem per cell).

To summarize, I have to decompose in singular values a 3x3 tensor, then to compute the minimum of a function of three variables depending on the singular values of the tensor.

To do so, I use classical python functions from Numpy and scipy.optimize
I suspect that this slows down strongly the code execution.

Is it possible to compute the same mathematical operations using Fenics built-in operators.

The script I have coded can be found below.

Thank you in advance for your kind assistance.

Best regards

Xavier  


# -*- coding: utf-8 -*-
from fenics import *
from mshr import *
import numpy as np
from scipy.optimize import fsolve

# Material coefficient (Neo Hookean and Money Rivlin)
C = Constant(1.)
C2 = Constant(1.)

# Algorithm parameters
r = Constant(100.)   # paramètre r
eta = 2*C            # fonction eta
kmax = 3             # itérations sur k
niter = 5            # itérations sur n
tol = 1e-3           # tolérance pour la convergence

# Mesh
mesh = Cylinder( Point(1,0), Point(0,0), 0.5, 20)

# Define geometry
cylinder = Cylinder(Point(0, 0, -1), Point(0, 0, 1), 1.0, 0.5)
# Generate mesh
mesh = generate_mesh(cylinder , 20)
d = 3
# Initial values for F and Lambda
F = Identity(d)
Lambda = Expression( (('1','2','3'),('4','5','6'),('7','8','9')) ,degree =1)
# Functional spaces
V = VectorFunctionSpace(mesh,"CG",1)
V0 = TensorFunctionSpace(mesh,"DG",0)
# Functions
v = TestFunction(V)
u = TrialFunction(V)
u_nouveau = Function(V)
u_ancien = Function(V)
F0 = Function(V0)
# bilinear form
a = (2*C+eta*r)*inner(grad(u),grad(v))*dx
#
# Boundary conditions
def top(x, on_boundary):
    return near(x[2],-1) and on_boundary
def bottom(x, on_boundary):
    return near(x[2],1) and on_boundary
delta = 0.3
bc0 =  DirichletBC(V, Constant((0,0,-delta )), top)
bc1 =  DirichletBC(V, Constant((0,0, delta)), bottom)
bc = [bc0, bc1]

# main loop
compteur_n=0
erreur = tol + 1
for n in range (niter):
    compteur_n += 1
    print "n =", compteur_n
    compteur_k=0
    # secondary loop
    for k in range( kmax ):
            compteur_k += 1
            print "k =", compteur_k
            # new linear form
            XC = (-2*C-eta*r)*Identity(d)+eta*r*F-eta*Lambda
            XC = project(XC,V0)
            L = inner(XC,grad(v))*dx
            # The linear systeme is solved
            solve (a == L , u_nouveau, bc)
            # Define A as an array
            A = eta*(r*(Identity(d)+grad(u_nouveau))+Lambda)
            A = project (A,V0)
            A_array = A.vector().array()
            F_Array = []
            tet_inds = [cell.index() for cell in cells(mesh)]
            # loop on elements
            for l in tet_inds:
                  # Singulars values decomposition on each element
                  Q,D,R = np.linalg.svd(np.reshape(A_array[(d**2)*l:(d**2)*(l+1)], (3,3)))
                  # non linear equation on each element
                  def equations(t):
                      Jc_1 = r * eta * (2 * t[0] - 2 / t[0] ** 3 / t[1] ** 2) / 2 - D[0] + D[2] / t[0] ** 2 / t[1]
                      Jc_2 = r * eta * (2 * t[1] - 2 / t[0] ** 2 / t[1] ** 3) / 2 - D[1] + D[2] / t[0] / t[1] ** 2
                      # The values of W2_1 and W2_2 depend on the particular behavior (o for neo Hookean material)
                      W2_1 = 0
                      W2_2 = 0
                      return [Jc_1 + W2_1 , Jc_2 + W2_2]
                  t1, t2 =  fsolve(equations, [1, 1])
                  # Compute T for each element
                  T = np.diag([t1,t2,1/(t1*t2)])
                  F_l = np.dot(np.dot(Q,T), R)
                  # store deformation gradient in an array
                  F_Array = np.append(F_Array, F_l )
            # project deformation gradient on V0
            F = Function (V0)
            F.vector()[:] = F_Array

    # new Lambda
    Lambda = project(Lambda + r * (Identity (d) + grad(u_nouveau) - F),V0)
    # Error for convergence
    num1 = norm(project(F - Identity(d) - grad(u_nouveau) , V0), 'L2')
    denum1 = norm(F,'L2')
    err1 = num1/denum1
    print "erreur1 = " , err1
    num2 = norm (project(u_nouveau-u_ancien, V), 'L2')
    denum2 = norm (u_nouveau, 'L2')
    err2 = num2/denum2
    print "erreur2 = " , err2
    erreur = err1 + err2
    u_ancien.assign(u_nouveau)
print "L'algo a converge en ", compteur_n, " iterations"
           



Community: FEniCS Project
1
I haven't quite parsed your algorithm from the code, but I once implemented the algorithm from this Wikipedia page

https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices

to get the eigenvalues of a 3x3 matrix.  It is non-iterative and can be implemented using the conditional() structure in UFL:

http://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#conditional-operators

However, my implementation always seemed a bit shaky, and needed some tuning with epsilons.  I even had to try a few different mathematically-equivalent expressions for the form to even compile (although, to be fair, I was pushing a function of the eigenvalues through derivative() twice, which was probably asking a bit much).  There is also a non-iterative algorithm for polar decomposition of a deformation gradient in Box 7.1 of the Simo and Hughes "Computational Inelasticity" book.  It looks like the authors have put some thought into numerical stability issues, although I haven't tried implementing that one myself in UFL.
written 6 months ago by David Kamensky  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »