/[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 3090 - (show annotations)
Wed Aug 11 00:51:28 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 17696 byte(s)
Removed:
DUDLEY_CONTACT_ELEMENTS_1
DUDLEY_CONTACT_ELEMENTS_2
DUDLEY_REDUCED_CONTACT_ELEMENTS_1
DUDLEY_REDUCED_CONTACT_ELEMENTS_2

escript tests now query if the domain supports contact elements
before trying to use them.


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, numSub=0, isub=0;
38 type_t data_type=getFunctionSpaceType(data);
39 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 index_t s_offset=0, *nodes_selector=NULL;
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 numSub=jac->numSub;
82 numQuad=jac->numQuadTotal/numSub;
83 s_offset=jac->offsets[0];
84 s_offset=jac->offsets[0];
85 localGradSize=sizeof(double)*numDim*numQuad*numSub*numComps;
86 if ( (data_type==DUDLEY_REDUCED_NODES) || (DUDLEY_REDUCED_DEGREES_OF_FREEDOM==data_type) ) {
87 nodes_selector=refElement->Type->linearNodes;
88 numShapesTotal2=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
89 } else {
90 nodes_selector=refElement->Type->subElementNodes;
91 numShapesTotal2=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
92 }
93 /* check the dimensions of data */
94
95 if (! numSamplesEqual(grad_data,numQuad*numSub,elements->numElements)) {
96 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples in gradient Data object");
97 } else if (! numSamplesEqual(data,1,numNodes)) {
98 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of samples of input Data object");
99 } else if (numDim*numComps!=getDataPointSize(grad_data)) {
100 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: illegal number of components in gradient data object.");
101 } else if (!isExpanded(grad_data)) {
102 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_gradient: expanded Data object is expected for output data.");
103 } else if (! (s_offset+numShapes <= numShapesTotal)) {
104 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
105 } else if (! (s_offset+numShapes <= numShapesTotal)) {
106 Dudley_setError(SYSTEM_ERROR,"Dudley_Assemble_gradient: offset test failed.");
107 }
108 }
109 /* now we can start */
110
111 if (Dudley_noError()) {
112 requireWrite(grad_data);
113 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e, isub)
114 {
115
116 if (data_type==DUDLEY_NODES) {
117 if (numDim==1) {
118 #define DIM 1
119 #pragma omp for schedule(static)
120 for (e=0;e<elements->numElements;e++) {
121 grad_data_e=getSampleDataRW(grad_data,e);
122 memset(grad_data_e,0, localGradSize);
123 for (isub=0; isub<numSub; isub++) {
124 for (s=0;s<numShapes;s++) {
125 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
126 data_array=getSampleDataRO(data,n);
127 for (q=0;q<numQuad;q++) {
128 #pragma ivdep
129 for (l=0;l<numComps;l++) {
130 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
131 }
132 }
133 }
134 }
135 }
136 #undef DIM
137 } else if (numDim==2) {
138 #define DIM 2
139 #pragma omp for schedule(static)
140 for (e=0;e<elements->numElements;e++) {
141 grad_data_e=getSampleDataRW(grad_data,e);
142 memset(grad_data_e,0, localGradSize);
143 for (isub=0; isub<numSub; isub++) {
144 for (s=0;s<numShapes;s++) {
145 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
146 data_array=getSampleDataRO(data,n);
147 for (q=0;q<numQuad;q++) {
148 #pragma ivdep
149 for (l=0;l<numComps;l++) {
150 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
151 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
152 /*printf("data_array of l=%d = %e\n",l,data_array[l]); */
153 }
154 }
155 }
156 }
157 }
158 #undef DIM
159 } else if (numDim==3) {
160 #define DIM 3
161 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
162 for (e=0;e<elements->numElements;e++) {
163 grad_data_e=getSampleDataRW(grad_data,e);
164 memset(grad_data_e,0, localGradSize);
165 for (isub=0; isub<numSub; isub++) {
166 for (s=0;s<numShapes;s++) {
167 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
168 data_array=getSampleDataRO(data,n);
169 for (q=0;q<numQuad;q++) {
170 #pragma ivdep
171 for (l=0;l<numComps;l++) {
172 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
173 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
174 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
175 }
176 }
177 }
178 }
179 }
180 #undef DIM
181 }
182 } else if (data_type==DUDLEY_REDUCED_NODES) {
183 if (numDim==1) {
184 #define DIM 1
185 #pragma omp for schedule(static)
186 for (e=0;e<elements->numElements;e++) {
187 grad_data_e=getSampleDataRW(grad_data,e);
188 memset(grad_data_e,0, localGradSize);
189 for (isub=0; isub<numSub; isub++) {
190 for (s=0;s<numShapes;s++) {
191 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
192 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
193 for (q=0;q<numQuad;q++) {
194 #pragma ivdep
195 for (l=0;l<numComps;l++) {
196 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
197 }
198 }
199 }
200 }
201 }
202 #undef DIM
203 } else if (numDim==2) {
204 #define DIM 2
205 #pragma omp for schedule(static)
206 for (e=0;e<elements->numElements;e++) {
207 grad_data_e=getSampleDataRW(grad_data,e);
208 memset(grad_data_e,0, localGradSize);
209 for (isub=0; isub<numSub; isub++) {
210 for (s=0;s<numShapes;s++) {
211 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
212 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
213 for (q=0;q<numQuad;q++) {
214 #pragma ivdep
215 for (l=0;l<numComps;l++) {
216 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
217 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
218 }
219 }
220 }
221 }
222 }
223 #undef DIM
224 } else if (numDim==3) {
225 #define DIM 3
226 #pragma omp for schedule(static)
227 for (e=0;e<elements->numElements;e++) {
228 grad_data_e=getSampleDataRW(grad_data,e);
229 memset(grad_data_e,0, localGradSize);
230 for (isub=0; isub<numSub; isub++) {
231 for (s=0;s<numShapes;s++) {
232 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
233 data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
234 for (q=0;q<numQuad;q++) {
235 #pragma ivdep
236 for (l=0;l<numComps;l++) {
237 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
238 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
239 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
240 }
241 }
242 }
243 }
244 }
245 #undef DIM
246 }
247 } else if (data_type==DUDLEY_DEGREES_OF_FREEDOM) {
248
249 if (numDim==1) {
250 #define DIM 1
251 #pragma omp for schedule(static)
252 for (e=0;e<elements->numElements;e++) {
253 grad_data_e=getSampleDataRW(grad_data,e);
254 memset(grad_data_e,0, localGradSize);
255 for (isub=0; isub<numSub; isub++) {
256 for (s=0;s<numShapes;s++) {
257 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
258 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
259 for (q=0;q<numQuad;q++) {
260 #pragma ivdep
261 for (l=0;l<numComps;l++) {
262 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
263 }
264 }
265 }
266 }
267 }
268 #undef DIM
269 } else if (numDim==2) {
270 #define DIM 2
271 #pragma omp for schedule(static)
272 for (e=0;e<elements->numElements;e++) {
273 grad_data_e=getSampleDataRW(grad_data,e);
274 memset(grad_data_e,0, localGradSize);
275 for (isub=0; isub<numSub; isub++) {
276 for (s=0;s<numShapes;s++) {
277 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
278 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
279 for (q=0;q<numQuad;q++) {
280 #pragma ivdep
281 for (l=0;l<numComps;l++) {
282 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
283 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
284 }
285 }
286 }
287 }
288 }
289 #undef DIM
290 } else if (numDim==3) {
291 #define DIM 3
292 #pragma omp for schedule(static)
293 for (e=0;e<elements->numElements;e++) {
294 grad_data_e=getSampleDataRW(grad_data,e);
295 memset(grad_data_e,0, localGradSize);
296 for (isub=0; isub<numSub; isub++) {
297 for (s=0;s<numShapes;s++) {
298 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
299 data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
300 for (q=0;q<numQuad;q++) {
301 #pragma ivdep
302 for (l=0;l<numComps;l++) {
303 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
304 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
305 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
306 }
307 }
308 }
309 }
310 }
311 #undef DIM
312 }
313 } else if (data_type==DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
314 if (numDim==1) {
315 #define DIM 1
316 #pragma omp for schedule(static)
317 for (e=0;e<elements->numElements;e++) {
318 grad_data_e=getSampleDataRW(grad_data,e);
319 memset(grad_data_e,0, localGradSize);
320 for (isub=0; isub<numSub; isub++) {
321 for (s=0;s<numShapes;s++) {
322 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
323 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
324 for (q=0;q<numQuad;q++) {
325 #pragma ivdep
326 for (l=0;l<numComps;l++) {
327 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
328 }
329 }
330 }
331 }
332 }
333 #undef DIM
334 } else if (numDim==2) {
335 #define DIM 2
336 #pragma omp for schedule(static)
337 for (e=0;e<elements->numElements;e++) {
338 grad_data_e=getSampleDataRW(grad_data,e);
339 memset(grad_data_e,0, localGradSize);
340 for (isub=0; isub<numSub; isub++) {
341 for (s=0;s<numShapes;s++) {
342 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
343 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
344 for (q=0;q<numQuad;q++) {
345 #pragma ivdep
346 for (l=0;l<numComps;l++) {
347 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
348 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
349 }
350 }
351 }
352 }
353 }
354 #undef DIM
355
356 } else if (numDim==3) {
357 #define DIM 3
358 #pragma omp for schedule(static)
359 for (e=0;e<elements->numElements;e++) {
360 grad_data_e=getSampleDataRW(grad_data,e);
361 memset(grad_data_e,0, localGradSize);
362 for (isub=0; isub<numSub; isub++) {
363 for (s=0;s<numShapes;s++) {
364 n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
365 data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
366 for (q=0;q<numQuad;q++) {
367 #pragma ivdep
368 for (l=0;l<numComps;l++) {
369 grad_data_e[INDEX4(l,0,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
370 grad_data_e[INDEX4(l,1,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
371 grad_data_e[INDEX4(l,2,q,isub, numComps,DIM,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e, numShapesTotal,DIM,numQuad,numSub)];
372 }
373 }
374 }
375 }
376 }
377 #undef DIM
378 }
379 }
380 } /* end parallel region */
381 }
382 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26