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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2770 - (show annotations)
Wed Nov 25 01:24:51 2009 UTC (10 years ago) by jfenwick
File MIME type: text/plain
File size: 17984 byte(s)
Removed buffer implementation of Lazy.
Removed supporting Alloc/Free Sample buffer calls.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,
30 escriptDataC* grad_data,escriptDataC* data) {
31
32 Finley_ReferenceElement* refElement=NULL;
33 size_t localGradSize=0;
34 register dim_t e,q,l,s,n;
35 register __const double *data_array;
36 register double *grad_data_e;
37 dim_t numNodes=0, numShapes=0, numShapesTotal=0, numComps, NN=0, numDim=0, numShapesTotal2=0, numQuad=0, numSub=0, isub=0;
38 type_t data_type=getFunctionSpaceType(data);
39 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 index_t s_offset=0, *nodes_selector=NULL;
41 Finley_ElementFile_Jacobeans* jac=NULL;
42 type_t grad_data_type=getFunctionSpaceType(grad_data);
43
44 Finley_resetError();
45 if (nodes==NULL || elements==NULL) return;
46 numComps=getDataPointSize(data);
47 NN=elements->numNodes;
48 reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
49
50 if (data_type==FINLEY_NODES) {
51 reducedShapefunction=FALSE;
52 numNodes=nodes->nodesMapping->numTargets;
53 } else if (data_type==FINLEY_REDUCED_NODES) {
54 reducedShapefunction=TRUE;
55 numNodes=nodes->reducedNodesMapping->numTargets;
56 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
57 if (elements->MPIInfo->size>1) {
58 Finley_setError(TYPE_ERROR,"Finley_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==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
64 if (elements->MPIInfo->size>1) {
65 Finley_setError(TYPE_ERROR,"Finley_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 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
72 }
73
74 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
75 refElement=Finley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, reducedIntegrationOrder);
76
77 if (Finley_noError()) {
78
79 numDim=jac->numDim;
80 numShapes=jac->BasisFunctions->Type->numShapes;
81 numShapesTotal=jac->numShapesTotal;
82 numSub=jac->numSub;
83 numQuad=jac->numQuadTotal/numSub;
84 if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2 || grad_data_type== FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
85 s_offset=jac->offsets[1];
86 s_offset=jac->offsets[1];
87 } else {
88 s_offset=jac->offsets[0];
89 s_offset=jac->offsets[0];
90 }
91 localGradSize=sizeof(double)*numDim*numQuad*numSub*numComps;
92 if ( (data_type==FINLEY_REDUCED_NODES) || (FINLEY_REDUCED_DEGREES_OF_FREEDOM==data_type) ) {
93 nodes_selector=refElement->Type->linearNodes;
94 numShapesTotal2=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
95 } else {
96 nodes_selector=refElement->Type->subElementNodes;
97 numShapesTotal2=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
98 }
99 /* check the dimensions of data */
100
101 if (! numSamplesEqual(grad_data,numQuad*numSub,elements->numElements)) {
102 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");
103 } else if (! numSamplesEqual(data,1,numNodes)) {
104 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");
105 } else if (numDim*numComps!=getDataPointSize(grad_data)) {
106 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");
107 } else if (!isExpanded(grad_data)) {
108 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");
109 } else if (! (s_offset+numShapes <= numShapesTotal)) {
110 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
111 } else if (! (s_offset+numShapes <= numShapesTotal)) {
112 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: offset test failed.");
113 }
114 }
115 /* now we can start */
116
117 if (Finley_noError()) {
118 requireWrite(grad_data);
119 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e, isub)
120 {
121
122 if (data_type==FINLEY_NODES) {
123 if (numDim==1) {
124 #define DIM 1
125 #pragma omp for schedule(static)
126 for (e=0;e<elements->numElements;e++) {
127 grad_data_e=getSampleDataRW(grad_data,e);
128 memset(grad_data_e,0, localGradSize);
129 for (isub=0; isub<numSub; isub++) {
130 for (s=0;s<numShapes;s++) {
131 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
132 data_array=getSampleDataRO(data,n);
133 for (q=0;q<numQuad;q++) {
134 #pragma ivdep
135 for (l=0;l<numComps;l++) {
136 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
137 }
138 }
139 }
140 }
141 }
142 #undef DIM
143 } else if (numDim==2) {
144 #define DIM 2
145 #pragma omp for schedule(static)
146 for (e=0;e<elements->numElements;e++) {
147 grad_data_e=getSampleDataRW(grad_data,e);
148 memset(grad_data_e,0, localGradSize);
149 for (isub=0; isub<numSub; isub++) {
150 for (s=0;s<numShapes;s++) {
151 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
152 data_array=getSampleDataRO(data,n);
153 for (q=0;q<numQuad;q++) {
154 #pragma ivdep
155 for (l=0;l<numComps;l++) {
156 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
157 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
158 /*printf("data_array of l=%d = %e\n",l,data_array[l]); */
159 }
160 }
161 }
162 }
163 }
164 #undef DIM
165 } else if (numDim==3) {
166 #define DIM 3
167 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
168 for (e=0;e<elements->numElements;e++) {
169 grad_data_e=getSampleDataRW(grad_data,e);
170 memset(grad_data_e,0, localGradSize);
171 for (isub=0; isub<numSub; isub++) {
172 for (s=0;s<numShapes;s++) {
173 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
174 data_array=getSampleDataRO(data,n);
175 for (q=0;q<numQuad;q++) {
176 #pragma ivdep
177 for (l=0;l<numComps;l++) {
178 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
179 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
180 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
181 }
182 }
183 }
184 }
185 }
186 #undef DIM
187 }
188 } else if (data_type==FINLEY_REDUCED_NODES) {
189 if (numDim==1) {
190 #define DIM 1
191 #pragma omp for schedule(static)
192 for (e=0;e<elements->numElements;e++) {
193 grad_data_e=getSampleDataRW(grad_data,e);
194 memset(grad_data_e,0, localGradSize);
195 for (isub=0; isub<numSub; isub++) {
196 for (s=0;s<numShapes;s++) {
197 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
198 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
199 for (q=0;q<numQuad;q++) {
200 #pragma ivdep
201 for (l=0;l<numComps;l++) {
202 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
203 }
204 }
205 }
206 }
207 }
208 #undef DIM
209 } else if (numDim==2) {
210 #define DIM 2
211 #pragma omp for schedule(static)
212 for (e=0;e<elements->numElements;e++) {
213 grad_data_e=getSampleDataRW(grad_data,e);
214 memset(grad_data_e,0, localGradSize);
215 for (isub=0; isub<numSub; isub++) {
216 for (s=0;s<numShapes;s++) {
217 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
218 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
219 for (q=0;q<numQuad;q++) {
220 #pragma ivdep
221 for (l=0;l<numComps;l++) {
222 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
223 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
224 }
225 }
226 }
227 }
228 }
229 #undef DIM
230 } else if (numDim==3) {
231 #define DIM 3
232 #pragma omp for schedule(static)
233 for (e=0;e<elements->numElements;e++) {
234 grad_data_e=getSampleDataRW(grad_data,e);
235 memset(grad_data_e,0, localGradSize);
236 for (isub=0; isub<numSub; isub++) {
237 for (s=0;s<numShapes;s++) {
238 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
239 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
240 for (q=0;q<numQuad;q++) {
241 #pragma ivdep
242 for (l=0;l<numComps;l++) {
243 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
244 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
245 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
246 }
247 }
248 }
249 }
250 }
251 #undef DIM
252 }
253 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
254
255 if (numDim==1) {
256 #define DIM 1
257 #pragma omp for schedule(static)
258 for (e=0;e<elements->numElements;e++) {
259 grad_data_e=getSampleDataRW(grad_data,e);
260 memset(grad_data_e,0, localGradSize);
261 for (isub=0; isub<numSub; isub++) {
262 for (s=0;s<numShapes;s++) {
263 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
264 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
265 for (q=0;q<numQuad;q++) {
266 #pragma ivdep
267 for (l=0;l<numComps;l++) {
268 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
269 }
270 }
271 }
272 }
273 }
274 #undef DIM
275 } else if (numDim==2) {
276 #define DIM 2
277 #pragma omp for schedule(static)
278 for (e=0;e<elements->numElements;e++) {
279 grad_data_e=getSampleDataRW(grad_data,e);
280 memset(grad_data_e,0, localGradSize);
281 for (isub=0; isub<numSub; isub++) {
282 for (s=0;s<numShapes;s++) {
283 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
284 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
285 for (q=0;q<numQuad;q++) {
286 #pragma ivdep
287 for (l=0;l<numComps;l++) {
288 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
289 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
290 }
291 }
292 }
293 }
294 }
295 #undef DIM
296 } else if (numDim==3) {
297 #define DIM 3
298 #pragma omp for schedule(static)
299 for (e=0;e<elements->numElements;e++) {
300 grad_data_e=getSampleDataRW(grad_data,e);
301 memset(grad_data_e,0, localGradSize);
302 for (isub=0; isub<numSub; isub++) {
303 for (s=0;s<numShapes;s++) {
304 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
305 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
306 for (q=0;q<numQuad;q++) {
307 #pragma ivdep
308 for (l=0;l<numComps;l++) {
309 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
310 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
311 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
312 }
313 }
314 }
315 }
316 }
317 #undef DIM
318 }
319 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
320 if (numDim==1) {
321 #define DIM 1
322 #pragma omp for schedule(static)
323 for (e=0;e<elements->numElements;e++) {
324 grad_data_e=getSampleDataRW(grad_data,e);
325 memset(grad_data_e,0, localGradSize);
326 for (isub=0; isub<numSub; isub++) {
327 for (s=0;s<numShapes;s++) {
328 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
329 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
330 for (q=0;q<numQuad;q++) {
331 #pragma ivdep
332 for (l=0;l<numComps;l++) {
333 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
334 }
335 }
336 }
337 }
338 }
339 #undef DIM
340 } else if (numDim==2) {
341 #define DIM 2
342 #pragma omp for schedule(static)
343 for (e=0;e<elements->numElements;e++) {
344 grad_data_e=getSampleDataRW(grad_data,e);
345 memset(grad_data_e,0, localGradSize);
346 for (isub=0; isub<numSub; isub++) {
347 for (s=0;s<numShapes;s++) {
348 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
349 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
350 for (q=0;q<numQuad;q++) {
351 #pragma ivdep
352 for (l=0;l<numComps;l++) {
353 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
354 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
355 }
356 }
357 }
358 }
359 }
360 #undef DIM
361
362 } else if (numDim==3) {
363 #define DIM 3
364 #pragma omp for schedule(static)
365 for (e=0;e<elements->numElements;e++) {
366 grad_data_e=getSampleDataRW(grad_data,e);
367 memset(grad_data_e,0, localGradSize);
368 for (isub=0; isub<numSub; isub++) {
369 for (s=0;s<numShapes;s++) {
370 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
371 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
372 for (q=0;q<numQuad;q++) {
373 #pragma ivdep
374 for (l=0;l<numComps;l++) {
375 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
376 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
377 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
378 }
379 }
380 }
381 }
382 }
383 #undef DIM
384 }
385 }
386 } /* end parallel region */
387 }
388 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26