/[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 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 14906 byte(s)
Assorted spelling fixes.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26