### 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.

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
#
# 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)
# The linear systeme is solved
solve (a == L , u_nouveau, bc)
# Define A as an array
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: