/[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 3206 - (show annotations)
Fri Sep 24 03:20:22 2010 UTC (9 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 13920 byte(s)
Remove references to reducedShapeFunction

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26