#
# ..    # gedit: set fileencoding=utf8 :
#
# .. raw:: html
#
#  <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><p align="center"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png"/></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>
#
# .. _ReissnerMindlinDG:
#
# ==============================================================
# Reissner-Mindlin plate with a Discontinuous-Galerkin approach
# ==============================================================
#
# -------------
# Introduction
# -------------
#
# This program solves the Reissner-Mindlin plate equations on the unit
# square with uniform transverse loading and clamped boundary conditions.
# The corresponding file can be obtained from :download:`reissner_mindlin_dg.py`.
#
# It uses a Discontinuous Galerkin interpolation for the rotation field to
# remove shear-locking issues in the thin plate limit. Details of the formulation
# can be found in [HAN2011]_.
#
# The solution for :math:`\theta_x` on the middle line of equation :math:`y=0.5`
# will look as follows for 10 elements and a stabilization parameter :math:`s=1`:
#
# .. image:: dg_rotation_N10_s1.png
#    :scale: 15%
#
#
#
# ---------------
# Implementation
# ---------------
#
#
# Material properties and loading are the same as in :ref:`ReissnerMindlinQuads`::

from dolfin import *

E = Constant(1e3)
nu = Constant(0.3)
thick = Constant(1e-2)
D = E*thick**3/(1-nu**2)/12.
F = E/2/(1+nu)*thick*5./6.
f = Constant(-thick**3)

# The unit square mesh is here divided in triangles and we get the facet MeshFunction for the integration measure :math:`\text{d}s`::

N = 40
mesh = UnitSquareMesh(N, N)
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
ds = Measure("ds", subdomain_data=facets)

# Continuous interpolation using of degree 2 is chosen for the deflection :math:`w`
# whereas the rotation field :math:`\underline{\theta}` is discretized using discontinuous linear polynomials::

We = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Te = VectorElement("DG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([We, Te]))

# Clamped boundary conditions on the lateral boundary are defined as::

def border(x, on_boundary):
    return on_boundary

bc =  [DirichletBC(V.sub(0), Constant(0.), border)]


# Standard part of the variational form is the same (without full integration)::

def strain2voigt(eps):
    return as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S):
    return as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u):
    (w, theta) = split(u)
    return sym(grad(theta))
def shear_strain(u):
    (w, theta) = split(u)
    return theta-grad(w)
def bending_moment(u):
    DD = as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
    return voigt2stress(dot(DD,strain2voigt(curv(u))))
def shear_force(u):
    return F*shear_strain(u)

u = Function(V)
u_ = TestFunction(V)
du = TrialFunction(V)


L = f*u_[0]*dx
a = inner(bending_moment(u_), curv(du))*dx + dot(shear_force(u_), shear_strain(du))*dx


# We then add the contribution of jumps in rotation across all internal facets plus
# a stabilization term involing a user-defined parameter :math:`s`::

n = FacetNormal(mesh)
h = CellVolume(mesh)
h_avg = (h('+')+h('-'))/2
stabilization = Constant(10.)

(dw, dtheta) = split(du)
(w_, theta_) = split(u_)

a -= dot(avg(dot(bending_moment(u_), n)), jump(dtheta))*dS + dot(avg(dot(bending_moment(du), n)), jump(theta_))*dS \
   - stabilization*D/h_avg*dot(jump(theta_), jump(dtheta))*dS

# Because of the clamped boundary conditions, we also need to add the corresponding
# contributions of the external facets (the imposed rotation is zero on the boundary
# so that no term arise in the linear functional)::

a -= dot(dot(bending_moment(u_), n), dtheta)*ds + dot(dot(bending_moment(du), n), theta_)*ds \
   - 2*stabilization*D/h*dot(theta_, dtheta)*ds

# We then solve for the solution and export the relevant fields to XDMF files ::

solve(a == L, u, bc)

(w, theta) = split(u)

Vw = FunctionSpace(mesh, We)
Vt = FunctionSpace(mesh, Te)
ww = u.sub(0, True)
ww.rename("Deflection", "")
tt = u.sub(1, True)
tt.rename("Rotation", "")

file_results = XDMFFile("RM_DG_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(ww, 0.)
file_results.write(tt, 0.)

# The solution is compared to the Kirchhoff analytical solution::

print("Kirchhoff deflection:", -1.265319087e-3*float(f/D))
print("Reissner-Mindlin FE deflection:", -ww(0.5, 0.5))

# For :math:`h=0.001` and 50 elements per side, one finds :math:`w_{FE} = 1.38322\text{e-5}`  against :math:`w_{\text{Kirchhoff}} = 1.38173\text{e-5}` for the thin plate solution.
#
# -----------
# References
# -----------
#
# .. [HAN2011] Peter Hansbo, David Heintz, Mats G. Larson, A finite element method with discontinuous rotations for the Mindlin-Reissner plate model, *Computer Methods in Applied Mechanics and Engineering*, 200, 5-8, 2011, pp. 638-648, https://doi.org/10.1016/j.cma.2010.09.009.