Error to assemble_system

Question

Dear Fenics users,

I am trying to assembly a linear system, but fenics is giving me the following error message:

Traceback (most recent call last):
File "", line 1, in
File "/usr/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 540, in runfile
execfile(filename, namespace)
File "/home/fenics-doll/Desktop/FENICs_RUN/PlaneStrain/Beam_test2.py", line 117, in
A, b = assemble_system(a, L, bcs);
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 243, in assemble_system
A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
raise TypeError("Invalid form type %s" % (type(form),))
TypeError: Invalid form type

Here is my code:

from dolfin import * 
import numpy as np

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True

ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


# Create mesh and define function space
w = 1000 # 1000 mm comprimento, x
h = 250  # 250 mm largura, y

# Element Average Size
el_size = 10;   

nx = int(w/el_size);
ny = int(h/el_size);

# Mesh for cantilever beam
p0 = Point(0,0);
p1 = Point(w, h);

p_load = Point(w, h/2);

model = RectangleMesh(p0, p1, nx, ny)

# Compute all entities and connectivity.
model.init()

V = VectorFunctionSpace(model, "Lagrange", 2)

# Define Dirichlet boundary
# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class load_point(SubDomain):
    def inside(self, x, on_boundary):
        delta = 2*el_size
        return near(x[0],w) and (x[1] > h/2 - delta) and (x[1] < h/2 + delta) 

left = Left()
pt_load = load_point()

dim = model.topology().dim()-1;

sub_domains = MeshFunction("size_t", model, dim)

# Set all sub_domains in 0
sub_domains.set_all(0)

# Mark the sub_domain selection
left.mark(sub_domains, 1)
pt_load.mark(sub_domains, 2)

#plot(sub_domains_double)

ds = Measure("ds")[sub_domains]

# Elasticity parameters
E, nu = 3000.0, 0.4
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

Load_Area = 2*el_size;

force_x = 0;
force_y = -500;

press_x = force_x/Load_Area;
press_y = force_y/Load_Area;

g_B = Constant((0.0, 0.0))            # Body force per unit volume
g_T = Constant((press_x, press_y))    # Traction force on the boundary

# Define functions
u = TrialFunction(V)                        # Incremental displacement
v = TestFunction(V)                        # Test function 

#Kinematics
I = Identity(u.geometric_dimension());  # Identity tensor
F = variable(I + grad(u));              # Deformation Gradient
invF = inv(F);
J = det(F);
C = variable(F.T*F);                # Right Cauchy-Green tensor
Ic = tr(C);     # Invariants of deformation tensors

# Stored strain energy density (compressible neo-Hookean)
phi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2;

Pd = diff(phi,F);

Pa = (lmbda/2)*(J**2-1)*invF.T + mu*(F-invF).T;

a = inner(Pd,grad(v));
L = inner(g_B,v)*dx + inner(g_T,v)*ds(2);

# Apply the boundary condition
left_end_displacement = Constant((0.0, 0.0));

bcs = [ DirichletBC(V, left_end_displacement, left) ];

# Assembly the system
A, b = assemble_system(a, L, bcs);

Any idea what's going on?

Best Regards,

Comments

Do you need a measure for 'a' - i.e. *dx or something

---

Yes, I am sorry, I forgot *dx,

But even including *dx another error appear:

raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs (Argument(VectorElement('Lagrange', Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_0', data=''), 2, dim=2, quad_scheme=None), 1, None),).

Here is the modified code:

# -*- coding: utf-8 -*-
"""
Created on Fri Nov 27 13:47:19 2015

@author: Ricardo Lahuerta
"""

from dolfin import * 
import numpy as np


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True

ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


# Create mesh and define function space
w = 1000 # 1000 mm comprimento, x
h = 250  # 250 mm largura, y

# Element Average Size
el_size = 10;   

nx = int(w/el_size);
ny = int(h/el_size);

# Mesh for cantilever beam
p0 = Point(0,0);
p1 = Point(w, h);

p_load = Point(w, h/2);

model = RectangleMesh(p0, p1, nx, ny)

# Compute all entities and connectivity.
model.init()

V = VectorFunctionSpace(model, "Lagrange", 2)

# Define Dirichlet boundary
# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class load_point(SubDomain):
    def inside(self, x, on_boundary):
        delta = 2*el_size
        return near(x[0],w) and (x[1] > h/2 - delta) and (x[1] < h/2 + delta) 

left = Left()
pt_load = load_point()

dim = model.topology().dim()-1;

sub_domains = MeshFunction("size_t", model, dim)

# Set all sub_domains in 0
sub_domains.set_all(0)

# Mark the sub_domain selection
left.mark(sub_domains, 1)
pt_load.mark(sub_domains, 2)

ds = Measure("ds")[sub_domains]

# Elasticity parameters
E, nu = 3000.0, 0.4
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

Load_Area = 2*el_size;

force_x = 0;
force_y = -500;

press_x = force_x/Load_Area;
press_y = force_y/Load_Area;

g_B = Constant((0.0, 0.0))            # Body force per unit volume
g_T = Constant((press_x, press_y))    # Traction force on the boundary

# Define functions
u = TrialFunction(V)                        # Incremental displacement
v = TestFunction(V)                        # Test function 

#Kinematics
I = Identity(u.geometric_dimension());  # Identity tensor
F = variable(I + grad(u));              # Deformation Gradient
J = det(F);
C = variable(F.T*F);                # Right Cauchy-Green tensor
Ic = tr(C);     # Invariants of deformation tensors

# Stored strain energy density (compressible neo-Hookean)
phi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2;

Pd = diff(phi,F);

a = inner(Pd,grad(v))*dx;
L = inner(g_B,v)*dx + inner(g_T,v)*ds(2);

# Apply the boundary condition
left_end_displacement = Constant((0.0, 0.0));

bcs = [ DirichletBC(V, left_end_displacement, left) ];

# Assembly the system
A, b = assemble_system(a, L, bcs);

Any idea?

Answer 1

Hello All,

A collegue of mine solve the problem, here is the solution:

# Stored strain energy density (compressible neo-Hookean)
phi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

P = diff(phi, F)

a = inner(P, grad(v))*dx
L = inner(g_B,v)*dx + inner(g_T,v)*ds(2)

jacobian = derivative(a, u, du)

# Apply the boundary condition
left_end_displacement = Constant((0.0, 0.0))

bcs = [ DirichletBC(V, left_end_displacement, left) ]

# Assembly the system
A, b = assemble_system(jacobian, L, bcs)


# Note that in this case the second argument has to be a generic
# vector instead of the raw Function:
solve(A, u.vector(), b)