How to interpolate or project a strain rate tensor into a TensorFunctionSpace?


80
views
0
29 days ago by
Hi,

I'm having trouble figuring out how to interpolate or project a strain rate tensor into a TensorFunctionSpace. I thought it would be very straightforward but I'm getting multiple errors with different things.

This gives me an error because as_tensor can only take 2 arguments

QTS = TensorFunctionSpace(mesh,'CG',3)
dudx, dudy, dudz = gradu[0].dx(0) , gradu[0].dx(1) , gradu[0].dx(2)
dvdx, dvdy, dvdz = gradu[1].dx(0) , gradu[1].dx(1) , gradu[1].dx(2)
dwdx, dwdy, dwdz = gradu[2].dx(0) , gradu[2].dx(1) , gradu[2].dx(2)
sij=project(as_tensor(   (dudx,0.5*(dudy+dvdx),0.5*(dudz+dwdx)) , (0.5*(dvdx+dudy),dvdy,0.5*(dvdz+dwdy)) , (0.5*(dwdx+dudz),0.5*(dwdy+dvdz),dwdz) ) ,QTS)​



And then this:

#### TENSORS #####
QTS = TensorFunctionSpace(mesh,'CG',3)
def Sij(u):
    return 0.5*(nabla_grad(u)+transpose(nabla_grad(u)))

sij = project(Sij(v0),QTS)


gives me this error:

File "HMGTENSOR.py", line 89, in Sij
return 0.5*(nabla_grad(u)+transpose(nabla_grad(u)))
File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/ufl/operators.py", line 119, in transpose
return Transposed(A)
File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/ufl/tensoralgebra.py", line 112, in __init__
error("Transposed is only defined for rank 2 tensors.")
File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/ufl/log.py", line 171, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Transposed is only defined for rank 2 tensors.
Aborted (core dumped)


Community: FEniCS Project
2
Can you include your definition of u, or rather of v0 which you supply to your function? If your argument comes from a 3D vector function space, this should work:
lagrangeElement = FiniteElement("P", tetrahedron, 1)
vectorElement = VectorElement(lagrangeElement)
tensorElement = TensorElement(lagrangeElement)

# Function spaces for u, sij
vectorSpace = FunctionSpace(mesh, vectorElement)
tensorSpace = FunctionSpace(mesh, tensorElement)

u = Function(vectorSpace)

Sij = sym(nabla_grad(u))
sij = project(Sij, tensorSpace)​
written 29 days ago by klunkean  
from dolfin import *

mesh = BoxMesh(Point(0, 0, 0),Point( 1, 1, 1), 10, 10, 10)
ff   = FacetFunction("size_t", mesh, 0) # face pointer

# DEFINE FACETS FOR INTEGRATION DIRICHLET BC
for fno in facets(mesh):
  n   = fno.normal()
  tol = 1e-3
  # right == 1 :
  if n.x() >= tol and fno.exterior():
    ff[fno] = 1
  # left == 2 :
  elif n.x() <= -tol and fno.exterior():
    ff[fno] = 2
  # back == 3 :
  elif n.y() >= tol and fno.exterior():
    ff[fno] = 3
  # front == 4 :
  elif n.y() <= -tol and fno.exterior():
    ff[fno] = 4

#### SCALARS #####
# define function space :
Q = FunctionSpace(mesh, "CG", 1)
# Any expression for now
u0 = Expression('x[2]+x[1]+x[0]', degree=0)
# Integrating in the x-direction
def integrate(u, b, d, ff, Q):
  phi    = TestFunction(Q)
  v      = TrialFunction(Q)
  bc     = DirichletBC(Q, 0, ff, b)
  a      = v.dx(d) * phi * dx
  L      = u * phi * dx
  v      = Function(Q)
  solve(a == L, v, bc)
  return v
# Extrude solution back from face to x-domain
def extrude(f, b, d, ff, Q):
  # define test and trial based on function space :
  phi = TestFunction(Q)
  v   = TrialFunction(Q)
  a  = v.dx(d) * phi * dx
  L  = DOLFIN_EPS * phi * dx  # really close to zero to fool FFC
  bc = DirichletBC(Q, f, ff, b)
  v  = Function(Q)
  solve(a == L, v, bc)
  return v
# plot solution :
v0 = integrate(u0, 2, 0, ff, Q)   # x- integral of u
w0 = extrude(v0, 1, 0, ff, Q) # x- extrusion
q0 = integrate(w0, 4, 1, ff, Q) # y- integral
p0 = extrude(q0, 3, 1, ff, Q) # y- extrusion


#### VECTORS #####
# Use the scalar field resulting from previous homogenization
# Create a gradient from it to test on vectors
QV = VectorFunctionSpace(mesh,'CG',2, dim=3)
gradu = project(nabla_grad(v0),QV)
# Integrating function in x-direction for vectors
def integrate_V(u, b, d, ff, Q):
  phi    = TestFunction(Q)
  v      = TrialFunction(Q)
  bcu     = DirichletBC(Q, (0,0,0), ff, b)
  a      = dot(v.dx(d), phi) * dx
  L      = dot(u, phi) * dx
  v      = Function(Q)
  solve(a == L, v, bcu)
  return v
def extrude_V(f, b, d, ff, Q):
  # define test and trial based on function space :
  pseudozero = as_vector( [DOLFIN_EPS]+[DOLFIN_EPS]+[DOLFIN_EPS])
  phi = TestFunction(Q)
  v   = TrialFunction(Q)
  bcu    = DirichletBC(Q, f, ff, b)
  a  = dot(v.dx(d), phi) * dx
  L  = dot(pseudozero, phi) * dx  # really close to zero to fool FFC
  v  = Function(Q)
  solve(a == L, v, bcu)
  return v

gradI = integrate_V(gradu,2,0,ff,QV)
gradE = extrude_V(gradI,1,0,ff,QV)

#### TENSORS #####
QTS = TensorFunctionSpace(mesh,'CG',3)
def Sij(u):
    return 0.5*(nabla_grad(u)+transpose(nabla_grad(u)))

sij = project(Sij(gradI),QTS)

# Integrating function in x-direction for tensors
def integrate_TS(u, b, d, ff, Q):
  phi    = TestFunction(Q)
  v      = TrialFunction(Q)
  bcu     = DirichletBC(Q, (0,0,0,0,0,0,0,0,0), ff, b)
  a      = inner(v.dx(d), phi) * dx
  L      = inner(u, phi) * dx
  v      = Function(Q)
  solve(a == L, v, bcu)
  return v

sijI = integrate_TS(sij,2,0,ff,QTS)
file1 = File('SIJ.pvd')
file1 << sijI
###########


This is what I was trying before but I get an error when I try to project Sij to the TensorFunctionSpace, is there any way around this?

UMFPACK V5.7.6 (May 4, 2016): ERROR: out of memory

Traceback (most recent call last):
  File "HMGTENSOR.py", line 91, in <module>
    sij = project(Sij(gradI),QTS)
  File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/fem/projection.py", line 147, in project
    cpp.la_solve(A, function.vector(), b, solver_type, preconditioner_type)
  File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/cpp/la.py", line 5274, in la_solve
    return _la.la_solve(*args)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /feedstock_root/build_artefacts/fenics_1489214170396/work/dolfin-2016.2.0/dolfin/la/PETScLUSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2016.2.0
*** Git changeset:  3a4c1d9243e205211b5a19200e755b30bd3c3b6c
*** -------------------------------------------------------------------------




written 29 days ago by Luz Imelda Pacheco  
1
Well maybe you really just are out of memory. If you put
set_log_level(DEBUG)​
somewhere, you can see the size of your linear systems.

I reduced the orders of all your function spaces to 1, i.e.

QV = VectorFunctionSpace(mesh,'CG',1, dim=3)
QTS = TensorFunctionSpace(mesh,'CG',1)​

and the code runs without errors.

Just for a quick check:
print QTS.dim()​
returns 268119.
That's a lot, I guess, since the corresponding mass matrix for projection is of size 2681192

If the results are acceptable with lower order function spaces is up to your judgment, I did not check the results.
written 29 days ago by klunkean  
Thank you very much!
written 29 days ago by Luz Imelda Pacheco  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »