Hi, I’m using Legacy fenics for solving a custom nonlinear elasticity problem with 8 total unknown variables
I’ve specified the quadrature degree using Measure when constructing the variational form which does fine on its own, but when trying to postprocess the solutions, the number of integration points shoots way up into an absurdly high number.
Removing the project() command causes the code to run just fine, so I’m sure it has to do with the function space used in the projection not using any quadrature methods. I’ve seen ways to create a “Quadrature” function space, but when implementing that (not in the code example below), it says there a mismatch of quadrature points. Is there any way to otherwise pass quadrature information when projecting to a function space?
Thanks in advance for any help! Here’s a relatively minimal working example:
from dolfin import *
import ufl
import time
import math
# Material parameters
mu = 20.0e6 # Matrix shear modulus [Pa]
nu = 0.33 # Matrix Poisson's ratio
kappa = 2*mu*(1+nu)/(3*(1-2*nu)) # Matrix bulk modulus
Lame = 2*mu*nu/(1 - 2*nu) # Matrix first Lame Parameter
r = 1e-3 # Fiber cross-section radius [m]
E_fiber = 23e9 # Fiber Young's modulus [Pa]
A = math.pi*r**2 # Fiber cross-section area [m^2]
I = math.pi*(r**4)/4 # Fiber moment of inertia [m^4]
J = math.pi*(r**4)/2 # Fiber polar moment of inertia [m^4]
k = E_fiber*A # Fiber axial stiffness [N]
F = E_fiber*I # Fiber flexural stiffness [Nm^2]
T_fiber= E_fiber*J
# ----------------------------------------------------------------
# Mesh
# ----------------------------------------------------------------
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 0.5
z_min, z_max = 0.0, 0.1
nx, ny, nz = 80, 40, 1
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
dim = mesh.geometric_dimension()
x = SpatialCoordinate(mesh)
# Define the fiber trajectory
alpha = 0.2
wavenumber = 10
beta = 2*math.pi*wavenumber
D1 = 1 / ((alpha*sin(beta*x[0]))**2 + 1)
D2 = alpha*sin(beta*x[0]) / ((alpha*sin(beta*x[0]))**2 + 1)
D3 = 0
quad_degree = 2
# ----------------------------------------------------------------
# Define function spaces
# ----------------------------------------------------------------
u_elem = VectorElement('CG', mesh.ufl_cell(), degree = 2, dim = 3) # displacement
theta_elem = VectorElement('CG', mesh.ufl_cell(), degree = 1, dim = 2) # rotation
lambda_elem = VectorElement('CG', mesh.ufl_cell(), degree = 1, dim = 2) # Force Lagrange multiplier
mixedUT = MixedElement([u_elem, theta_elem, lambda_elem])
V = FunctionSpace(mesh, mixedUT)
U, T, Lam = V.split()
U_0, U_1, U_2 = U.split()
T_0, T_1 = T.split()
# ----------------------------------------------------------------
# Define boundary conditions
# ----------------------------------------------------------------
left = CompiledSubDomain("near(x[0], mesh_xmin) && on_boundary", mesh_xmin = x_min)
right = CompiledSubDomain("near(x[0], mesh_xmax) && on_boundary", mesh_xmax = x_max)
# Dirichlet boundaries
BC_left_x = DirichletBC(U_0, Constant(0.0), left)
BC_right_x = DirichletBC(U_0,Constant(0.1), right)
BC = [BC_left_x, BC_right_x]
# ----------------------------------------------------------------
# Kinematic and kinetic terms
# ----------------------------------------------------------------
def mynorm(aa):
return sqrt(inner(aa,aa))
Dvec = as_vector([D1,D2,D3])
Dvec = Dvec / mynorm(Dvec)
def prime(aa):
return dot(grad(aa),Dvec)
def tangent_1(d):
e2 = ufl.as_vector([0, 1, 0])
t1 = ufl.cross(e2, d)
t1 = t1 / mynorm(t1)
return t1
def tangent_2(d, t1):
t2 = ufl.cross(d, t1)
t2 = t2 / mynorm(t2)
return t2
D1vec = tangent_1(Dvec)
D2vec = tangent_2(Dvec,D1vec)
R0_ufl = ufl.as_matrix([[D1vec[0], D2vec[0], Dvec[0]], [D1vec[1], D2vec[1], Dvec[1]], [D1vec[2], D2vec[2], Dvec[2]]])
def my_rotation_matrix(theta,R0):
Lambda_matrix = ufl.as_matrix([[ufl.cos(theta[1]), ufl.sin(theta[0])*ufl.sin(theta[1]), ufl.sin(theta[1]) * ufl.cos(theta[0])],
[0, ufl.cos(theta[0]), -ufl.sin(theta[0])],
[-ufl.sin(theta[1]), ufl.cos(theta[1])*ufl.sin(theta[0]), ufl.cos(theta[1]) * ufl.cos(theta[0])]])
R_hale = dot(R0,Lambda_matrix)
return dot(R_hale, R0.T)
def gamma(R):
result = R.T * prime(R)
return ufl.as_vector([-result[1,2], result[0,2], -result[0,1]])
def Ften(u):
return Identity(dim) + grad(u)
def Eten(u,R):
return R.T * Ften(u)
def J(u):
return det(Ften(u))
def fibep(u,R):
Lambda2 = inner(Eten(u,R)*Dvec,Eten(u,R)*Dvec)
return (Lambda2-1)/2
DoD = outer(Dvec,Dvec)
D1oD1 = outer(D1vec,D1vec)
D2oD2 = outer(D2vec,D2vec)
def sigma(u,theta):
R = my_rotation_matrix(theta,R0_ufl)
return mu*Eten(u,R) + (Lame*ln(J(u))-mu)*inv(Eten(u,R)).T + k*fibep(u,R)*dot(Eten(u,R),DoD)
# ----------------------------------------------------------------
# Variational problem
# ----------------------------------------------------------------
xh, x_test, delx = Function(V), TestFunction(V), TrialFunction(V)
u, theta, lam = split(xh)
eta, zeta, tau = split(x_test)
du, dtheta, dlam = split(delx)
# Relevant kinematic terms
R_ufl = my_rotation_matrix(theta,R0_ufl)
E = Eten(u , R_ufl)
epsilon = fibep(u , R_ufl)
gam = gamma(R_ufl)
# Strain energy density
psi_matrix = (mu/2)*(inner(E,E) - 3) - mu*ln(det(E)) + (Lame/2)*(ln(det(E)))**2
psi_fiber = (k/2)*(epsilon)**2 + (F/2)*(dot(gam,D1vec)**2 + dot(gam,D2vec)**2) + (T_fiber/2)*(dot(gam,Dvec)**2)
psi_constraint = lam[0]*inner(D1vec,E*Dvec) + lam[1]*inner(D2vec,E*Dvec)
# Quadrature degree and total strain energy
dx = Measure("dx", domain=mesh, metadata={"quadrature_degree": quad_degree})
Pi = (psi_matrix + psi_fiber + psi_constraint) * dx
# Define residual and Jacobian
Fres = derivative(Pi, xh, x_test)
Jac = derivative(Fres, xh, delx)
# Define variational problem & solver & output results
solve(Fres == 0, xh, BC, J=Jac,solver_parameters={'newton_solver': {'linear_solver': 'mumps', "absolute_tolerance":1e-7}})
u, theta, lam = xh.split()
# Generate space for stress and project
T_space = TensorFunctionSpace(mesh, 'CG', 1)
s1 = project(sigma(u,theta), T_space)
and here’s the output:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.100e+06 (tol = 1.000e-07) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.316e+03 (tol = 1.000e-07) r (rel) = 3.796e-04 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 4.793e+00 (tol = 1.000e-07) r (rel) = 7.856e-07 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 2.062e-05 (tol = 1.000e-07) r (rel) = 3.380e-12 (tol = 1.000e-09)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Calling FFC just-in-time (JIT) compiler, this may take some time.
WARNING: The number of integration points for each cell will be: 5841725401
Consider using the option 'quadrature_degree' to reduce the number of points