defining DirichletBC of a Tensor function space

Hello Veena,

does something like this work?

zero = Constant(0)

bc_2 = DirichletBC(VV.sub(2), zero, BOUNDARY)
bc_5 = DirichletBC(VV.sub(5), zero, BOUNDARY)
bc_6 = DirichletBC(VV.sub(6), zero, BOUNDARY)
 ... and so on

Cheers,

Jan

Thank you Jan

I have tried to used that but it seems that I need the bc_2, bc_5, bc_6, bc_7 and bc_8 written in a tensor form - that one look like a vector form.

Actually, my problem is a mixed formulation. The vector function space represents displacements and the tensor function space is a strain tensor

   V = VectorFunctionSpace(mesh,'CG',1) # displacement vector
   V_x, V_y, V_z= V.split()
   VV = TensorFunctionSpace(mesh,'CG',1,shape=(3,3)) # strain tensor
   VV_0, VV_1, VV_2, VV_3, VV_4, VV_5, VV_6, VV_7, VV_8= VV.split()
   MX = MixedFunctionSpace([V,VV])

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

So, when I define the dirichlet boundary condition of displacements, V_z = 0.0 - it is OK

   bc1=[DirichletBC(V.sub(2),(0.0), side1)]

But, when I do the same approach on the strains, VV_2, VV_5. VV_6, VV_7,VV_8

   bc2=[DirichletBC(VV.sub(2),(0.0), side2)]
   bc3=[DirichletBC(VV.sub(5),(0.0), side2)]
   .... and so on

and the variational form I have, the VV variables are in a tensor form u_Strain [i,j]

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

When I run the code, It shows error which says;

ERROR: Unable to solve variational problem
REASON: Unable to extract boundary condition arguments

Hi,

> I have tried to used that but it seems that I need the bc_2, bc_5, bc_6, bc_7
> and bc_8 written in a tensor form - that one look like a vector form.

I'm not sure what you mean here, but AFAIK the tensor is somewhat linearized in
a row-major style. So the tensor T

    / T_11 T_12 T_13 \
T = | T_21 T_22 T_23 |
    \ T_31 T_32 T_33 /

is indexed as

T_linearized = [T_11 T_12 T_13 T_21 T_22 T_23 T_31 T_32 T_33]

this way you only need one index to access the elements of the tensor.

Maybe someone how knows can answer this better.

>
> Actually, my problem is a mixed formulation. The vector function space
> represents displacements and the tensor function space is a strain tensor
>
> ...
>
> MX = MixedFunctionSpace([V,VV])
>

Maybe this is the problem(?)
In a mixed formulation you're supposed to apply the boundary conditions on the
mixed function space (see for example:
http://fenicsproject.org/documentation/dolfin/1.1.0/python/demo/pde/mixedpoisson/python/documentation.html#implementation
somewhere in the middle of the section "5.2. Implementation").

So instead of

> bc1=[DirichletBC(V.sub(2),(0.0), side1)]

you should write

   bc1=[DirichletBC(MX.sub(0).sub(2),(0.0), side1)]

and for the tensor

   bc2=[DirichletBC(MX.sub(1).sub(2),(0.0), side2)]
   bc3=[DirichletBC(MX.sub(1).sub(5),(0.0), side2)]
   .... and so on

Also, maybe (I'm guessing here), if you split your function spaces in a mixed
form you might need to write something like this(??):

V_x, V_y, V_z= MX.sub(0).split()
VV_0, VV_1, VV_2, VV_3, VV_4, VV_5, VV_6, VV_7, VV_8= MX.sub(1).split()

Anyhow, I don't really understand why you do this? Sorry!

Cheers,

Jan

I am trying to do this

from dolfin import *
set_log_level(DEBUG)

#geometry is a 3D beam
xlength=100.0#[m]
ylength=1.0#[m]
zlength=1.0#[m]

mesh=Box(0,0,0,xlength,ylength,zlength,100,1,1)
coor = mesh.coordinates()

V = VectorFunctionSpace(mesh,'CG',1) #displacement
VV = TensorFunctionSpace(mesh,'CG',1,shape=(3,3)) #strain
MX = MixedFunctionSpace([V,VV])
V_x, V_y, V_z= MX.sub(0).split() #u_x, u_y, u_z
VV_0, VV_1, VV_2, VV_3, VV_4, VV_5, VV_6, VV_7, VV_8= MX.sub(1).split() #e11 e12 e13, e21 e22 e23, e31 e32 e33
(u, u_Strain) = TrialFunctions(MX)
(w, w_Strain) = TestFunctions(MX)

#for integer marking of the surface(surface is one dimensionless than volume)
#from the domain V find the parts of the surface to clamp and tract
D = mesh.topology().dim()
left = compile_subdomains('x[0]==0.0')
right = compile_subdomains('x[0]==length')
right.length = xlength
bottom = compile_subdomains('x[1]==0.0')
top = compile_subdomains('x[1]==length')
top.length = ylength
back =compile_subdomains('x[2]==0.0')
front = compile_subdomains('x[2]==length')
front.length = zlength
neumann_domains = MeshFunction("uint",mesh,D-1)

#marking the parts of the surface
neumann_domains.set_all(0)
right.mark(neumann_domains,1) #for traction

# For the ds(1) symbols to work we must obviously connect them (or a and L)
# to the mesh function marking parts of the boundary. This is done by a
# certain keyword argument to the assemble function:
ds = Measure("ds")[neumann_domains]

#loading in the surface normal direction, the traction vector in MPa
tr=Expression(('0.0','0.0','0.0'))
#Body Force, the gravity accerlation in m/s^2
bf=Expression(('10.0','0.0','0.0'))

#a zero for the left clamped in the wall
null=Constant((0.0,0.0,0.0))

#stating the boundary values,which will be set after the assembly
#bc=[DirichletBC(V,null,left)] # This is working
bc1=[DirichletBC(MX.sub(0).sub(0),(0.0),left)]
bc2=[DirichletBC(MX.sub(0).sub(1),(0.0),left)]
bc3=[DirichletBC(MX.sub(0).sub(2),(0.0),left)]
bc=[bc1, bc2, bc3]

#material coeff.s of a stainlesss tell
nu=0.3
Rho=10 #kg/m^3
E=210e6 #MPa
G=E/(2.0*(1.0+nu))
Lambda=2.0*G*nu/(1.0-2.0*nu)
mu=G

#metric
delta=Identity(V.cell().d)#Identity(3)

#index notation
i,j,k,l,A,B,K,L,M=indices(9)

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

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

#Stress tensor based on Strain (from extra field, strain)
def StressTS(u_Strain):
 return as_tensor((Lambda*u_Strain[k,k]*delta[i,j] + 2*mu*u_Strain[i,j]), [i,j])

#variational form
F1 = StressT(u)[i,j]*w[i].dx(j)*dx - w[i]*tr[i]*ds(1) - w[i]*Rho*bf[i]*dx
F2 = (StressTS(u_Strain)[i,j] - StressT(u)[i,j])*w_Strain[i,j]*dx

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

mx = Function(MX)

# Solving linear variational problem
solve(a == L, mx, bc)
u, u_Strain = mx.split()

I have used "bc1=[DirichletBC(MX.sub(0).sub(0),(0.0),left)]" to define u_x is equal to 0.0 on the left surface.
and used "bc2=[DirichletBC(MX.sub(0).sub(1),(0.0),left)]" and "bc3=[DirichletBC(MX.sub(0).sub(2),(0.0),left)]"
for u_y = 0 and u_z =0 but it is not working.

Hi,

I think I spotted the problem.

> bc1=[DirichletBC(MX.sub(0).sub(0),(0.0),left)]
> bc2=[DirichletBC(MX.sub(0).sub(1),(0.0),left)]
> bc3=[DirichletBC(MX.sub(0).sub(2),(0.0),left)]
> bc=[bc1, bc2, bc3]

with this code you're defining a list of lists (of DirichletBC's), meaning bc
looks like

bc = [[DirichletBC(MX.sub(0).sub(0),(0.0),left)],
      [DirichletBC(MX.sub(0).sub(1),(0.0),left)],
      [DirichletBC(MX.sub(0).sub(2),(0.0),left)]]

but solve expects a flat list of DirichletBC's.

Just remove the square brackets around the DirichletBC(MX...

like:

bc1=DirichletBC(MX.sub(0).sub(0),(0.0),left)
bc2=DirichletBC(MX.sub(0).sub(1),(0.0),left)
bc3=DirichletBC(MX.sub(0).sub(2),(0.0),left)
bc=[bc1, bc2, bc3]

this way your code works for me. (I assumed F = F1 + F2)

I also tried a DirichletBC for the strain

bc4=DirichletBC(MX.sub(1).sub(0),(0.0),top)
bc.append(bc4)

and it works too.

Cheers,

Jan

the link in answer #3 is wrong

http://fenicsproject.org/documentation/dolfin/1.1.0/python/demo/pde/mixed-poisson/python/documentation.html#implementation

sorry!

Just for completeness

you could also concatenate the lists of DirichletBC's to one flat list

bc1=[DirichletBC(MX.sub(0).sub(0),(0.0),left)]
bc2=[DirichletBC(MX.sub(0).sub(1),(0.0),left)]
bc3=[DirichletBC(MX.sub(0).sub(2),(0.0),left)]
bc = bc1 + bc2 + bc3

Jan

Thank you so much Jan that solved my problem.

However, according to what you said that the tensor is somewhat linearized in a row-major style

       / T_11 T_12 T_13 \
T = | T_21 T_22 T_23 |
       \ T_31 T_32 T_33 /

T is indexed as
T_linearized = [T_11 T_12 T_13 T_21 T_22 T_23 T_31 T_32 T_33]

this way I only need one index to access the elements of the tensor but the variational form I am using is in the form of T[i,j], (i and j are indices) is that represent the same thing, using two indices in the variational form and using one index in the function space?

and also when define the stress tensor StressT(u) [i,j] as

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

Is the StressT(u)[i,j] linearized as well?

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

Do I have to use one index to access the elements of StressT(u)?

Hi,

sorry for the late reply.

> However, according to what you said that the tensor is somewhat linearized
> in a row-major style
>
> / T_11 T_12 T_13 \
> T = | T_21 T_22 T_23 |
> \ T_31 T_32 T_33 /
>
> T is indexed as
> T_linearized = [T_11 T_12 T_13 T_21 T_22 T_23 T_31 T_32 T_33]

here I meant only the TensorFunctionSpace. Sorry, I wasn't clear about this.

> this way I only need one index to access the elements of the tensor but the
> variational form I am using is in the form of T[i,j], (i and j are indices)
> is that represent the same thing, using two indices in the variational form
> and using one index in the function space?

The function space and the form are not the same thing.

The function space basically defines properties of the solution, like how the
solution can vary over an element and some other properties, e.g. if the
derivative of the solution is continuous or not.

The form specifies what equation (the weak form) the solution should fulfill.

(For more please look in some book on finite elements. A good start, with
references to some books, may be
http://fenicsproject.org/documentation/tutorial/fundamentals.html)

> and also when define the stress tensor StressT(u) [i,j] as
>
> def StressT(u):
> return as_tensor((Lambda*StrainT(u)[k,k]*delta[i,j]+2*mu*StrainT(u)[i,j]),[i,j])
>
> Is the StressT(u)[i,j] linearized as well?
>
> #variational form
> F1 = StressT(u)[i,j]*w[i].dx(j)*dx - w[i]*tr[i]*ds(1) - w[i]*Rho*bf[i]*dx
> F2 = (StressTS(u_Strain)[i,j] - StressT(u)[i,j])*w_Strain[i,j]*dx
> F=F1+F2
>
> Do I have to use one index to access the elements of StressT(u)?

No, in the form you must use two indices (StressT(u)[i,j]) for 2nd order
tensors, one index for vectors and no index for scalars. This is like in math
and like you did it. For indexing of the function space on the other hand you
must use one index. That is somewhat strange but the way it is and I cannot
tell you why, sorry.

I hope this helps you.

Jan

Yes, Thank you so much.