/[escript]/branches/domexper/dudley/src/Assemble_gradient.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_gradient.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3204 - (show annotations)
Thu Sep 23 23:59:39 2010 UTC (8 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 14112 byte(s)
More params collapsed

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
17 /* assemblage jacobiens: calculate the gradient of nodal data at quadrature points */
18
19 /**************************************************************/
20
21 #include "Assemble.h"
22 #include "Util.h"
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26 /*****************************************************************/
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 Dudley_ReferenceElement* refElement=NULL;
35 size_t localGradSize=0;
36 register dim_t e,q,l,s,n;
37 register __const double *data_array;
38 register double *grad_data_e;
39 dim_t numNodes=0, numShapesTotal=0, numComps, NN=0, numDim=0, numShapesTotal2=0, numQuad=0;
40 type_t data_type=getFunctionSpaceType(data);
41 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
42 Dudley_ElementFile_Jacobeans* jac=NULL;
43
44 Dudley_resetError();
45 if (nodes==NULL || elements==NULL) return;
46 numComps=getDataPointSize(data);
47 NN=elements->numNodes;
48 reducedIntegrationOrder=Dudley_Assemble_reducedIntegrationOrder(grad_data);
49
50 if (data_type==DUDLEY_NODES) {
51 reducedShapefunction=FALSE;
52 numNodes=nodes->nodesMapping->numTargets;
53 } else if (data_type==DUDLEY_REDUCED_NODES) {
54 reducedShapefunction=TRUE;
55 numNodes=nodes->reducedNodesMapping->numTargets;
56 } else if (data_type==DUDLEY_DEGREES_OF_FREEDOM) {
57 if (elements->MPIInfo->size>1) {
58 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
59 return;
60 }
61 reducedShapefunction=FALSE;
62 numNodes=nodes->degreesOfFreedomMapping->numTargets;
63 } else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
64 if (elements->MPIInfo->size>1) {
65 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
66 return;
67 }
68 reducedShapefunction=TRUE;
69 numNodes=nodes->reducedDegreesOfFreedomMapping->numTargets;
70 } else {
71 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
72 }
73
74 jac=Dudley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
75 refElement=Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, reducedIntegrationOrder);
76
77 if (Dudley_noError())
78 {
79 numDim=jac->numDim;
80 // numShapes=jac->BasisFunctions->Type->numShapes;
81 numShapesTotal=jac->numShapes;
82 numQuad=jac->numQuadTotal;
83 localGradSize=sizeof(double)*numDim*numQuad*numComps;
84 numShapesTotal2=refElement->BasisFunctions->Type->numShapes;
85 /* check the dimensions of data */
86
87 if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) {
88 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
89 } else if (! numSamplesEqual(data,1,numNodes)) {
90 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples of input Data object");
91 } else if (numDim*numComps!=getDataPointSize(grad_data)) {
92 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of components in gradient data object.");
93 } else if (!isExpanded(grad_data)) {
94 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: expanded Data object is expected for output data.");
95 }
96 }
97 /* now we can start */
98
99 if (Dudley_noError())
100 {
101 requireWrite(grad_data);
102 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
103 {
104 if (data_type==DUDLEY_NODES)
105 {
106 if (numDim==1)
107 {
108 const dim_t numShapes=2;
109 #define DIM 1
110 #pragma omp for schedule(static)
111 for (e=0;e<elements->numElements;e++)
112 {
113 grad_data_e=getSampleDataRW(grad_data,e);
114 memset(grad_data_e,0, localGradSize);
115 for (s=0;s<numShapes;s++)
116 {
117 n=elements->Nodes[INDEX2(s,e, NN)];
118 data_array=getSampleDataRO(data,n);
119 for (q=0;q<numQuad;q++)
120 {
121 #pragma ivdep
122 for (l=0;l<numComps;l++)
123 {
124 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
125 }
126 }
127 }
128 }
129 #undef DIM
130 }
131 else if (numDim==2)
132 {
133 const dim_t numShapes=3;
134 #define DIM 2
135 #pragma omp for schedule(static)
136 for (e=0;e<elements->numElements;e++)
137 {
138 grad_data_e=getSampleDataRW(grad_data,e);
139 memset(grad_data_e,0, localGradSize);
140 for (s=0;s<numShapes;s++)
141 {
142 n=elements->Nodes[INDEX2(s,e, NN)];
143 data_array=getSampleDataRO(data,n);
144 for (q=0;q<numQuad;q++)
145 {
146 #pragma ivdep
147 for (l=0;l<numComps;l++)
148 {
149 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
150 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
151 }
152 }
153 }
154 }
155 #undef DIM
156 }
157 else if (numDim==3)
158 {
159 const dim_t numShapes=4;
160 #define DIM 3
161 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
162 for (e=0;e<elements->numElements;e++)
163 {
164 grad_data_e=getSampleDataRW(grad_data,e);
165 memset(grad_data_e,0, localGradSize);
166 for (s=0;s<numShapes;s++)
167 {
168 n=elements->Nodes[INDEX2(s,e, NN)];
169 data_array=getSampleDataRO(data,n);
170 for (q=0;q<numQuad;q++)
171 {
172 #pragma ivdep
173 for (l=0;l<numComps;l++)
174 {
175 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
176 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
177 grad_data_e[INDEX4(l,2,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,2,q,0,e, numShapesTotal,DIM,numQuad,1)];
178 }
179 }
180 }
181 }
182 #undef DIM
183 }
184 }
185 else if (data_type==DUDLEY_REDUCED_NODES)
186 {
187 if (numDim==1)
188 {
189 const dim_t numShapes=2;
190 #define DIM 1
191 #pragma omp for schedule(static)
192 for (e=0;e<elements->numElements;e++)
193 {
194 grad_data_e=getSampleDataRW(grad_data,e);
195 memset(grad_data_e,0, localGradSize);
196 for (s=0;s<numShapes;s++)
197 {
198 n=elements->Nodes[INDEX2(s,e, NN)];
199 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
200 for (q=0;q<numQuad;q++)
201 {
202 #pragma ivdep
203 for (l=0;l<numComps;l++)
204 {
205 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
206 }
207 }
208 }
209 }
210 #undef DIM
211 }
212 else if (numDim==2)
213 {
214 const dim_t numShapes=3;
215 #define DIM 2
216 #pragma omp for schedule(static)
217 for (e=0;e<elements->numElements;e++)
218 {
219 grad_data_e=getSampleDataRW(grad_data,e);
220 memset(grad_data_e,0, localGradSize);
221 for (s=0;s<numShapes;s++)
222 {
223 n=elements->Nodes[INDEX2(s,e, NN)];
224 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
225 for (q=0;q<numQuad;q++)
226 {
227 #pragma ivdep
228 for (l=0;l<numComps;l++)
229 {
230 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
231 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
232 }
233 }
234 }
235 }
236 #undef DIM
237 }
238 else if (numDim==3)
239 {
240 const dim_t numShapes=4;
241 #define DIM 3
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)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
257 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
258 grad_data_e[INDEX4(l,2,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,2,q,0,e, numShapesTotal,DIM,numQuad,1)];
259 }
260 }
261 }
262 }
263 #undef DIM
264 }
265 }
266 else if (data_type==DUDLEY_DEGREES_OF_FREEDOM)
267 {
268 if (numDim==1)
269 {
270 const dim_t numShapes=2;
271 #define DIM 1
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->degreesOfFreedomMapping->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)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
287 }
288 }
289 }
290 }
291 #undef DIM
292 }
293 else if (numDim==2)
294 {
295 const dim_t numShapes=3;
296 #define DIM 2
297 #pragma omp for schedule(static)
298 for (e=0;e<elements->numElements;e++)
299 {
300 grad_data_e=getSampleDataRW(grad_data,e);
301 memset(grad_data_e,0, localGradSize);
302 for (s=0;s<numShapes;s++)
303 {
304 n=elements->Nodes[INDEX2(s,e, NN)];
305 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
306 for (q=0;q<numQuad;q++)
307 {
308 #pragma ivdep
309 for (l=0;l<numComps;l++)
310 {
311 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
312 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
313 }
314 }
315 }
316 }
317 #undef DIM
318 }
319 else if (numDim==3)
320 {
321 const dim_t numShapes=4;
322 #define DIM 3
323 #pragma omp for schedule(static)
324 for (e=0;e<elements->numElements;e++)
325 {
326 grad_data_e=getSampleDataRW(grad_data,e);
327 memset(grad_data_e,0, localGradSize);
328 for (s=0;s<numShapes;s++)
329 {
330 n=elements->Nodes[INDEX2(s,e, NN)];
331 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
332 for (q=0;q<numQuad;q++)
333 {
334 #pragma ivdep
335 for (l=0;l<numComps;l++)
336 {
337 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
338 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
339 grad_data_e[INDEX4(l,2,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,2,q,0,e, numShapesTotal,DIM,numQuad,1)];
340 }
341 }
342 }
343 }
344 #undef DIM
345 }
346 }
347 else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM)
348 {
349 if (numDim==1)
350 {
351 const dim_t numShapes=2;
352 #define DIM 1
353 #pragma omp for schedule(static)
354 for (e=0;e<elements->numElements;e++)
355 {
356 grad_data_e=getSampleDataRW(grad_data,e);
357 memset(grad_data_e,0, localGradSize);
358 for (s=0;s<numShapes;s++)
359 {
360 n=elements->Nodes[INDEX2(s,e, NN)];
361 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
362 for (q=0;q<numQuad;q++)
363 {
364 #pragma ivdep
365 for (l=0;l<numComps;l++)
366 {
367 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
368 }
369 }
370 }
371 }
372 #undef DIM
373 }
374 else if (numDim==2)
375 {
376 const dim_t numShapes=3;
377 #define DIM 2
378 #pragma omp for schedule(static)
379 for (e=0;e<elements->numElements;e++)
380 {
381 grad_data_e=getSampleDataRW(grad_data,e);
382 memset(grad_data_e,0, localGradSize);
383 for (s=0;s<numShapes;s++)
384 {
385 n=elements->Nodes[INDEX2(s,e, NN)];
386 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
387 for (q=0;q<numQuad;q++)
388 {
389 #pragma ivdep
390 for (l=0;l<numComps;l++)
391 {
392 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
393 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
394 }
395 }
396 }
397 }
398 #undef DIM
399 }
400 else if (numDim==3)
401 {
402 const dim_t numShapes=4;
403 #define DIM 3
404 #pragma omp for schedule(static)
405 for (e=0;e<elements->numElements;e++)
406 {
407 grad_data_e=getSampleDataRW(grad_data,e);
408 memset(grad_data_e,0, localGradSize);
409 for (s=0;s<numShapes;s++)
410 {
411 n=elements->Nodes[INDEX2(s,e, NN)];
412 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
413 for (q=0;q<numQuad;q++)
414 {
415 #pragma ivdep
416 for (l=0;l<numComps;l++)
417 {
418 grad_data_e[INDEX4(l,0,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
419 grad_data_e[INDEX4(l,1,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,1,q,0,e, numShapesTotal,DIM,numQuad,1)];
420 grad_data_e[INDEX4(l,2,q,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,2,q,0,e, numShapesTotal,DIM,numQuad,1)];
421 }
422 }
423 }
424 }
425 #undef DIM
426 }
427 }
428 } /* end parallel region */
429 }
430 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26