Alright. I just created a MWE. The same geometry. There are 2 cubes. The left cube should be solved using linear elasticity and the right one should be solved using hyperelasticity. The left face of the left cube is fixed (Marked as 1) and a load is applied the right face of the right cube (Marked as 2). The mesh is structured (6 elements on each cube). I created the mesh in GMSH and imported into the FEniCS. Here is the GMSH file:
Point(1) = {0, 0, 0,1};
Point(2) = {1, 0, 0,1};
Point(3) = {2, 0, 0,1};
Point(4) = {0, 1, 0,1};
Point(5) = {1, 1, 0,1};
Point(6) = {2, 1, 0,1};
Point(7) = {0, 0, 1,1};
Point(8) = {1, 0, 1,1};
Point(9) = {2, 0, 1,1};
Point(10) = {0, 1, 1,1};
Point(11) = {1, 1, 1,1};
Point(12) = {2, 1, 1,1};
Line(1) = {10, 11};
Line(2) = {12, 11};
Line(3) = {10, 4};
Line(4) = {4, 1};
Line(5) = {1, 7};
Line(6) = {7, 10};
Line(7) = {5, 11};
Line(8) = {6, 12};
Line(9) = {6, 5};
Line(10) = {5, 4};
Line(11) = {1, 2};
Line(12) = {2, 3};
Line(13) = {3, 9};
Line(14) = {9, 12};
Line(15) = {9, 8};
Line(16) = {8, 7};
Line(17) = {8, 2};
Line(18) = {8, 11};
Line(19) = {6, 3};
Line(20) = {2, 5};
Line Loop(21) = {3, -10, 7, -1};
Plane Surface(22) = {21};
Line Loop(23) = {3, 4, 5, 6};
Plane Surface(24) = {23};
Line Loop(25) = {1, -18, 16, 6};
Plane Surface(26) = {25};
Line Loop(27) = {4, 11, 20, 10};
Plane Surface(28) = {27};
Line Loop(29) = {11, -17, 16, -5};
Plane Surface(30) = {29};
Line Loop(31) = {18, -7, -20, -17};
Plane Surface(32) = {31};
Line Loop(33) = {18, -2, -14, 15};
Plane Surface(34) = {33};
Line Loop(35) = {7, -2, -8, 9};
Plane Surface(36) = {35};
Line Loop(37) = {20, -9, 19, -12};
Plane Surface(38) = {37};
Line Loop(39) = {15, 17, 12, 13};
Plane Surface(40) = {39};
Line Loop(41) = {14, -8, 19, 13};
Plane Surface(42) = {41};
Surface Loop(43) = {22, 24, 28, 30, 26, 32};
Volume(44) = {43};
Surface Loop(45) = {36, 34, 42, 38, 40, 32};
Volume(46) = {45};
Transfinite Surface {22,24,26,28,30,32,34,36,38,40,42};
Transfinite Line {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20} = 2;
Physical Volume(1) = {44};
Physical Volume(2) = {46};
Physical Surface(1) = {24};
Physical Surface(2) = {42};
Now the parameter \alpha ONLY exists on the variational form on the right subdomain and I want to assign 6 different values of \alpha to these elements. The global indexes of the elements on the right cube are: [6,7,8,9,10,11]
Here is the implementation based on previous suggestions:
from dolfin import *
mesh = Mesh("MESH.xml")
boundaries = MeshFunction("size_t", mesh, "MESH_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "MESH_physical_region.xml")
n = FacetNormal(mesh)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})
V = VectorFunctionSpace(mesh,"Lagrange",degree=2)
u = Function(V)
v = TestFunction(V)
nu = 0.3
force = 1000.
E = 1E6
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
return sym(grad(v))
def sigma(v):
return lmbda * tr(eps(v)) * Identity(3) + 2.0 * mu * eps(v)
bc = DirichletBC(V, Constant((0.,0.,0)), boundaries , 1)
F_lin = inner(sigma(u), eps(v))*dx(1) # left subdomain is solved using linear elastcity
# Right subdomain is solved using nonlinear hyperelasticity
I = Identity(3)
F = I + grad(u)
C = F.T*F
Ic = tr(C)
J = det(F)
DG = FunctionSpace(mesh, "DG", 0) # The parameter alpha ONLY exists on the right subdomain
alpha = Function(DG)
alpha.vector()[:] = [1000 , 200 , 700, 100, 2000 , 300] # Assigning different values of alpha to 6 elements on the right subdomain
alpha = project(alpha, DG)
lmbda_2 = 5E6
psi = (alpha/2)*(Ic - 3) - alpha*ln(J) + (lmbda_2/2)*(ln(J))**2
Pi = psi*dx(2) - inner(-force*n, u)*ds(2)
F_hyp = derivative(Pi, u, v)
F = F_lin + F_hyp
J = derivative(F, u)
solve(F==0, u, bc, J=J)
This returns this error:
*** Error: Unable to set local values of PETSc vector.
*** Reason: Size of values array is not equal to local vector size.
It seems like FEniCS expects the alpha to exist in all elements (Both subdomains) , while it only appears on the variational form on the right subdomain.
I hope I was able to show what I am trying to do. Do you see any solution to handle such problem?
Thanks!
