Thanks for your reply Nate! So I prepared a minimal code that still fails in parallel, albeit with a different error. For reference I follow a similar approach as mentioned in COMET for creating a periodic map between the dofs, and the hyperelasticity example in demos.
from dolfin import *
import numpy as np
import os
comm = MPI.comm_world
set_log_level(20)
workdir = os.path.join(os.getcwd(), "DiscourseTrial", "CR")
os.makedirs(workdir, exist_ok=True)
parameters["ghost_mode"] = "shared_facet"
parameters["form_compiler"]["quadrature_degree"] = 5
parameters["form_compiler"]["representation"]="uflacs"
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "lu", #lu or gmres or cg or bicgstab 'preconditioner: ilu, amg, jacobi'
"preconditioner": "petsc_amg",
"maximum_iterations": 50,
"report": True,
"line_search": 'basic',
"error_on_nonconvergence": False,
"relative_tolerance":1e-7,
"absolute_tolerance":1e-8}}
msh = UnitCubeMesh(6, 6, 6) # tried with a periodic mesh
class PeriodicBoundary(SubDomain):
def __init__(self, vertices, tolerance=DOLFIN_EPS):
""" vertices stores the coordinates of the 4 unit cell corners"""
SubDomain.__init__(self, tolerance)
self.tol = tolerance
self.vv = vertices
self.a1 = self.vv[1,:]-self.vv[0,:] # first vector generating periodicity
self.a2 = self.vv[2,:]-self.vv[0,:] # second vector generating periodicity
self.a3 = self.vv[3,:]-self.vv[0,:] # third vector generating periodicity
def inside(self, x, on_boundary):
"""
return True if on left, bottom, or back faces
and not on one of the top, front or right faces
or associate edges (and vertices) as defined below
"""
# faces
left = near(x[0],self.vv[0,0])
bottom = near(x[1],self.vv[0,1])
back = near(x[2],self.vv[0,2])
right = near(x[0],self.vv[1,0])
top = near(x[1],self.vv[2,1])
front = near(x[2],self.vv[3,2])
# line-segments (bottom 4; top 4; vertical 4)
bottom_front = bottom and front
bottom_right = bottom and right
top_left = top and left
top_back = top and back
left_front = left and front
right_back = right and back
return bool((left or back or bottom) and
(not( (top_left) or (left_front) or (top_back) or (right_back) or (bottom_right) or (bottom_front))) and on_boundary)
# return bool( bool(left and not((top_left) or (left_front))) or bool(back and not((top_back) or (right_back))) or
# bool(bottom and not((bottom_right) or (bottom_front))) and on_boundary )
def map(self, x, y):
""" Mapping the right boundary to left and top to bottom"""
# faces
right = near(x[0],self.vv[1,0])
top = near(x[1],self.vv[2,1])
front = near(x[2],self.vv[3,2])
# line-segments
top_right = top and right
top_front = top and front
right_front = right and front
point_6 = right and front and top
if point_6:
y[0] = x[0] - (self.a1[0] + self.a2[0] + self.a3[0])
y[1] = x[1] - (self.a1[1] + self.a2[1] + self.a3[1])
y[2] = x[2] - (self.a1[2] + self.a2[2] + self.a3[2])
elif top_right:
y[0] = x[0] - (self.a1[0] + self.a2[0])
y[1] = x[1] - (self.a1[1] + self.a2[1])
y[2] = x[2] - (self.a1[2] + self.a2[2])
elif top_front:
y[0] = x[0] - (self.a2[0] + self.a3[0])
y[1] = x[1] - (self.a2[1] + self.a3[1])
y[2] = x[2] - (self.a2[2] + self.a3[2])
elif right_front:
y[0] = x[0] - (self.a1[0] + self.a3[0])
y[1] = x[1] - (self.a1[1] + self.a3[1])
y[2] = x[2] - (self.a1[2] + self.a3[2])
elif right:
y[0] = x[0] - (self.a1[0])
y[1] = x[1] - (self.a1[1])
y[2] = x[2] - (self.a1[2])
elif front:
y[0] = x[0] - (self.a3[0])
y[1] = x[1] - (self.a3[1])
y[2] = x[2] - (self.a3[2])
elif top:
y[0] = x[0] - (self.a2[0])
y[1] = x[1] - (self.a2[1])
y[2] = x[2] - (self.a2[2])
else:
y[0] = -1.
y[1] = -1.
y[2] = -1.
def freeEnergyNG(u, params):
I = Identity(3) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
J = det(F)
I1 = tr(C)
mu1, mu2, alph1, alph2, mu_pr = params
psiEq = 3**(1-alph1)/(2.*alph1) * mu1 * (I1**alph1 - 3**alph1) + \
3**(1-alph2)/(2.*alph2) * mu2 * (I1**alph2 - 3**alph2) - \
(mu1 + mu2) * ln(J) + mu_pr/2*(J-1)**2
# psiEq = (mu1+mu2)/2.*(I1-3.) - (mu1+mu2)*ln(J) + mu_pr/2.*(J-1)**2
return psiEq
def assignLoading(i, fVal):
Gamm_Voight = np.zeros(9)
Gamm_Voight[0] = 1. * fVal
Gamm_Voight[4] = Gamm_Voight[0]**(-0.5)
Gamm_Voight[8] = Gamm_Voight[0]**(-0.5)
return np.array([[Gamm_Voight[0]-1, Gamm_Voight[1], Gamm_Voight[2]],
[Gamm_Voight[3], Gamm_Voight[4]-1 , Gamm_Voight[5]],
[Gamm_Voight[6], Gamm_Voight[7], Gamm_Voight[8]-1 ]])
deg = 1
Lx = Ly = Lz = 1
vertices = np.array([[0., 0, 0], # 0: Origin
[Lx, 0, 0], # 1: Right
[0, Ly, 0], # 2: Top
[0, 0, Lz]]) # 3
kap_by_mu = Constant(1000)
rad = (0.2/(4./3*np.pi))**(1./3)
origin = CompiledSubDomain("near(x[0],0) && near(x[1],0) && near(x[2],0) && on_boundary")
rightPoint = CompiledSubDomain("near(x[0], Lx) && near(x[1],0) && near(x[2],0) && on_boundary", Lx = Lx)
topPoint = CompiledSubDomain("near(x[0],0) && near(x[1], Ly) && near(x[2],0) && on_boundary", Ly = Ly)
frontPoint = CompiledSubDomain("near(x[0],0) && near(x[1],0) && near(x[2],Lz) && on_boundary", Lz = Lz)
Wu = VectorElement("CR", msh.ufl_cell(), deg) # for displacement
Wlam = VectorElement("R", msh.ufl_cell(), 0)
V = FunctionSpace(msh, MixedElement([Wu, Wlam]), constrained_domain=PeriodicBoundary(vertices))
wper = Function(V)
uper, lamper = split(wper)
wper_ = TrialFunction(V)
uper_, lamper_ = split(wper)
dw_per = TestFunction(V)
du_per, dlam_per = split(dw_per)
Vplot = VectorFunctionSpace(msh, "CG", deg)
VQeplot = TensorFunctionSpace(msh, "DG", deg-1) # at least one order lower
uplot = Function(Vplot, name="disp")
# *********************Non Gaussian Parameters ***********************
mu1 = Constant(1)
mu2 = Constant(1)
mu_pr = Constant(kap_by_mu*2)
alph1 = Constant(1)
alph2 = Constant(1)
facMul = Constant(1000)
matParamsMatrix = [mu1, mu2, alph1, alph2, mu_pr]
matParamsParts = [mu1*facMul, mu2*facMul, alph1, alph2, mu_pr*facMul]
matParams = np.vstack((matParamsMatrix, matParamsParts))
qvals = Constant(20*(mu1 + mu2))
h = FacetArea(msh)
hAvg = 1./2*(h('+') + h('-'))
Hbar = Constant( ( (0,0,0), (0,0,0), (0,0,0) ))
x = SpatialCoordinate(msh)
bcOrigin = DirichletBC(V.sub(0), Constant((0., 0., 0.)), origin, method="pointwise")
bcAll = [bcOrigin]
amuv = derivative(freeEnergyNG(uper + dot(Hbar, x), matParams[0]), wper, dw_per) * dx
# amuv += inner(lamper_, du_per) * dx + inner(dlam_per, uper_) * dx
amuv += qvals/hAvg*inner(jump(uper), jump(du_per))*dS # penalize the jump across inner facets
Jac = derivative(amuv, wper, wper_)
prblm = NonlinearVariationalProblem(amuv, wper, bcAll, J=Jac)
solver = NonlinearVariationalSolver(prblm)
solver.parameters.update(snes_solver_parameters)
wfil = XDMFFile(comm, os.path.join(workdir, 'dispTrialSphereTension.xdmf'))
wfil.parameters["flush_output"] = True
wfil.parameters["functions_share_mesh"] = True
wfil.parameters["rewrite_function_mesh"] = False
Hbar.assign(Constant(assignLoading(0, 1.5)))
solver.solve()
uper, lamper = wper.split(True)
uplot.assign(project(uper + dot(Hbar, x) , Vplot))
wfil.write(uplot, 0)
gives the following error when run in parallel
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 1
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
Edit: If requred I can furnish the entire code (with the mesh) that should reproduce the exact error that I posted in the first place. However I feel that the issue is with using “CR” elements on a periodic domain, and running the code in parallel