/[escript]/trunk/dudley/src/Assemble_gradient.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_gradient.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3490 - (show annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 14718 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 /**************************************************************/
15
16 /* assemblage jacobiens: calculate the gradient of nodal data at quadrature points */
17
18 /**************************************************************/
19
20 #include "Assemble.h"
21 #include "Util.h"
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25 /*****************************************************************/
26
27 /* Unless the loops in here get complicated again this file should be compiled for loop unrolling */
28
29 void Dudley_Assemble_gradient(Dudley_NodeFile * nodes, Dudley_ElementFile * elements,
30 escriptDataC * grad_data, escriptDataC * data)
31 {
32 size_t localGradSize = 0;
33 register dim_t e, q, l, s, n;
34 register __const double *data_array;
35 register double *grad_data_e;
36 dim_t numNodes = 0, numShapesTotal = 0, numComps, NN = 0, numDim = 0, numQuad = 0;
37 type_t data_type = getFunctionSpaceType(data);
38 bool_t reducedIntegrationOrder = FALSE;
39 Dudley_ElementFile_Jacobeans *jac = NULL;
40
41 Dudley_resetError();
42 if (nodes == NULL || elements == NULL)
43 return;
44 numComps = getDataPointSize(data);
45 NN = elements->numNodes;
46 reducedIntegrationOrder = Dudley_Assemble_reducedIntegrationOrder(grad_data);
47
48 if (data_type == DUDLEY_NODES)
49 {
50 numNodes = nodes->nodesMapping->numTargets;
51 }
52 else if (data_type == DUDLEY_REDUCED_NODES)
53 {
54 numNodes = nodes->reducedNodesMapping->numTargets;
55 }
56 else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)
57 {
58 if (elements->MPIInfo->size > 1)
59 {
60 Dudley_setError(TYPE_ERROR,
61 "Dudley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
62 return;
63 }
64 numNodes = nodes->degreesOfFreedomMapping->numTargets;
65 }
66 else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)
67 {
68 if (elements->MPIInfo->size > 1)
69 {
70 Dudley_setError(TYPE_ERROR,
71 "Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
72 return;
73 }
74 numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
75 }
76 else
77 {
78 Dudley_setError(TYPE_ERROR,
79 "Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
80 }
81
82 jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);
83 if (Dudley_noError())
84 {
85 numDim = jac->numDim;
86 numShapesTotal = jac->numShapes;
87 numQuad = jac->numQuad;
88 localGradSize = sizeof(double) * numDim * numQuad * numComps;
89 /* check the dimensions of data */
90
91 if (!numSamplesEqual(grad_data, numQuad, elements->numElements))
92 {
93 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
94 }
95 else if (!numSamplesEqual(data, 1, numNodes))
96 {
97 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: illegal number of samples of input Data object");
98 }
99 else if (numDim * numComps != getDataPointSize(grad_data))
100 {
101 Dudley_setError(TYPE_ERROR,
102 "Dudley_Assemble_gradient: illegal number of components in gradient data object.");
103 }
104 else if (!isExpanded(grad_data))
105 {
106 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_gradient: expanded Data object is expected for output data.");
107 }
108 }
109 /* now we can start */
110
111 if (Dudley_noError())
112 {
113 requireWrite(grad_data);
114 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
115 {
116 if (data_type == DUDLEY_NODES)
117 {
118 if (numDim == 1)
119 {
120 const dim_t numShapes = 2;
121 #define DIM 1
122 #pragma omp for schedule(static)
123 for (e = 0; e < elements->numElements; e++)
124 {
125 grad_data_e = getSampleDataRW(grad_data, e);
126 memset(grad_data_e, 0, localGradSize);
127 for (s = 0; s < numShapes; s++)
128 {
129 n = elements->Nodes[INDEX2(s, e, NN)];
130 data_array = getSampleDataRO(data, n);
131 for (q = 0; q < numQuad; q++)
132 {
133 #pragma ivdep
134 for (l = 0; l < numComps; l++)
135 {
136 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
137 data_array[l] *
138 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
139 }
140 }
141 }
142 }
143 #undef DIM
144 }
145 else if (numDim == 2)
146 {
147 const dim_t numShapes = 3;
148 #define DIM 2
149 #pragma omp for schedule(static)
150 for (e = 0; e < elements->numElements; e++)
151 {
152 grad_data_e = getSampleDataRW(grad_data, e);
153 memset(grad_data_e, 0, localGradSize);
154 for (s = 0; s < numShapes; s++)
155 {
156 n = elements->Nodes[INDEX2(s, e, NN)];
157 data_array = getSampleDataRO(data, n);
158 for (q = 0; q < numQuad; q++)
159 {
160 #pragma ivdep
161 for (l = 0; l < numComps; l++)
162 {
163 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
164 data_array[l] *
165 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
166 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
167 data_array[l] *
168 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
169 }
170 }
171 }
172 }
173 #undef DIM
174 }
175 else if (numDim == 3)
176 {
177 const dim_t numShapes = 4;
178 #define DIM 3
179 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
180 for (e = 0; e < elements->numElements; e++)
181 {
182 grad_data_e = getSampleDataRW(grad_data, e);
183 memset(grad_data_e, 0, localGradSize);
184 for (s = 0; s < numShapes; s++)
185 {
186 n = elements->Nodes[INDEX2(s, e, NN)];
187 data_array = getSampleDataRO(data, n);
188 for (q = 0; q < numQuad; q++)
189 {
190 #pragma ivdep
191 for (l = 0; l < numComps; l++)
192 {
193 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
194 data_array[l] *
195 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
196 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
197 data_array[l] *
198 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
199 grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
200 data_array[l] *
201 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
202 }
203 }
204 }
205 }
206 #undef DIM
207 }
208 }
209 else if (data_type == DUDLEY_REDUCED_NODES)
210 {
211 if (numDim == 1)
212 {
213 const dim_t numShapes = 2;
214 #define DIM 1
215 #pragma omp for schedule(static)
216 for (e = 0; e < elements->numElements; e++)
217 {
218 grad_data_e = getSampleDataRW(grad_data, e);
219 memset(grad_data_e, 0, localGradSize);
220 for (s = 0; s < numShapes; s++)
221 {
222 n = elements->Nodes[INDEX2(s, e, NN)];
223 data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
224 for (q = 0; q < numQuad; q++)
225 {
226 #pragma ivdep
227 for (l = 0; l < numComps; l++)
228 {
229 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
230 data_array[l] *
231 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
232 }
233 }
234 }
235 }
236 #undef DIM
237 }
238 else if (numDim == 2)
239 {
240 const dim_t numShapes = 3;
241 #define DIM 2
242 #pragma omp for schedule(static)
243 for (e = 0; e < elements->numElements; e++)
244 {
245 grad_data_e = getSampleDataRW(grad_data, e);
246 memset(grad_data_e, 0, localGradSize);
247 for (s = 0; s < numShapes; s++)
248 {
249 n = elements->Nodes[INDEX2(s, e, NN)];
250 data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
251 for (q = 0; q < numQuad; q++)
252 {
253 #pragma ivdep
254 for (l = 0; l < numComps; l++)
255 {
256 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
257 data_array[l] *
258 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
259 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
260 data_array[l] *
261 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
262 }
263 }
264 }
265 }
266 #undef DIM
267 }
268 else if (numDim == 3)
269 {
270 const dim_t numShapes = 4;
271 #define DIM 3
272 #pragma omp for schedule(static)
273 for (e = 0; e < elements->numElements; e++)
274 {
275 grad_data_e = getSampleDataRW(grad_data, e);
276 memset(grad_data_e, 0, localGradSize);
277 for (s = 0; s < numShapes; s++)
278 {
279 n = elements->Nodes[INDEX2(s, e, NN)];
280 data_array = getSampleDataRO(data, nodes->reducedNodesMapping->target[n]);
281 for (q = 0; q < numQuad; q++)
282 {
283 #pragma ivdep
284 for (l = 0; l < numComps; l++)
285 {
286 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
287 data_array[l] *
288 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
289 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
290 data_array[l] *
291 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
292 grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
293 data_array[l] *
294 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
295 }
296 }
297 }
298 }
299 #undef DIM
300 }
301 }
302 else if (data_type == DUDLEY_DEGREES_OF_FREEDOM)
303 {
304 if (numDim == 1)
305 {
306 const dim_t numShapes = 2;
307 #define DIM 1
308 #pragma omp for schedule(static)
309 for (e = 0; e < elements->numElements; e++)
310 {
311 grad_data_e = getSampleDataRW(grad_data, e);
312 memset(grad_data_e, 0, localGradSize);
313 for (s = 0; s < numShapes; s++)
314 {
315 n = elements->Nodes[INDEX2(s, e, NN)];
316 data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
317 for (q = 0; q < numQuad; q++)
318 {
319 #pragma ivdep
320 for (l = 0; l < numComps; l++)
321 {
322 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
323 data_array[l] *
324 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
325 }
326 }
327 }
328 }
329 #undef DIM
330 }
331 else if (numDim == 2)
332 {
333 const dim_t numShapes = 3;
334 #define DIM 2
335 #pragma omp for schedule(static)
336 for (e = 0; e < elements->numElements; e++)
337 {
338 grad_data_e = getSampleDataRW(grad_data, e);
339 memset(grad_data_e, 0, localGradSize);
340 for (s = 0; s < numShapes; s++)
341 {
342 n = elements->Nodes[INDEX2(s, e, NN)];
343 data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
344 for (q = 0; q < numQuad; q++)
345 {
346 #pragma ivdep
347 for (l = 0; l < numComps; l++)
348 {
349 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
350 data_array[l] *
351 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
352 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
353 data_array[l] *
354 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
355 }
356 }
357 }
358 }
359 #undef DIM
360 }
361 else if (numDim == 3)
362 {
363 const dim_t numShapes = 4;
364 #define DIM 3
365 #pragma omp for schedule(static)
366 for (e = 0; e < elements->numElements; e++)
367 {
368 grad_data_e = getSampleDataRW(grad_data, e);
369 memset(grad_data_e, 0, localGradSize);
370 for (s = 0; s < numShapes; s++)
371 {
372 n = elements->Nodes[INDEX2(s, e, NN)];
373 data_array = getSampleDataRO(data, nodes->degreesOfFreedomMapping->target[n]);
374 for (q = 0; q < numQuad; q++)
375 {
376 #pragma ivdep
377 for (l = 0; l < numComps; l++)
378 {
379 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
380 data_array[l] *
381 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
382 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
383 data_array[l] *
384 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
385 grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
386 data_array[l] *
387 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
388 }
389 }
390 }
391 }
392 #undef DIM
393 }
394 }
395 else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM)
396 {
397 if (numDim == 1)
398 {
399 const dim_t numShapes = 2;
400 #define DIM 1
401 #pragma omp for schedule(static)
402 for (e = 0; e < elements->numElements; e++)
403 {
404 grad_data_e = getSampleDataRW(grad_data, e);
405 memset(grad_data_e, 0, localGradSize);
406 for (s = 0; s < numShapes; s++)
407 {
408 n = elements->Nodes[INDEX2(s, e, NN)];
409 data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
410 for (q = 0; q < numQuad; q++)
411 {
412 #pragma ivdep
413 for (l = 0; l < numComps; l++)
414 {
415 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
416 data_array[l] *
417 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
418 }
419 }
420 }
421 }
422 #undef DIM
423 }
424 else if (numDim == 2)
425 {
426 const dim_t numShapes = 3;
427 #define DIM 2
428 #pragma omp for schedule(static)
429 for (e = 0; e < elements->numElements; e++)
430 {
431 grad_data_e = getSampleDataRW(grad_data, e);
432 memset(grad_data_e, 0, localGradSize);
433 for (s = 0; s < numShapes; s++)
434 {
435 n = elements->Nodes[INDEX2(s, e, NN)];
436 data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
437 for (q = 0; q < numQuad; q++)
438 {
439 #pragma ivdep
440 for (l = 0; l < numComps; l++)
441 {
442 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
443 data_array[l] *
444 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
445 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
446 data_array[l] *
447 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
448 }
449 }
450 }
451 }
452 #undef DIM
453 }
454 else if (numDim == 3)
455 {
456 const dim_t numShapes = 4;
457 #define DIM 3
458 #pragma omp for schedule(static)
459 for (e = 0; e < elements->numElements; e++)
460 {
461 grad_data_e = getSampleDataRW(grad_data, e);
462 memset(grad_data_e, 0, localGradSize);
463 for (s = 0; s < numShapes; s++)
464 {
465 n = elements->Nodes[INDEX2(s, e, NN)];
466 data_array = getSampleDataRO(data, nodes->reducedDegreesOfFreedomMapping->target[n]);
467 for (q = 0; q < numQuad; q++)
468 {
469 #pragma ivdep
470 for (l = 0; l < numComps; l++)
471 {
472 grad_data_e[INDEX4(l, 0, q, 0, numComps, DIM, numQuad)] +=
473 data_array[l] *
474 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
475 grad_data_e[INDEX4(l, 1, q, 0, numComps, DIM, numQuad)] +=
476 data_array[l] *
477 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
478 grad_data_e[INDEX4(l, 2, q, 0, numComps, DIM, numQuad)] +=
479 data_array[l] *
480 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, DIM, numQuad, 1)];
481 }
482 }
483 }
484 }
485 #undef DIM
486 }
487 }
488 } /* end parallel region */
489 }
490 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26