FEniCS Project

Problem on using a=lhs(F)

Asked by Veena P on 2012-09-05

I have tries to solve weak form of elastic theories and strain gradient elastic theories using mixed formulation (displacement and strain variables)

For mixed formulation for elasticity problem, I defined

V=VectorFunctionSpace(mesh,'CG',1)
VV=TensorFunctionSpace(mesh,'CG',1,shape=(3,3))
MX=MixedFunctionSpace([V,VV])

(u,u_Strain)=TrialFunctions(MX)
(w,w_Strain)=TestFunctions(MX)

delta=Identity(V.cell().d)
i,j,k,l,p,A,B,K,L,M=indices(10)

def StrainT(u):
return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])

def StressT(u):
return as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])

def StressTS(u_Strain):
return as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])

F1=w[i].dx(j)*StressT(u)[j,i]*dx - w[i]*tr[i]*ds(1)
F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
F=F1+F2

a=lhs(F)
L=rhs(F)

it is OK..., but when solving mixed formulation for strain gradient elasticity problem

def StrainT(u):
return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])

def StressT(u):
return as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])

def StressTS(u_Strain):
return as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])

def DBStressT(u):
return as_tensor(nonlocal*nonlocal*(Lambda*StrainT(u)[p,p].dx(i)*delta[j,k]+2*mu*StrainT(u)[j,k].dx(i)),[i,j,k])

F1=w[i].dx(j)*StressT(u)[j,i]*dx + w_Strain[j,k].dx(i)*DBStressT(u)[i,j,k]*dx - w[i]*tr[i]*ds(1)
F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
F=F1+F2

a=lhs(F)
L=rhs(F)

It shows error on "a=lhs(F)". I need some help!

Question information

Language:: English Edit question

Status:: Answered

For:: FEniCS Project Edit question

Assignee:: No assignee Edit question

Last query:: 2012-09-05

Last reply:: 2018-01-25

Link existing bug

Revision history for this message

Kent-Andre Mardal (kent-and) said on 2012-09-11:

Could you provide the error message? It is a lot easier for us to
understand the
problem if we know what the error actually is :)

BTW: Using CG1 for both variables in a mixed elasticity problem will give
you an
unstable problem.

Kent

On 5 September 2012 18:06, Veena P <email address hidden>wrote:

> New question #207760 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/207760
>
> I have tries to solve weak form of elastic theories and strain gradient
> elastic theories using mixed formulation (displacement and strain variables)
>
> For mixed formulation for elasticity problem, I defined
>
> V=VectorFunctionSpace(mesh,'CG',1)
> VV=TensorFunctionSpace(mesh,'CG',1,shape=(3,3))
> MX=MixedFunctionSpace([V,VV])
>
> (u,u_Strain)=TrialFunctions(MX)
> (w,w_Strain)=TestFunctions(MX)
>
> delta=Identity(V.cell().d)
> i,j,k,l,p,A,B,K,L,M=indices(10)
>
> def StrainT(u):
> return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])
>
> def StressT(u):
> return
> as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])
>
> def StressTS(u_Strain):
> return
> as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])
>
> F1=w[i].dx(j)*StressT(u)[j,i]*dx - w[i]*tr[i]*ds(1)
> F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
> F=F1+F2
>
> a=lhs(F)
> L=rhs(F)
>
> it is OK..., but when solving mixed formulation for strain gradient
> elasticity problem
>
> def StrainT(u):
> return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])
>
> def StressT(u):
> return
> as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])
>
> def StressTS(u_Strain):
> return
> as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])
>
> def DBStressT(u):
> return
> as_tensor(nonlocal*nonlocal*(Lambda*StrainT(u)[p,p].dx(i)*delta[j,k]+2*mu*StrainT(u)[j,k].dx(i)),[i,j,k])
>
> F1=w[i].dx(j)*StressT(u)[j,i]*dx +
> w_Strain[j,k].dx(i)*DBStressT(u)[i,j,k]*dx - w[i]*tr[i]*ds(1)
> F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
> F=F1+F2
>
> a=lhs(F)
> L=rhs(F)
>
> It shows error on "a=lhs(F)". I need some help!
>
>
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

Could you provide the error message? It is a lot easier for us to
understand the
problem if we know what the error actually is :)

BTW: Using CG1 for both variables in a mixed elasticity problem will give
you an
unstable problem.

Kent

On 5 September 2012 18:06, Veena P <question207760@answers.launchpad.net>wrote:

> New question #207760 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/207760
>
> I have tries to solve weak form of elastic theories and strain gradient
> elastic theories using mixed formulation (displacement and strain variables)
>
> For mixed formulation for elasticity problem, I defined
>
> V=VectorFunctionSpace(mesh,'CG',1)
> VV=TensorFunctionSpace(mesh,'CG',1,shape=(3,3))
> MX=MixedFunctionSpace([V,VV])
>
> (u,u_Strain)=TrialFunctions(MX)
> (w,w_Strain)=TestFunctions(MX)
>
> delta=Identity(V.cell().d)
> i,j,k,l,p,A,B,K,L,M=indices(10)
>
> def StrainT(u):
>     return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])
>
> def StressT(u):
>     return
> as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])
>
> def StressTS(u_Strain):
>     return
> as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])
>
> F1=w[i].dx(j)*StressT(u)[j,i]*dx - w[i]*tr[i]*ds(1)
> F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
> F=F1+F2
>
> a=lhs(F)
> L=rhs(F)
>
> it is OK..., but when solving mixed formulation for strain gradient
> elasticity problem
>
> def StrainT(u):
>     return as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),[i,j])
>
> def StressT(u):
>     return
> as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])
>
> def StressTS(u_Strain):
>     return
> as_tensor((Lambda*u_Strain[k,k]*delta[i,j]+2*mu*u_Strain[i,j]),[i,j])
>
> def DBStressT(u):
>     return
> as_tensor(nonlocal*nonlocal*(Lambda*StrainT(u)[p,p].dx(i)*delta[j,k]+2*mu*StrainT(u)[j,k].dx(i)),[i,j,k])
>
> F1=w[i].dx(j)*StressT(u)[j,i]*dx +
> w_Strain[j,k].dx(i)*DBStressT(u)[i,j,k]*dx - w[i]*tr[i]*ds(1)
> F2=w_Strain[i,j]*(StressTS(u_Strain)[i,j]-StressT(u)[i,j])*dx
> F=F1+F2
>
> a=lhs(F)
> L=rhs(F)
>
> It shows error on "a=lhs(F)". I need some help!
>
>
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

Revision history for this message

arnaud lejeune (arnaud-lejeune) said on 2018-01-25:

UP !!!
Hello,

This is the full error message :
*************
Missing indices set([Index(9)]) in expression 0.01 * (800.0 * (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9
(dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{10}, i_{11}]) + ((I)[i_{10}, i_{11}]) * 266.666666667 * (sum_{i_{13}} (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10}
(dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{13}, i_{13}]) )).

FAILURE in reuse_if_possible:
type(o) = <class 'ufl.tensors.ComponentTensor'>
operands =

0.01 * (800.0 * (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_
9))[1]])[i_{10}])) })[i_{10}, i_{11}]) + ((I)[i_{10}, i_{11}]) * 266.666666667 * (sum_{i_{13}} (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx
_i_9))[1]])[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{13}, i_{13}]) ))

i_9, i_{10}, i_{11}

stack =
////////////////////////////////////////////////////////////////////////////////
Visit stack in Transformer:
<class 'ufl.algebra.Sum'> ; sum_{i_{10}} sum_{i_9} (([
[(v_{-2})[2], (v_{-2})[
<class 'ufl.algebra.Sum'> ; sum_{i_{11}} sum_{i_{10}} sum_{i_9} (({ A | A_{i_9,
<class 'ufl.indexsum.IndexSum'> ; sum_{i_{11}} sum_{i_{10}} sum_{i_9} (({ A | A_
<class 'ufl.indexsum.IndexSum'> ; sum_{i_{10}} sum_{i_9} (({ A | A_{i_9, i_{10},
<class 'ufl.indexsum.IndexSum'> ; sum_{i_9} (({ A | A_{i_9, i_{10}, i_{11}} = 0.
<class 'ufl.algebra.Product'> ; (({ A | A_{i_9, i_{10}, i_{11}} = 0.01 * (800.0
<class 'ufl.indexed.Indexed'> ; ({ A | A_{i_9, i_{10}, i_{11}} = 0.01 * (800.0 *
<class 'ufl.tensors.ComponentTensor'> ; { A | A_{i_9, i_{10}, i_{11}} = 0.01 * (
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

Traceback (most recent call last):
  File "expleFaq2.py", line 55, in <module>
    a=lhs(F)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\formoperators.py", line 65, in lhs
    form = expand_derivatives(form)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\ad.py", line 95, in expand_derivatives
    return transform_integrands(form, _expand_derivatives)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 837, in transform_integrands
    integrand = transform(itg.integrand())
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\ad.py", line 93, in _expand_derivatives
    return aa.visit(expression)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 156, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\ad.py", line 46, in expr
    e = Transformer.reuse_if_possible(self, e, *ops)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\algorithms\transformations.py", line 182, in reuse_if_possible
    r = o.reconstruct(*operands)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\expr.py", line 214, in reconstruct
    return self.__class__._uflclass(*operands)
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\tensors.py", line 182, in __init__
    error("Missing indices %s in expression %s." % (missingset, expression))
  File "G:\APPLICATIONS-FEMTO-ST\FEniCS\lib\site-packages\ufl\log.py", line 148, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Missing indices set([Index(9)]) in expression 0.01 * (800.0 * (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])
[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{10}, i_{11}]) + ((I)[i_{10}, i_{11}]) * 266.666666667 * (sum_{i_{13}} (({ A | A_{i_9, i_{10}} = 0
.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{13}, i_{13}]) )).

****************************

The error comes from the term
w_Strain[j,k].dx(i)*DBStressT(u)[i,j,k]*dx
in F1

I am not sure to well understand the behavior of function "i,j,k,l,p,A,B,K,L,M = indices(10)".

Second question, is there a difference between
as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)), [i,j] )
and
as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)), (i,j) )
??

Thanks for your help and/or ideas

Arnaud

UP !!!
Hello,

Sorry to reopen this question which seems not close.
 I try to solve quite the same problem and I have an error which looks like the one mentioned by "Veena P". It seems to be due to the "free indices". So to convince me I run the same example that was posted. 
 
This is the full error message  :
*************
Missing indices set([Index(9)]) in expression 0.01 * (800.0 * (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10} (dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9
(dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{10}, i_{11}]) + ((I)[i_{10}, i_{11}]) * 266.666666667 * (sum_{i_{13}} (({ A | A_{i_9, i_{10}} = 0.5 * ((([(d/dx_i_{10}
(dv_{-1}/dx_i_9))[0], (d/dx_i_{10} (dv_{-1}/dx_i_9))[1]])[i_9]) + (([(d/dx_i_9 (dv_{-1}/dx_i_9))[0], (d/dx_i_9 (dv_{-1}/dx_i_9))[1]])[i_{10}])) })[i_{13}, i_{13}]) )).

FAILURE in reuse_if_possible:
type(o) = <class 'ufl.tensors.ComponentTensor'>
operands =

i_9, i_{10}, i_{11}

stack =
////////////////////////////////////////////////////////////////////////////////
Visit stack in Transformer:
<class 'ufl.algebra.Sum'> ; sum_{i_{10}} sum_{i_9} (([
  [(v_{-2})[2], (v_{-2})[
<class 'ufl.algebra.Sum'> ; sum_{i_{11}} sum_{i_{10}} sum_{i_9} (({ A | A_{i_9,
<class 'ufl.indexsum.IndexSum'> ; sum_{i_{11}} sum_{i_{10}} sum_{i_9} (({ A | A_
<class 'ufl.indexsum.IndexSum'> ; sum_{i_{10}} sum_{i_9} (({ A | A_{i_9, i_{10},
<class 'ufl.indexsum.IndexSum'> ; sum_{i_9} (({ A | A_{i_9, i_{10}, i_{11}} = 0.
<class 'ufl.algebra.Product'> ; (({ A | A_{i_9, i_{10}, i_{11}} = 0.01 * (800.0
<class 'ufl.indexed.Indexed'> ; ({ A | A_{i_9, i_{10}, i_{11}} = 0.01 * (800.0 *
<class 'ufl.tensors.ComponentTensor'> ; { A | A_{i_9, i_{10}, i_{11}} = 0.01 * (
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

****************************

The error comes from the term
w_Strain[j,k].dx(i)*DBStressT(u)[i,j,k]*dx 
in F1

I am not sure to well understand  the behavior of function "i,j,k,l,p,A,B,K,L,M = indices(10)".

Second question, is there a difference between
 as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),     [i,j]     )
and
 as_tensor(1./2.*(u[i].dx(j)+u[j].dx(i)),     (i,j)     )
??

Thanks for your help and/or ideas

Arnaud

Can you help with this problem?

Provide an answer of your own, or ask Veena P for more information if necessary.

To post a message you must log in.

Ask a question

Edit question

FEniCS Project

Problem on using a=lhs(F)

Question information

Related bugs

Related FAQ:

Can you help with this problem?

Subscribers