/[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 3152 - (show annotations)
Fri Sep 3 05:48:31 2010 UTC (9 years, 2 months ago) by jfenwick
File MIME type: text/plain
File size: 15355 byte(s)
Collapsed duplicate fields
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 numShapesTotal2=refElement->BasisFunctions->Type->numShapes;
83 /* check the dimensions of data */
84
85 if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) {
86 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
87 } else if (! numSamplesEqual(data,1,numNodes)) {
88 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples of input Data object");
89 } else if (numDim*numComps!=getDataPointSize(grad_data)) {
90 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of components in gradient data object.");
91 } else if (!isExpanded(grad_data)) {
92 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: expanded Data object is expected for output data.");
93 } else if (! (numShapes <= numShapesTotal)) {
94 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
95 } else if (! (numShapes <= numShapesTotal)) {
96 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: offset test failed.");
97 }
98 }
99 /* now we can start */
100
101 if (Dudley_noError()) {
102 requireWrite(grad_data);
103 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
104 {
105
106 if (data_type==DUDLEY_NODES) {
107 if (numDim==1) {
108 #define DIM 1
109 #pragma omp for schedule(static)
110 for (e=0;e<elements->numElements;e++) {
111 grad_data_e=getSampleDataRW(grad_data,e);
112 memset(grad_data_e,0, localGradSize);
113 for (s=0;s<numShapes;s++) {
114 n=elements->Nodes[INDEX2(s,e, NN)];
115 data_array=getSampleDataRO(data,n);
116 for (q=0;q<numQuad;q++) {
117 #pragma ivdep
118 for (l=0;l<numComps;l++) {
119 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)];
120 }
121 }
122 }
123 }
124 #undef DIM
125 } else if (numDim==2) {
126 #define DIM 2
127 #pragma omp for schedule(static)
128 for (e=0;e<elements->numElements;e++) {
129 grad_data_e=getSampleDataRW(grad_data,e);
130 memset(grad_data_e,0, localGradSize);
131 for (s=0;s<numShapes;s++) {
132 n=elements->Nodes[INDEX2(s,e, NN)];
133 data_array=getSampleDataRO(data,n);
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,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
138 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)];
139 /*printf("data_array of l=%d = %e\n",l,data_array[l]); */
140 }
141 }
142 }
143 }
144 #undef DIM
145 } else if (numDim==3) {
146 #define DIM 3
147 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
148 for (e=0;e<elements->numElements;e++) {
149 grad_data_e=getSampleDataRW(grad_data,e);
150 memset(grad_data_e,0, localGradSize);
151 for (s=0;s<numShapes;s++) {
152 n=elements->Nodes[INDEX2(s,e, NN)];
153 data_array=getSampleDataRO(data,n);
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,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
158 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)];
159 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)];
160 }
161 }
162 }
163 }
164 #undef DIM
165 }
166 } else if (data_type==DUDLEY_REDUCED_NODES) {
167 if (numDim==1) {
168 #define DIM 1
169 #pragma omp for schedule(static)
170 for (e=0;e<elements->numElements;e++) {
171 grad_data_e=getSampleDataRW(grad_data,e);
172 memset(grad_data_e,0, localGradSize);
173 for (s=0;s<numShapes;s++) {
174 n=elements->Nodes[INDEX2(s,e, NN)];
175 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
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,0, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s,0,q,0,e, numShapesTotal,DIM,numQuad,1)];
180 }
181 }
182 }
183 }
184 #undef DIM
185 } else if (numDim==2) {
186 #define DIM 2
187 #pragma omp for schedule(static)
188 for (e=0;e<elements->numElements;e++) {
189 grad_data_e=getSampleDataRW(grad_data,e);
190 memset(grad_data_e,0, localGradSize);
191 for (s=0;s<numShapes;s++) {
192 n=elements->Nodes[INDEX2(s,e, NN)];
193 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
194 for (q=0;q<numQuad;q++) {
195 #pragma ivdep
196 for (l=0;l<numComps;l++) {
197 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)];
198 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)];
199 }
200 }
201 }
202 }
203 #undef DIM
204 } else if (numDim==3) {
205 #define DIM 3
206 #pragma omp for schedule(static)
207 for (e=0;e<elements->numElements;e++) {
208 grad_data_e=getSampleDataRW(grad_data,e);
209 memset(grad_data_e,0, localGradSize);
210 for (s=0;s<numShapes;s++) {
211 n=elements->Nodes[INDEX2(s,e, NN)];
212 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
213 for (q=0;q<numQuad;q++) {
214 #pragma ivdep
215 for (l=0;l<numComps;l++) {
216 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)];
217 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)];
218 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)];
219 }
220 }
221 }
222 }
223 #undef DIM
224 }
225 } else if (data_type==DUDLEY_DEGREES_OF_FREEDOM) {
226
227 if (numDim==1) {
228 #define DIM 1
229 #pragma omp for schedule(static)
230 for (e=0;e<elements->numElements;e++) {
231 grad_data_e=getSampleDataRW(grad_data,e);
232 memset(grad_data_e,0, localGradSize);
233 for (s=0;s<numShapes;s++) {
234 n=elements->Nodes[INDEX2(s,e, NN)];
235 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
236 for (q=0;q<numQuad;q++) {
237 #pragma ivdep
238 for (l=0;l<numComps;l++) {
239 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)];
240 }
241 }
242 }
243 }
244 #undef DIM
245 } else if (numDim==2) {
246 #define DIM 2
247 #pragma omp for schedule(static)
248 for (e=0;e<elements->numElements;e++) {
249 grad_data_e=getSampleDataRW(grad_data,e);
250 memset(grad_data_e,0, localGradSize);
251 for (s=0;s<numShapes;s++) {
252 n=elements->Nodes[INDEX2(s,e, NN)];
253 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
254 for (q=0;q<numQuad;q++) {
255 #pragma ivdep
256 for (l=0;l<numComps;l++) {
257 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)];
258 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)];
259 }
260 }
261 }
262 }
263 #undef DIM
264 } else if (numDim==3) {
265 #define DIM 3
266 #pragma omp for schedule(static)
267 for (e=0;e<elements->numElements;e++) {
268 grad_data_e=getSampleDataRW(grad_data,e);
269 memset(grad_data_e,0, localGradSize);
270 for (s=0;s<numShapes;s++) {
271 n=elements->Nodes[INDEX2(s,e, NN)];
272 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
273 for (q=0;q<numQuad;q++) {
274 #pragma ivdep
275 for (l=0;l<numComps;l++) {
276 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)];
277 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)];
278 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)];
279 }
280 }
281 }
282 }
283 #undef DIM
284 }
285 } else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
286 if (numDim==1) {
287 #define DIM 1
288 #pragma omp for schedule(static)
289 for (e=0;e<elements->numElements;e++) {
290 grad_data_e=getSampleDataRW(grad_data,e);
291 memset(grad_data_e,0, localGradSize);
292 for (s=0;s<numShapes;s++) {
293 n=elements->Nodes[INDEX2(s,e, NN)];
294 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
295 for (q=0;q<numQuad;q++) {
296 #pragma ivdep
297 for (l=0;l<numComps;l++) {
298 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)];
299 }
300 }
301 }
302 }
303 #undef DIM
304 } else if (numDim==2) {
305 #define DIM 2
306 #pragma omp for schedule(static)
307 for (e=0;e<elements->numElements;e++) {
308 grad_data_e=getSampleDataRW(grad_data,e);
309 memset(grad_data_e,0, localGradSize);
310 for (s=0;s<numShapes;s++) {
311 n=elements->Nodes[INDEX2(s,e, NN)];
312 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
313 for (q=0;q<numQuad;q++) {
314 #pragma ivdep
315 for (l=0;l<numComps;l++) {
316 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)];
317 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)];
318 }
319 }
320 }
321 }
322 #undef DIM
323
324 } else if (numDim==3) {
325 #define DIM 3
326 #pragma omp for schedule(static)
327 for (e=0;e<elements->numElements;e++) {
328 grad_data_e=getSampleDataRW(grad_data,e);
329 memset(grad_data_e,0, localGradSize);
330 for (s=0;s<numShapes;s++) {
331 n=elements->Nodes[INDEX2(s,e, NN)];
332 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
333 for (q=0;q<numQuad;q++) {
334 #pragma ivdep
335 for (l=0;l<numComps;l++) {
336 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)];
337 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)];
338 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)];
339 }
340 }
341 }
342 }
343 #undef DIM
344 }
345 }
346 } /* end parallel region */
347 }
348 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26