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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
| /** # This file include mathematical operator for coord struct */
typedef struct ten{
coord x;
coord y;
#if dimension == 3
coord z;
#endif
}ten;
/** declaration of the null coord*/
const coord Coord_nul = {0};
const tens tens_nul = {{0}};
#define Identity\
{{1,0,0},\
{0,1,0},\
{0,0,1}}
const tens I = Identity;
/** and operator on coord */
coord and_coord(coord a,coord b){
coord c = {(int)a.x && (int)b.x,(int)a.y &&(int) b.y,(int)a.z &&(int) b.z};
return c;
}
coord add_coord(coord a,coord b){
#if dimension == 2
coord c = {a.x+b.x, a.y+b.y};
#elif dimension == 3
coord c = {a.x+b.x, a.y+b.y, a.z+b.z};
#endif
return c;
}
coord diff_coord(coord a,coord b){
#if dimension == 2
coord c = {a.x-b.x,a.y-b.y};
#elif dimension == 3
coord c = {a.x-b.x, a.y-b.y, a.z-b.z};
#endif
return c;
}
coord mult_coord(coord a,double b){
#if dimension == 2
coord c = {a.x*b, a.y*b};
#elif dimension == 3
coord c = {a.x*b, a.y*b, a.z*b};
#endif
return c;
}
coord div_coord(coord a,double b){
#if dimension == 2
coord c = {a.x/b, a.y/b};
#elif dimension == 3
coord c = {a.x/b, a.y/b, a.z/b};
#endif
return c;
}
tens div_tens(tens a,double b){
einstein_sum(i,j){
a.i.j /= b;
}
return a;
}
tens mult_tens(tens a,double b){
einstein_sum(i,j){
a.i.j *= b;
}
return a;
}
tens add_tens(tens A,tens B){
einstein_sum(i,j){
A.i.j += B.i.j;
}
return A;
}
double tens_norm(tens A){
double norm = 0;
einstein_sum(k,l){
norm = A.k.l * A.k.l;
}
norm = sqrt(norm);
return norm;
}
double normL2_coord(coord a){
#if dimension == 2
double c = sqrt(pow(a.x,2)+pow(a.y,2));
#elif dimension == 3
double c = sqrt(pow(a.x,2)+pow(a.y,2)+pow(a.z,2));
#endif
return c;
}
double normL2_coord_sq(coord a){
#if dimension == 2
double c = sq(a.x)+sq(a.y);
#elif dimension == 3
double c = sq(a.x)+sq(a.y)+sq(a.z);
#endif
return c;
}
double normL2_ten(ten a){
#if dimension == 2
double c = sqrt(pow(a.x.x,2)+pow(a.y.y,2)+pow(a.x.y,2)+pow(a.y.x,2));
#elif dimension == 3
double c = sqrt(pow(a.x.x,2)+pow(a.y.y,2)+pow(a.z.z,2)+pow(a.x.y,2)+pow(a.y.x,2)+pow(a.x.z,2)+pow(a.z.x,2)+pow(a.z.y,2)+pow(a.y.z,2));
#endif
return c;
}
//only in 3D for the cross product
#if dimension == 3
coord cross_coord(coord a,coord b){
coord c = {a.y*b.z-a.z*b.y,
a.z*b.x-a.x*b.z,
a.x*b.y-a.y*b.x};
#elif dimension == 2
double cross_coord(coord a,coord b){
double c = a.x*b.y-a.y*b.x;
#endif
return c;
}
coord EigenValue(ten A){
double Delta = pow(A.x.x,2.) - 2.*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x;
coord Eig = {0.,0.};
if(Delta > 0){
Eig.x = 1./2.*( A.x.x + A.y.y - sqrt( Delta ));
Eig.y = 1./2.*( A.x.x + A.y.y + sqrt( Delta ));
}
return Eig;
}
ten EigenVectors(ten A){
ten V;
if(A.y.x != 0){
V.x.x = -1./(2.*A.y.x)*( - A.x.x + A.y.y + sqrt( pow(A.x.x,2.) - 2.*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x ) );
V.x.y = 1.;
V.y.x = -1./(2.*A.y.x)*( - A.x.x + A.y.y - sqrt( pow(A.x.x,2.) - 2*A.x.x*A.y.y+pow(A.y.y,2.) +4.*A.x.y*A.y.x ) );
V.y.y = 1.;
}else if((A.x.x-A.y.y)!=0){
V.x.x = 1.;
V.x.y = 0.;
V.y.x = -A.x.y/(A.x.x-A.y.y);
V.y.y = 1.;
}else{
V.x.x = 1.;
V.x.y = 0.;
V.y.x = 0.;
V.y.y = 1.;
}
V.x = div_coord(V.x,normL2_coord(V.x));
V.y = div_coord(V.y,normL2_coord(V.y));
return V;
}
/** Set the coordinate of a coord at the positive side of the domain */
coord POS_perio(coord pos,coord per){
foreach_dimension(){
if(per.x) pos.x = pos.x>0?pos.x:pos.x+Ls;
}
return pos;
}
/** compute the periodic distance between two coord */
coord dist_perio_coord(coord a,coord b){
foreach_dimension(){
if(a.x>0) b.x = (b.x<a.x-Ls/2)*Ls+b.x;
else b.x = (b.x>a.x+Ls/2)*(-Ls)+b.x;
}
return diff_coord(a,b);
}
double dist_perio(coord a,coord b){
foreach_dimension(){
if(a.x>0) b.x = (b.x<a.x-Ls/2)*Ls+b.x;
else b.x = (b.x>a.x+Ls/2)*(-Ls)+b.x;
}
return normL2_coord(diff_coord(a,b));
}
/** compute the tranpsoe of a tensor */
tens transpsoe(tens A){
tens B;
einstein_sum(i,j){
B.i.j = A.j.i;
}
return B;
}
tens sym(tens A){
tens B;
einstein_sum(i,j){
B. i.j = 0.5*(A.i.j + A.j.i);
}
return B;
}
tens coord_prod(coord A,coord B){
tens C;
einstein_sum(i,j){
C.i.j = A.i*B.j;
}
return C;
}
/** compute the symetric part of a tensor */
|
