/[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 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 8 months ago) by trankine
File MIME type: text/plain
File size: 20894 byte(s)
And get the *(&(*&(* name right
1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* assemblage jacobiens: calculate the gradient of nodal data at quadrature points */
19
20 /**************************************************************/
21
22 #include "Assemble.h"
23 #include "Util.h"
24 #ifdef _OPENMP
25 #include <omp.h>
26 #endif
27 /*****************************************************************/
28
29
30 void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,
31 escriptDataC* grad_data,escriptDataC* data) {
32
33 register dim_t e,q,l,s,n;
34 register double* data_array, *grad_data_e;
35 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
36 type_t data_type=getFunctionSpaceType(data);
37 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
38 index_t dof_offset, s_offset;
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 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
103 {
104 if (data_type==FINLEY_NODES) {
105 if (jac->numDim==1) {
106 #define DIM 1
107 #pragma omp for schedule(static)
108 for (e=0;e<elements->numElements;e++) {
109 grad_data_e=getSampleData(grad_data,e);
110 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
111 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
112 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
113 data_array=getSampleData(data,n);
114 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
115 for (l=0;l<numComps;l++) {
116 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
117 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
118 }
119 }
120 }
121 }
122
123 #undef DIM
124 } else if (jac->numDim==2) {
125 #define DIM 2
126 #pragma omp for schedule(static)
127 for (e=0;e<elements->numElements;e++) {
128 grad_data_e=getSampleData(grad_data,e);
129 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
130 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
131 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
132 data_array=getSampleData(data,n);
133 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
134 for (l=0;l<numComps;l++) {
135 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
136 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
137 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
138 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
139 }
140 }
141 }
142 }
143 #undef DIM
144 } else if (jac->numDim==3) {
145 #define DIM 3
146 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
147 for (e=0;e<elements->numElements;e++) {
148 grad_data_e=getSampleData(grad_data,e);
149 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
150 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
151 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
152 data_array=getSampleData(data,n);
153 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
154 for (l=0;l<numComps;l++) {
155 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
156 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
157 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
158 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
159 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
160 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
161 }
162 }
163 }
164 }
165 #undef DIM
166 }
167
168 } else if (data_type==FINLEY_REDUCED_NODES) {
169 if (jac->numDim==1) {
170 #define DIM 1
171 #pragma omp for schedule(static)
172 for (e=0;e<elements->numElements;e++) {
173 grad_data_e=getSampleData(grad_data,e);
174 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
175 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
176 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
177 data_array=getSampleData(data,nodes->reducedNodesMapping->target[n]);
178 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
179 for (l=0;l<numComps;l++) {
180 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
181 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
182 }
183 }
184 }
185 }
186
187 #undef DIM
188 } else if (jac->numDim==2) {
189 #define DIM 2
190 #pragma omp for schedule(static)
191 for (e=0;e<elements->numElements;e++) {
192 grad_data_e=getSampleData(grad_data,e);
193 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
194 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
195 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
196 data_array=getSampleData(data,nodes->reducedNodesMapping->target[n]);
197 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
198 for (l=0;l<numComps;l++) {
199 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
200 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
201 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
202 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
203 }
204 }
205 }
206 }
207
208 #undef DIM
209
210 } else if (jac->numDim==3) {
211 #define DIM 3
212 #pragma omp for schedule(static)
213 for (e=0;e<elements->numElements;e++) {
214 grad_data_e=getSampleData(grad_data,e);
215 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
216 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
217 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
218 data_array=getSampleData(data,nodes->reducedNodesMapping->target[n]);
219 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
220 for (l=0;l<numComps;l++) {
221 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
222 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
223 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
224 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
225 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
226 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
227 }
228 }
229 }
230 }
231 #undef DIM
232 }
233 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
234
235 if (jac->numDim==1) {
236 #define DIM 1
237 #pragma omp for schedule(static)
238 for (e=0;e<elements->numElements;e++) {
239 grad_data_e=getSampleData(grad_data,e);
240 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
241 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
242 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
243 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
244 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
245 for (l=0;l<numComps;l++) {
246 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
247 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
248 }
249 }
250 }
251 }
252
253 #undef DIM
254 } else if (jac->numDim==2) {
255 #define DIM 2
256 #pragma omp for schedule(static)
257 for (e=0;e<elements->numElements;e++) {
258 grad_data_e=getSampleData(grad_data,e);
259 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
260 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
261 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
262 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
263 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
264 for (l=0;l<numComps;l++) {
265 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
266 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
267 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
268 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
269 }
270 }
271 }
272 }
273 #undef DIM
274 } else if (jac->numDim==3) {
275 #define DIM 3
276 #pragma omp for schedule(static)
277 for (e=0;e<elements->numElements;e++) {
278 grad_data_e=getSampleData(grad_data,e);
279 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
280 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
281 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
282 data_array=getSampleData(data,nodes->degreesOfFreedomMapping->target[n]);
283 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
284 for (l=0;l<numComps;l++) {
285 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
286 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
287 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
288 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
289 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
290 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
291 }
292 }
293 }
294 }
295 #undef DIM
296 }
297 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
298 if (jac->numDim==1) {
299 #define DIM 1
300 #pragma omp for schedule(static)
301 for (e=0;e<elements->numElements;e++) {
302 grad_data_e=getSampleData(grad_data,e);
303 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
304 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
305 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
306 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
307 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
308 for (l=0;l<numComps;l++) {
309 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
310 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
311 }
312 }
313 }
314 }
315
316 #undef DIM
317 } else if (jac->numDim==2) {
318 #define DIM 2
319 #pragma omp for schedule(static)
320 for (e=0;e<elements->numElements;e++) {
321 grad_data_e=getSampleData(grad_data,e);
322 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
323 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
324 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
325 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
326 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
327 for (l=0;l<numComps;l++) {
328 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
329 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
330 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
331 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
332 }
333 }
334 }
335 }
336
337 #undef DIM
338
339 } else if (jac->numDim==3) {
340 #define DIM 3
341 #pragma omp for schedule(static)
342 for (e=0;e<elements->numElements;e++) {
343 grad_data_e=getSampleData(grad_data,e);
344 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
345 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
346 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
347 data_array=getSampleData(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
348 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
349 for (l=0;l<numComps;l++) {
350 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
351 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
352 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
353 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
354 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
355 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
356 }
357 }
358 }
359 }
360 #undef DIM
361 }
362 }
363 } /* end parallel region */
364 }
365 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26