/[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 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 20165 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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 namespace finley {
28
29 void Assemble_gradient(const NodeFile* nodes, const ElementFile* elements,
30 escript::Data& grad_data, const escript::Data& data)
31 {
32 Finley_resetError();
33 if (!nodes || !elements)
34 return;
35
36 const int numComps=data.getDataPointSize();
37 const int NN=elements->numNodes;
38 const bool reducedOrder=util::hasReducedIntegrationOrder(grad_data);
39 const int data_type=data.getFunctionSpace().getTypeCode();
40
41 bool reducedShapefunction = false;
42 int numNodes = 0;
43 if (data_type == FINLEY_NODES) {
44 numNodes = nodes->getNumNodes();
45 } else if (data_type==FINLEY_REDUCED_NODES) {
46 reducedShapefunction = true;
47 numNodes = nodes->getNumReducedNodes();
48 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
49 if (elements->MPIInfo->size > 1) {
50 Finley_setError(TYPE_ERROR, "Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
51 return;
52 }
53 numNodes = nodes->getNumDegreesOfFreedom();
54 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
55 if (elements->MPIInfo->size > 1) {
56 Finley_setError(TYPE_ERROR, "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->getNumReducedDegreesOfFreedom();
61 } else {
62 Finley_setError(TYPE_ERROR, "Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
63 return;
64 }
65
66 ElementFile_Jacobians *jac = elements->borrowJacobians(nodes,
67 reducedShapefunction, reducedOrder);
68 ReferenceElement *refElement = ReferenceElementSet_borrowReferenceElement(
69 elements->referenceElementSet, reducedOrder);
70 const int numDim=jac->numDim;
71 const int numShapes=jac->BasisFunctions->Type->numShapes;
72 const int numShapesTotal=jac->numShapesTotal;
73 const int numSub=jac->numSub;
74 const int numQuad=jac->numQuadTotal/numSub;
75 int numShapesTotal2=0;
76 int s_offset=0, *nodes_selector=NULL;
77
78 if (Finley_noError()) {
79 const int grad_data_type=grad_data.getFunctionSpace().getTypeCode();
80 if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2 || grad_data_type==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
81 s_offset=jac->offsets[1];
82 s_offset=jac->offsets[1];
83 } else {
84 s_offset=jac->offsets[0];
85 s_offset=jac->offsets[0];
86 }
87 if (data_type==FINLEY_REDUCED_NODES || data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
88 nodes_selector=refElement->Type->linearNodes;
89 numShapesTotal2=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
90 } else {
91 nodes_selector=refElement->Type->subElementNodes;
92 numShapesTotal2=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
93 }
94
95 // check the dimensions of data
96 if (!grad_data.numSamplesEqual(numQuad*numSub, elements->numElements)) {
97 Finley_setError(TYPE_ERROR, "Assemble_gradient: illegal number of samples in gradient Data object");
98 } else if (!data.numSamplesEqual(1, numNodes)) {
99 Finley_setError(TYPE_ERROR, "Assemble_gradient: illegal number of samples of input Data object");
100 } else if (numDim*numComps != grad_data.getDataPointSize()) {
101 Finley_setError(TYPE_ERROR, "Assemble_gradient: illegal number of components in gradient data object.");
102 } else if (!grad_data.actsExpanded()) {
103 Finley_setError(TYPE_ERROR, "Assemble_gradient: expanded Data object is expected for output data.");
104 } else if (!(s_offset+numShapes <= numShapesTotal)) {
105 Finley_setError(SYSTEM_ERROR, "Assemble_gradient: nodes per element is inconsistent with number of jacobians.");
106 }
107 }
108
109 if (!Finley_noError())
110 return;
111
112 const size_t localGradSize=sizeof(double)*numDim*numQuad*numSub*numComps;
113 escript::Data& in(*const_cast<escript::Data*>(&data));
114 grad_data.requireWrite();
115 #pragma omp parallel
116 {
117 if (data_type==FINLEY_NODES) {
118 if (numDim==1) {
119 #pragma omp for
120 for (int e=0; e<elements->numElements; e++) {
121 double *grad_data_e=grad_data.getSampleDataRW(e);
122 memset(grad_data_e, 0, localGradSize);
123 for (int isub=0; isub<numSub; isub++) {
124 for (int s=0; s<numShapes; s++) {
125 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
126 const double *data_array=in.getSampleDataRO(n);
127 for (int q=0; q<numQuad; q++) {
128 #pragma ivdep
129 for (int l=0; l<numComps; l++) {
130 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)];
131 }
132 }
133 }
134 }
135 }
136 } else if (numDim==2) {
137 #pragma omp for
138 for (int e=0; e<elements->numElements; e++) {
139 double *grad_data_e=grad_data.getSampleDataRW(e);
140 memset(grad_data_e, 0, localGradSize);
141 for (int isub=0; isub<numSub; isub++) {
142 for (int s=0; s<numShapes; s++) {
143 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
144 const double *data_array=in.getSampleDataRO(n);
145 for (int q=0; q<numQuad; q++) {
146 #pragma ivdep
147 for (int l=0; l<numComps; l++) {
148 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)];
149 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)];
150 }
151 }
152 }
153 }
154 }
155 } else if (numDim==3) {
156 #pragma omp for
157 for (int e=0; e<elements->numElements; e++) {
158 double *grad_data_e=grad_data.getSampleDataRW(e);
159 memset(grad_data_e,0, localGradSize);
160 for (int isub=0; isub<numSub; isub++) {
161 for (int s=0; s<numShapes; s++) {
162 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
163 const double *data_array=in.getSampleDataRO(n);
164 for (int q=0; q<numQuad; q++) {
165 #pragma ivdep
166 for (int l=0; l<numComps; l++) {
167 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)];
168 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)];
169 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)];
170 }
171 }
172 }
173 }
174 }
175 }
176 } else if (data_type==FINLEY_REDUCED_NODES) {
177 const int* target = nodes->borrowTargetReducedNodes();
178 if (numDim==1) {
179 #pragma omp for
180 for (int e=0; e<elements->numElements; e++) {
181 double *grad_data_e=grad_data.getSampleDataRW(e);
182 memset(grad_data_e, 0, localGradSize);
183 for (int isub=0; isub<numSub; isub++) {
184 for (int s=0;s<numShapes;s++) {
185 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
186 const double *data_array=in.getSampleDataRO(target[n]);
187 for (int q=0; q<numQuad; q++) {
188 #pragma ivdep
189 for (int 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
198 for (int e=0; e<elements->numElements; e++) {
199 double *grad_data_e=grad_data.getSampleDataRW(e);
200 memset(grad_data_e, 0, localGradSize);
201 for (int isub=0; isub<numSub; isub++) {
202 for (int s=0; s<numShapes; s++) {
203 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
204 const double *data_array=in.getSampleDataRO(target[n]);
205 for (int q=0; q<numQuad; q++) {
206 #pragma ivdep
207 for (int 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
217 for (int e=0;e<elements->numElements;e++) {
218 double *grad_data_e=grad_data.getSampleDataRW(e);
219 memset(grad_data_e, 0, localGradSize);
220 for (int isub=0; isub<numSub; isub++) {
221 for (int s=0; s<numShapes; s++) {
222 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
223 const double *data_array=in.getSampleDataRO(target[n]);
224 for (int q=0; q<numQuad; q++) {
225 #pragma ivdep
226 for (int 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 const int* target = nodes->borrowTargetDegreesOfFreedom();
238 if (numDim==1) {
239 #pragma omp for
240 for (int e=0; e<elements->numElements; e++) {
241 double *grad_data_e=grad_data.getSampleDataRW(e);
242 memset(grad_data_e, 0, localGradSize);
243 for (int isub=0; isub<numSub; isub++) {
244 for (int s=0; s<numShapes; s++) {
245 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
246 const double *data_array=in.getSampleDataRO(target[n]);
247 for (int q=0; q<numQuad; q++) {
248 #pragma ivdep
249 for (int 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
258 for (int e=0; e<elements->numElements; e++) {
259 double *grad_data_e=grad_data.getSampleDataRW(e);
260 memset(grad_data_e, 0, localGradSize);
261 for (int isub=0; isub<numSub; isub++) {
262 for (int s=0; s<numShapes; s++) {
263 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
264 const double *data_array=in.getSampleDataRO(target[n]);
265 for (int q=0; q<numQuad; q++) {
266 #pragma ivdep
267 for (int 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
277 for (int e=0; e<elements->numElements; e++) {
278 double *grad_data_e=grad_data.getSampleDataRW(e);
279 memset(grad_data_e, 0, localGradSize);
280 for (int isub=0; isub<numSub; isub++) {
281 for (int s=0; s<numShapes; s++) {
282 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
283 const double *data_array=in.getSampleDataRO(target[n]);
284 for (int q=0; q<numQuad; q++) {
285 #pragma ivdep
286 for (int 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 const int* target = nodes->borrowTargetReducedDegreesOfFreedom();
298 if (numDim==1) {
299 #pragma omp for
300 for (int e=0; e<elements->numElements; e++) {
301 double *grad_data_e=grad_data.getSampleDataRW(e);
302 memset(grad_data_e,0, localGradSize);
303 for (int isub=0; isub<numSub; isub++) {
304 for (int s=0; s<numShapes; s++) {
305 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
306 const double *data_array=in.getSampleDataRO(target[n]);
307 for (int q=0; q<numQuad; q++) {
308 #pragma ivdep
309 for (int l=0; l<numComps; l++) {
310 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)];
311 }
312 }
313 }
314 }
315 }
316 } else if (numDim==2) {
317 #pragma omp for
318 for (int e=0; e<elements->numElements; e++) {
319 double *grad_data_e=grad_data.getSampleDataRW(e);
320 memset(grad_data_e, 0, localGradSize);
321 for (int isub=0; isub<numSub; isub++) {
322 for (int s=0; s<numShapes; s++) {
323 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
324 const double *data_array=in.getSampleDataRO(target[n]);
325 for (int q=0; q<numQuad; q++) {
326 #pragma ivdep
327 for (int l=0; l<numComps; l++) {
328 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)];
329 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)];
330 }
331 }
332 }
333 }
334 }
335
336 } else if (numDim==3) {
337 #pragma omp for
338 for (int e=0; e<elements->numElements; e++) {
339 double *grad_data_e=grad_data.getSampleDataRW(e);
340 memset(grad_data_e,0, localGradSize);
341 for (int isub=0; isub<numSub; isub++) {
342 for (int s=0; s<numShapes; s++) {
343 const int n=elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
344 const double *data_array=in.getSampleDataRO(target[n]);
345 for (int q=0; q<numQuad; q++) {
346 #pragma ivdep
347 for (int l=0; l<numComps; l++) {
348 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)];
349 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)];
350 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)];
351 }
352 }
353 }
354 }
355 }
356 } // numDim
357 } // data_type
358 } // end parallel region
359 }
360
361 } // namespace finley
362

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