/[escript]/trunk/ripley/src/Assemble_Interpolation_2D.c
ViewVC logotype

Contents of /trunk/ripley/src/Assemble_Interpolation_2D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3641 - (show annotations)
Thu Oct 27 02:16:12 2011 UTC (7 years, 5 months ago) by gross
File MIME type: text/plain
File size: 8902 byte(s)
more work on the assemblage
1 Ripley_Assemble_Interpolation_2D(Ripley_Grid *grid, Escript in, Escript out)
2 {
3 /* GENERATOR SNIP TOP */
4 if (out_data_type==RIPLEY_ELEMENTS) {
5 const double tmp0_0 = 0.10000000000000000000;
6 const double tmp0_4 = 0.056350832689629155741;
7 const double tmp0_1 = 0.012701665379258311482;
8 const double tmp0_5 = 0.25000000000000000000;
9 const double tmp0_2 = 0.78729833462074168852;
10 const double tmp0_3 = 0.44364916731037084426;
11 #pragma omp parallel for private(i,k1,k0)
12 for (k1 =0; k1 < N1; ++k1) {
13 for (k0 =0; k0 < N0; ++k0) {
14 for (i =0; i < NCOMP; ++i) {
15 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,k1, M0),NCOMP)];
16 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,k1+1, M0),NCOMP)];
17 const register double f_01 = in[INDEX2(i,INDEX2(k0,k1+1, M0),NCOMP)];
18 const register double f_00 = in[INDEX2(i,INDEX2(k0,k1, M0),NCOMP)];
19 out[INDEX3(i,0,INDEX2(k0,k1,N0),NCOMP,9)] = f_00*tmp0_2 + f_11*tmp0_1 + tmp0_0*(f_01 + f_10);
20 out[INDEX3(i,1,INDEX2(k0,k1,N0),NCOMP,9)] = tmp0_3*(f_00 + f_10) + tmp0_4*(f_01 + f_11);
21 out[INDEX3(i,2,INDEX2(k0,k1,N0),NCOMP,9)] = f_01*tmp0_1 + f_10*tmp0_2 + tmp0_0*(f_00 + f_11);
22 out[INDEX3(i,3,INDEX2(k0,k1,N0),NCOMP,9)] = tmp0_3*(f_00 + f_01) + tmp0_4*(f_10 + f_11);
23 out[INDEX3(i,4,INDEX2(k0,k1,N0),NCOMP,9)] = tmp0_5*(f_00 + f_01 + f_10 + f_11);
24 out[INDEX3(i,5,INDEX2(k0,k1,N0),NCOMP,9)] = tmp0_3*(f_10 + f_11) + tmp0_4*(f_00 + f_01);
25 out[INDEX3(i,6,INDEX2(k0,k1,N0),NCOMP,9)] = f_01*tmp0_2 + f_10*tmp0_1 + tmp0_0*(f_00 + f_11);
26 out[INDEX3(i,7,INDEX2(k0,k1,N0),NCOMP,9)] = tmp0_3*(f_01 + f_11) + tmp0_4*(f_00 + f_10);
27 out[INDEX3(i,8,INDEX2(k0,k1,N0),NCOMP,9)] = f_00*tmp0_1 + f_11*tmp0_2 + tmp0_0*(f_01 + f_10);
28 } /* close component loop i */
29 } /* close k0 loop */
30 } /* close k1 loop */
31 } else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {
32 const double tmp0_0 = 0.25000000000000000000;
33 #pragma omp parallel for private(i,k1,k0)
34 for (k1 =0; k1 < N1; ++k1) {
35 for (k0 =0; k0 < N0; ++k0) {
36 for (i =0; i < NCOMP; ++i) {
37 const register double f_01 = in[INDEX2(i,INDEX2(k0,k1+1, M0),NCOMP)];
38 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,k1, M0),NCOMP)];
39 const register double f_00 = in[INDEX2(i,INDEX2(k0,k1, M0),NCOMP)];
40 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,k1+1, M0),NCOMP)];
41 out[INDEX3(i,0,INDEX2(k0,k1,N0),NCOMP,1)] = tmp0_0*(f_00 + f_01 + f_10 + f_11);
42 } /* close component loop i */
43 } /* close k0 loop */
44 } /* close k1 loop */
45 } else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {
46 if (face_offset(0)>-1) {
47 const double tmp0_2 = 0.50000000000000000000;
48 const double tmp0_1 = 0.88729833462074168852;
49 const double tmp0_0 = 0.11270166537925831148;
50 #pragma omp parallel for private(i,k1)
51 for (k1 =0; k1 < N1; ++k1) {
52 for (i =0; i < NCOMP; ++i) {
53 const register double f_01 = in[INDEX2(i,INDEX2(0,k1+1, M0),NCOMP)];
54 const register double f_00 = in[INDEX2(i,INDEX2(0,k1, M0),NCOMP)];
55 out[INDEX3(i,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,3)] = f_00*tmp0_1 + f_01*tmp0_0;
56 out[INDEX3(i,1,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,3)] = tmp0_2*(f_00 + f_01);
57 out[INDEX3(i,2,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,3)] = f_00*tmp0_0 + f_01*tmp0_1;
58 } /* close component loop i */
59 } /* close k1 loop */
60 } /* end of face 0 */
61 if (face_offset(1)>-1) {
62 const double tmp0_2 = 0.50000000000000000000;
63 const double tmp0_1 = 0.11270166537925831148;
64 const double tmp0_0 = 0.88729833462074168852;
65 #pragma omp parallel for private(i,k1)
66 for (k1 =0; k1 < N1; ++k1) {
67 for (i =0; i < NCOMP; ++i) {
68 const register double f_10 = in[INDEX2(i,INDEX2(M0-1,k1, M0),NCOMP)];
69 const register double f_11 = in[INDEX2(i,INDEX2(M0-1,k1+1, M0),NCOMP)];
70 out[INDEX3(i,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,3)] = f_10*tmp0_0 + f_11*tmp0_1;
71 out[INDEX3(i,1,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,3)] = tmp0_2*(f_10 + f_11);
72 out[INDEX3(i,2,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,3)] = f_10*tmp0_1 + f_11*tmp0_0;
73 } /* close component loop i */
74 } /* close k1 loop */
75 } /* end of face 1 */
76 if (face_offset(2)>-1) {
77 const double tmp0_2 = 0.50000000000000000000;
78 const double tmp0_1 = 0.88729833462074168852;
79 const double tmp0_0 = 0.11270166537925831148;
80 #pragma omp parallel for private(i,k0)
81 for (k0 =0; k0 < N0; ++k0) {
82 for (i =0; i < NCOMP; ++i) {
83 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,0, M0),NCOMP)];
84 const register double f_00 = in[INDEX2(i,INDEX2(k0,0, M0),NCOMP)];
85 out[INDEX3(i,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,3)] = f_00*tmp0_1 + f_10*tmp0_0;
86 out[INDEX3(i,1,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,3)] = tmp0_2*(f_00 + f_10);
87 out[INDEX3(i,2,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,3)] = f_00*tmp0_0 + f_10*tmp0_1;
88 } /* close component loop i */
89 } /* close k0 loop */
90 } /* end of face 2 */
91 if (face_offset(3)>-1) {
92 const double tmp0_2 = 0.50000000000000000000;
93 const double tmp0_1 = 0.88729833462074168852;
94 const double tmp0_0 = 0.11270166537925831148;
95 #pragma omp parallel for private(i,k0)
96 for (k0 =0; k0 < N0; ++k0) {
97 for (i =0; i < NCOMP; ++i) {
98 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,M1-1, M0),NCOMP)];
99 const register double f_01 = in[INDEX2(i,INDEX2(k0,M1-1, M0),NCOMP)];
100 out[INDEX3(i,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,3)] = f_01*tmp0_1 + f_11*tmp0_0;
101 out[INDEX3(i,1,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,3)] = tmp0_2*(f_01 + f_11);
102 out[INDEX3(i,2,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,3)] = f_01*tmp0_0 + f_11*tmp0_1;
103 } /* close component loop i */
104 } /* close k0 loop */
105 } /* end of face 3 */
106 } else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {
107 if (face_offset(0)>-1) {
108 const double tmp0_0 = 0.50000000000000000000;
109 #pragma omp parallel for private(i,k1)
110 for (k1 =0; k1 < N1; ++k1) {
111 for (i =0; i < NCOMP; ++i) {
112 const register double f_00 = in[INDEX2(i,INDEX2(0,k1, M0),NCOMP)];
113 const register double f_01 = in[INDEX2(i,INDEX2(0,k1+1, M0),NCOMP)];
114 out[INDEX3(i,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,1)] = tmp0_0*(f_00 + f_01);
115 } /* close component loop i */
116 } /* close k1 loop */
117 } /* end of face 0 */
118 if (face_offset(1)>-1) {
119 const double tmp0_0 = 0.50000000000000000000;
120 #pragma omp parallel for private(i,k1)
121 for (k1 =0; k1 < N1; ++k1) {
122 for (i =0; i < NCOMP; ++i) {
123 const register double f_10 = in[INDEX2(i,INDEX2(M0-1,k1, M0),NCOMP)];
124 const register double f_11 = in[INDEX2(i,INDEX2(M0-1,k1+1, M0),NCOMP)];
125 out[INDEX3(i,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,1)] = tmp0_0*(f_10 + f_11);
126 } /* close component loop i */
127 } /* close k1 loop */
128 } /* end of face 1 */
129 if (face_offset(2)>-1) {
130 const double tmp0_0 = 0.50000000000000000000;
131 #pragma omp parallel for private(i,k0)
132 for (k0 =0; k0 < N0; ++k0) {
133 for (i =0; i < NCOMP; ++i) {
134 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,0, M0),NCOMP)];
135 const register double f_00 = in[INDEX2(i,INDEX2(k0,0, M0),NCOMP)];
136 out[INDEX3(i,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,1)] = tmp0_0*(f_00 + f_10);
137 } /* close component loop i */
138 } /* close k0 loop */
139 } /* end of face 2 */
140 if (face_offset(3)>-1) {
141 const double tmp0_0 = 0.50000000000000000000;
142 #pragma omp parallel for private(i,k0)
143 for (k0 =0; k0 < N0; ++k0) {
144 for (i =0; i < NCOMP; ++i) {
145 const register double f_01 = in[INDEX2(i,INDEX2(k0,M1-1, M0),NCOMP)];
146 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,M1-1, M0),NCOMP)];
147 out[INDEX3(i,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,1)] = tmp0_0*(f_01 + f_11);
148 } /* close component loop i */
149 } /* close k0 loop */
150 } /* end of face 3 */
151 } /* end of out_data_type branching
152 /* GENERATOR SNIP BOTTOM */
153 }

  ViewVC Help
Powered by ViewVC 1.1.26