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

Contents of /trunk/finley/src/Assemble_gradient.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4421 - (show annotations)
Tue May 21 05:16:39 2013 UTC (6 years ago) by caltinay
File size: 20801 byte(s)
Finley... there were actually some long-standing errors in the checks.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /****************************************************************************
18
19 Assemblage of jacobians: calculates the gradient of nodal data at
20 quadrature points
21
22 *****************************************************************************/
23
24 #include "Assemble.h"
25 #include "Util.h"
26
27 void Finley_Assemble_gradient(Finley_NodeFile* nodes,
28 Finley_ElementFile* elements,
29 escriptDataC* grad_data,
30 escriptDataC* data)
31 {
32 Finley_resetError();
33 if (!nodes || !elements)
34 return;
35
36 const dim_t numComps=getDataPointSize(data);
37 const dim_t NN=elements->numNodes;
38 const bool_t reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
39
40 const type_t data_type=getFunctionSpaceType(data);
41 bool_t reducedShapefunction = FALSE;
42 dim_t numNodes = 0;
43 if (data_type == FINLEY_NODES) {
44 numNodes = nodes->nodesMapping->numTargets;
45 } else if (data_type==FINLEY_REDUCED_NODES) {
46 reducedShapefunction = TRUE;
47 numNodes = nodes->reducedNodesMapping->numTargets;
48 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
49 if (elements->MPIInfo->size > 1) {
50 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
51 return;
52 }
53 numNodes = nodes->degreesOfFreedomMapping->numTargets;
54 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
55 if (elements->MPIInfo->size > 1) {
56 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
57 return;
58 }
59 reducedShapefunction = TRUE;
60 numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
61 } else {
62 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
63 return;
64 }
65
66 Finley_ElementFile_Jacobeans *jac = Finley_ElementFile_borrowJacobeans(
67 elements, nodes, reducedShapefunction, reducedIntegrationOrder);
68 Finley_ReferenceElement *refElement =
69 Finley_ReferenceElementSet_borrowReferenceElement(
70 elements->referenceElementSet, reducedIntegrationOrder);
71 const dim_t numDim=jac->numDim;
72 const dim_t numShapes=jac->BasisFunctions->Type->numShapes;
73 const dim_t numShapesTotal=jac->numShapesTotal;
74 const dim_t numSub=jac->numSub;
75 const dim_t numQuad=jac->numQuadTotal/numSub;
76 dim_t numShapesTotal2=0;
77 index_t s_offset=0, *nodes_selector=NULL;
78
79 if (Finley_noError()) {
80 const type_t grad_data_type=getFunctionSpaceType(grad_data);
81 if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2 || grad_data_type==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
82 s_offset=jac->offsets[1];
83 s_offset=jac->offsets[1];
84 } else {
85 s_offset=jac->offsets[0];
86 s_offset=jac->offsets[0];
87 }
88 if (data_type==FINLEY_REDUCED_NODES || data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
89 nodes_selector=refElement->Type->linearNodes;
90 numShapesTotal2=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
91 } else {
92 nodes_selector=refElement->Type->subElementNodes;
93 numShapesTotal2=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
94 }
95
96 // check the dimensions of data
97 if (!numSamplesEqual(grad_data, numQuad*numSub, elements->numElements)) {
98 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: illegal number of samples in gradient Data object");
99 } else if (!numSamplesEqual(data, 1, numNodes)) {
100 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: illegal number of samples of input Data object");
101 } else if (numDim*numComps != getDataPointSize(grad_data)) {
102 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: illegal number of components in gradient data object.");
103 } else if (!isExpanded(grad_data)) {
104 Finley_setError(TYPE_ERROR, "Finley_Assemble_gradient: expanded Data object is expected for output data.");
105 } else if (!(s_offset+numShapes <= numShapesTotal)) {
106 Finley_setError(SYSTEM_ERROR, "Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobians.");
107 }
108 }
109
110 // now we can start
111 if (!Finley_noError())
112 return;
113
114 const size_t localGradSize=sizeof(double)*numDim*numQuad*numSub*numComps;
115 requireWrite(grad_data);
116 #pragma omp parallel
117 {
118 if (data_type==FINLEY_NODES) {
119 if (numDim==1) {
120 #pragma omp for schedule(static)
121 for (dim_t e=0; e<elements->numElements; e++) {
122 double *grad_data_e=getSampleDataRW(grad_data,e);
123 memset(grad_data_e, 0, localGradSize);
124 for (dim_t isub=0; isub<numSub; isub++) {
125 for (dim_t s=0; s<numShapes; s++) {
126 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
127 const double *data_array=getSampleDataRO(data,n);
128 for (dim_t q=0; q<numQuad; q++) {
129 #pragma ivdep
130 for (dim_t l=0; l<numComps; l++) {
131 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
132 }
133 }
134 }
135 }
136 }
137 } else if (numDim==2) {
138 #pragma omp for schedule(static)
139 for (dim_t e=0; e<elements->numElements; e++) {
140 double *grad_data_e=getSampleDataRW(grad_data,e);
141 memset(grad_data_e, 0, localGradSize);
142 for (dim_t isub=0; isub<numSub; isub++) {
143 for (dim_t s=0; s<numShapes; s++) {
144 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
145 const double *data_array=getSampleDataRO(data,n);
146 for (dim_t q=0; q<numQuad; q++) {
147 #pragma ivdep
148 for (dim_t l=0; l<numComps; l++) {
149 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
150 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
151 }
152 }
153 }
154 }
155 }
156 } else if (numDim==3) {
157 #pragma omp for schedule(static)
158 for (dim_t e=0; e<elements->numElements; e++) {
159 double *grad_data_e=getSampleDataRW(grad_data,e);
160 memset(grad_data_e,0, localGradSize);
161 for (dim_t isub=0; isub<numSub; isub++) {
162 for (dim_t s=0; s<numShapes; s++) {
163 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
164 const double *data_array=getSampleDataRO(data,n);
165 for (dim_t q=0; q<numQuad; q++) {
166 #pragma ivdep
167 for (dim_t l=0; l<numComps; l++) {
168 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
169 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
170 grad_data_e[INDEX4(l,2,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
171 }
172 }
173 }
174 }
175 }
176 }
177 } else if (data_type==FINLEY_REDUCED_NODES) {
178 if (numDim==1) {
179 #pragma omp for schedule(static)
180 for (dim_t e=0; e<elements->numElements; e++) {
181 double *grad_data_e=getSampleDataRW(grad_data,e);
182 memset(grad_data_e, 0, localGradSize);
183 for (dim_t isub=0; isub<numSub; isub++) {
184 for (dim_t s=0;s<numShapes;s++) {
185 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
186 const double *data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
187 for (dim_t q=0; q<numQuad; q++) {
188 #pragma ivdep
189 for (dim_t l=0; l<numComps; l++) {
190 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
191 }
192 }
193 }
194 }
195 }
196 } else if (numDim==2) {
197 #pragma omp for schedule(static)
198 for (dim_t e=0; e<elements->numElements; e++) {
199 double *grad_data_e=getSampleDataRW(grad_data,e);
200 memset(grad_data_e, 0, localGradSize);
201 for (dim_t isub=0; isub<numSub; isub++) {
202 for (dim_t s=0; s<numShapes; s++) {
203 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
204 const double *data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
205 for (dim_t q=0; q<numQuad; q++) {
206 #pragma ivdep
207 for (dim_t l=0; l<numComps; l++) {
208 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
209 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
210 }
211 }
212 }
213 }
214 }
215 } else if (numDim==3) {
216 #pragma omp for schedule(static)
217 for (dim_t e=0;e<elements->numElements;e++) {
218 double *grad_data_e=getSampleDataRW(grad_data,e);
219 memset(grad_data_e, 0, localGradSize);
220 for (dim_t isub=0; isub<numSub; isub++) {
221 for (dim_t s=0; s<numShapes; s++) {
222 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
223 const double *data_array=getSampleDataRO(data,nodes->reducedNodesMapping->target[n]);
224 for (dim_t q=0; q<numQuad; q++) {
225 #pragma ivdep
226 for (dim_t l=0;l<numComps;l++) {
227 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
228 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
229 grad_data_e[INDEX4(l,2,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
230 }
231 }
232 }
233 }
234 }
235 }
236 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
237
238 if (numDim==1) {
239 #pragma omp for schedule(static)
240 for (dim_t e=0; e<elements->numElements; e++) {
241 double *grad_data_e=getSampleDataRW(grad_data,e);
242 memset(grad_data_e, 0, localGradSize);
243 for (dim_t isub=0; isub<numSub; isub++) {
244 for (dim_t s=0; s<numShapes; s++) {
245 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
246 const double *data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
247 for (dim_t q=0; q<numQuad; q++) {
248 #pragma ivdep
249 for (dim_t l=0; l<numComps; l++) {
250 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
251 }
252 }
253 }
254 }
255 }
256 } else if (numDim==2) {
257 #pragma omp for schedule(static)
258 for (dim_t e=0; e<elements->numElements; e++) {
259 double *grad_data_e=getSampleDataRW(grad_data,e);
260 memset(grad_data_e, 0, localGradSize);
261 for (dim_t isub=0; isub<numSub; isub++) {
262 for (dim_t s=0; s<numShapes; s++) {
263 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
264 const double *data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
265 for (dim_t q=0; q<numQuad; q++) {
266 #pragma ivdep
267 for (dim_t l=0; l<numComps; l++) {
268 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
269 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
270 }
271 }
272 }
273 }
274 }
275 } else if (numDim==3) {
276 #pragma omp for schedule(static)
277 for (dim_t e=0; e<elements->numElements; e++) {
278 double *grad_data_e=getSampleDataRW(grad_data,e);
279 memset(grad_data_e, 0, localGradSize);
280 for (dim_t isub=0; isub<numSub; isub++) {
281 for (dim_t s=0; s<numShapes; s++) {
282 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
283 const double *data_array=getSampleDataRO(data,nodes->degreesOfFreedomMapping->target[n]);
284 for (dim_t q=0; q<numQuad; q++) {
285 #pragma ivdep
286 for (dim_t l=0; l<numComps; l++) {
287 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
288 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
289 grad_data_e[INDEX4(l,2,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
290 }
291 }
292 }
293 }
294 }
295 }
296 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
297 if (numDim==1) {
298 #pragma omp for schedule(static)
299 for (dim_t e=0; e<elements->numElements; e++) {
300 double *grad_data_e=getSampleDataRW(grad_data,e);
301 memset(grad_data_e,0, localGradSize);
302 for (dim_t isub=0; isub<numSub; isub++) {
303 for (dim_t s=0; s<numShapes; s++) {
304 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
305 const double *data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
306 for (dim_t q=0; q<numQuad; q++) {
307 #pragma ivdep
308 for (dim_t l=0; l<numComps; l++) {
309 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
310 }
311 }
312 }
313 }
314 }
315 } else if (numDim==2) {
316 #pragma omp for schedule(static)
317 for (dim_t e=0;e<elements->numElements;e++) {
318 double *grad_data_e=getSampleDataRW(grad_data,e);
319 memset(grad_data_e, 0, localGradSize);
320 for (dim_t isub=0; isub<numSub; isub++) {
321 for (dim_t s=0; s<numShapes; s++) {
322 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
323 const double *data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
324 for (dim_t q=0; q<numQuad; q++) {
325 #pragma ivdep
326 for (dim_t l=0; l<numComps; l++) {
327 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
328 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
329 }
330 }
331 }
332 }
333 }
334
335 } else if (numDim==3) {
336 #pragma omp for schedule(static)
337 for (dim_t e=0; e<elements->numElements; e++) {
338 double *grad_data_e=getSampleDataRW(grad_data,e);
339 memset(grad_data_e,0, localGradSize);
340 for (dim_t isub=0; isub<numSub; isub++) {
341 for (dim_t s=0; s<numShapes; s++) {
342 const dim_t n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
343 const double *data_array=getSampleDataRO(data,nodes->reducedDegreesOfFreedomMapping->target[n]);
344 for (dim_t q=0; q<numQuad; q++) {
345 #pragma ivdep
346 for (dim_t l=0; l<numComps; l++) {
347 grad_data_e[INDEX4(l,0,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,0,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
348 grad_data_e[INDEX4(l,1,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,1,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
349 grad_data_e[INDEX4(l,2,q,isub,numComps,numDim,numQuad)]+=data_array[l]*jac->DSDX[INDEX5(s_offset+s,2,q,isub,e,numShapesTotal,numDim,numQuad,numSub)];
350 }
351 }
352 }
353 }
354 }
355 } // numDim
356 } // data_type
357 } // end parallel region
358 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/Assemble_gradient.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_gradient.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_gradient.cpp:3871-3891 /branches/restext/finley/src/Assemble_gradient.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_gradient.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_gradient.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_gradient.cpp:3471-3974 /release/3.0/finley/src/Assemble_gradient.cpp:2591-2601 /trunk/finley/src/Assemble_gradient.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26