/[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 2271 - (show annotations)
Mon Feb 16 05:08:29 2009 UTC (10 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 21135 byte(s)
Merging version 2269 to trunk

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26