/[escript]/trunk/finley/src/Assemble_gradient.c
ViewVC logotype

Annotation of /trunk/finley/src/Assemble_gradient.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 20859 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26