Mesh refinement with an adaptiv solution in linear elasticity

Asked by Roman Moritz on 2013-04-25

Hello Community
I try to solve the contactproblem between a wheel and an a surface. Therefor i want a meshrefinement at the contactwidth to. I tried to implement the automated adaptivity and error control. My Goalfunction is the contactpressure given by n * sigma(u) * -n.
Unfortunately i get following Error:
Error: Unable to assemble form.
*** Reason: Invalid value rank for coefficient 1 (got 1 but expecting 0). You might have forgotten to specify the value rank correctly in an Expression subclass.

But the belongs to the goalfunction which i defined as :
n = FacetNormal(mesh)
F = sigma(u1)*(-n)
Fn = dot(F, n)
scaling = 1.0/FacetArea(mesh)

Ln = scaling*Fn*ds

Fn should be rank 0 since the contacpressure is a scalar but i don't know where the failure in the implementation is.

The whole "minimal" running example is:

from dolfin import *

mesh= Mesh("Reifen2Dlow.xml");
file=File("Boundary2.pvd")
file << mesh

V=VectorFunctionSpace(mesh,"CG",2)
scalar=FunctionSpace(mesh, "CG",1)

nu = 0.3
F=4000.0
R=317.0
bfk=220.0
E=15.0
pi=3.141592654
lam=E*nu/((1 + nu)*(1 - 2*nu))
mu=E/(2*(1 + nu))

def Restboundary2(x,on_boundary):
 tol = 1E-14
 return on_boundary and sqrt(x[0]**2+x[1]**2)< sqrt((0.22+tol)**2)

d=F/(bfk/1000.0)/(E*1.0E6)*(1.0-nu*nu)*4.0/pi

u_i=Constant((0.0,-d))
bci=DirichletBC(V,u_i,Restboundary2)

b=(2.0*sqrt(F/bfk*R/pi/E*(1.0-nu*nu)))/1000.0

def Restboundary1(x,on_boundary):
 tol = 1E-14
 return on_boundary and x[1] < -0.3+tol and x[0]< b and x[0] > -b

r=0.317
class newexpression(Expression):
 def eval(self, values, x):
  values[0]=0
  values[1]=-(r-sqrt((r*r)-(x[0]*x[0])))
 def value_shape(self):
  return (2,)

u_a=newexpression()
bca=DirichletBC(V,u_a,Restboundary1)
bcall=[bca,bci]

u=TrialFunction(V)
v=TestFunction(V)

def sigma(v):
    return 2.0*mu*sym(grad(v)) +lam*tr(sym(grad(v)))*Identity(v.cell().d)

f=Constant((0,-(9.81*(20/0.3))*(1E-6)))
a = inner(sigma(u), sym(grad(v)))*dx
L=dot(f,v)*dx

u1=Function(V)

n = FacetNormal(mesh)
F = sigma(u1)*(-n)
Fn = dot(F, n)
scaling = 1.0/FacetArea(mesh)

Ln = scaling*Fn*ds

solve(a==L,u1,bcall,tol =0.01, M=Ln)

some Help with this or maybe an other solution would be really nice
thx for help

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-04-25
Last reply:
2013-05-09
Marie Rognes (meg-simula) said : #1

Hi Roman,

Could you provide the mesh you are using or reduce the example to a problem on a non-custom mesh (say a UnitSquare)?

Anders Logg (logg) said : #3

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

Anders Logg (logg) said : #4

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

Can you help with this problem?

Provide an answer of your own, or ask Roman Moritz for more information if necessary.

To post a message you must log in.