/[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 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 2 weeks ago) by uqaeller
File size: 21247 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18
19 /****************************************************************************
20
21 Assemblage of jacobians: calculates the gradient of nodal data at
22 quadrature points
23
24 *****************************************************************************/
25
26 #include "Assemble.h"
27 #include "Util.h"
28
29 #include <escript/index.h>
30
31 namespace finley {
32
33 template<typename Scalar>
34 void Assemble_gradient(const NodeFile* nodes, const ElementFile* elements,
35 escript::Data& out, const escript::Data& data)
36 {
37 if (!nodes || !elements)
38 return;
39
40 const int numComps = data.getDataPointSize();
41 const int NN = elements->numNodes;
42 const bool reducedOrder = util::hasReducedIntegrationOrder(out);
43 const int dataType = data.getFunctionSpace().getTypeCode();
44
45 bool reducedShapefunction = false;
46 dim_t numNodes = 0;
47 if (dataType == FINLEY_NODES) {
48 numNodes = nodes->getNumNodes();
49 } else if (dataType == FINLEY_REDUCED_NODES) {
50 reducedShapefunction = true;
51 numNodes = nodes->getNumReducedNodes();
52 } else if (dataType == FINLEY_DEGREES_OF_FREEDOM) {
53 if (elements->MPIInfo->size > 1) {
54 throw escript::ValueError("Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
55 }
56 numNodes = nodes->getNumDegreesOfFreedom();
57 } else if (dataType == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
58 if (elements->MPIInfo->size > 1) {
59 throw escript::ValueError("Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
60 }
61 reducedShapefunction = true;
62 numNodes = nodes->getNumReducedDegreesOfFreedom();
63 } else {
64 throw escript::ValueError("Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
65 }
66
67 ElementFile_Jacobians* jac = elements->borrowJacobians(nodes,
68 reducedShapefunction, reducedOrder);
69 const_ReferenceElement_ptr refElement(elements->referenceElementSet->
70 borrowReferenceElement(reducedOrder));
71 const int numDim = jac->numDim;
72 const int numShapes = jac->BasisFunctions->Type->numShapes;
73 const int numShapesTotal = jac->numShapesTotal;
74 const int numSub = jac->numSub;
75 const int numQuad = jac->numQuadTotal/numSub;
76 int numShapesTotal2 = 0;
77 int s_offset = 0;
78 const int* nodes_selector = NULL;
79
80 const int gradDataType = out.getFunctionSpace().getTypeCode();
81 if (gradDataType==FINLEY_CONTACT_ELEMENTS_2 || gradDataType==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
82 s_offset = jac->offsets[1];
83 } else {
84 s_offset = jac->offsets[0];
85 }
86 if (dataType==FINLEY_REDUCED_NODES || dataType==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
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
94 // check the dimensions of data
95 if (!out.numSamplesEqual(numQuad*numSub, elements->numElements)) {
96 throw escript::ValueError("Assemble_gradient: illegal number of samples in gradient Data object");
97 } else if (!data.numSamplesEqual(1, numNodes)) {
98 throw escript::ValueError("Assemble_gradient: illegal number of samples of input Data object");
99 } else if (numDim * numComps != out.getDataPointSize()) {
100 throw escript::ValueError("Assemble_gradient: illegal number of components in gradient data object.");
101 } else if (!out.actsExpanded()) {
102 throw escript::ValueError("Assemble_gradient: expanded Data object is expected for output data.");
103 } else if (!(s_offset+numShapes <= numShapesTotal)) {
104 throw escript::ValueError("Assemble_gradient: nodes per element is inconsistent with number of jacobians.");
105 }
106
107 const Scalar zero = static_cast<Scalar>(0);
108 const size_t localGradSize = numDim*numQuad*numSub*numComps;
109 out.requireWrite();
110 #pragma omp parallel
111 {
112 if (dataType == FINLEY_NODES) {
113 if (numDim == 1) {
114 #pragma omp for
115 for (index_t e = 0; e < elements->numElements; e++) {
116 Scalar* gradData_e = out.getSampleDataRW(e, zero);
117 std::fill(gradData_e, gradData_e+localGradSize, zero);
118 for (int isub = 0; isub < numSub; isub++) {
119 for (int s = 0; s < numShapes; s++) {
120 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
121 const Scalar* data_array = data.getSampleDataRO(n, zero);
122 for (int q = 0; q < numQuad; q++) {
123 #pragma ivdep
124 for (int l = 0; l < numComps; l++) {
125 gradData_e[INDEX4(l, 0, q, isub, numComps, numDim, numQuad)] +=
126 data_array[l] * jac->DSDX[INDEX5
127 (s_offset+s, 0, q, isub, e, numShapesTotal, numDim, numQuad, numSub)];
128 }
129 }
130 }
131 }
132 }
133 } else if (numDim == 2) {
134 #pragma omp for
135 for (index_t e = 0; e < elements->numElements; e++) {
136 Scalar* gradData_e = out.getSampleDataRW(e, zero);
137 std::fill(gradData_e, gradData_e+localGradSize, zero);
138 for (int isub = 0; isub < numSub; isub++) {
139 for (int s = 0; s < numShapes; s++) {
140 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
141 const Scalar* data_array = data.getSampleDataRO(n, zero);
142 for (int q = 0; q < numQuad; q++) {
143 #pragma ivdep
144 for (int l = 0; l < numComps; l++) {
145 gradData_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)];
146 gradData_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)];
147 }
148 }
149 }
150 }
151 }
152 } else if (numDim == 3) {
153 #pragma omp for
154 for (index_t e = 0; e < elements->numElements; e++) {
155 Scalar* gradData_e = out.getSampleDataRW(e, zero);
156 std::fill(gradData_e, gradData_e+localGradSize, zero);
157 for (int isub = 0; isub < numSub; isub++) {
158 for (int s = 0; s < numShapes; s++) {
159 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
160 const Scalar* data_array = data.getSampleDataRO(n, zero);
161 for (int q = 0; q < numQuad; q++) {
162 #pragma ivdep
163 for (int l = 0; l < numComps; l++) {
164 gradData_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)];
165 gradData_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)];
166 gradData_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)];
167 }
168 }
169 }
170 }
171 }
172 }
173 } else if (dataType == FINLEY_REDUCED_NODES) {
174 const index_t* target = nodes->borrowTargetReducedNodes();
175 if (numDim == 1) {
176 #pragma omp for
177 for (index_t e = 0; e < elements->numElements; e++) {
178 Scalar* gradData_e = out.getSampleDataRW(e, zero);
179 std::fill(gradData_e, gradData_e+localGradSize, zero);
180 for (int isub = 0; isub < numSub; isub++) {
181 for (int s = 0; s < numShapes; s++) {
182 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
183 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
184 for (int q = 0; q < numQuad; q++) {
185 #pragma ivdep
186 for (int l = 0; l < numComps; l++) {
187 gradData_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)];
188 }
189 }
190 }
191 }
192 }
193 } else if (numDim==2) {
194 #pragma omp for
195 for (index_t e = 0; e < elements->numElements; e++) {
196 Scalar* gradData_e = out.getSampleDataRW(e, zero);
197 std::fill(gradData_e, gradData_e+localGradSize, zero);
198 for (int isub = 0; isub < numSub; isub++) {
199 for (int s = 0; s < numShapes; s++) {
200 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
201 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
202 for (int q = 0; q < numQuad; q++) {
203 #pragma ivdep
204 for (int l = 0; l < numComps; l++) {
205 gradData_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)];
206 gradData_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)];
207 }
208 }
209 }
210 }
211 }
212 } else if (numDim==3) {
213 #pragma omp for
214 for (index_t e = 0; e < elements->numElements; e++) {
215 Scalar* gradData_e = out.getSampleDataRW(e, zero);
216 std::fill(gradData_e, gradData_e+localGradSize, zero);
217 for (int isub = 0; isub < numSub; isub++) {
218 for (int s = 0; s < numShapes; s++) {
219 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
220 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
221 for (int q = 0; q < numQuad; q++) {
222 #pragma ivdep
223 for (int l = 0; l < numComps;l++) {
224 gradData_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)];
225 gradData_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)];
226 gradData_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)];
227 }
228 }
229 }
230 }
231 }
232 }
233 } else if (dataType==FINLEY_DEGREES_OF_FREEDOM) {
234 const index_t* target = nodes->borrowTargetDegreesOfFreedom();
235 if (numDim==1) {
236 #pragma omp for
237 for (index_t e = 0; e < elements->numElements; e++) {
238 Scalar* gradData_e = out.getSampleDataRW(e, zero);
239 std::fill(gradData_e, gradData_e+localGradSize, zero);
240 for (int isub = 0; isub < numSub; isub++) {
241 for (int s = 0; s < numShapes; s++) {
242 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
243 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
244 for (int q = 0; q < numQuad; q++) {
245 #pragma ivdep
246 for (int l = 0; l < numComps; l++) {
247 gradData_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)];
248 }
249 }
250 }
251 }
252 }
253 } else if (numDim==2) {
254 #pragma omp for
255 for (index_t e = 0; e < elements->numElements; e++) {
256 Scalar* gradData_e = out.getSampleDataRW(e, zero);
257 std::fill(gradData_e, gradData_e+localGradSize, zero);
258 for (int isub = 0; isub < numSub; isub++) {
259 for (int s = 0; s < numShapes; s++) {
260 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
261 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
262 for (int q = 0; q < numQuad; q++) {
263 #pragma ivdep
264 for (int l = 0; l < numComps; l++) {
265 gradData_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)];
266 gradData_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)];
267 }
268 }
269 }
270 }
271 }
272 } else if (numDim==3) {
273 #pragma omp for
274 for (index_t e = 0; e < elements->numElements; e++) {
275 Scalar* gradData_e = out.getSampleDataRW(e, zero);
276 std::fill(gradData_e, gradData_e+localGradSize, zero);
277 for (int isub = 0; isub < numSub; isub++) {
278 for (int s = 0; s < numShapes; s++) {
279 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
280 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
281 for (int q = 0; q < numQuad; q++) {
282 #pragma ivdep
283 for (int l = 0; l < numComps; l++) {
284 gradData_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)];
285 gradData_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)];
286 gradData_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)];
287 }
288 }
289 }
290 }
291 }
292 }
293 } else if (dataType==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
294 const index_t* target = nodes->borrowTargetReducedDegreesOfFreedom();
295 if (numDim==1) {
296 #pragma omp for
297 for (index_t e = 0; e < elements->numElements; e++) {
298 Scalar* gradData_e = out.getSampleDataRW(e, zero);
299 std::fill(gradData_e, gradData_e+localGradSize, zero);
300 for (int isub = 0; isub < numSub; isub++) {
301 for (int s = 0; s < numShapes; s++) {
302 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
303 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
304 for (int q = 0; q < numQuad; q++) {
305 #pragma ivdep
306 for (int l = 0; l < numComps; l++) {
307 gradData_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)];
308 }
309 }
310 }
311 }
312 }
313 } else if (numDim==2) {
314 #pragma omp for
315 for (index_t e = 0; e < elements->numElements; e++) {
316 Scalar* gradData_e = out.getSampleDataRW(e, zero);
317 std::fill(gradData_e, gradData_e+localGradSize, zero);
318 for (int isub = 0; isub < numSub; isub++) {
319 for (int s = 0; s < numShapes; s++) {
320 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
321 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
322 for (int q = 0; q < numQuad; q++) {
323 #pragma ivdep
324 for (int l = 0; l < numComps; l++) {
325 gradData_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)];
326 gradData_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)];
327 }
328 }
329 }
330 }
331 }
332
333 } else if (numDim==3) {
334 #pragma omp for
335 for (index_t e = 0; e < elements->numElements; e++) {
336 Scalar* gradData_e = out.getSampleDataRW(e, zero);
337 std::fill(gradData_e, gradData_e+localGradSize, zero);
338 for (int isub = 0; isub < numSub; isub++) {
339 for (int s = 0; s < numShapes; s++) {
340 const index_t n = elements->Nodes[INDEX2(nodes_selector[INDEX2(s_offset+s,isub,numShapesTotal2)],e, NN)];
341 const Scalar* data_array = data.getSampleDataRO(target[n], zero);
342 for (int q = 0; q < numQuad; q++) {
343 #pragma ivdep
344 for (int l = 0; l < numComps; l++) {
345 gradData_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)];
346 gradData_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)];
347 gradData_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)];
348 }
349 }
350 }
351 }
352 }
353 } // numDim
354 } // dataType
355 } // end parallel region
356 }
357
358 // instantiate our two supported versions
359 template void Assemble_gradient<escript::DataTypes::real_t>(
360 const NodeFile* nodes, const ElementFile* elements,
361 escript::Data& out, const escript::Data& data);
362 template void Assemble_gradient<escript::DataTypes::cplx_t>(
363 const NodeFile* nodes, const ElementFile* elements,
364 escript::Data& out, const escript::Data& data);
365
366 } // namespace finley
367

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_gradient.cpp:5567-5588 /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 /branches/trilinos_from_5897/finley/src/Assemble_gradient.cpp:5898-6118 /release/3.0/finley/src/Assemble_gradient.cpp:2591-2601 /release/4.0/finley/src/Assemble_gradient.cpp:5380-5406 /trunk/finley/src/Assemble_gradient.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26