Skip to content

Instantly share code, notes, and snippets.

@genkuroki

genkuroki/Falling 2D Membranes Part 2.ipynb

Last active April 4, 2020 03:47
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save genkuroki/56d3063f90f3e3a16d529f8feca2ba5b to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save genkuroki/56d3063f90f3e3a16d529f8feca2ba5b to your computer and use it in GitHub Desktop.
Download ZIP
Falling 2D Membranes Part 2 rev.3
Display the source blob
Display the rendered blob
Raw
Falling 2D Membranes Part 2.ipynb
{
"cells": [
{
"metadata": {},
"cell_type": "markdown",
"source": "# 落下する膜 Part 2\n\n黒木玄\n\n2019-05-16, 2020-04-04\n\n$h, \\delta > 0$ であるとし, \n\n$$\nu_{tt} = u_{xx} + u_{yy} - g(u)\n$$\n\nを次のように差分化する:\n\n$$\n\\begin{aligned}\n&\n\\frac{u(t-\\delta h,x,y) + u(t+\\delta h,x,y) - 2u(t,x,y)}{\\delta^2 h^2} \n\\\\ &= \n\\frac{u(t,x-h,y) + u(t,x+h,y) + u(t,x,y-h) + u(t,x,y+h) - 4u(t,x,y)}{h^2} - g(u(x,y)).\n\\end{aligned}\n$$\n\nこれは次のように書き直される:\n\n$$\n\\begin{aligned}\nu(t+\\delta h,x,y) &= \n2u(t,x,y) - u(t-\\delta h,x,y) \n\\\\ & + \n\\delta^2 (u(t,x-h,y) + u(t,x+h,y) + u(t,x,y-h) + u(t,x,y+h) - 4u(t,x,y))\n\\\\ &\\, -\ng(u(t,x,y))\\delta^2 h^2.\n\\end{aligned}\n$$\n\nこれを任意の領域 $D$ 上で扱う. \n\n**注意:** $\\delta > 0$ をある程度小さくしないと数値計算に失敗する."
},
{
"metadata": {
"toc": true
},
"cell_type": "markdown",
"source": "<h1>目次<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#諸定義\" data-toc-modified-id=\"諸定義-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>諸定義</a></span><ul class=\"toc-item\"><li><span><a href=\"#領域の取り扱い\" data-toc-modified-id=\"領域の取り扱い-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>領域の取り扱い</a></span></li><li><span><a href=\"#Falling-Membrane-Problem-型の定義\" data-toc-modified-id=\"Falling-Membrane-Problem-型の定義-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Falling Membrane Problem 型の定義</a></span></li><li><span><a href=\"#境界条件をみたす状態を作る函数達\" data-toc-modified-id=\"境界条件をみたす状態を作る函数達-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>境界条件をみたす状態を作る函数達</a></span></li><li><span><a href=\"#時間発展\" data-toc-modified-id=\"時間発展-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>時間発展</a></span></li><li><span><a href=\"#数値解を求める函数\" data-toc-modified-id=\"数値解を求める函数-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>数値解を求める函数</a></span></li><li><span><a href=\"#プロット時に数値解の値の制限を計算する函数\" data-toc-modified-id=\"プロット時に数値解の値の制限を計算する函数-1.6\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>プロット時に数値解の値の制限を計算する函数</a></span></li><li><span><a href=\"#2次元プロット\" data-toc-modified-id=\"2次元プロット-1.7\"><span class=\"toc-item-num\">1.7&nbsp;&nbsp;</span>2次元プロット</a></span></li><li><span><a href=\"#3次元プロット\" data-toc-modified-id=\"3次元プロット-1.8\"><span class=\"toc-item-num\">1.8&nbsp;&nbsp;</span>3次元プロット</a></span></li></ul></li><li><span><a href=\"#使用例\" data-toc-modified-id=\"使用例-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>使用例</a></span><ul class=\"toc-item\"><li><span><a href=\"#波動方程式-(Neumann-(3))\" data-toc-modified-id=\"波動方程式-(Neumann-(3))-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>波動方程式 (Neumann (3))</a></span></li><li><span><a href=\"#波動方程式-(Neumann-(4))\" data-toc-modified-id=\"波動方程式-(Neumann-(4))-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>波動方程式 (Neumann (4))</a></span></li><li><span><a href=\"#波動方程式-(Dirichlet-(2))\" data-toc-modified-id=\"波動方程式-(Dirichlet-(2))-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>波動方程式 (Dirichlet (2))</a></span></li><li><span><a href=\"#波動方程式-(Neumann-(5))\" data-toc-modified-id=\"波動方程式-(Neumann-(5))-2.4\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>波動方程式 (Neumann (5))</a></span></li><li><span><a href=\"#波動方程式-(Neumann-(6))\" data-toc-modified-id=\"波動方程式-(Neumann-(6))-2.5\"><span class=\"toc-item-num\">2.5&nbsp;&nbsp;</span>波動方程式 (Neumann (6))</a></span></li></ul></li></ul></div>"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## 諸定義"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "using Base64\nusing Statistics\n\nusing Plots\nclibrary(:colorcet)\ndefault(fmt=:svg, size=(240, 240))\ngr()\n\nusing Measures: cm\n\nshowimg(mime, fn; scale=\"\") = open(fn) do f\n base64 = base64encode(f)\n if scale == \"\"\n display(\"text/html\", \"\"\"<img src=\"data:$mime;base64,$base64\">\"\"\")\n else\n display(\"text/html\", \"\"\"<img src=\"data:$mime;base64,$base64\" width=\"$scale\">\"\"\")\n end\nend\n\nshowvideo(mime, fn; scale=\"\") = open(fn) do f\n base64 = base64encode(f)\n display(\"text/html\", \"\"\"<video controls><source src=\"data:$mime;base64,$base64\" type=\"$mime\"></video>\"\"\")\nend",
"execution_count": 1,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 1,
"data": {
"text/plain": "showvideo (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 領域の取り扱い"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function domain(isinDomain, x, y)\n Domain = isinDomain.(x, y')\n Domain[1:end, 1] .= false\n Domain[1:end, end] .= false\n Domain[1, 1:end] .= false\n Domain[end, 1:end] .= false\n Domain\nend\n\ninteriorof(Domain) = [\n (i,j) for i in 1:size(Domain,1) for j in 1:size(Domain,2) if Domain[i,j]\n]\n\ncomplementof(Domain) = [\n (i,j) for i in 1:size(Domain,1) for j in 1:size(Domain,2) if !Domain[i,j]\n]\n\nisinBoundary(Domain, i, j) = (\n (i != 1 && Domain[i-1,j]) ||\n (i != size(Domain,1) && Domain[i+1,j]) || \n (j != 1 && Domain[i,j-1]) || \n (j != size(Domain,2) && Domain[i,j+1])\n)\n\nboundaryof(Domain, Complement) = [(i,j) for (i,j) in Complement if isinBoundary(Domain, i, j)]\n\nexteriorof(Domain, Complement) = [(i,j) for (i,j) in Complement if !isinBoundary(Domain, i, j)]",
"execution_count": 2,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 2,
"data": {
"text/plain": "exteriorof (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function boundaryformof(Domain, Boundary)\n BF = Array{Tuple{eltype(Boundary), Array{eltype(Boundary),1}}}(undef, length(Boundary))\n V = [(-1,0), (1,0), (0,-1), (0,1)]\n for (k, (i,j)) in enumerate(Boundary)\n A = [(i+p, j+q) for (p,q) in V if (1 ≤ i+p ≤ size(Domain,1)) && (1 ≤ j+q ≤ size(Domain,2)) && Domain[i+p,j+q]]\n BF[k] = ((i,j), A)\n end\n BF\nend",
"execution_count": 3,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 3,
"data": {
"text/plain": "boundaryformof (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "### test\nx = -0.6:0.2:0.6\ny = -0.5:0.2:0.5\nisinDomain(x,y) = (x^2+y^2 < 0.5^2)\nDomain = domain(isinDomain, x, y)\nInterior = interiorof(Domain)\nComplement = complementof(Domain)\nBoundary = boundaryof(Domain, Complement)\nExterior = exteriorof(Domain, Complement)\n\nheatmap(Domain', colorbar=false)\nscatter!(Interior, legend=false, color=:blue)\nscatter!(Complement, legend=false, color=:cyan, markersize=8)\nscatter!(Boundary, legend=false, color=:red)\nscatter!(Exterior, legend=false, color=:yellow)",
"execution_count": 4,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 4,
"data": {
"image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"240\" height=\"240\" viewBox=\"0 0 960 960\">\n<defs>\n <clipPath id=\"clip4700\">\n <rect x=\"0\" y=\"0\" width=\"960\" height=\"960\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip4700)\" d=\"\nM0 960 L960 960 L960 0 L0 0 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip4701\">\n <rect x=\"192\" y=\"96\" width=\"673\" height=\"673\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip4700)\" d=\"\nM113.754 847.475 L912.756 847.475 L912.756 47.2441 L113.754 47.2441 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip4702\">\n <rect x=\"113\" y=\"47\" width=\"800\" height=\"801\"/>\n </clipPath>\n</defs>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 190.209,847.475 190.209,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 297.891,847.475 297.891,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 405.573,847.475 405.573,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 513.255,847.475 513.255,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 620.937,847.475 620.937,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 728.619,847.475 728.619,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 836.302,847.475 836.302,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,761.915 912.756,761.915 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,636.093 912.756,636.093 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,510.271 912.756,510.271 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,384.448 912.756,384.448 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,258.626 912.756,258.626 \n \"/>\n<polyline clip-path=\"url(#clip4702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 113.754,132.803 912.756,132.803 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,847.475 912.756,847.475 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,847.475 113.754,47.2441 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 190.209,847.475 190.209,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 297.891,847.475 297.891,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 405.573,847.475 405.573,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 513.255,847.475 513.255,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 620.937,847.475 620.937,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 728.619,847.475 728.619,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 836.302,847.475 836.302,837.872 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,761.915 123.342,761.915 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,636.093 123.342,636.093 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,510.271 123.342,510.271 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,384.448 123.342,384.448 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,258.626 123.342,258.626 \n \"/>\n<polyline clip-path=\"url(#clip4700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 113.754,132.803 123.342,132.803 \n \"/>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 190.209, 894.275)\" x=\"190.209\" y=\"894.275\">1</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 297.891, 894.275)\" x=\"297.891\" y=\"894.275\">2</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 405.573, 894.275)\" x=\"405.573\" y=\"894.275\">3</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 513.255, 894.275)\" x=\"513.255\" y=\"894.275\">4</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 620.937, 894.275)\" x=\"620.937\" y=\"894.275\">5</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 728.619, 894.275)\" x=\"728.619\" y=\"894.275\">6</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 836.302, 894.275)\" x=\"836.302\" y=\"894.275\">7</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 779.415)\" x=\"104.154\" y=\"779.415\">1</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 653.593)\" x=\"104.154\" y=\"653.593\">2</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 527.771)\" x=\"104.154\" y=\"527.771\">3</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 401.948)\" x=\"104.154\" y=\"401.948\">4</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 276.126)\" x=\"104.154\" y=\"276.126\">5</text>\n</g>\n<g clip-path=\"url(#clip4700)\">\n<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 104.154, 150.303)\" x=\"104.154\" y=\"150.303\">6</text>\n</g>\n<g clip-path=\"url(#clip4702)\">\n<image width=\"754\" height=\"755\" xlink:href=\"data:image/png;base64,\niVBORw0KGgoAAAANSUhEUgAAAvIAAALzCAYAAABzx+PpAAAPbElEQVR4nO3YwQlDMQwFQSmk/5ad\nFpJDMPuZqeCBMSzamTkDAACkvG4PAAAAfifkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJ\neQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjI\nAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIe\nAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIA\nABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcA\ngCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAA\nBAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAg\nSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABB\nQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS\n8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQ\nBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8\nAAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQe/bA+DJzjm3JwBct7u3J8AjucgDAECQ\nkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKE\nPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTk\nAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEP\nAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkA\nAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMA\nQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAA\ngoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQ\nJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAg\nIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJ\neQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjI\nAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIe\nAACChDwAAATtzJzbI/jeOZ4LAPif3b09gS+5yAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAI\nEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQ\nkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKE\nPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTk\nAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEP\nAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkA\nAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMA\nQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAA\ngoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQ\nJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAg\nIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJ\neQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjI\nAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIe\nAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIA\nABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcA\ngCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAA\nBAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAg\nSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABB\nQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS\n8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQ\nBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8\nAAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQB\nACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8A\nAEE7M+f2CHiqc3wvgN29PQEeyUUeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEP\nAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkA\nAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMA\nQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAA\ngoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQ\nJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAg\nIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJ\neQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjI\nAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIe\nAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIA\nABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcA\ngCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBoZ+bcHgEAAPzGRR4AAIKEPAAABAl5AAAI\nEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQ\nkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKE\nPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTk\nAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEP\nAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAAgoQ8AAAECXkA\nAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQJOQBACBIyAMA\nQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAgIQ8AAEFCHgAA\ngoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJeQAACBLyAAAQ\nJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABAk5AEAIEjIAwBAkJAHAIAg\nIQ8AAEFCHgAAgoQ8AAAECXkAAAgS8gAAECTkAQAgSMgDAECQkAcAgCAhDwAAQUIeAACChDwAAAQJ\neQAACBLyAAAQJOQBACBIyAMAQJCQBwCAICEPAABBQh4AAIKEPAAABAl5AAAIEvIAABD0AQQBEuHd\njAIkAAAAAElFTkSuQmCC\n\" transform=\"translate(136, 70)\"/>\n</g>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"510.271\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"384.448\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"636.093\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"510.271\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"384.448\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"258.626\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"636.093\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"510.271\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"384.448\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"258.626\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"636.093\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"510.271\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"384.448\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"258.626\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"510.271\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"384.448\" r=\"14\" fill=\"#0000ff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"636.093\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"510.271\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"384.448\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"258.626\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"636.093\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"258.626\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"636.093\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"258.626\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"761.915\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"636.093\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"510.271\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"384.448\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"258.626\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"132.803\" r=\"28\" fill=\"#00ffff\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"510.271\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"384.448\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"636.093\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"258.626\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"761.915\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"405.573\" cy=\"132.803\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"761.915\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"513.255\" cy=\"132.803\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"761.915\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"620.937\" cy=\"132.803\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"636.093\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"258.626\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"510.271\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"384.448\" r=\"14\" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"761.915\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"636.093\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"258.626\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"190.209\" cy=\"132.803\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"761.915\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"297.891\" cy=\"132.803\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"761.915\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"728.619\" cy=\"132.803\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"761.915\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"636.093\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"258.626\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n<circle clip-path=\"url(#clip4702)\" cx=\"836.302\" cy=\"132.803\" r=\"14\" fill=\"#ffff00\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"1.92\"/>\n</svg>\n"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "### test\nBoundaryForm = boundaryformof(Domain, Boundary)",
"execution_count": 5,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 5,
"data": {
"text/plain": "14-element Array{Tuple{Tuple{Int64,Int64},Array{Tuple{Int64,Int64},1}},1}:\n ((1, 3), [(2, 3)])\n ((1, 4), [(2, 4)])\n ((2, 2), [(3, 2), (2, 3)])\n ((2, 5), [(3, 5), (2, 4)])\n ((3, 1), [(3, 2)])\n ((3, 6), [(3, 5)])\n ((4, 1), [(4, 2)])\n ((4, 6), [(4, 5)])\n ((5, 1), [(5, 2)])\n ((5, 6), [(5, 5)])\n ((6, 2), [(5, 2), (6, 3)])\n ((6, 5), [(5, 5), (6, 4)])\n ((7, 3), [(6, 3)])\n ((7, 4), [(6, 4)])"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### Falling Membrane Problem 型の定義"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "struct FMProblem{\n Q<:Integer,\n R<:AbstractFloat,\n Sx<:AbstractArray{R,1}, Sy<:AbstractArray{R,1}, St<:AbstractArray{R,1},\n D<:AbstractArray{Bool,2}, BC, \n TI<:AbstractArray, TB<:AbstractArray, TE<:AbstractArray, TBF<:AbstractArray,\n U1<:AbstractArray{R,2}, U0<:AbstractArray{R,2},\n Fg, Fprepare, Fupdate\n }\n nskips::Q; nsteps::Q\n Δx::R; Δy::R; Δt::R\n δ::R; δ²::R; Δt²::R\n x::Sx; y::Sy; t::St\n Domain::D; bc::BC;\n Interior::TI; Boundary::TB; Exterior::TE; BoundaryForm::TBF\n u1::U1; u0::U0\n g::Fg; prepare!::Fprepare; update!::Fupdate\nend",
"execution_count": 6,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function spanstointervals(Δx, Δt, xspan, yspan, tspan)\n nskips = round(Int, Δx/Δt)\n δ = 1/nskips\n h, Δy, Δt = Δx, Δx, δ*Δx\n δ², Δt² = δ^2, Δt^2\n \n xmax = xspan[1] + ceil((xspan[2] - xspan[1])/h) * h\n ymax = yspan[1] + ceil((yspan[2] - yspan[1])/h) * h\n tmax = tspan[1] + ceil((tspan[2] - tspan[1])/h) * h\n x = range(xspan[1], xmax, step=h)\n y = range(yspan[1], ymax, step=h)\n t = range(tspan[1], tmax, step=h)\n nsteps = length(t) - 1\n\n nskips, δ, h, Δy, Δt, δ², Δt², x, y, t, nsteps\nend\n\nfunction domainsetting(Domain)\n Complement = complementof(Domain)\n Interior = interiorof(Domain)\n Boundary = boundaryof(Domain, Complement)\n Exterior = exteriorof(Domain, Complement)\n BoundaryForm = boundaryformof(Domain, Boundary)\n\n Complement, Interior, Boundary, Exterior, BoundaryForm\nend\n\nfunction solversetting(bc)\n if lowercase(string(bc)) == \"dirichlet\"\n prepare! = prepare_dirichlet!\n update! = update_dirichlet!\n elseif lowercase(string(bc)) == \"neumann\"\n prepare! = prepare_neumann!\n update! = update_neumann!\n else\n error(\"bc must be Dirichlet or Neumann\")\n end\n prepare!, update!\nend\n\nfunction FMProblem(\n F1, F0, \n Δx::R, Δt::R,\n xspan::Tuple{R,R}, yspan::Tuple{R,R}, tspan::Tuple{R,R},\n isinDomain, bc, g\n ) where R<:AbstractFloat\n nskips, δ, h, Δy, Δt, δ², Δt², x, y, t, nsteps = spanstointervals(Δx, Δt, xspan, yspan, tspan)\n Domain = domain(isinDomain, x, y)\n Complement, Interior, Boundary, Exterior, BoundaryForm = domainsetting(Domain)\n u1 = F1.(x,y')\n u0 = F0.(x,y') \n prepare!, update! = solversetting(bc)\n \n FMProblem(\n nskips, nsteps,\n Δx, Δy, Δt,\n δ, δ², Δt²,\n x, y, t,\n Domain, bc, \n Interior, Boundary, Exterior, BoundaryForm,\n u1, u0,\n g, prepare!, update!\n )\nend\n\nfunction FMProblem(\n F1, F0, \n Δx, Δt,\n xspan, yspan, tspan,\n isinDomain, bc, g\n )\n FMProblem(\n F1, F0, \n float(Δx), float(Δt),\n float.(xspan), float.(yspan), float.(tspan),\n isinDomain, bc, g\n )\nend",
"execution_count": 7,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 7,
"data": {
"text/plain": "FMProblem"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 境界条件をみたす状態を作る函数達"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function fillexteriorwithnans!(s::FMProblem, u)\n for (i,j) in s.Exterior\n u[i,j] = NaN\n end\nend\n\nfunction fillboundarywithzeros!(s::FMProblem, u)\n for (i,j) in s.Boundary\n u[i,j] = 0\n end\nend\n\nfunction makeboundaryfree!(s::FMProblem, u)\n for ((i,j), A) in s.BoundaryForm\n u[i,j] = mean(u[k,l] for (k,l) in A)\n end\nend\n\nfunction prepare_dirichlet!(s::FMProblem, u)\n fillboundarywithzeros!(s, u)\n fillexteriorwithnans!(s, u)\nend\n\nfunction prepare_neumann!(s::FMProblem, u)\n makeboundaryfree!(s, u)\n fillexteriorwithnans!(s, u)\nend",
"execution_count": 8,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 8,
"data": {
"text/plain": "prepare_neumann! (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 時間発展"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function update_interior!(s::FMProblem, u, u1, u0)\n for (i,j) in s.Interior\n u[i,j] = 2u1[i,j] - u0[i,j] + s.δ²*(-4u1[i,j] + u1[i-1,j] + u1[i+1,j] + u1[i,j-1] + u1[i,j+1]) - s.Δt²*s.g(u1[i,j])\n end\nend\n\nfunction update_dirichlet!(s::FMProblem, u, u1, u0)\n update_interior!(s::FMProblem, u, u1, u0)\nend\n \nfunction update_neumann!(s::FMProblem, u, u1, u0)\n update_interior!(s, u, u1, u0)\n makeboundaryfree!(s, u)\nend",
"execution_count": 9,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 9,
"data": {
"text/plain": "update_neumann! (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 数値解を求める函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function solve(s::FMProblem)\n nx, ny = size(s.Domain)\n u = Array{eltype(s.u1), 3}(undef, nx, ny, s.nsteps+1)\n u[:,:,1] .= s.u1\n for j in 1:s.nsteps+1\n @views s.prepare!(s, u[:,:,j])\n end\n\n utmp = Array{eltype(s.u1), 3}(undef, nx, ny, s.nskips+2)\n utmp[:,:,end] .= s.u1\n utmp[:,:,end-1] .= s.u0\n for i in 1:s.nskips+2\n @views s.prepare!(s, utmp[:,:,i])\n end\n \n for j in 1:s.nsteps\n @. @views utmp[:,:,1] = utmp[:,:,end-1]\n @. @views utmp[:,:,2] = utmp[:,:,end]\n for i in 2:s.nskips+1\n @views s.update!(s, utmp[:,:,i+1], utmp[:,:,i], utmp[:,:,i-1])\n end\n @. @views u[:,:,j+1] = utmp[:,:,end]\n end\n u\nend",
"execution_count": 10,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 10,
"data": {
"text/plain": "solve (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### プロット時に数値解の値の制限を計算する函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "rmnans(u) = u[.!isnan.(u)]\n\nulim_auto(ut) = :auto\n\nfunction ulim_dirichlet(ut)\n umax = min(eltype(ut)(10.0), maximum(abs.(rmnans(ut))))\n umin = -umax \n umin, umax\nend\n\nfunction make_ulim_fix(u)\n uminmax = minimum(maximum(rmnans(@view u[:,:,j])) for j in 1:size(u,3))\n umaxmin = maximum(minimum(rmnans(@view u[:,:,j])) for j in 1:size(u,3))\n @show uminmax, umaxmin\n umin1 = uminmax - 0.15*(umaxmin - uminmax)\n umax1 = umaxmin + 0.15*(umaxmin - uminmax)\n @show umin1, umax1\n umin_all = minimum(rmnans(u))\n umax_all = maximum(rmnans(u))\n @show umin_all, umax_all\n umax = min(umax1, umax_all)\n umin = max(umin1, umin_all)\n @show umin, umax\n ulim_neumann(ut) = umin, umax\n ulim_neumann\nend",
"execution_count": 11,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 11,
"data": {
"text/plain": "make_ulim_fix (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 2次元プロット"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function plot_sol(u, ut, clim, color; kwargs...)\n plot(; legend=false, showaxis=false, grid=false, kwargs...)\n heatmap!(ut', color=color, clim=clim)\nend\n\nfunction plot_sol25(s::FMProblem, u; color=:RdBu, ulim=ulim_dirichlet, kwargs...)\n @views P = plot_sol(u, u[:,:,1], ulim(u[:,:,1]), color; kwargs...)\n PP = [P]\n for k in 1:24\n j = k*size(u,3) ÷ 24\n @views P = plot_sol(u, u[:,:,j], ulim(u[:,:,j]), color; kwargs...)\n push!(PP, P)\n end\n plot(PP..., size=(800, 700), layout=(5,5))\nend\n\nfunction animate_sol(s::FMProblem, u; fps=10, figsize=(220, 220), color=:RdBu, ulim=ulim_dirichlet)\n M = size(u,3)\n tskip = round(Int, 1/(fps*s.Δx))\n anim = @animate for j in [fill(1,10); 1:tskip:M; fill(M, 10)]\n plot(; size=figsize, legend=false, showaxis=false, grid=false)\n plot!(left_margin=-1cm, right_margin=-1cm, top_margin=-1cm, bottom_margin=-1cm)\n @views heatmap!(u[:,:,j]', color=color, clim=ulim(u[:,:,j]))\n end\n anim\nend\n\nfunction gif_sol(s::FMProblem, u; fps=10, figsize=(220, 220), color=:RdBu, gifname=\"test.gif\", ulim=ulim_dirichlet)\n @time anim = animate_sol(s, u; fps=fps, figsize=figsize, color=color, ulim=ulim)\n gif(anim, gifname, fps=fps)\n sleep(0.1)\n showimg(\"image/gif\", gifname)\nend",
"execution_count": 12,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 12,
"data": {
"text/plain": "gif_sol (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 3次元プロット"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function plot_sol3d(s::FMProblem, u, j; zlim=extrema(rmnans(u)), color=:RdYlBu, kwargs...)\n plot(; colorbar=false, gridalpha=0.3, tickfontsize=8, kwargs...)\n @views surface!(s.x, s.y, u[:,:,j]'; zlim=zlim, color=color)\nend\n\nfunction plot_sol3d25(s::FMProblem, u; zlim=extrema(rmnans(u)), color=:RdYlBu, kwargs...)\n P = plot_sol3d(s, u, 1; zlim=zlim, color=color, tickfontsize=4, kwargs...)\n PP = [P]\n for k in 2:25\n j = k*size(u,3) ÷ 25\n P = plot_sol3d(s, u, j; zlim=zlim, color=color, tickfontsize=4, kwargs...)\n push!(PP, P)\n end\n plot(PP..., layout=(5,5), size=(900, 800))\nend\n\nfunction animate_sol3d(s::FMProblem, u; zlim=extrema(rmnans(u)), fps=10,\n figsize=(400, 400), color=:RdYlBu, titlestr=\"\",\n kwargs...\n )\n M = size(u,3)\n tskip = max(1, round(Int, 1/(3.33fps*s.Δx)))\n @show tskip, M\n anim = @animate for j in [fill(1,10); 1:tskip:M; fill(M, 10)]\n plot_sol3d(s, u, j; zlim=zlim, color=color)\n plot!(; size=figsize, kwargs...)\n if titlestr != \"\"\n title!(titlestr, titlefontsize=10)\n end\n end\n anim\nend\n\nfunction gif_sol3d(s::FMProblem, u; zlim=extrema(rmnans(u)), fps=10,\n figsize=(400, 400), color=:RdYlBu, titlestr=\"\", gifname=\"test3d.gif\", kwargs...\n )\n @time anim = animate_sol3d(s, u; zlim=zlim, fps=fps,\n figsize=figsize, color=color, titlestr=titlestr,\n kwargs...\n )\n gif(anim, gifname, fps=fps)\n sleep(0.1)\n showimg(\"image/gif\", gifname)\nend",
"execution_count": 13,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 13,
"data": {
"text/plain": "gif_sol3d (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## 使用例"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 波動方程式 (Neumann (3))"
},
{
"metadata": {
"scrolled": false,
"trusted": true
},
"cell_type": "code",
"source": "### x軸とt軸の刻み幅\n###\n### Δt は Δx の正の整数分の1にする. 以下では10分の1にしてある.\n###\nΔx, Δt = 1/100, 1/2000\n\n### x軸とy軸とt軸の範囲\nxspan = (-1,1)\nyspan = (-1,1)\ntspan = (0,10)\n\n### 領域内で真になり, 領域外で偽になる函数\nisinDomain(x,y) = true\n\n### 境界条件 :Neumann または :Dirichlet\nbc = :Neumann\n\n### 高さ u での重力加速度\ng(u) = 0.0\n\n### 初期値\n###\n### F1(x, y) = u(t₀, x, y)\n### F0(x, y) = u(t₀ - Δt, x, y)\n###\nF1(x,y) = exp(-((x+y+1)/√2)^2/(2*0.02^2))\nF0(x,y) = F1(x + Δt/√2, y + Δt/√2)\n\n### 以上のデータで決まる偏微分方程式の問題を作成\n@time s6 = FMProblem(F1, F0, Δx, Δt, xspan, yspan, tspan, isinDomain, bc, g)\n\n### その問題の数値解を計算\n@time u6 = solve(s6)\n\n### プロット\n@time plot_sol25(s6, u6)",
"execution_count": 14,
"outputs": [
{
"output_type": "stream",
"text": " 12.979496 seconds (30.36 M allocations: 8.266 GiB, 9.83% gc time)\n 11.463166 seconds (35.38 M allocations: 1.143 GiB, 3.79% gc time)\n 2.515727 seconds (6.24 M allocations: 316.824 MiB, 3.44% gc time)\n",
"name": "stdout"
},
{
"output_type": "execute_result",
"execution_count": 14,
"data": {
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment