/[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 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (10 years ago) by gross
File MIME type: text/plain
File size: 18149 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
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 void* buffer=allocSampleBuffer(data);
119 requireWrite(grad_data);
120 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e, isub)
121 {
122
123 if (data_type==FINLEY_NODES) {
124 if (numDim==1) {
125 #define DIM 1
126 #pragma omp for schedule(static)
127 for (e=0;e<elements->numElements;e++) {
128 grad_data_e=getSampleDataRW(grad_data,e);
129 memset(grad_data_e,0, localGradSize);
130 for (isub=0; isub<numSub; isub++) {
131 for (s=0;s<numShapes;s++) {
132 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
133 data_array=getSampleDataRO(data,n,buffer);
134 for (q=0;q<numQuad;q++) {
135 #pragma ivdep
136 for (l=0;l<numComps;l++) {
137 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)];
138 }
139 }
140 }
141 }
142 }
143 #undef DIM
144 } else if (numDim==2) {
145 #define DIM 2
146 #pragma omp for schedule(static)
147 for (e=0;e<elements->numElements;e++) {
148 grad_data_e=getSampleDataRW(grad_data,e);
149 memset(grad_data_e,0, localGradSize);
150 for (isub=0; isub<numSub; isub++) {
151 for (s=0;s<numShapes;s++) {
152 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
153 data_array=getSampleDataRO(data,n,buffer);
154 for (q=0;q<numQuad;q++) {
155 #pragma ivdep
156 for (l=0;l<numComps;l++) {
157 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)];
158 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)];
159 /*printf("data_array of l=%d = %e\n",l,data_array[l]); */
160 }
161 }
162 }
163 }
164 }
165 #undef DIM
166 } else if (numDim==3) {
167 #define DIM 3
168 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
169 for (e=0;e<elements->numElements;e++) {
170 grad_data_e=getSampleDataRW(grad_data,e);
171 memset(grad_data_e,0, localGradSize);
172 for (isub=0; isub<numSub; isub++) {
173 for (s=0;s<numShapes;s++) {
174 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
175 data_array=getSampleDataRO(data,n,buffer);
176 for (q=0;q<numQuad;q++) {
177 #pragma ivdep
178 for (l=0;l<numComps;l++) {
179 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)];
180 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)];
181 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)];
182 }
183 }
184 }
185 }
186 }
187 #undef DIM
188 }
189 } else if (data_type==FINLEY_REDUCED_NODES) {
190 if (numDim==1) {
191 #define DIM 1
192 #pragma omp for schedule(static)
193 for (e=0;e<elements->numElements;e++) {
194 grad_data_e=getSampleDataRW(grad_data,e);
195 memset(grad_data_e,0, localGradSize);
196 for (isub=0; isub<numSub; isub++) {
197 for (s=0;s<numShapes;s++) {
198 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
199 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n],buffer);
200 for (q=0;q<numQuad;q++) {
201 #pragma ivdep
202 for (l=0;l<numComps;l++) {
203 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)];
204 }
205 }
206 }
207 }
208 }
209 #undef DIM
210 } else if (numDim==2) {
211 #define DIM 2
212 #pragma omp for schedule(static)
213 for (e=0;e<elements->numElements;e++) {
214 grad_data_e=getSampleDataRW(grad_data,e);
215 memset(grad_data_e,0, localGradSize);
216 for (isub=0; isub<numSub; isub++) {
217 for (s=0;s<numShapes;s++) {
218 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
219 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n],buffer);
220 for (q=0;q<numQuad;q++) {
221 #pragma ivdep
222 for (l=0;l<numComps;l++) {
223 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)];
224 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)];
225 }
226 }
227 }
228 }
229 }
230 #undef DIM
231 } else if (numDim==3) {
232 #define DIM 3
233 #pragma omp for schedule(static)
234 for (e=0;e<elements->numElements;e++) {
235 grad_data_e=getSampleDataRW(grad_data,e);
236 memset(grad_data_e,0, localGradSize);
237 for (isub=0; isub<numSub; isub++) {
238 for (s=0;s<numShapes;s++) {
239 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
240 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n],buffer);
241 for (q=0;q<numQuad;q++) {
242 #pragma ivdep
243 for (l=0;l<numComps;l++) {
244 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)];
245 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)];
246 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)];
247 }
248 }
249 }
250 }
251 }
252 #undef DIM
253 }
254 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
255
256 if (numDim==1) {
257 #define DIM 1
258 #pragma omp for schedule(static)
259 for (e=0;e<elements->numElements;e++) {
260 grad_data_e=getSampleDataRW(grad_data,e);
261 memset(grad_data_e,0, localGradSize);
262 for (isub=0; isub<numSub; isub++) {
263 for (s=0;s<numShapes;s++) {
264 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
265 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n],buffer);
266 for (q=0;q<numQuad;q++) {
267 #pragma ivdep
268 for (l=0;l<numComps;l++) {
269 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)];
270 }
271 }
272 }
273 }
274 }
275 #undef DIM
276 } else if (numDim==2) {
277 #define DIM 2
278 #pragma omp for schedule(static)
279 for (e=0;e<elements->numElements;e++) {
280 grad_data_e=getSampleDataRW(grad_data,e);
281 memset(grad_data_e,0, localGradSize);
282 for (isub=0; isub<numSub; isub++) {
283 for (s=0;s<numShapes;s++) {
284 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
285 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n],buffer);
286 for (q=0;q<numQuad;q++) {
287 #pragma ivdep
288 for (l=0;l<numComps;l++) {
289 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)];
290 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)];
291 }
292 }
293 }
294 }
295 }
296 #undef DIM
297 } else if (numDim==3) {
298 #define DIM 3
299 #pragma omp for schedule(static)
300 for (e=0;e<elements->numElements;e++) {
301 grad_data_e=getSampleDataRW(grad_data,e);
302 memset(grad_data_e,0, localGradSize);
303 for (isub=0; isub<numSub; isub++) {
304 for (s=0;s<numShapes;s++) {
305 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
306 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n],buffer);
307 for (q=0;q<numQuad;q++) {
308 #pragma ivdep
309 for (l=0;l<numComps;l++) {
310 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)];
311 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)];
312 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)];
313 }
314 }
315 }
316 }
317 }
318 #undef DIM
319 }
320 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
321 if (numDim==1) {
322 #define DIM 1
323 #pragma omp for schedule(static)
324 for (e=0;e<elements->numElements;e++) {
325 grad_data_e=getSampleDataRW(grad_data,e);
326 memset(grad_data_e,0, localGradSize);
327 for (isub=0; isub<numSub; isub++) {
328 for (s=0;s<numShapes;s++) {
329 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
330 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n],buffer);
331 for (q=0;q<numQuad;q++) {
332 #pragma ivdep
333 for (l=0;l<numComps;l++) {
334 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)];
335 }
336 }
337 }
338 }
339 }
340 #undef DIM
341 } else if (numDim==2) {
342 #define DIM 2
343 #pragma omp for schedule(static)
344 for (e=0;e<elements->numElements;e++) {
345 grad_data_e=getSampleDataRW(grad_data,e);
346 memset(grad_data_e,0, localGradSize);
347 for (isub=0; isub<numSub; isub++) {
348 for (s=0;s<numShapes;s++) {
349 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
350 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n],buffer);
351 for (q=0;q<numQuad;q++) {
352 #pragma ivdep
353 for (l=0;l<numComps;l++) {
354 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)];
355 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)];
356 }
357 }
358 }
359 }
360 }
361 #undef DIM
362
363 } else if (numDim==3) {
364 #define DIM 3
365 #pragma omp for schedule(static)
366 for (e=0;e<elements->numElements;e++) {
367 grad_data_e=getSampleDataRW(grad_data,e);
368 memset(grad_data_e,0, localGradSize);
369 for (isub=0; isub<numSub; isub++) {
370 for (s=0;s<numShapes;s++) {
371 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
372 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n],buffer);
373 for (q=0;q<numQuad;q++) {
374 #pragma ivdep
375 for (l=0;l<numComps;l++) {
376 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)];
377 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)];
378 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)];
379 }
380 }
381 }
382 }
383 }
384 #undef DIM
385 }
386 }
387 } /* end parallel region */
388 freeSampleBuffer(buffer);
389 }
390 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26