v0.9.2
acoustic.cpp
Go to the documentation of this file.
1 #include <cmath>
2 #include <cstdlib>
3 #include <iostream>
4 #include <vector>
5 using namespace std;
6 
7 #include "../../FTensor.hpp"
8 
9 void initial(double P1[], double P2[], double P3[], double c[], const int N)
10 {
11  int cr = N / 2 - 1;
12  int cc = 7.0 * N / 8.0 - 1;
13  float s2 = 64.0 * 9.0 / pow(N / 2.0, 2.0);
14 
15  for(int i = 0; i < N; ++i)
16  for(int j = 0; j < N; ++j)
17  {
18  c[i + N * j] = 0.2;
19  P1[i + N * j] = 0.0;
20  P2[i + N * j] = exp(-(pow(i - cr, 2.0) + pow(j - cc, 2.0)) * s2);
21  P3[i + N * j] = 0.0;
22  }
23 
24  const int blockLeft = 0;
25  const int blockRight = 2 * N / 5.0;
26  const int blockTop = N / 3.0;
27  const int blockBottom = 2 * N / 3.0;
28 
29  for(int i = blockTop; i < blockBottom; ++i)
30  for(int j = blockLeft; j < blockRight; ++j)
31  c[i + N * j] = 0.5;
32 
33  int channelLeft = 4 * N / 5.0;
34  int channelRight = N;
35  int channel1Height = 3 * N / 8.0;
36  int channel2Height = 5 * N / 8.0;
37  for(int i = channelLeft; i < channelRight; ++i)
38  for(int j = channel1Height; j < channel2Height; ++j)
39  c[i + N * j] = 0;
40 }
41 
42 void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
43 {
44  FTensor::Tensor0<double *> P1(P1_), P2(P2_), P3(P3_), c(c_);
45 
46  FTensor::Tensor1<int, 2> d_ijk(1, N);
47  FTensor::Tensor1<double, 2> d_xyz(1, 1);
50  for(int j = 0; j < N; ++j)
51  for(int i = 0; i < N; ++i)
52  {
53  if(i != 0 && i != N - 1 && j != 0 && j != N - 1)
54  {
55  // P3 = 2*P2 - P1;
56  P3 = 2 * P2 - P1
57  + c
58  * (dd(P2, l, m, d_ijk, d_xyz)(0, 0)
59  + dd(P2, l, m, d_ijk, d_xyz)(1, 1));
60  }
61  ++P1;
62  ++P2;
63  ++P3;
64  ++c;
65  }
66  // cout << P3_[(N/2-1)+N*((7*N)/8-1)] << endl;
67 }
68 
69 void simple(double P1[], double P2[], double P3[], double c[], const int N)
70 {
71  for(int j = 1; j < N - 1; ++j)
72  for(int i = 1; i < N - 1; ++i)
73  {
74  // P3[i+N*j] = 2*P2[i+N*j] - P1[i+N*j];
75  P3[i + N * j] = (2 - 4 * c[i + N * j]) * P2[i + N * j] - P1[i + N * j]
76  + c[i + N * j]
77  * (P2[i - 1 + N * j] + P2[i + 1 + N * j]
78  + P2[i + N * (j - 1)] + P2[i + N * (j + 1)]);
79  }
80  // cout << P3[(N/2-1)+N*((7*N)/8-1)] << endl;
81 }
82 
83 int main(int argc, char *argv[])
84 {
85  const int N = atoi(argv[1]);
86 
87  vector<double> P1(N * N), P2(N * N), P3(N * N), c(N * N);
88 
89  initial(&c[0], &P1[0], &P2[0], &P3[0], N);
90 
91  const int niters = atoi(argv[2]);
92 
93  for(int i = 0; i < niters; ++i)
94  {
95  // simple(&P1[0],&P2[0],&P3[0],&c[0],N);
96  // simple(&P2[0],&P3[0],&P1[0],&c[0],N);
97  // simple(&P3[0],&P1[0],&P2[0],&c[0],N);
98  FTen(&P1[0], &P2[0], &P3[0], &c[0], N);
99  FTen(&P2[0], &P3[0], &P1[0], &c[0], N);
100  FTen(&P3[0], &P1[0], &P2[0], &c[0], N);
101  }
102 }
void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
Definition: acoustic.cpp:42
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
void initial(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:9
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
int main(int argc, char *argv[])
Definition: acoustic.cpp:83
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
const int N
Definition: speed_test.cpp:3
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26