as_vector & as_matrix notations in plasticity


178
views
0
4 months ago by
Hello all,

I am trying to expand a plasticity code and started with a Python script that has been shown in a paper.

from dolfin import *
from mshr import *
import numpy as np
L = 1; W = 0.2
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 1, 1, 1)

element = FunctionSpace(mesh,"Lagrange",2)
elementT = VectorElement("Quadrature", tetrahedron, 2, 36)
elementS = VectorElement("Quadrature", tetrahedron, 2, 6)
u, w = TrialFunction(element), TestFunction(element)
b, h = Coefficient(element), Coefficient(element)
t, s = Coefficient(elementT), Coefficient(elementS)

class DirichletBoundary(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] <= DOLFIN_EPS and on_boundary
    
def eps(u):
    return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1,2)]])
def sigma(s):
    return as_matrix([[s[0], s[3], s[4]],[s[3], s[1], s[5]],[s[4], s[5], s[2]]])
def tangent(t):
    return as_matrix([[t[i*6 + j] for j in range(6)] for i in range(6)])
# Incremental stresses
def dsigma(t, u):
    return inner(tangent(t), eps(u))
a = inner(dot(tangent(t), eps(u)), eps(w))*dx
L = inner(sigma(s), grad(w))*dx - dot(b, w)*dx - dot(h, w)*ds​

However I got an error like :

IndexError                                Traceback (most recent call last)
<ipython-input-1-12d3443a9faf> in <module>()
     26 def dsigma(t, u):
     27     return inner(tangent(t), eps(u))
---> 28 a = inner(dot(tangent(t), eps(u)), eps(w))*dx
     29 L = inner(sigma(s), grad(w))*dx - dot(b, w)*dx - dot(h, w)*ds

<ipython-input-1-12d3443a9faf> in eps(u)
     18 
     19 def eps(u):
---> 20     return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1,2)]])
     21 def sigma(s):
     22     return as_matrix([[s[0], s[3], s[4]],[s[3], s[1], s[5]],[s[4], s[5], s[2]]])

<ipython-input-1-12d3443a9faf> in <listcomp>(.0)
     18 
     19 def eps(u):
---> 20     return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1,2)]])
     21 def sigma(s):
     22     return as_matrix([[s[0], s[3], s[4]],[s[3], s[1], s[5]],[s[4], s[5], s[2]]])

/usr/local/lib/python3.5/dist-packages/ufl/exproperators.py in _getitem(self, component)
    447 
    448     # Analyse slices (:) and Ellipsis (...)
--> 449     all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
    450 
    451     # Check that we have the right number of indices for a tensor with

/usr/local/lib/python3.5/dist-packages/ufl/index_combination_utils.py in create_slice_indices(component, shape, fi)
    163             all_indices.append(ind)
    164         elif isinstance(ind, int):
--> 165             if int(ind) >= shape[len(all_indices)]:
    166                 error("Index out of bounds.")
    167             all_indices.append(FixedIndex(ind))

IndexError: tuple index out of range

I could not find any index error in tuples in the code. Can it be about the usage of as_vector and as_matrix attributes ?

Regards

Community: FEniCS Project

2 Answers


3
4 months ago by
The initial issue was that `u` was in a scalar space, so `u[i]` didn't make sense.  I got some other weird errors after that, but worked on the function space definitions a bit until the code ran:
from dolfin import *
from mshr import *
import numpy as np

# I've had strange problems in the past with quadrature elements; you may
# need these options, especially if solving problems in quadrature spaces:
#import warnings
#from ffc.quadrature.deprecation \
#    import QuadratureRepresentationDeprecationWarning
#warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
#
#parameters["form_compiler"]["representation"] = "quadrature" 

# I find this to be the most convenient way to avoid quadrature mismatches
# when using quadrature elements:
QUAD_DEG = 2
dx = dx(metadata={'quadrature_degree': QUAD_DEG})

L = 1; W = 0.2
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 1, 1, 1)

# I modified some definitions here:s
V = VectorFunctionSpace(mesh,"Lagrange",2)
elementT = VectorElement("Quadrature", mesh.ufl_cell(), degree=QUAD_DEG,
                         quad_scheme="default",dim=36)
VT = FunctionSpace(mesh,elementT)
elementS = VectorElement("Quadrature", mesh.ufl_cell(), degree=QUAD_DEG,
                         quad_scheme="default",dim=6)
VS = FunctionSpace(mesh,elementS)
u, w = TrialFunction(V), TestFunction(V)
b, h = Coefficient(V), Coefficient(V)
t, s = Coefficient(VT), Coefficient(VS)

class DirichletBoundary(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] <= DOLFIN_EPS and on_boundary


def eps(u):    
    return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i,j in [(0, 1), (0, 2), (1,2)]])
def sigma(s):
    return as_matrix([[s[0], s[3], s[4]],[s[3], s[1], s[5]],[s[4], s[5], s[2]]])
def tangent(t):
    return as_matrix([[t[i*6 + j] for j in range(6)] for i in range(6)])
# Incremental stresses
def dsigma(t, u):
    return inner(tangent(t), eps(u))
a = inner(dot(tangent(t), eps(u)), eps(w))*dx
L = inner(sigma(s), grad(w))*dx - dot(b, w)*dx - dot(h, w)*ds
​

There is also some nice support for index notation in UFL, with abstract index objects,

https://fenicsproject.org/pub/documents/ufl/ufl-user-manual/ufl-user-manual.pdf

(see Section 2.5) but, if a bunch of formulas from the paper you're reproducing are already in Voigt notation, I can see how the integer indexing approach might be more convenient.  

0
4 months ago by
Hello David,

I am really appreciate that you answered my problem very clearly. But since I am not experienced user, I need to ask you some questions related your answer.
  • You said that you faced some strange kinds of problems by using Quadrature. Were they about the compiling difficulties or that causes wrong results ? If so, is there any alternative function that can be used instead of Quadrature ?
  • You suggested to use:
#from ffc.quadrature.deprecation \
#    import QuadratureRepresentationDeprecationWarning
#warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
#
#parameters["form_compiler"]["representation"] = "quadrature" ​
However, "from ffc.quadrature.deprecation" does not work on me. May it be about the version ?
  • You used 
    dx = dx(metadata={'quadrature_degree': QUAD_DEG})​
in the code. What if you dont use that command ? Does it take the QUAD_DEG as 1 in default ?

  • What the "mesh.ufl_cell()" represent for in the VectorElement definitions ?
By the way I will not use the Voigt notation initially. Lastly, if you have some advices for me to work on plasticity, I will be glad for that.

Thanks !

Regards,


1
> You said that you faced some strange kinds of problems by using Quadrature. Were they about the compiling difficulties or that causes wrong results ?

It was mainly to do with compiling forms. 

> If so, is there any alternative function that can be used instead of Quadrature ?

In some papers on plasticity, things usually presented as being history variables at quadrature points are considered to be \(L^2\), and the "point values" (which would be undefined in that functional setting) are understood to be coefficients of constant basis functions supported on little regions around the quadrature points.  For instance, that is the interpretation advocated in Remark 2.1 of

https://doi.org/10.1016/j.cma.2009.06.021

However, as in the linked paper, you could also take some other subspace, e.g., constants over elements, which may be more robust in FEniCS. 

> However, "from ffc.quadrature.deprecation" does not work on me. May it be about the version ?

I'm using 2017.2.  The first three lines are just to suppress an annoying warning message, since the quadrature representation is being removed in the next version.  However, in my experience, solving problems over quadrature spaces is not robust in the new UFLACS representation yet.  The fourth line should work in the past few versions.  

> What if you dont use that command ? Does it take the QUAD_DEG as 1 in default ?

The default behavior is to automatically estimate an appropriate quadrature degree, based on the expression that you're trying to integrate.  See the discussion here

https://www.allanswered.com/post/reggq/how-to-find-quadrature-used/

for more information on that.  This can be a problem if the automatically-estimated degree does not line up with the degree of the quadrature function space, since the quadrature function will not be defined at the points of the automatically-generated quadrature rule. 

> What the "mesh.ufl_cell()" represent for in the VectorElement definitions ?

It just says what type of cell (e.g., tetrahedron) is associated with the mesh.
written 4 months ago by David Kamensky  
Thanks a lot David !
written 4 months ago by Christian  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »