v0.13.1
Functions
acoustic.cpp File Reference
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>
#include "../../FTensor.hpp"

Go to the source code of this file.

Functions

void initial (double P1[], double P2[], double P3[], double c[], const int N)
 
void FTen (double P1_[], double P2_[], double P3_[], double c_[], const int N)
 
void simple (double P1[], double P2[], double P3[], double c[], const int N)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ FTen()

void FTen ( double  P1_[],
double  P2_[],
double  P3_[],
double  c_[],
const int  N 
)

Definition at line 42 of file acoustic.cpp.

43{
44 FTensor::Tensor0<double *> P1(P1_), P2(P2_), P3(P3_), c(c_);
45
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}
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
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
const int N
Definition: speed_test.cpp:3

◆ initial()

void initial ( double  P1[],
double  P2[],
double  P3[],
double  c[],
const int  N 
)

Definition at line 9 of file acoustic.cpp.

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}

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 83 of file acoustic.cpp.

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 initial(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:9
void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
Definition: acoustic.cpp:42

◆ simple()

void simple ( double  P1[],
double  P2[],
double  P3[],
double  c[],
const int  N 
)
Examples
NonlinearElasticElementInterface.hpp, approx_sphere.cpp, contact.cpp, eigen_elastic.cpp, free_surface.cpp, heat_equation.cpp, heat_method.cpp, hello_world.cpp, nonlinear_elastic.cpp, photon_diffusion.cpp, plastic.cpp, plate.cpp, scalar_check_approximation.cpp, shallow_wave.cpp, thermo_elastic.cpp, and wave_equation.cpp.

Definition at line 69 of file acoustic.cpp.

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}