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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 781 - (show annotations)
Fri Jul 14 08:47:38 2006 UTC (13 years, 3 months ago) by gross
File MIME type: text/plain
File size: 18035 byte(s)
grad functions linked into the persistent jacobean scheme
1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13 /**************************************************************/
14
15 /* assemblage rjacines: calculate the gradient of nodal data at quadrature points */
16
17 /**************************************************************/
18
19 /* Copyrights by ACcESS Australia, 2003,2004,2005 */
20 /* author: gross@access.edu.au */
21 /* version: $Id$ */
22
23 /**************************************************************/
24
25 #include "Assemble.h"
26 #include "Util.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30 /*****************************************************************/
31
32
33 void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,
34 escriptDataC* grad_data,escriptDataC* data) {
35
36 Finley_resetError();
37 if (nodes==NULL || elements==NULL) return;
38 dim_t numComps=getDataPointSize(data);
39 dim_t NN=elements->ReferenceElement->Type->numNodes;
40 dim_t numNodes, numShapes, numLocalNodes;
41 index_t dof_offset, s_offset;
42 type_t data_type=getFunctionSpaceType(data);
43 type_t grad_data_type=getFunctionSpaceType(grad_data);
44 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
45
46 if (data_type==FINLEY_NODES) {
47 reducedShapefunction=FALSE;
48 numNodes=nodes->numNodes;
49 numShapes=elements->ReferenceElement->Type->numShapes;
50 numLocalNodes=elements->ReferenceElement->Type->numNodes;
51 }
52 /* lock these two options jac for the MPI version */
53 #ifndef PASO_MPI
54 else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
55 reducedShapefunction=FALSE;
56 numNodes=nodes->numDegreesOfFreedom;
57 numShapes=elements->ReferenceElement->Type->numShapes;
58 numLocalNodes=elements->ReferenceElement->Type->numNodes;
59 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
60 reducedShapefunction=TRUE;
61 numNodes=nodes->reducedNumDegreesOfFreedom;
62 numShapes=elements->LinearReferenceElement->Type->numShapes;
63 numLocalNodes=elements->LinearReferenceElement->Type->numNodes;
64 }
65 #endif
66 else {
67 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
68 }
69 if (grad_data_type==FINLEY_ELEMENTS) {
70 reducedIntegrationOrder=FALSE;
71 dof_offset=0;
72 } else if (grad_data_type==FINLEY_FACE_ELEMENTS) {
73 reducedIntegrationOrder=FALSE;
74 dof_offset=0;
75 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_1) {
76 reducedIntegrationOrder=FALSE;
77 dof_offset=0;
78 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2) {
79 reducedIntegrationOrder=FALSE;
80 /* don't reset offset */
81 } else {
82 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculated for requested location.");
83 }
84
85
86 Finley_ElementFile_Jacobeans* jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
87
88 if (Finley_noError()) {
89
90 if (grad_data_type==FINLEY_ELEMENTS) {
91 dof_offset=0;
92 s_offset=0;
93 } else if (grad_data_type==FINLEY_FACE_ELEMENTS) {
94 dof_offset=0;
95 s_offset=0;
96 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_1) {
97 dof_offset=0;
98 s_offset=0;
99 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2) {
100 dof_offset=numShapes;
101 s_offset=jac->ReferenceElement->Type->numShapes;
102 }
103
104 /* check the dimensions of data */
105
106 if (! numSamplesEqual(grad_data,jac->ReferenceElement->numQuadNodes,elements->numElements)) {
107 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");
108 } else if (! numSamplesEqual(data,1,numNodes)) {
109 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");
110 } else if (jac->numDim*numComps!=getDataPointSize(grad_data)) {
111 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");
112 } else if (!isExpanded(grad_data)) {
113 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");
114 } else if (! (dof_offset+jac->ReferenceElement->Type->numShapes <= numLocalNodes)) {
115 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
116 } else if (! (s_offset+jac->ReferenceElement->Type->numShapes <= jac->ReferenceElement->Type->numNodes)) {
117 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: offset test failed.");
118 }
119 }
120 /* now we can start */
121
122 if (Finley_noError()) {
123 #pragma omp parallel
124 {
125 register dim_t e,q,l,s,n;
126 register double* data_array, *grad_data_e;
127 if (data_type==FINLEY_NODES) {
128 if (jac->numDim==1) {
129 #define DIM 1
130 #pragma omp for schedule(static)
131 for (e=0;e<elements->numElements;e++) {
132 grad_data_e=getSampleData(grad_data,e);
133 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
134 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
135 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
136 data_array=getSampleData(data,n);
137 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
138 for (l=0;l<numComps;l++) {
139 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
140 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
141 }
142 }
143 }
144 }
145
146 #undef DIM
147 } else if (jac->numDim==2) {
148 #define DIM 2
149 #pragma omp for schedule(static)
150 for (e=0;e<elements->numElements;e++) {
151 grad_data_e=getSampleData(grad_data,e);
152 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
153 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
154 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
155 data_array=getSampleData(data,n);
156 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
157 for (l=0;l<numComps;l++) {
158 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
159 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
160 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
161 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
162 }
163 }
164 }
165 }
166 #undef DIM
167 } else if (jac->numDim==3) {
168 #define DIM 3
169 #pragma omp for schedule(static)
170 for (e=0;e<elements->numElements;e++) {
171 grad_data_e=getSampleData(grad_data,e);
172 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
173 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
174 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
175 data_array=getSampleData(data,n);
176 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
177 for (l=0;l<numComps;l++) {
178 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
179 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
180 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
181 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
182 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
183 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
184 }
185 }
186 }
187 }
188 #undef DIM
189 }
190
191 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
192 if (jac->numDim==1) {
193 #define DIM 1
194 #pragma omp for schedule(static)
195 for (e=0;e<elements->numElements;e++) {
196 grad_data_e=getSampleData(grad_data,e);
197 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
198 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
199 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
200 data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
201 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
202 for (l=0;l<numComps;l++) {
203 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
204 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
205 }
206 }
207 }
208 }
209
210 #undef DIM
211 } else if (jac->numDim==2) {
212 #define DIM 2
213 #pragma omp for schedule(static)
214 for (e=0;e<elements->numElements;e++) {
215 grad_data_e=getSampleData(grad_data,e);
216 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
217 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
218 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
219 data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
220 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
221 for (l=0;l<numComps;l++) {
222 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
223 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
224 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
225 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
226 }
227 }
228 }
229 }
230 #undef DIM
231 } else if (jac->numDim==3) {
232 #define DIM 3
233 #pragma omp for schedule(static)
234 for (e=0;e<elements->numElements;e++) {
235 grad_data_e=getSampleData(grad_data,e);
236 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
237 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
238 n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
239 data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
240 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
241 for (l=0;l<numComps;l++) {
242 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
243 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
244 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
245 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
246 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
247 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
248 }
249 }
250 }
251 }
252 #undef DIM
253 }
254 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
255 if (jac->numDim==1) {
256 #define DIM 1
257 #pragma omp for schedule(static)
258 for (e=0;e<elements->numElements;e++) {
259 grad_data_e=getSampleData(grad_data,e);
260 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
261 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
262 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
263 data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
264 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
265 for (l=0;l<numComps;l++) {
266 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
267 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
268 }
269 }
270 }
271 }
272
273 #undef DIM
274 } else if (jac->numDim==2) {
275 #define DIM 2
276 #pragma omp for schedule(static)
277 for (e=0;e<elements->numElements;e++) {
278 grad_data_e=getSampleData(grad_data,e);
279 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
280 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
281 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
282 data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
283 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
284 for (l=0;l<numComps;l++) {
285 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
286 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
287 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
288 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
289 }
290 }
291 }
292 }
293
294 #undef DIM
295
296 } else if (jac->numDim==3) {
297 #define DIM 3
298 #pragma omp for schedule(static)
299 for (e=0;e<elements->numElements;e++) {
300 grad_data_e=getSampleData(grad_data,e);
301 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
302 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
303 n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
304 data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
305 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
306 for (l=0;l<numComps;l++) {
307 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
308 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
309 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
310 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
311 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
312 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
313 }
314 }
315 }
316 }
317
318 #undef DIM
319 }
320 }
321 } /* end parallel region */
322 }
323 }
324 /*
325 * $Log$
326 * Revision 1.6 2005/09/15 03:44:21 jgs
327 * Merge of development branch dev-02 back to main trunk on 2005-09-15
328 *
329 * Revision 1.5.2.1 2005/09/07 06:26:17 gross
330 * the solver from finley are put into the standalone package paso now
331 *
332 * Revision 1.5 2005/07/08 04:07:47 jgs
333 * Merge of development branch back to main trunk on 2005-07-08
334 *
335 * Revision 1.4 2004/12/15 07:08:32 jgs
336 * *** empty log message ***
337 * Revision 1.1.1.1.2.2 2005/06/29 02:34:48 gross
338 * some changes towards 64 integers in finley
339 *
340 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
341 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
342 *
343 *
344 *
345 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26