/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_gradient.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_gradient.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 19320 byte(s)
Much needed sync with trunk...

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 /****************************************************************************
18
19 Assemblage of Jacobians: calculate the gradient of nodal data at quadrature
20 points
21
22 *****************************************************************************/
23
24 #include "Assemble.h"
25 #include "Util.h"
26
27 // Unless the loops in here get complicated again this file should be compiled
28 // with loop unrolling
29
30 namespace dudley {
31
32 void Assemble_gradient(Dudley_NodeFile* nodes, Dudley_ElementFile* elements,
33 escript::Data* grad_data, const escript::Data* data)
34 {
35 if (!nodes || !elements)
36 return;
37
38 const int numComps = data->getDataPointSize();
39 const int NN = elements->numNodes;
40 const bool reducedIntegrationOrder = Assemble_reducedIntegrationOrder(grad_data);
41 const int data_type = data->getFunctionSpace().getTypeCode();
42
43 dim_t numNodes = 0;
44 if (data_type == DUDLEY_NODES) {
45 numNodes = nodes->nodesMapping->numTargets;
46 } else if (data_type == DUDLEY_REDUCED_NODES) {
47 numNodes = nodes->reducedNodesMapping->numTargets;
48 } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
49 if (elements->MPIInfo->size > 1) {
50 throw DudleyException("Assemble_gradient: for more than one "
51 "processor DEGREES_OF_FREEDOM data are not accepted as input.");
52 }
53 numNodes = nodes->degreesOfFreedomMapping->numTargets;
54 } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
55 if (elements->MPIInfo->size > 1) {
56 throw DudleyException("Assemble_gradient: for more than one "
57 "processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
58 }
59 numNodes = nodes->reducedDegreesOfFreedomMapping->numTargets;
60 } else {
61 throw DudleyException("Assemble_gradient: Cannot calculate gradient "
62 "of data because of unsuitable input data representation.");
63 }
64
65 Dudley_ElementFile_Jacobians* jac = Dudley_ElementFile_borrowJacobians(
66 elements, nodes, reducedIntegrationOrder);
67 const int numDim = jac->numDim;
68 const int numShapesTotal = jac->numShapes;
69 const int numQuad = jac->numQuad;
70 const size_t localGradSize = sizeof(double) * numDim * numQuad * numComps;
71
72 // check the dimensions of data
73 if (!grad_data->numSamplesEqual(numQuad, elements->numElements)) {
74 throw DudleyException("Assemble_gradient: illegal number of samples in gradient Data object");
75 } else if (!data->numSamplesEqual(1, numNodes)) {
76 throw DudleyException("Assemble_gradient: illegal number of samples of input Data object");
77 } else if (numDim * numComps != grad_data->getDataPointSize()) {
78 throw DudleyException("Assemble_gradient: illegal number of components in gradient data object.");
79 } else if (!grad_data->actsExpanded()) {
80 throw DudleyException("Assemble_gradient: expanded Data object is expected for output data.");
81 }
82
83 grad_data->requireWrite();
84 #pragma omp parallel
85 {
86 if (data_type == DUDLEY_NODES) {
87 if (numDim == 1) {
88 const dim_t numShapes = 2;
89 #pragma omp for
90 for (index_t e = 0; e < elements->numElements; e++) {
91 double* grad_data_e = grad_data->getSampleDataRW(e);
92 memset(grad_data_e, 0, localGradSize);
93 for (int s = 0; s < numShapes; s++) {
94 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
95 const double* data_array = data->getSampleDataRO(n);
96 for (int q = 0; q < numQuad; q++) {
97 #pragma ivdep
98 for (int l = 0; l < numComps; l++) {
99 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
100 data_array[l] *
101 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
102 }
103 }
104 }
105 }
106 } else if (numDim == 2) {
107 const dim_t numShapes = 3;
108 #pragma omp for
109 for (index_t e = 0; e < elements->numElements; e++) {
110 double* grad_data_e = grad_data->getSampleDataRW(e);
111 memset(grad_data_e, 0, localGradSize);
112 for (int s = 0; s < numShapes; s++) {
113 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
114 const double* data_array = data->getSampleDataRO(n);
115 for (int q = 0; q < numQuad; q++) {
116 #pragma ivdep
117 for (int l = 0; l < numComps; l++) {
118 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
119 data_array[l] *
120 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
121 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
122 data_array[l] *
123 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
124 }
125 }
126 }
127 }
128 } else if (numDim == 3) {
129 const dim_t numShapes = 4;
130 #pragma omp for
131 for (index_t e = 0; e < elements->numElements; e++) {
132 double* grad_data_e = grad_data->getSampleDataRW(e);
133 memset(grad_data_e, 0, localGradSize);
134 for (int s = 0; s < numShapes; s++) {
135 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
136 const double* data_array = data->getSampleDataRO(n);
137 for (int q = 0; q < numQuad; q++) {
138 #pragma ivdep
139 for (int l = 0; l < numComps; l++) {
140 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
141 data_array[l] *
142 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
143 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
144 data_array[l] *
145 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
146 grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
147 data_array[l] *
148 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
149 }
150 }
151 }
152 }
153 }
154 } else if (data_type == DUDLEY_REDUCED_NODES) {
155 if (numDim == 1) {
156 const dim_t numShapes = 2;
157 #pragma omp for
158 for (index_t e = 0; e < elements->numElements; e++) {
159 double* grad_data_e = grad_data->getSampleDataRW(e);
160 memset(grad_data_e, 0, localGradSize);
161 for (int s = 0; s < numShapes; s++) {
162 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
163 const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[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, 0, numComps, numDim, numQuad)] +=
168 data_array[l] *
169 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
170 }
171 }
172 }
173 }
174 } else if (numDim == 2) {
175 const dim_t numShapes = 3;
176 #pragma omp for
177 for (index_t e = 0; e < elements->numElements; e++) {
178 double* grad_data_e = grad_data->getSampleDataRW(e);
179 memset(grad_data_e, 0, localGradSize);
180 for (int s = 0; s < numShapes; s++) {
181 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
182 const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->target[n]);
183 for (int q = 0; q < numQuad; q++) {
184 #pragma ivdep
185 for (int l = 0; l < numComps; l++) {
186 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
187 data_array[l] *
188 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
189 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
190 data_array[l] *
191 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
192 }
193 }
194 }
195 }
196 } else if (numDim == 3) {
197 const dim_t numShapes = 4;
198 #pragma omp for
199 for (index_t e = 0; e < elements->numElements; e++) {
200 double* grad_data_e = grad_data->getSampleDataRW(e);
201 memset(grad_data_e, 0, localGradSize);
202 for (int s = 0; s < numShapes; s++) {
203 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
204 const double* data_array = data->getSampleDataRO(nodes->reducedNodesMapping->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, 0, numComps, numDim, numQuad)] +=
209 data_array[l] *
210 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
211 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
212 data_array[l] *
213 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
214 grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
215 data_array[l] *
216 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
217 }
218 }
219 }
220 }
221 }
222 } else if (data_type == DUDLEY_DEGREES_OF_FREEDOM) {
223 if (numDim == 1)
224 {
225 const dim_t numShapes = 2;
226 #pragma omp for
227 for (index_t e = 0; e < elements->numElements; e++) {
228 double* grad_data_e = grad_data->getSampleDataRW(e);
229 memset(grad_data_e, 0, localGradSize);
230 for (int s = 0; s < numShapes; s++) {
231 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
232 const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
233 for (int q = 0; q < numQuad; q++) {
234 #pragma ivdep
235 for (int l = 0; l < numComps; l++) {
236 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
237 data_array[l] *
238 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
239 }
240 }
241 }
242 }
243 } else if (numDim == 2) {
244 const dim_t numShapes = 3;
245 #pragma omp for
246 for (index_t e = 0; e < elements->numElements; e++) {
247 double* grad_data_e = grad_data->getSampleDataRW(e);
248 memset(grad_data_e, 0, localGradSize);
249 for (int s = 0; s < numShapes; s++) {
250 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
251 const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
252 for (int q = 0; q < numQuad; q++) {
253 #pragma ivdep
254 for (int l = 0; l < numComps; l++) {
255 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
256 data_array[l] *
257 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
258 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
259 data_array[l] *
260 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
261 }
262 }
263 }
264 }
265 } else if (numDim == 3) {
266 const dim_t numShapes = 4;
267 #pragma omp for
268 for (index_t e = 0; e < elements->numElements; e++) {
269 double* grad_data_e = grad_data->getSampleDataRW(e);
270 memset(grad_data_e, 0, localGradSize);
271 for (int s = 0; s < numShapes; s++) {
272 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
273 const double* data_array = data->getSampleDataRO(nodes->degreesOfFreedomMapping->target[n]);
274 for (int q = 0; q < numQuad; q++) {
275 #pragma ivdep
276 for (int l = 0; l < numComps; l++) {
277 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
278 data_array[l] *
279 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
280 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
281 data_array[l] *
282 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
283 grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
284 data_array[l] *
285 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
286 }
287 }
288 }
289 }
290 }
291 } else if (data_type == DUDLEY_REDUCED_DEGREES_OF_FREEDOM) {
292 if (numDim == 1) {
293 const dim_t numShapes = 2;
294 #pragma omp for
295 for (index_t e = 0; e < elements->numElements; e++) {
296 double* grad_data_e = grad_data->getSampleDataRW(e);
297 memset(grad_data_e, 0, localGradSize);
298 for (int s = 0; s < numShapes; s++) {
299 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
300 const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
301 for (int q = 0; q < numQuad; q++) {
302 #pragma ivdep
303 for (int l = 0; l < numComps; l++) {
304 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
305 data_array[l] *
306 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
307 }
308 }
309 }
310 }
311 } else if (numDim == 2) {
312 const dim_t numShapes = 3;
313 #pragma omp for
314 for (index_t e = 0; e < elements->numElements; e++) {
315 double* grad_data_e = grad_data->getSampleDataRW(e);
316 memset(grad_data_e, 0, localGradSize);
317 for (int s = 0; s < numShapes; s++) {
318 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
319 const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
320 for (int q = 0; q < numQuad; q++) {
321 #pragma ivdep
322 for (int l = 0; l < numComps; l++) {
323 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
324 data_array[l] *
325 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
326 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
327 data_array[l] *
328 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
329 }
330 }
331 }
332 }
333 } else if (numDim == 3) {
334 const dim_t numShapes = 4;
335 #pragma omp for
336 for (index_t e = 0; e < elements->numElements; e++) {
337 double* grad_data_e = grad_data->getSampleDataRW(e);
338 memset(grad_data_e, 0, localGradSize);
339 for (int s = 0; s < numShapes; s++) {
340 const index_t n = elements->Nodes[INDEX2(s, e, NN)];
341 const double* data_array = data->getSampleDataRO(nodes->reducedDegreesOfFreedomMapping->target[n]);
342 for (int q = 0; q < numQuad; q++) {
343 #pragma ivdep
344 for (int l = 0; l < numComps; l++) {
345 grad_data_e[INDEX4(l, 0, q, 0, numComps, numDim, numQuad)] +=
346 data_array[l] *
347 jac->DSDX[INDEX5(s, 0, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
348 grad_data_e[INDEX4(l, 1, q, 0, numComps, numDim, numQuad)] +=
349 data_array[l] *
350 jac->DSDX[INDEX5(s, 1, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
351 grad_data_e[INDEX4(l, 2, q, 0, numComps, numDim, numQuad)] +=
352 data_array[l] *
353 jac->DSDX[INDEX5(s, 2, q, 0, e, numShapesTotal, numDim, numQuad, 1)];
354 }
355 }
356 }
357 }
358 }
359 }
360 } // end parallel region
361 }
362
363 } // namespace dudley
364

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Assemble_gradient.cpp:5567-5588 /branches/complex/dudley/src/Assemble_gradient.cpp:5866-5937 /branches/lapack2681/finley/src/Assemble_gradient.cpp:2682-2741 /branches/pasowrap/dudley/src/Assemble_gradient.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Assemble_gradient.cpp:3871-3891 /branches/restext/finley/src/Assemble_gradient.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Assemble_gradient.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_gradient.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Assemble_gradient.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_gradient.cpp:3517-3974 /release/3.0/finley/src/Assemble_gradient.cpp:2591-2601 /release/4.0/dudley/src/Assemble_gradient.cpp:5380-5406 /trunk/dudley/src/Assemble_gradient.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Assemble_gradient.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26