/[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 1811 - (show 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
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 double* data_array, *grad_data_e;
34 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
35 type_t data_type=getFunctionSpaceType(data);
36 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
37 index_t dof_offset, s_offset;
38 Finley_ElementFile_Jacobeans* jac=NULL;
39 type_t grad_data_type=getFunctionSpaceType(grad_data);
40 Finley_resetError();
41 if (nodes==NULL || elements==NULL) return;
42 numComps=getDataPointSize(data);
43 NN=elements->numNodes;
44 reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
45
46 if (data_type==FINLEY_NODES) {
47 reducedShapefunction=FALSE;
48 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 reducedShapefunction=FALSE;
58 numNodes=nodes->degreesOfFreedomMapping->numTargets;
59 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
60 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 reducedShapefunction=TRUE;
65 numNodes=nodes->reducedDegreesOfFreedomMapping->numTargets;
66 } else {
67 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
68 }
69
70 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
71 if (Finley_noError()) {
72
73 numShapes=jac->ReferenceElement->Type->numShapes;
74 numLocalNodes=jac->ReferenceElement->Type->numNodes;
75 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 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 }
99 /* now we can start */
100 if (Finley_noError()) {
101 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
102 {
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
122 #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 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
146 for (e=0;e<elements->numElements;e++) {
147 grad_data_e=getSampleData(grad_data,e);
148 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 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
155 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
156 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
157 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
158 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
159 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
167 } 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 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
233
234 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 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
243 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 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
262 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 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
282 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 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
306 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 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
325 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 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
347 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 }
364 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26