API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 23rd, 2017
66
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 14.63 KB | None | 0 0
raw download clone embed print report
  1. #define N 25
  2. #include <stdio.h>
  3. #include <math.h>
  4. #include <Core>
  5. #include <Dense>
  6. #include <IterativeLinearSolvers>
  7.  
  8. using namespace Eigen;
  9.  
  10. double b[N*N*N] = {0};
  11. double u[N*N*N] = {0};
  12. double p[N*N*N] = {0};
  13. int EqNum;
  14.  
  15.  
  16.  
  17. double F (int i, int j, int k){
  18. double H = 2.0/(N+1);
  19. double x = -1 + (i+1)*H;
  20. double y = -1 + (j+1)*H;
  21. double z = -1 + (k+1)*H;
  22.  
  23. return 100*(1-z*z)*(1-z*z)*(1-y*y)*(1-x*x);
  24. }
  25. double f(int i, int j, int k) {
  26. double H = 2.0/(N+1);
  27. double x = -1 + (i+1)*H;
  28. double y = -1 + (j+1)*H;
  29. double z = -1 + (k+1)*H;
  30.  
  31. return 100*(4 *(x*x - 1)*(y*y - 1)*(3*z*z - 1) + 2*(x*x - 1)*(1 - z*z)*(1 - z*z)+ 2*(y*y - 1)*(1 - z*z)*(1 - z*z));
  32. }
  33.  
  34. double CoefA(int EqNum, int CoefNum){
  35. int i = CoefNum / (N*N);
  36. int j = CoefNum % (N*N) / N;
  37. int k = CoefNum % N;
  38. if ((EqNum != CoefNum) && (EqNum != CoefNum + 1) && (EqNum != CoefNum - 1) && (EqNum != CoefNum + N) && (EqNum != CoefNum - N)&& (EqNum != CoefNum + N*N) && (EqNum != CoefNum - N*N)) return 0;
  39. //Внутренние узлы
  40. if ((i != 0) && (j != 0) && (k != 0) && (i != N - 1) && (j != N - 1) && (k != N - 1)) {
  41. if (EqNum == CoefNum) return -6;
  42. if (EqNum == CoefNum + N*N) return 1;
  43. if (EqNum == CoefNum - N*N) return 1;
  44. if (EqNum == CoefNum + N) return 1;
  45. if (EqNum == CoefNum - N) return 1;
  46. if (EqNum == CoefNum + 1) return 1;
  47. if (EqNum == CoefNum - 1) return 1;
  48. }
  49. //Грани
  50. if ((i == 0) && (j != 0) && (k != 0) && (j != N - 1) && (k != N - 1)) {
  51. if (EqNum == CoefNum) return -6;
  52. if (EqNum == CoefNum + N*N) return 1;
  53. if (EqNum == CoefNum - N*N) return 0;
  54. if (EqNum == CoefNum + N) return 1;
  55. if (EqNum == CoefNum - N) return 1;
  56. if (EqNum == CoefNum + 1) return 1;
  57. if (EqNum == CoefNum - 1) return 1;
  58. }
  59.  
  60. if ((i == N - 1) && (j != 0) && (k != 0) && (j != N - 1) && (k != N - 1)) {
  61.  
  62. if (EqNum == CoefNum) return -6;
  63. if (EqNum == CoefNum + N*N) return 0;
  64. if (EqNum == CoefNum - N*N) return 1;
  65. if (EqNum == CoefNum + N) return 1;
  66. if (EqNum == CoefNum - N) return 1;
  67. if (EqNum == CoefNum + 1) return 1;
  68. if (EqNum == CoefNum - 1) return 1;
  69. }
  70. if ((j == 0) && (i != 0) && (k != 0) && (i != N - 1) && (k != N - 1)) {
  71.  
  72. if (EqNum == CoefNum) return -6;
  73. if (EqNum == CoefNum + N*N) return 1;
  74. if (EqNum == CoefNum - N*N) return 1;
  75. if (EqNum == CoefNum + N) return 1;
  76. if (EqNum == CoefNum - N) return 0;
  77. if (EqNum == CoefNum + 1) return 1;
  78. if (EqNum == CoefNum - 1) return 1;
  79. }
  80.  
  81. if ((j == N - 1) && (i != 0) && (k != 0) && (i != N - 1) && (k != N - 1)) {
  82. if (EqNum == CoefNum) return -6;
  83. if (EqNum == CoefNum + N*N) return 1;
  84. if (EqNum == CoefNum - N*N) return 1;
  85. if (EqNum == CoefNum + N) return 0;
  86. if (EqNum == CoefNum - N) return 1;
  87. if (EqNum == CoefNum + 1) return 1;
  88. if (EqNum == CoefNum - 1) return 1;
  89.  
  90. }
  91.  
  92. if ((k == 0) && (i != 0) && (j != 0) && (i != N - 1) && (j != N - 1)) {
  93. if (EqNum == CoefNum) return -5;
  94. if (EqNum == CoefNum + N*N) return 1;
  95. if (EqNum == CoefNum - N*N) return 1;
  96. if (EqNum == CoefNum + N) return 1;
  97. if (EqNum == CoefNum - N) return 1;
  98. if (EqNum == CoefNum + 1) return 1;
  99. if (EqNum == CoefNum - 1) return 0;
  100. }
  101.  
  102. if ((k == N - 1) && (i != 0) && (j != 0) && (i != N - 1) && (j != N - 1)) {
  103. if (EqNum == CoefNum) return -5;
  104. if (EqNum == CoefNum + N*N) return 1;
  105. if (EqNum == CoefNum - N*N) return 1;
  106. if (EqNum == CoefNum + N) return 1;
  107. if (EqNum == CoefNum - N) return 1;
  108. if (EqNum == CoefNum + 1) return 0;
  109. if (EqNum == CoefNum - 1) return 1;
  110. }
  111.  
  112.  
  113. //Рёбра
  114.  
  115. if ((i == 0) && (j == 0) && (k != 0) && (k != N - 1)) {
  116. if (EqNum == CoefNum) return -6;
  117. if (EqNum == CoefNum + N*N) return 1;
  118. if (EqNum == CoefNum - N*N) return 0;
  119. if (EqNum == CoefNum + N) return 1;
  120. if (EqNum == CoefNum - N) return 0;
  121. if (EqNum == CoefNum + 1) return 1;
  122. if (EqNum == CoefNum - 1) return 1;
  123.  
  124. }
  125. if ((i == 0) && (j == N - 1) && (k != 0) && (k != N - 1)) {
  126. if (EqNum == CoefNum) return -6;
  127. if (EqNum == CoefNum + N*N) return 1;
  128. if (EqNum == CoefNum - N*N) return 0;
  129. if (EqNum == CoefNum + N) return 0;
  130. if (EqNum == CoefNum - N) return 1;
  131. if (EqNum == CoefNum + 1) return 1;
  132. if (EqNum == CoefNum - 1) return 1;
  133.  
  134. }
  135. if ((i == N - 1) && (j == 0) && (k != 0) && (k != N - 1)) {
  136. if (EqNum == CoefNum) return -6;
  137. if (EqNum == CoefNum + N*N) return 0;
  138. if (EqNum == CoefNum - N*N) return 1;
  139. if (EqNum == CoefNum + N) return 1;
  140. if (EqNum == CoefNum - N) return 0;
  141. if (EqNum == CoefNum + 1) return 1;
  142. if (EqNum == CoefNum - 1) return 1;
  143. }
  144. if ((i == N - 1) && (j == N - 1) && (k != 0) && (k != N - 1)) {
  145. if (EqNum == CoefNum) return -6;
  146. if (EqNum == CoefNum + N*N) return 0;
  147. if (EqNum == CoefNum - N*N) return 1;
  148. if (EqNum == CoefNum + N) return 0;
  149. if (EqNum == CoefNum - N) return 1;
  150. if (EqNum == CoefNum + 1) return 1;
  151. if (EqNum == CoefNum - 1) return 1;
  152.  
  153. }
  154.  
  155.  
  156. if ((i == 0) && (k == 0) && (j != 0) && (j != N - 1)) {
  157. if (EqNum == CoefNum) return -5;
  158. if (EqNum == CoefNum + N*N) return 1;
  159. if (EqNum == CoefNum - N*N) return 0;
  160. if (EqNum == CoefNum + N) return 1;
  161. if (EqNum == CoefNum - N) return 1;
  162. if (EqNum == CoefNum + 1) return 1;
  163. if (EqNum == CoefNum - 1) return 0;
  164. }
  165. if ((i == 0) && (k == N - 1) && (j != 0) && (j != N - 1)) {
  166. if (EqNum == CoefNum) return -5;
  167. if (EqNum == CoefNum + N*N) return 1;
  168. if (EqNum == CoefNum - N*N) return 0;
  169. if (EqNum == CoefNum + N) return 1;
  170. if (EqNum == CoefNum - N) return 1;
  171. if (EqNum == CoefNum + 1) return 0;
  172. if (EqNum == CoefNum - 1) return 1;
  173. }
  174. if ((i == N - 1) && (k == 0) && (j != 0) && (j != N - 1)) {
  175. if (EqNum == CoefNum) return -5;
  176. if (EqNum == CoefNum + N*N) return 0;
  177. if (EqNum == CoefNum - N*N) return 1;
  178. if (EqNum == CoefNum + N) return 1;
  179. if (EqNum == CoefNum - N) return 1;
  180. if (EqNum == CoefNum + 1) return 1;
  181. if (EqNum == CoefNum - 1) return 0;
  182. }
  183. if ((i == N - 1) && (k == N - 1) && (j != 0) && (j != N - 1)) {
  184. if (EqNum == CoefNum) return -5;
  185. if (EqNum == CoefNum + N*N) return 0;
  186. if (EqNum == CoefNum - N*N) return 1;
  187. if (EqNum == CoefNum + N) return 1;
  188. if (EqNum == CoefNum - N) return 1;
  189. if (EqNum == CoefNum + 1) return 0;
  190. if (EqNum == CoefNum - 1) return 1;
  191. }
  192.  
  193. if ((j == 0) && (k == 0) && (i != 0) && (i != N - 1)) {
  194. if (EqNum == CoefNum) return -5;
  195. if (EqNum == CoefNum + N*N) return 1;
  196. if (EqNum == CoefNum - N*N) return 1;
  197. if (EqNum == CoefNum + N) return 1;
  198. if (EqNum == CoefNum - N) return 0;
  199. if (EqNum == CoefNum + 1) return 1;
  200. if (EqNum == CoefNum - 1) return 0;
  201.  
  202. }
  203. if ((j == 0) && (k == N - 1) && (i != 0) && (i != N - 1)) {
  204. if (EqNum == CoefNum) return -5;
  205. if (EqNum == CoefNum + N*N) return 1;
  206. if (EqNum == CoefNum - N*N) return 1;
  207. if (EqNum == CoefNum + N) return 1;
  208. if (EqNum == CoefNum - N) return 0;
  209. if (EqNum == CoefNum + 1) return 0;
  210. if (EqNum == CoefNum - 1) return 1;
  211. }
  212. if ((j == N - 1) && (k == 0) && (i != 0) && (i != N - 1)) {
  213. if (EqNum == CoefNum) return -5;
  214. if (EqNum == CoefNum + N*N) return 1;
  215. if (EqNum == CoefNum - N*N) return 1;
  216. if (EqNum == CoefNum + N) return 0;
  217. if (EqNum == CoefNum - N) return 1;
  218. if (EqNum == CoefNum + 1) return 1;
  219. if (EqNum == CoefNum - 1) return 0;
  220. }
  221. if ((j == N - 1) && (k == N - 1) && (i != 0) && (i != N - 1)) {
  222. if (EqNum == CoefNum) return -5;
  223. if (EqNum == CoefNum + N*N) return 1;
  224. if (EqNum == CoefNum - N*N) return 1;
  225. if (EqNum == CoefNum + N) return 0;
  226. if (EqNum == CoefNum - N) return 1;
  227. if (EqNum == CoefNum + 1) return 0;
  228. if (EqNum == CoefNum - 1) return 1;
  229. }
  230. //Узлы, примыкающие к вершинам
  231.  
  232. if ((i == 0) && (j == 0) && (k == 0)) {
  233. if (EqNum == CoefNum) return -5;
  234. if (EqNum == CoefNum + N*N) return 1;
  235. if (EqNum == CoefNum - N*N) return 0;
  236. if (EqNum == CoefNum + N) return 1;
  237. if (EqNum == CoefNum - N) return 0;
  238. if (EqNum == CoefNum + 1) return 1;
  239. if (EqNum == CoefNum - 1) return 0;
  240. }
  241. if ((i == N - 1) && (j == 0) && (k == 0)) {
  242. if (EqNum == CoefNum) return -5;
  243. if (EqNum == CoefNum + N*N) return 0;
  244. if (EqNum == CoefNum - N*N) return 1;
  245. if (EqNum == CoefNum + N) return 1;
  246. if (EqNum == CoefNum - N) return 0;
  247. if (EqNum == CoefNum + 1) return 1;
  248. if (EqNum == CoefNum - 1) return 0;
  249. }
  250. if ((i == 0) && (j == N - 1) && (k == 0)) {
  251. if (EqNum == CoefNum) return -5;
  252. if (EqNum == CoefNum + N*N) return 1;
  253. if (EqNum == CoefNum - N*N) return 0;
  254. if (EqNum == CoefNum + N) return 0;
  255. if (EqNum == CoefNum - N) return 1;
  256. if (EqNum == CoefNum + 1) return 1;
  257. if (EqNum == CoefNum - 1) return 0;
  258. }
  259. if ((i == 0) && (j == 0) && (k == N - 1)) {
  260. if (EqNum == CoefNum) return -5;
  261. if (EqNum == CoefNum + N*N) return 1;
  262. if (EqNum == CoefNum - N*N) return 0;
  263. if (EqNum == CoefNum + N) return 1;
  264. if (EqNum == CoefNum - N) return 0;
  265. if (EqNum == CoefNum + 1) return 0;
  266. if (EqNum == CoefNum - 1) return 1;
  267. }
  268. if ((i == N - 1) && (j == N - 1) && (k == 0)) {
  269. if (EqNum == CoefNum) return -5;
  270. if (EqNum == CoefNum + N*N) return 0;
  271. if (EqNum == CoefNum - N*N) return 1;
  272. if (EqNum == CoefNum + N) return 0;
  273. if (EqNum == CoefNum - N) return 1;
  274. if (EqNum == CoefNum + 1) return 1;
  275. if (EqNum == CoefNum - 1) return 0;
  276. }
  277. if ((i == 0) && (j == N - 1) && (k == N - 1)) {
  278. if (EqNum == CoefNum) return -5;
  279. if (EqNum == CoefNum + N*N) return 1;
  280. if (EqNum == CoefNum - N*N) return 0;
  281. if (EqNum == CoefNum + N) return 0;
  282. if (EqNum == CoefNum - N) return 1;
  283. if (EqNum == CoefNum + 1) return 0;
  284. if (EqNum == CoefNum - 1) return 1;
  285. }
  286. if ((i == N - 1) && (j == 0) && (k == N - 1)) {
  287. if (EqNum == CoefNum) return -5;
  288. if (EqNum == CoefNum + N*N) return 0;
  289. if (EqNum == CoefNum - N*N) return 1;
  290. if (EqNum == CoefNum + N) return 1;
  291. if (EqNum == CoefNum - N) return 0;
  292. if (EqNum == CoefNum + 1) return 0;
  293. if (EqNum == CoefNum - 1) return 1;
  294. }
  295.  
  296. if (EqNum == CoefNum) return -5;
  297. if (EqNum == CoefNum + N*N) return 0;
  298. if (EqNum == CoefNum - N*N) return 1;
  299. if (EqNum == CoefNum + N) return 0;
  300. if (EqNum == CoefNum - N) return 1;
  301. if (EqNum == CoefNum + 1) return 0;
  302. if (EqNum == CoefNum - 1) return 1;
  303. }
  304.  
  305. int main(){
  306. EqNum = N*N*N;
  307. VectorXd u(EqNum), b(EqNum);
  308. SparseMatrix<double> A(EqNum,EqNum);
  309. for (int i = 0; i < EqNum; i++){
  310. b(i) = f(i / (N*N), i % (N*N) / N, i % N)*(2.0/(N+1))*(2.0/(N+1);
  311. }
  312. for (int i = 0; i < EqNum; i++){
  313. for (int j = 0; j < EqNum; j++){
  314. A(i,j) = CoefA(i,j);
  315. }
  316. BiCGSTAB<SparseMatrix<double> > solver;
  317. solver.compute(A);
  318. x = solver.solve(b);
  319.  
  320.  
  321. for (int i = 0; i < EqNum; i++) printf("%d %d %d %lf %lf %lf\n",i / (N*N),i % (N*N) / N, i % N,u[i], F(i / (N*N), i % (N*N) / N, i % N), fabs(u[i]-F(i / (N*N), i % (N*N) / N, i % N)));
  322. scanf("%d",&EqNum);
  323. return 0;
  324. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!