{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda config --add channels conda-forge\n",
    "!conda install -y --update-dependencies fenics mshr matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --upgrade sympy==1.1.1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def press_bottom(x, on_bound):\n",
    "    return on_bound and near(x[1], -pressHeight)\n",
    " \n",
    "def press_top(x, on_bound):\n",
    "    return on_bound and near(x[1], blockHeight + pressHeight)\n",
    " \n",
    "def border_all(x, on_bound):\n",
    "    return on_bound\n",
    "\n",
    "class CoefExp(UserExpression):\n",
    "    def value_shape(s):\n",
    "        return []\n",
    "    def eval_cell(s, v, x, c):\n",
    "        if s.d[c.index] == 0:\n",
    "            v[0] = s.v1\n",
    "        else:\n",
    "            v[0] = s.v2\n",
    "\n",
    "def Coef(domains, v1, v2):\n",
    "    W = FunctionSpace(domains.mesh(), \"DG\", 0)\n",
    "    C = CoefExp(degree=0)\n",
    "    C.d = domains\n",
    "    C.v1 = v1\n",
    "    C.v2 = v2\n",
    "    return project(C, W)\n",
    "    \n",
    "def epsilon(v):\n",
    "    return 0.5*(grad(v) + grad(v).T)\n",
    " \n",
    "def sigma(lmd, mu, v):\n",
    "    return lmd*tr(epsilon(v))*Identity(2)+2*mu*epsilon(v)\n",
    "            \n",
    "def problem(name):\n",
    "    mesh = Mesh(\"meshes/%s.xml\" % name)\n",
    "    domains = MeshFunction(\"size_t\", mesh, \"meshes/%s_domains.xml\" % name)\n",
    "    bounds = MeshFunction(\"size_t\", mesh, 1, 0)\n",
    "    AutoSubDomain(press_bottom).mark(bounds, 1)\n",
    "    AutoSubDomain(press_top).mark(bounds, 2)\n",
    "    \n",
    "    E   = Coef(domains, E1, E2)\n",
    "    nu  = Coef(domains, nu1, nu2)\n",
    "    lmd = nu * E / (1 + nu) / (1 - 2 * nu)\n",
    "    mu  = E / 2 / (1 + nu)\n",
    "    \n",
    "    ds = Measure(\"ds\", subdomain_data=bounds)\n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    bc = DirichletBC(V, Constant((0, 0)), bounds, 1)\n",
    "    a = inner(sigma(lmd, mu, u), epsilon(v)) * dx\n",
    "    L = inner(f, v) * ds(2)\n",
    "    y = Function(V)\n",
    "    solve(a == L, y, bc)\n",
    "    \n",
    "    File(\"%s/u.pvd\" % name) << y\n",
    "    File(\"%s/s.pvd\" % name) << project(sigma(lmd, mu, y), TensorFunctionSpace(mesh, 'CG', 1)) \n",
    "    plt.figure(figsize=(16, 12))\n",
    "    plot(y)\n",
    "    return y\n",
    "\n",
    "def hom_coef(lmd, mu, y, i):\n",
    "    eps11 =     assemble(epsilon(y)[0, 0]*dx)\n",
    "    eps12 = 2 * assemble(epsilon(y)[0, 1]*dx)\n",
    "    eps22 =     assemble(epsilon(y)[1, 1]*dx)\n",
    "    print(\"epsilon: \", eps11, eps12, eps22)\n",
    "    if i == 0:\n",
    "        epsj = eps11\n",
    "    elif i == 1:\n",
    "        epsj = eps12\n",
    "    else:\n",
    "        epsj = eps22\n",
    "    Ej1 = assemble(sigma(lmd, mu, y)[0, 0]*dx)/epsj\n",
    "    Ej2 = assemble(sigma(lmd, mu, y)[0, 1]*dx)/epsj\n",
    "    Ej3 = assemble(sigma(lmd, mu, y)[1, 1]*dx)/epsj\n",
    "    return Ej1, Ej2, Ej3\n",
    "    \n",
    "def hom_local(name):\n",
    "    mesh = Mesh(\"meshes/%s.xml\" % name)\n",
    "    domains = MeshFunction(\"size_t\", mesh, \"meshes/%s_domains.xml\" % name)\n",
    "    \n",
    "    E   = Coef(domains, E1, E2)\n",
    "    nu  = Coef(domains, nu1, nu2)\n",
    "    lmd = nu * E / (1 + nu) / (1 - 2 * nu)\n",
    "    mu  = E / 2 / (1 + nu)\n",
    "    \n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    bcs = [DirichletBC(V, Expression(('x[0]', '0'), degree=1), border_all),\n",
    "           DirichletBC(V, Expression(('x[1]/2', 'x[0]/2'), degree=1), border_all),\n",
    "           DirichletBC(V, Expression(('0', 'x[1]'), degree=1), border_all)]\n",
    "\n",
    "    a = inner(sigma(lmd, mu, u), epsilon(v)) * dx\n",
    "    L = inner(Constant([0, 0]), v) * ds\n",
    "    y = Function(V)\n",
    "    Eh = []\n",
    "    for i in range(len(bcs)):\n",
    "        solve(a == L, y, bcs[i])\n",
    "        File(\"%s/u%d.pvd\" % (name, i)) << y\n",
    "        File(\"%s/s%d.pvd\" % (name, i)) << project(sigma(lmd, mu, y), TensorFunctionSpace(mesh, 'CG', 1)) \n",
    "        plt.figure(figsize=(16, 12))\n",
    "        plot(y)\n",
    "        Eh.append(hom_coef(lmd, mu, y, i))\n",
    "    return Eh\n",
    "   \n",
    "class ElasticityCoefExp(UserExpression):\n",
    "    def value_shape(s):\n",
    "        return (3, 3)\n",
    "    def eval_cell(s, v, x, c):\n",
    "        i = s.d[c.index]\n",
    "        EC = s.EC2\n",
    "        if i == 0:\n",
    "            EC = s.EC0\n",
    "        elif i == 1:\n",
    "            EC = s.EC1\n",
    "        for j in range(len(v)):\n",
    "            v[j] = EC[j]\n",
    "\n",
    "def ElasticityCoef(domains, EC0, EC1, EC2):\n",
    "    C = ElasticityCoefExp(degree=0)\n",
    "    C.d = domains\n",
    "    C.EC0 = EC0\n",
    "    C.EC1 = EC1\n",
    "    C.EC2 = EC2\n",
    "    return C\n",
    "\n",
    "def hom_epsilon(v):\n",
    "    return as_vector((v[0].dx(0), (v[0].dx(1)+v[1].dx(0)), v[1].dx(1)))\n",
    "    \n",
    "def hom(name, Eh):\n",
    "    mesh = Mesh(\"meshes/%s.xml\" % name)\n",
    "    domains = MeshFunction(\"size_t\", mesh, \"meshes/%s_domains.xml\" % name)\n",
    "    bounds = MeshFunction(\"size_t\", mesh, 1, 0)\n",
    "    AutoSubDomain(press_bottom).mark(bounds, 1)\n",
    "    AutoSubDomain(press_top).mark(bounds, 2)\n",
    "    \n",
    "    lmd1 = nu1 * E1 / (1 + nu1) / (1 - 2 * nu1)\n",
    "    mu1  = E1 / 2 / (1 + nu1)\n",
    "    lmd2 = nu2 * E2 / (1 + nu2) / (1 - 2 * nu2)\n",
    "    mu2  = E2 / 2 / (1 + nu2)\n",
    "\n",
    "    EC0 = [Eh[0][0], Eh[1][0], Eh[2][0], Eh[0][1], Eh[1][1], Eh[2][1], Eh[0][2], Eh[1][2], Eh[2][2]]\n",
    "    EC1 = [lmd1+2*mu1, 0, lmd1, 0, mu1, 0, lmd1, 0, lmd1+2*mu1]\n",
    "    EC2 = [lmd2+2*mu2, 0, lmd2, 0, mu2, 0, lmd2, 0, lmd2+2*mu2]\n",
    "    print(EC0)\n",
    "    print(EC1)\n",
    "    print(EC2)\n",
    "\n",
    "    E = ElasticityCoef(domains, EC1, EC1, EC1)\n",
    "    \n",
    "    ds = Measure(\"ds\", subdomain_data=bounds)\n",
    "    V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    bc = DirichletBC(V, Constant((0, 0)), bounds, 1)\n",
    "    a = inner(E * hom_epsilon(u), hom_epsilon(v)) * dx\n",
    "    L = inner(f, v) * ds(2)\n",
    "    y = Function(V)\n",
    "    solve(a == L, y, bc)\n",
    "    \n",
    "    File(\"%s/u.pvd\" % name) << y\n",
    "    s = E * hom_epsilon(y)\n",
    "    File(\"%s/s.pvd\" % name) << project(as_matrix(((s[0], s[1]/2), (s[1]/2, s[2]))), TensorFunctionSpace(mesh, 'CG', 1)) \n",
    "    plt.figure(figsize=(16, 12))\n",
    "    plot(y)\n",
    "    return y\n",
    "    \n",
    "blockHeight = 0.5\n",
    "pressHeight = 0.1\n",
    "E1, E2 = 4e10, 2e11\n",
    "nu1, nu2 = 0.15, 0.3\n",
    "f = Constant((0, -1e5))\n",
    "\n",
    "Eh = hom_local(\"rve\")\n",
    "y = problem(\"block\")\n",
    "ys = hom(\"sparse\", Eh)\n",
    "yc = hom(\"coarse\", Eh)\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}