Mixed problem with tensor function subspaces without solution

dolfinx
1

Hello FENICSx users,
I need to solve this elementary problem from strain gradient elasticity using mixed formulation

\begin{equation} \begin{split} &\int_{\Omega}\big[\boldsymbol{\sigma}:\delta\boldsymbol{\varepsilon} +\boldsymbol{\tau}:\cdotp\delta\boldsymbol{\eta} +\boldsymbol{\rho}:\big(\delta\boldsymbol{\psi} -\nabla\delta\boldsymbol{u}\big)\big]\mathrm{d}\Omega\ &\thickapprox \int_{\Omega}\boldsymbol{f}\cdotp\delta\boldsymbol{u}\mathrm{d}\Omega, \end{split} \end{equation}
\int_{\Omega}\big(\psi_{jk}-\partial_j u_k\big)\delta\rho_{jk}\mathrm{d}\Omega=0,

where

\begin{equation} \begin{split} \varepsilon_{ij} =& \frac{1}{2}\big(\partial_j u_i+\partial_i u_j\big), \\ \eta_{ijk}=& \frac{1}{2}\big(\partial_i\psi_{jk}+\partial_i \psi_{kj}\big), \\ \end{split} \end{equation}
\begin{equation} \begin{split} \sigma_{ij} =& \lambda\delta_{ij}\varepsilon_{ll} +2\mu\varepsilon_{ij}, \\ \tau_{ijk} =& l^2\partial_i\big[\lambda\delta_{jk}\psi_{ll} +\mu(\psi_{jk}+\psi_{kj})\big]. \end{split} \end{equation}

with Dirichlet boundary conditions

\boldsymbol{u}=0,\ \boldsymbol{\psi}=0\ \mathrm{on}\ y=0.

I think the problem is in the tensor function subspaces, because I used a very similar script to solve a mixed problem with scalar and vector function subspaces without any problems. Below is the script, as concise as possible, which does not give any error, but unfortunately no result either. I am afraid that there is some fundamental error that I cannot detect.

Thank you very much for any advice.

from dolfinx import mesh,fem,plot,io,default_scalar_type,default_real_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl,basix
import numpy as np

#-----[ elastic constants ]-----
#-----[ >MPa< and >nm< ]-----

mu=2.0e+04
nu=3.0e-1
nonloc=2.5e-01

#-----[ >MPa< to >N/nm**2< ]-----

mu=mu*1.e-12
lambda_=2.*mu*nu/(1.-2.*nu)

#-----[ FENICSx ]-----

def Strain(u):
  Du=ufl.grad(u)
  return 0.5*(Du+Du.T)

def Stress(u):
  return lambda_*ufl.tr(Strain(u))*ufl.Identity(2)+2*mu*Strain(u)

def HStrain(Du):
  DDu=ufl.grad(Du)
  DDuT=ufl.grad(Du.T)
  return 0.5*(DDu+DDuT)

def HStress(Du):
  DDu=ufl.grad(Du)
  DDuT=ufl.grad(Du.T)
  trDu=ufl.tr(Du)*ufl.Identity(2)
  DtrDu=ufl.grad(trDu)
  return nonloc*nonloc*(lambda_*DtrDu+mu*(DDu+DDuT))

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

f=fem.Constant(domain,default_scalar_type((0,1.e0)))

dim=domain.geometry.dim
V_el=basix.ufl.element("Lagrange",domain.topology.cell_name(),2,shape=(dim,))
VV_el=basix.ufl.element("Lagrange",domain.topology.cell_name(),1,shape=(dim,dim))
VVV_el=basix.ufl.element("DG",domain.topology.cell_name(),0,shape=(dim,dim))
M_el=basix.ufl.mixed_element([V_el,VV_el,VVV_el])
V=fem.functionspace(domain,M_el)

V0,V0_dofs=V.sub(0).collapse()
V1,V1_dofs=V.sub(1).collapse()

(u_trial,psi_trial,rho_trial)=ufl.TrialFunctions(V) 
(u_test,psi_test,rho_test)=ufl.TestFunctions(V)

facets_dim=domain.topology.dim-1

def pin_zero(x):
  return np.isclose(x[0],0.)

def displ_zero(x):
  return np.vstack((0.,0.))

def psi_zero(x):
  return np.vstack((0.,0.,0.,0.))

dofs_V0=fem.locate_dofs_geometrical(V=(V.sub(0),V0),marker=pin_zero)
dofs_V1=fem.locate_dofs_geometrical(V=(V.sub(1),V1),marker=pin_zero)
fixed_displ=fem.Function(V0)
fixed_displ.interpolate(displ_zero)
fixed_displ=fem.Function(V1)
fixed_displ.interpolate(psi_zero)
bc_displ=fem.dirichletbc(value=fixed_displ,dofs=dofs_V0,V=V.sub(0))
bc_psi=fem.dirichletbc(value=fixed_displ,dofs=dofs_V1,V=V.sub(1))

bcs=[bc_displ,bc_psi]

F1=ufl.inner(Stress(u_trial),Strain(u_test))*ufl.dx-ufl.inner(rho_trial,ufl.grad(u_test))*ufl.dx
F2=ufl.inner(HStress(psi_trial),HStrain(psi_test))*ufl.dx+ufl.inner(rho_trial,psi_test)*ufl.dx
F30=(psi_trial[0,0]-ufl.Dx(u_trial[0],0))*rho_test[0,0]*ufl.dx
F31=(psi_trial[0,1]-ufl.Dx(u_trial[0],1))*rho_test[0,1]*ufl.dx
F32=(psi_trial[1,0]-ufl.Dx(u_trial[1],0))*rho_test[1,0]*ufl.dx
F33=(psi_trial[1,1]-ufl.Dx(u_trial[1],1))*rho_test[1,1]*ufl.dx
a=F1+F2+F30+F31+F32+F33

l=ufl.dot(f,u_test)*ufl.dx

solution=LinearProblem(a=a,L=l,bcs=bcs,
                       petsc_options_prefix="grad_",
                       petsc_options={"ksp_type":"preonly","pc_type":"lu"})

uh=solution.solve()
u,psi,rho=uh.sub(0).collapse(),uh.sub(1).collapse(),uh.sub(2).collapse()
print(u.x.array)
print(psi.x.array)
1 Like