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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1064 - (hide annotations)
Tue Mar 27 06:21:02 2007 UTC (16 years ago) by gross
File MIME type: text/plain
File size: 17355 byte(s) 
test for reduced integration order for grad, interpolate and integrate added.
Bug shown by the tests have been fixed.
1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
13     /**************************************************************/
14    
15 gross 781 /* assemblage rjacines: calculate the gradient of nodal data at quadrature points */
16 jgs 82
17     /**************************************************************/
18    
19 jgs 150 /* Copyrights by ACcESS Australia, 2003,2004,2005 */
20 jgs 82 /* author: gross@access.edu.au */
21 jgs 150 /* version: $Id$ */
22 jgs 82
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 gross 781 escriptDataC* grad_data,escriptDataC* data) {
35 jgs 82
36 gross 1028 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
37 jgs 123 type_t data_type=getFunctionSpaceType(data);
38 gross 781 type_t grad_data_type=getFunctionSpaceType(grad_data);
39     bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 gross 1028 index_t dof_offset, s_offset;
41     Finley_ElementFile_Jacobeans* jac=NULL;
42     Finley_resetError();
43     if (nodes==NULL || elements==NULL) return;
44     numComps=getDataPointSize(data);
45     NN=elements->ReferenceElement->Type->numNodes;
46 gross 1062 reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
47 jgs 82
48     if (data_type==FINLEY_NODES) {
49 gross 781 reducedShapefunction=FALSE;
50 jgs 82 numNodes=nodes->numNodes;
51 gross 1062 } else if (data_type==FINLEY_REDUCED_NODES) { /* TODO */
52     reducedShapefunction=FALSE;
53     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: FINLEY_REDUCED_NODES is not supported yet.");
54 bcumming 751 }
55 gross 781 /* lock these two options jac for the MPI version */
56     #ifndef PASO_MPI
57     else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
58     reducedShapefunction=FALSE;
59 jgs 82 numNodes=nodes->numDegreesOfFreedom;
60     } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
61 gross 781 reducedShapefunction=TRUE;
62 jgs 82 numNodes=nodes->reducedNumDegreesOfFreedom;
63 gross 781 }
64     #endif
65     else {
66     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
67 jgs 82 }
68    
69 gross 1028 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
70 gross 781 if (Finley_noError()) {
71    
72 gross 1064 numShapes=jac->ReferenceElement->Type->numShapes;
73     numLocalNodes=jac->ReferenceElement->Type->numNodes;
74 gross 1062 if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2 || grad_data_type== FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
75     dof_offset=numShapes;
76     s_offset=jac->ReferenceElement->Type->numShapes;
77     } else {
78 gross 781 dof_offset=0;
79     s_offset=0;
80     }
81    
82     /* check the dimensions of data */
83    
84     if (! numSamplesEqual(grad_data,jac->ReferenceElement->numQuadNodes,elements->numElements)) {
85     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");
86     } else if (! numSamplesEqual(data,1,numNodes)) {
87     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");
88     } else if (jac->numDim*numComps!=getDataPointSize(grad_data)) {
89     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");
90     } else if (!isExpanded(grad_data)) {
91     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");
92     } else if (! (dof_offset+jac->ReferenceElement->Type->numShapes <= numLocalNodes)) {
93     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
94     } else if (! (s_offset+jac->ReferenceElement->Type->numShapes <= jac->ReferenceElement->Type->numNodes)) {
95     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: offset test failed.");
96     }
97 jgs 82 }
98     /* now we can start */
99    
100 jgs 150 if (Finley_noError()) {
101 gross 853 register dim_t e,q,l,s,n;
102     register double* data_array, *grad_data_e;
103     #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
104 gross 781 {
105     if (data_type==FINLEY_NODES) {
106     if (jac->numDim==1) {
107     #define DIM 1
108     #pragma omp for schedule(static)
109     for (e=0;e<elements->numElements;e++) {
110     grad_data_e=getSampleData(grad_data,e);
111     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
112     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
113     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
114     data_array=getSampleData(data,n);
115     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
116     for (l=0;l<numComps;l++) {
117     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
118     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
119     }
120     }
121     }
122     }
123 jgs 82
124 gross 781 #undef DIM
125     } else if (jac->numDim==2) {
126     #define DIM 2
127     #pragma omp for schedule(static)
128     for (e=0;e<elements->numElements;e++) {
129     grad_data_e=getSampleData(grad_data,e);
130     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
131     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
132     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
133     data_array=getSampleData(data,n);
134     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
135     for (l=0;l<numComps;l++) {
136     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
137     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
138     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
139     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
140     }
141     }
142     }
143     }
144     #undef DIM
145     } else if (jac->numDim==3) {
146     #define DIM 3
147 gross 853 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
148 gross 781 for (e=0;e<elements->numElements;e++) {
149 gross 853 grad_data_e=getSampleData(grad_data,e);
150 gross 781 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
151     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
152     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
153     data_array=getSampleData(data,n);
154     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
155     for (l=0;l<numComps;l++) {
156 gross 853 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
157 gross 781 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
158 gross 853 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
159 gross 781 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
160 gross 853 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
161 gross 781 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
162     }
163     }
164     }
165     }
166     #undef DIM
167     }
168 gross 532
169 gross 781 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
170     if (jac->numDim==1) {
171     #define DIM 1
172     #pragma omp for schedule(static)
173     for (e=0;e<elements->numElements;e++) {
174     grad_data_e=getSampleData(grad_data,e);
175     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
176     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
177     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
178     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
179     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
180     for (l=0;l<numComps;l++) {
181     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
182     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
183     }
184     }
185     }
186     }
187    
188     #undef DIM
189     } else if (jac->numDim==2) {
190     #define DIM 2
191     #pragma omp for schedule(static)
192     for (e=0;e<elements->numElements;e++) {
193     grad_data_e=getSampleData(grad_data,e);
194     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
195     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
196     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
197     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
198     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
199     for (l=0;l<numComps;l++) {
200     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
201     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
202     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
203     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
204     }
205     }
206     }
207     }
208     #undef DIM
209     } else if (jac->numDim==3) {
210     #define DIM 3
211     #pragma omp for schedule(static)
212     for (e=0;e<elements->numElements;e++) {
213     grad_data_e=getSampleData(grad_data,e);
214     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
215     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
216     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
217     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
218     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
219     for (l=0;l<numComps;l++) {
220     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
221     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
222     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
223     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
224     grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
225     jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
226     }
227     }
228     }
229     }
230     #undef DIM
231     }
232     } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
233     if (jac->numDim==1) {
234     #define DIM 1
235     #pragma omp for schedule(static)
236     for (e=0;e<elements->numElements;e++) {
237     grad_data_e=getSampleData(grad_data,e);
238     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
239     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
240     n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
241     data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
242     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
243     for (l=0;l<numComps;l++) {
244     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
245     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
246     }
247     }
248     }
249     }
250    
251     #undef DIM
252     } else if (jac->numDim==2) {
253     #define DIM 2
254     #pragma omp for schedule(static)
255     for (e=0;e<elements->numElements;e++) {
256     grad_data_e=getSampleData(grad_data,e);
257     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
258     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
259     n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
260     data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
261     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
262     for (l=0;l<numComps;l++) {
263     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
264     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
265     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
266     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
267     }
268     }
269     }
270     }
271    
272     #undef DIM
273    
274     } else if (jac->numDim==3) {
275     #define DIM 3
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     grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
290     jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
291     }
292     }
293     }
294     }
295     #undef DIM
296     }
297     }
298     } /* end parallel region */
299 jgs 82 }
300     }
301     /*
302     * $Log$
303 jgs 150 * Revision 1.6 2005/09/15 03:44:21 jgs
304     * Merge of development branch dev-02 back to main trunk on 2005-09-15
305     *
306     * Revision 1.5.2.1 2005/09/07 06:26:17 gross
307     * the solver from finley are put into the standalone package paso now
308     *
309 jgs 123 * Revision 1.5 2005/07/08 04:07:47 jgs
310     * Merge of development branch back to main trunk on 2005-07-08
311     *
312 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
313 jgs 97 * *** empty log message ***
314 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:48 gross
315     * some changes towards 64 integers in finley
316 jgs 82 *
317 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
318     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
319 jgs 97 *
320 jgs 82 *
321 jgs 123 *
322 jgs 82 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26