This repository has been archived by the owner on Oct 9, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexample4-2.py
101 lines (74 loc) · 2.69 KB
/
example4-2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#!/usr/bin/env python3
# @author patrickdiehl@lsu.edu
# @author serge.prudhomme@polymtl.ca
# @date 10/05/2019
import numpy as np
import sys
import matplotlib.pyplot as plt
from decimal import Decimal
import matplotlib
matplotlib.rcParams["text.usetex"] = True
matplotlib.rcParams["font.size"] = 12
matplotlib.rcParams['text.latex.preamble'] = [
r'\usepackage{xfrac}']
eps=0.01
#############################################################################
# Solve the system
#############################################################################
def solve(M,f):
return np.linalg.solve(M,f)
#############################################################################
# Exact solution
#############################################################################
def exactSolution(x):
return x-((np.exp(-(1-x)/eps)-np.exp(-1/eps))/(1-np.exp(-1/eps)))
#############################################################################
# Loading
#############################################################################
def force(x):
return np.exp((x-1)/eps)/(eps*eps*(1-np.exp(-1/eps)))
#############################################################################
# Assemble the stiffness matrix for the varibale horizon model (VHM)
#############################################################################
def VHM(n,h):
MVHM = np.zeros([n,n])
MVHM[0][0] = 1.
MVHM[1][0] = -8.
MVHM[1][1] = 16.
MVHM[1][2] = -8.
for i in range(2,n-2):
MVHM[i][i-2] = -1.
MVHM[i][i-1] = -4.
MVHM[i][i] = 10.
MVHM[i][i+1] = -4.
MVHM[i][i+2] = -1.
MVHM[n-2][n-1] = -8.
MVHM[n-2][n-2] = 16.
MVHM[n-2][n-3] = -8.
MVHM[n-1][n-1] = 1
MVHM *= 1./(8.*h*h)
return MVHM
#############################################################################
#Computation
#############################################################################
print("n,h,VHM")
markers = ['s','o','x','.']
size= [2,4,8,10]
for i in range(6,10):
n = np.power(2,i)
h = 1./n
nodes = n+1
x = np.linspace(0,1.,nodes)
load = force(x)
load[len(x)-1] = 0
u=np.linalg.solve(VHM(nodes,h),load)
print(str(n)+","+str(h)+","+str(abs(max(abs((u[1:len(u)-2]-exactSolution(x)[1:len(u)-2])/exactSolution(x)[1:len(u)-2])))))
plt.plot(x,u-exactSolution(x),label="$\delta=\sfrac{1}{2^{"+str(int(n/2))+"}}$", marker=markers[i-6], c="black",markevery=size[i-6],ls='')
plt.title(r"Example with $\epsilon=$"+str(eps)+" Solution using VHM")
plt.xlabel("x")
plt.ylabel('Error in displacement')
plt.xscale('linear')
plt.yscale('linear')
plt.grid()
plt.legend()
plt.savefig("VHM-Error-eps-"+str(eps)+"Solution.pdf",bbox_inches='tight')