Created
December 22, 2018 04:40
-
-
Save kailaix/7807f444f30518588a0dcc407d96c8cb to your computer and use it in GitHub Desktop.
Plane Stress: An Example of Assembler Programming
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # This example solves the plane stress problem | |
| # E/(2(1+ν)) ∇^2 u_1 + E/(2(1-ν)) * ∂/∂x ( ∂u_1/∂x + ∂u_2/∂y ) = f1 | |
| # E/(2(1+ν)) ∇^2 u_2 + E/(2(1-ν)) * ∂/∂y ( ∂u_1/∂x + ∂u_2/∂y ) = f2 | |
| # Exact solution = sin(πx)sin(πy) | |
| # f1, f2, σ can be evaluated analytically | |
| using LinearAlgebra | |
| using PyPlot | |
| import Plots | |
| using PyCall | |
| using SparseArrays | |
| using LinearAlgebra | |
| @pyimport numpy as np | |
| # parameters | |
| for N = [20,40,80,160,320] | |
| x = LinRange(-1.,1., N+1) | |
| y = LinRange(-1.,1., N+1) | |
| h = x[2]-x[1] | |
| E = 1. | |
| ν = 0.3 | |
| a = E/(2(1+ν)) | |
| b = E/(2(1-ν)) | |
| Ni = N+1 | |
| h = 2/N | |
| # construct geometry | |
| Y, X = np.meshgrid(x, y) | |
| function marktraction(f) | |
| dmap = zeros(Int64, Ni, Ni) | |
| for i = 1:N+1 | |
| for j = 1:N+1 | |
| if i in [1,N+1] || j in [1,N+1] | |
| if f(X[i,j], Y[i,j]) | |
| dmap[i, j] = 2 | |
| else | |
| dmap[i, j] = 1 | |
| end | |
| end | |
| end | |
| end | |
| return dmap | |
| end | |
| function fmarker(x,y) | |
| if x<-1+h/2 | |
| return false | |
| else | |
| return true | |
| end | |
| # return false | |
| end | |
| dmap = marktraction(fmarker) | |
| function f1(x,y) | |
| return -2*E*pi^2*sin(pi*x)*sin(pi*y)/(2*ν + 2) + E*(-pi^2*sin(pi*x)*sin(pi*y) + pi^2*cos(pi*x)*cos(pi*y))/(-2*ν + 2) | |
| end | |
| function f2(x,y) | |
| return -2*E*pi^2*sin(pi*x)*sin(pi*y)/(2*ν + 2) + E*(-pi^2*sin(pi*x)*sin(pi*y) + pi^2*cos(pi*x)*cos(pi*y))/(-2*ν + 2) | |
| end | |
| function σ(x, y) | |
| return [pi*sin(pi*y)*cos(pi*x) pi*sin(pi*x)*cos(pi*y)/2+pi*sin(pi*y)*cos(pi*x)/2; | |
| pi*sin(pi*x)*cos(pi*y)/2+pi*sin(pi*y)*cos(pi*x)/2 pi*sin(pi*x)*cos(pi*y)] | |
| end | |
| # code assembler | |
| ii = Int64[] | |
| jj = Int64[] | |
| vv = Float64[] | |
| function add(i,j,k,l,v,o1=0,o2=0) | |
| if i<1 || i>Ni || j<1 || j>Ni || k<1 || k>Ni || l<1 || l>Ni | |
| return | |
| end | |
| idx = i + (j-1)*Ni + o1 * Ni^2 | |
| idy = k + (l-1)*Ni + o2 * Ni^2 | |
| push!(ii, idx) | |
| push!(jj, idy) | |
| push!(vv, v) | |
| end | |
| function assemble() | |
| return sparse(ii,jj,vv,2Ni*Ni,2Ni*Ni) | |
| end | |
| function add_rhs(i,j,f1::Function,f2::Function) | |
| idx = i + (j-1)*Ni | |
| rhs[idx] += f1(X[i,j], Y[i,j]) | |
| idx = i + (j-1)*Ni + Ni^2 | |
| rhs[idx] += f2(X[i,j], Y[i,j]) | |
| end | |
| function add_rhs(i,j,f1::Array{Float64}) | |
| idx = i + (j-1)*Ni | |
| rhs[idx] += f1[1] | |
| idx = i + (j-1)*Ni + Ni^2 | |
| rhs[idx] += f1[2] | |
| end | |
| function add_derive(i,j,c,order,o1,o2) | |
| if order=="0" | |
| add(i,j,i,j,c,o1,o2) | |
| elseif order=="1" | |
| if 1<i && i<Ni | |
| add(i,j,i+1,j,c/2h,o1,o2) | |
| add(i,j,i-1,j,-c/2h,o1,o2) | |
| elseif i==1 | |
| add(i,j,i+1,j,c/h,o1,o2) | |
| add(i,j,i, j,-c/h,o1,o2) | |
| elseif i==Ni | |
| add(i,j,i, j, c/h,o1,o2) | |
| add(i,j,i-1,j,-c/h,o1,o2) | |
| end | |
| elseif order=="2" | |
| if 1<j && j<Ni | |
| add(i,j,i,j+1,c/2h,o1,o2) | |
| add(i,j,i,j-1,-c/2h,o1,o2) | |
| elseif j==1 | |
| add(i,j,i,j+1,c/h,o1,o2) | |
| add(i,j,i,j,-c/h,o1,o2) | |
| elseif j==Ni | |
| add(i,j,i, j, c/h,o1,o2) | |
| add(i,j,i, j-1,-c/h,o1,o2) | |
| end | |
| elseif order=="11" | |
| add(i,j,i+1,j,c/h^2,o1,o2) | |
| add(i,j,i-1,j,c/h^2,o1,o2) | |
| add(i,j,i,j,-2c/h^2,o1,o2) | |
| elseif order=="22" | |
| add(i,j,i,j+1,c/h^2,o1,o2) | |
| add(i,j,i,j-1,c/h^2,o1,o2) | |
| add(i,j,i,j,-2c/h^2,o1,o2) | |
| elseif order=="12" | |
| add(i,j,i-1,j+1,-c/4h^2,o1,o2) | |
| add(i,j,i-1,j-1,c/4h^2,o1,o2) | |
| add(i,j,i+1,j+1,c/4h^2,o1,o2) | |
| add(i,j,i+1,j-1,-c/4h^2,o1,o2) | |
| elseif order=="Δ" | |
| add_derive(i,j,c,"11",o1,o2) | |
| add_derive(i,j,c,"22",o1,o2) | |
| else | |
| error("not implemented yet") | |
| end | |
| end | |
| function findnormal(i,j) | |
| n1 = 0. | |
| n2 = 0. | |
| i==1 && (n1 = -1.) | |
| i==Ni && (n1 = 1.) | |
| j==1 && (n2 = -1.) | |
| j==Ni && (n2 = 1.) | |
| n = norm([n1;n2]) | |
| [n1/n; n2/n] | |
| end | |
| rhs = zeros(2Ni^2) | |
| for i = 1:Ni | |
| for j = 1:Ni | |
| if dmap[i,j]==0 | |
| add_derive(i,j,a,"Δ",0,0) | |
| add_derive(i,j,b,"11",0,0) | |
| add_derive(i,j,b,"12",0,1) | |
| add_derive(i,j,a,"Δ",1,1) | |
| add_derive(i,j,b,"12",1,0) | |
| add_derive(i,j,b,"22",1,1) | |
| add_rhs(i,j,f1,f2) | |
| elseif dmap[i,j]==2 | |
| n = findnormal(i,j) | |
| add_derive(i,j,n[1],"1",0,0) | |
| add_derive(i,j,n[2]/2,"2",0,0) | |
| add_derive(i,j,n[2]/2,"1",0,1) | |
| add_derive(i,j,n[1]/2,"1",1,1) | |
| add_derive(i,j,n[1]/2,"2",1,0) | |
| add_derive(i,j,n[2],"2",1,1) | |
| Σ = σ(X[i,j], Y[i,j]) | |
| add_rhs(i, j, Σ*n) | |
| # @show Σ*n | |
| end | |
| end | |
| end | |
| A = assemble() | |
| # # Solve the equation | |
| I = vec(dmap) .!= 1 | |
| I = [I;I] | |
| u = zeros(2(N+1)^2) | |
| u[I] = A[I,I] \ rhs[I] | |
| ue = (x,y)->sin(pi*x)*sin(pi*y) | |
| Z = ue.(X,Y) | |
| err = maximum(abs.(Z-reshape(u[1:(N+1)^2], N+1,N+1))) | |
| @show N, err | |
| end | |
Author
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This script elaborates the idea of assembler programming. We solve a plane stress problem on [-1,1]^2
E/(2(1+ν)) ∇^2 u_1 + E/(2(1-ν)) * ∂/∂x ( ∂u_1/∂x + ∂u_2/∂y ) = f1
E/(2(1+ν)) ∇^2 u_2 + E/(2(1-ν)) * ∂/∂y ( ∂u_1/∂x + ∂u_2/∂y ) = f2
The analytical solution is given by sin(πx)sin(πy).