/[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 3149 - (show annotations)
Fri Sep 3 03:49:32 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 15546 byte(s)


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 void Dudley_Assemble_gradient(Dudley_NodeFile* nodes, Dudley_ElementFile* elements,
30 escriptDataC* grad_data,escriptDataC* data) {
31
32 Dudley_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;
38 type_t data_type=getFunctionSpaceType(data);
39 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 Dudley_ElementFile_Jacobeans* jac=NULL;
41
42 Dudley_resetError();
43 if (nodes==NULL || elements==NULL) return;
44 numComps=getDataPointSize(data);
45 NN=elements->numNodes;
46 reducedIntegrationOrder=Dudley_Assemble_reducedIntegrationOrder(grad_data);
47
48 if (data_type==DUDLEY_NODES) {
49 reducedShapefunction=FALSE;
50 numNodes=nodes->nodesMapping->numTargets;
51 } else if (data_type==DUDLEY_REDUCED_NODES) {
52 reducedShapefunction=TRUE;
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 reducedShapefunction=FALSE;
60 numNodes=nodes->degreesOfFreedomMapping->numTargets;
61 } else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
62 if (elements->MPIInfo->size>1) {
63 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
64 return;
65 }
66 reducedShapefunction=TRUE;
67 numNodes=nodes->reducedDegreesOfFreedomMapping->numTargets;
68 } else {
69 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
70 }
71
72 jac=Dudley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
73 refElement=Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, reducedIntegrationOrder);
74
75 if (Dudley_noError()) {
76
77 numDim=jac->numDim;
78 numShapes=jac->BasisFunctions->Type->numShapes;
79 numShapesTotal=jac->numShapesTotal;
80 numQuad=jac->numQuadTotal;
81 localGradSize=sizeof(double)*numDim*numQuad*numComps;
82 if ( (data_type==DUDLEY_REDUCED_NODES) || (DUDLEY_REDUCED_DEGREES_OF_FREEDOM==data_type) ) {
83 numShapesTotal2=refElement->LinearBasisFunctions->Type->numShapes;
84 } else {
85 numShapesTotal2=refElement->BasisFunctions->Type->numShapes;
86 }
87 /* check the dimensions of data */
88
89 if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) {
90 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
91 } else if (! numSamplesEqual(data,1,numNodes)) {
92 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples of input Data object");
93 } else if (numDim*numComps!=getDataPointSize(grad_data)) {
94 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of components in gradient data object.");
95 } else if (!isExpanded(grad_data)) {
96 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: expanded Data object is expected for output data.");
97 } else if (! (numShapes <= numShapesTotal)) {
98 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
99 } else if (! (numShapes <= numShapesTotal)) {
100 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: offset test failed.");
101 }
102 }
103 /* now we can start */
104
105 if (Dudley_noError()) {
106 requireWrite(grad_data);
107 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
108 {
109
110 if (data_type==DUDLEY_NODES) {
111 if (numDim==1) {
112 #define DIM 1
113 #pragma omp for schedule(static)
114 for (e=0;e<elements->numElements;e++) {
115 grad_data_e=getSampleDataRW(grad_data,e);
116 memset(grad_data_e,0, localGradSize);
117 for (s=0;s<numShapes;s++) {
118 n=elements->Nodes[INDEX2(s,e, NN)];
119 data_array=getSampleDataRO(data,n);
120 for (q=0;q<numQuad;q++) {
121 #pragma ivdep
122 for (l=0;l<numComps;l++) {
123 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)];
124 }
125 }
126 }
127 }
128 #undef DIM
129 } else if (numDim==2) {
130 #define DIM 2
131 #pragma omp for schedule(static)
132 for (e=0;e<elements->numElements;e++) {
133 grad_data_e=getSampleDataRW(grad_data,e);
134 memset(grad_data_e,0, localGradSize);
135 for (s=0;s<numShapes;s++) {
136 n=elements->Nodes[INDEX2(s,e, NN)];
137 data_array=getSampleDataRO(data,n);
138 for (q=0;q<numQuad;q++) {
139 #pragma ivdep
140 for (l=0;l<numComps;l++) {
141 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)];
142 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)];
143 /*printf("data_array of l=%d = %e\n",l,data_array[l]); */
144 }
145 }
146 }
147 }
148 #undef DIM
149 } else if (numDim==3) {
150 #define DIM 3
151 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
152 for (e=0;e<elements->numElements;e++) {
153 grad_data_e=getSampleDataRW(grad_data,e);
154 memset(grad_data_e,0, localGradSize);
155 for (s=0;s<numShapes;s++) {
156 n=elements->Nodes[INDEX2(s,e, NN)];
157 data_array=getSampleDataRO(data,n);
158 for (q=0;q<numQuad;q++) {
159 #pragma ivdep
160 for (l=0;l<numComps;l++) {
161 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)];
162 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)];
163 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)];
164 }
165 }
166 }
167 }
168 #undef DIM
169 }
170 } else if (data_type==DUDLEY_REDUCED_NODES) {
171 if (numDim==1) {
172 #define DIM 1
173 #pragma omp for schedule(static)
174 for (e=0;e<elements->numElements;e++) {
175 grad_data_e=getSampleDataRW(grad_data,e);
176 memset(grad_data_e,0, localGradSize);
177 for (s=0;s<numShapes;s++) {
178 n=elements->Nodes[INDEX2(s,e, NN)];
179 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
180 for (q=0;q<numQuad;q++) {
181 #pragma ivdep
182 for (l=0;l<numComps;l++) {
183 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)];
184 }
185 }
186 }
187 }
188 #undef DIM
189 } else if (numDim==2) {
190 #define DIM 2
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 (s=0;s<numShapes;s++) {
196 n=elements->Nodes[INDEX2(s,e, NN)];
197 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
198 for (q=0;q<numQuad;q++) {
199 #pragma ivdep
200 for (l=0;l<numComps;l++) {
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 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)];
203 }
204 }
205 }
206 }
207 #undef DIM
208 } else if (numDim==3) {
209 #define DIM 3
210 #pragma omp for schedule(static)
211 for (e=0;e<elements->numElements;e++) {
212 grad_data_e=getSampleDataRW(grad_data,e);
213 memset(grad_data_e,0, localGradSize);
214 for (s=0;s<numShapes;s++) {
215 n=elements->Nodes[INDEX2(s,e, NN)];
216 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
217 for (q=0;q<numQuad;q++) {
218 #pragma ivdep
219 for (l=0;l<numComps;l++) {
220 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)];
221 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)];
222 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)];
223 }
224 }
225 }
226 }
227 #undef DIM
228 }
229 } else if (data_type==DUDLEY_DEGREES_OF_FREEDOM) {
230
231 if (numDim==1) {
232 #define DIM 1
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 (s=0;s<numShapes;s++) {
238 n=elements->Nodes[INDEX2(s,e, NN)];
239 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->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,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
244 }
245 }
246 }
247 }
248 #undef DIM
249 } else if (numDim==2) {
250 #define DIM 2
251 #pragma omp for schedule(static)
252 for (e=0;e<elements->numElements;e++) {
253 grad_data_e=getSampleDataRW(grad_data,e);
254 memset(grad_data_e,0, localGradSize);
255 for (s=0;s<numShapes;s++) {
256 n=elements->Nodes[INDEX2(s,e, NN)];
257 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
258 for (q=0;q<numQuad;q++) {
259 #pragma ivdep
260 for (l=0;l<numComps;l++) {
261 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)];
262 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)];
263 }
264 }
265 }
266 }
267 #undef DIM
268 } else if (numDim==3) {
269 #define DIM 3
270 #pragma omp for schedule(static)
271 for (e=0;e<elements->numElements;e++) {
272 grad_data_e=getSampleDataRW(grad_data,e);
273 memset(grad_data_e,0, localGradSize);
274 for (s=0;s<numShapes;s++) {
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 #pragma ivdep
279 for (l=0;l<numComps;l++) {
280 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)];
281 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)];
282 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)];
283 }
284 }
285 }
286 }
287 #undef DIM
288 }
289 } else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
290 if (numDim==1) {
291 #define DIM 1
292 #pragma omp for schedule(static)
293 for (e=0;e<elements->numElements;e++) {
294 grad_data_e=getSampleDataRW(grad_data,e);
295 memset(grad_data_e,0, localGradSize);
296 for (s=0;s<numShapes;s++) {
297 n=elements->Nodes[INDEX2(s,e, NN)];
298 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
299 for (q=0;q<numQuad;q++) {
300 #pragma ivdep
301 for (l=0;l<numComps;l++) {
302 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)];
303 }
304 }
305 }
306 }
307 #undef DIM
308 } else if (numDim==2) {
309 #define DIM 2
310 #pragma omp for schedule(static)
311 for (e=0;e<elements->numElements;e++) {
312 grad_data_e=getSampleDataRW(grad_data,e);
313 memset(grad_data_e,0, localGradSize);
314 for (s=0;s<numShapes;s++) {
315 n=elements->Nodes[INDEX2(s,e, NN)];
316 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
317 for (q=0;q<numQuad;q++) {
318 #pragma ivdep
319 for (l=0;l<numComps;l++) {
320 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)];
321 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)];
322 }
323 }
324 }
325 }
326 #undef DIM
327
328 } else if (numDim==3) {
329 #define DIM 3
330 #pragma omp for schedule(static)
331 for (e=0;e<elements->numElements;e++) {
332 grad_data_e=getSampleDataRW(grad_data,e);
333 memset(grad_data_e,0, localGradSize);
334 for (s=0;s<numShapes;s++) {
335 n=elements->Nodes[INDEX2(s,e, NN)];
336 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
337 for (q=0;q<numQuad;q++) {
338 #pragma ivdep
339 for (l=0;l<numComps;l++) {
340 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)];
341 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)];
342 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)];
343 }
344 }
345 }
346 }
347 #undef DIM
348 }
349 }
350 } /* end parallel region */
351 }
352 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26