This repository has been archived by the owner on Dec 25, 2021. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
boundaries.cpp
123 lines (98 loc) · 4.53 KB
/
boundaries.cpp
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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#include "common.h"
void boundaryCond()
{
int edgeNum;
readEdges(); // reading condition types on each edge
// natural conditions
for (edgeNum = 0; edgeNum < 4; edgeNum++) // iterate over the edges: 0 - lower, 1 - right, 2 - upper, 3 - left
{
if (edges[edgeNum] == 2) // if second type
secondCond(edgeNum);
else if (edges[edgeNum] == 3) // if third type
thirdCond(edgeNum);
}
// first type condition
for (edgeNum = 0; edgeNum < 4; edgeNum++)
if (edges[edgeNum] == 1)
firstCond(edgeNum);
}
void firstCond(int edgeNum)
{
switch(edgeNum)
{
case 0: for (i = 0; i < nX; i++)
diOrig[i] = 1.0e+100, f[i] = 1.0e+100 * getAnalytical(gridX[i], gridY[0]);
break;
case 1: for (i = 0; i < nY; i++)
diOrig[nX * (i + 1) - 1] = 1.0e+100, f[nX * (i + 1) - 1] = 1.0e+100 * getAnalytical(gridX[nX - 1], gridY[i]);
break;
case 2: for (i = 0; i < nX; i++)
diOrig[nX * (nY - 1) + i] = 1.0e+100, f[nX * (nY - 1) + i] = 1.0e+100 * getAnalytical(gridX[i], gridY[nY - 1]);
break;
case 3: for (i = 0; i < nY; i++)
diOrig[nX * i] = 1.0e+100, f[nX * i] = 1.0e+100 * getAnalytical(gridX[0], gridY[i]);
break;
}
}
void secondCond(int edgeNum)
{
double b[2];
int edgeNode1 = 0, edgeNode2 = 0;
double theta1 = 0, theta2 = 0;
double h = 0;
switch(edgeNum)
{
case 0: edgeNode1 = 0, edgeNode2 = nX - 1, h = gridX[nX - 1] - gridX[0];
theta1 = - getLambda(gridX[0], gridY[0]), theta2 = - getLambda(gridX[nX - 1], gridY[0]);
break;
case 1: edgeNode1 = nX - 1, edgeNode2 = nX * nY - 1, h = gridY[nY - 1] - gridY[0];
theta1 = getLambda(gridX[nX - 1], gridY[0]), theta2 = getLambda(gridX[nX - 1], gridY[nY - 1]);
break;
case 2: edgeNode1 = nX * (nY - 1), edgeNode2 = nX * (nY - 1) + nX - 1, h = gridX[nX - 1] - gridX[0];
theta1 = getLambda(gridX[0], gridY[nY - 1]), theta2 = getLambda(gridX[nX - 1], gridY[nY - 1]);
break;
case 3: edgeNode1 = 0, edgeNode2 = nX * (nY - 1), h = gridY[nY - 1] - gridY[0];
theta1 = - getLambda(gridX[0], gridY[0]), theta2 = - getLambda(gridX[0], gridY[nY - 1]);
break;
}
b[0] = (h / 6) * (2 * theta1 + theta2);
b[1] = (h / 6) * (theta1 + 2 * theta2);
f[edgeNode1] += b[0];
f[edgeNode2] += b[1];
}
void thirdCond(int edgeNum)
{
double b[2], A[2][2];
int edgeNode1 = 0, edgeNode2 = 0;
double uBeta1 = 0, uBeta2 = 0;
double h = 0;
switch(edgeNum)
{
case 0: edgeNode1 = 0, edgeNode2 = nX - 1, h = gridX[nX - 1] - gridX[0];
uBeta1 = - getLambda(gridX[0], gridY[0]) + getBeta() * getAnalytical(gridX[0], gridY[0]);
uBeta2 = - getLambda(gridX[nX - 1], gridY[0]) + getBeta() * getAnalytical(gridX[nX - 1], gridY[0]);
break;
case 1: edgeNode1 = nX - 1, edgeNode2 = nX * nY - 1, h = gridY[nY - 1] - gridY[0];
uBeta1 = getLambda(gridX[nX - 1], gridY[0]) + getBeta() * getAnalytical(gridX[nX - 1], gridY[0]);
uBeta2 = getLambda(gridX[nX - 1], gridY[nY - 1]) + getBeta() * getAnalytical(gridX[nX - 1], gridY[nY - 1]);
break;
case 2: edgeNode1 = nX * (nY - 1), edgeNode2 = nX * (nY - 1) + nX - 1, h = gridX[nX - 1] - gridX[0];
uBeta1 = getLambda(gridX[0], gridY[nY - 1]) + getBeta() * getAnalytical(gridX[0], gridY[nY - 1]);
uBeta2 = getLambda(gridX[nX - 1], gridY[nY - 1]) + getBeta() * getAnalytical(gridX[nX - 1], gridY[nY - 1]);
break;
case 3: edgeNode1 = 0, edgeNode2 = nX * (nY - 1), h = gridY[nY - 1] - gridY[0];
uBeta1 = - getLambda(gridX[0], gridY[0]) + getBeta() * getAnalytical(gridX[0], gridY[0]);
uBeta2 = - getLambda(gridX[0], gridY[nY - 1]) + getBeta() * getAnalytical(gridX[0], gridY[nY - 1]);
break;
}
b[0] = getBeta() * (h / 6) * (2 * uBeta1 + uBeta2);
b[1] = getBeta() * (h / 6) * (uBeta1 + 2 * uBeta2);
f[edgeNode1] += b[0];
f[edgeNode2] += b[1];
A[0][0] = getBeta() * h / 3, A[0][1] = getBeta() * h / 6;
A[1][0] = getBeta() * h / 6, A[1][1] = getBeta() * h / 3;
addGlobal(A[0][0], edgeNode1, edgeNode1);
addGlobal(A[0][1], edgeNode1, edgeNode2);
addGlobal(A[1][0], edgeNode2, edgeNode1);
addGlobal(A[1][1], edgeNode2, edgeNode2);
}