/[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 1064 - (show annotations)
Tue Mar 27 06:21:02 2007 UTC (12 years, 7 months 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 /*
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 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
37 type_t data_type=getFunctionSpaceType(data);
38 type_t grad_data_type=getFunctionSpaceType(grad_data);
39 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 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 reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
47
48 if (data_type==FINLEY_NODES) {
49 reducedShapefunction=FALSE;
50 numNodes=nodes->numNodes;
51 } 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 }
55 /* 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 numNodes=nodes->numDegreesOfFreedom;
60 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
61 reducedShapefunction=TRUE;
62 numNodes=nodes->reducedNumDegreesOfFreedom;
63 }
64 #endif
65 else {
66 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
67 }
68
69 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
70 if (Finley_noError()) {
71
72 numShapes=jac->ReferenceElement->Type->numShapes;
73 numLocalNodes=jac->ReferenceElement->Type->numNodes;
74 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 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 }
98 /* now we can start */
99
100 if (Finley_noError()) {
101 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 {
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
124 #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 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
148 for (e=0;e<elements->numElements;e++) {
149 grad_data_e=getSampleData(grad_data,e);
150 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 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
157 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
158 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
159 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
160 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
161 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
169 } 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 }
300 }
301 /*
302 * $Log$
303 * 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 * 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 * Revision 1.4 2004/12/15 07:08:32 jgs
313 * *** empty log message ***
314 * Revision 1.1.1.1.2.2 2005/06/29 02:34:48 gross
315 * some changes towards 64 integers in finley
316 *
317 * 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 *
320 *
321 *
322 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26