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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26