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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 780 by bcumming, Mon Jun 26 01:46:34 2006 UTC revision 781 by gross, Fri Jul 14 08:47:38 2006 UTC
# Line 12  Line 12 
12    
13  /**************************************************************/  /**************************************************************/
14    
15  /*    assemblage routines: calculate the gradient of nodal data at quadrature points */  /*    assemblage rjacines: calculate the gradient of nodal data at quadrature points */
16    
17  /**************************************************************/  /**************************************************************/
18    
# Line 30  Line 30 
30  /*****************************************************************/  /*****************************************************************/
31    
32    
 #define NODES 0  
 #define DOF 1  
 #define REDUCED_DOF 2  
   
33  void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,  void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,
34               escriptDataC* grad_data,escriptDataC* data) {                                escriptDataC* grad_data,escriptDataC* data) {
35    
36    double *local_X=NULL, *local_data=NULL, *dVdv=NULL, *dvdV=NULL, *Vol=NULL, *d_datadv=NULL, *gradS=NULL,*data_array;    Finley_resetError();
   index_t node_offset,*resort_nodes=FALSE,dof_offset;  
   dim_t numNodes=0,e,i,q,NS_DOF=0,NN_DOF=0;  
   type_t type=DOF;  
37    if (nodes==NULL || elements==NULL) return;    if (nodes==NULL || elements==NULL) return;
38      dim_t numComps=getDataPointSize(data);
39    dim_t NN=elements->ReferenceElement->Type->numNodes;    dim_t NN=elements->ReferenceElement->Type->numNodes;
40    dim_t NS=elements->ReferenceElement->Type->numShapes;    dim_t numNodes, numShapes, numLocalNodes;
41    index_t id[NN];    index_t dof_offset, s_offset;
   dim_t numDim=nodes->numDim;  
42    type_t data_type=getFunctionSpaceType(data);    type_t data_type=getFunctionSpaceType(data);
43    dim_t numComps=getDataPointSize(data);    type_t grad_data_type=getFunctionSpaceType(grad_data);
44    dim_t numQuad=elements->ReferenceElement->numQuadNodes;    bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
   for (i=0;i<NN;i++) id[i]=i;  
   Finley_resetError();  
   
     /* set some parameter */  
45    
46    if (data_type==FINLEY_NODES) {    if (data_type==FINLEY_NODES) {
47         type=NODES;         reducedShapefunction=FALSE;
        resort_nodes=id;  
        NN_DOF=elements->ReferenceElement->Type->numNodes;  
        NS_DOF=elements->ReferenceElement->Type->numShapes;  
        gradS=elements->ReferenceElement->dSdv;  
48         numNodes=nodes->numNodes;         numNodes=nodes->numNodes;
49           numShapes=elements->ReferenceElement->Type->numShapes;
50           numLocalNodes=elements->ReferenceElement->Type->numNodes;
51    }    }
52  /* lock these two options out for the MPI version */    /* lock these two options jac for the MPI version */
53  #ifndef PASO_MPI    #ifndef PASO_MPI
54   else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {    else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
55         type=DOF;         reducedShapefunction=FALSE;
        resort_nodes=id;  
        NN_DOF=elements->ReferenceElement->Type->numNodes;  
        NS_DOF=elements->ReferenceElement->Type->numShapes;  
        gradS=elements->ReferenceElement->dSdv;  
56         numNodes=nodes->numDegreesOfFreedom;         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) {    } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
60         type=REDUCED_DOF;         reducedShapefunction=TRUE;
        resort_nodes=elements->ReferenceElement->Type->linearNodes;  
        NN_DOF=elements->LinearReferenceElement->Type->numNodes;  
        NS_DOF=elements->LinearReferenceElement->Type->numShapes;  
        gradS=elements->LinearReferenceElement->dSdv;  
61         numNodes=nodes->reducedNumDegreesOfFreedom;         numNodes=nodes->reducedNumDegreesOfFreedom;
62     }         numShapes=elements->LinearReferenceElement->Type->numShapes;
63  #endif         numLocalNodes=elements->LinearReferenceElement->Type->numNodes;
64  else {    }
65         Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data");    #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 (getFunctionSpaceType(grad_data)==FINLEY_CONTACT_ELEMENTS_2) {    if (grad_data_type==FINLEY_ELEMENTS) {
70         node_offset=NN-NS;         reducedIntegrationOrder=FALSE;
71         dof_offset=NN_DOF-NS_DOF;         dof_offset=0;
72    } else {    } else if (grad_data_type==FINLEY_FACE_ELEMENTS)  {
73         node_offset=0;         reducedIntegrationOrder=FALSE;
74         dof_offset=0;         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    }    }
                                                                                                                                                           
   /* check the dimensions of interpolated_data and data */  
84    
85    if (numDim!=elements->ReferenceElement->Type->numDim) {  
86       Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Spatial and element dimension must match.");    Finley_ElementFile_Jacobeans* jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
87    } else if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) {  
88         Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");    if (Finley_noError()) {
89    } else if (! numSamplesEqual(data,1,numNodes)) {  
90         Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");        if (grad_data_type==FINLEY_ELEMENTS) {
91    } else if (numDim*numComps!=getDataPointSize(grad_data)) {          dof_offset=0;
92         Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");          s_offset=0;
93    }  else if (!isExpanded(grad_data)) {        } else if (grad_data_type==FINLEY_FACE_ELEMENTS)  {
94         Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");         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 */    /* now we can start */
121    
122    if (Finley_noError()) {    if (Finley_noError()) {
123            #pragma omp parallel private(local_X,local_data,dvdV,dVdv,Vol,d_datadv)        #pragma omp parallel
124            {        {
125                 local_X=local_data=dVdv=dvdV=Vol=d_datadv=NULL;           register dim_t e,q,l,s,n;
126            /* allocation of work arrays */           register double* data_array,  *grad_data_e;
127            local_X=THREAD_MEMALLOC(NS*numDim,double);           if (data_type==FINLEY_NODES) {
128            local_data=THREAD_MEMALLOC(NS*numComps,double);              if (jac->numDim==1) {
129            dVdv=THREAD_MEMALLOC(numQuad*numDim*numDim,double);                  #define DIM 1
130            dvdV=THREAD_MEMALLOC(numQuad*numDim*numDim,double);                  #pragma omp for schedule(static)
131            Vol=THREAD_MEMALLOC(numQuad,double);              for (e=0;e<elements->numElements;e++) {
132            d_datadv=THREAD_MEMALLOC(numQuad*numComps*numDim,double);                      grad_data_e=getSampleData(grad_data,e);
133            if (!(Finley_checkPtr(local_X) || Finley_checkPtr(dVdv) || Finley_checkPtr(dvdV) || Finley_checkPtr(Vol) || Finley_checkPtr(d_datadv) || Finley_checkPtr(local_data) ))  {                      for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
134              /* open the element loop */                      for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
135                  #pragma omp for private(e,i,q,data_array) schedule(static)                         n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
136              for(e=0;e<elements->numElements;e++) {                         data_array=getSampleData(data,n);
137                /* gather local coordinates of nodes into local_X: */                         for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
138                Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);                             for (l=0;l<numComps;l++) {
139                /*  calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */                                  grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
140                Finley_Util_SmallMatMult(numDim,numDim*numQuad,dVdv,NS,local_X,elements->ReferenceElement->dSdv);                                        jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
141                /*  dvdV=invert(dVdv) */                             }
142                Finley_Util_InvertSmallMat(numQuad,numDim,dVdv,dvdV,Vol);                         }
143                /* gather local data into local_data(numComps,NS_DOF): */                      }
144                    switch (type) {                  }
145                       case NODES:  
146                          for (q=0;q<NS_DOF;q++) {                  #undef DIM
147                             i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)];              } else if (jac->numDim==2) {
148                             data_array=getSampleData(data,i);                  #define DIM 2
149                             Finley_copyDouble(numComps,data_array,local_data+q*numComps);                  #pragma omp for schedule(static)
150                          }              for (e=0;e<elements->numElements;e++) {
151                          break;                      grad_data_e=getSampleData(grad_data,e);
152                       case DOF:                      for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
153                          for (q=0;q<NS_DOF;q++) {                      for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
154                             i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)];                         n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
155                             data_array=getSampleData(data,nodes->degreeOfFreedom[i]);                         data_array=getSampleData(data,n);
156                             Finley_copyDouble(numComps,data_array,local_data+q*numComps);                         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                          break;                                          jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
160                       case REDUCED_DOF:                                 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
161                          for (q=0;q<NS_DOF;q++) {                                          jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
162                             i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)];                             }
163                             data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[i]);                         }
164                             Finley_copyDouble(numComps,data_array,local_data+q*numComps);                      }
165                          }                  }
166                          break;                  #undef DIM
167                    }              } else if (jac->numDim==3) {
168                /*  calculate grad_data(l,i,q)=local_data(l,n)* DSDV(n,i,q) */                  #define DIM 3
169                // Finley_Util_SmallMatMult(numQuad,numComps,numQuad*numDim,getSampleData(grad_data,e),NS_DOF,local_data,dSdV);                  #pragma omp for schedule(static)
170                for (e=0;e<elements->numElements;e++) {
171                /*  calculate d_datadv(l,i,q)=local_data(l,n)*DSDv(n,i,q) */                      grad_data_e=getSampleData(grad_data,e);
172                Finley_Util_SmallMatMult(numComps,numDim*numQuad,d_datadv,NS_DOF,local_data,gradS);                      for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
173                /*  calculate grad_data(l,i,q)=d_datadv(l,k,q)*dvdV(k,i,q) */                      for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
174                Finley_Util_SmallMatSetMult(numQuad,numComps,numDim,getSampleData(grad_data,e),numDim,d_datadv,dvdV);                         n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
175              } /* for */                         data_array=getSampleData(data,n);
176            }                         for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
177            THREAD_MEMFREE(local_X);                             for (l=0;l<numComps;l++) {
178            THREAD_MEMFREE(dVdv);                                 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
179            THREAD_MEMFREE(dvdV);                                      jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
180            THREAD_MEMFREE(Vol);                                 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
181            THREAD_MEMFREE(local_data);                                      jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
182            THREAD_MEMFREE(d_datadv);                                 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  }  }
 #undef NODES  
 #undef DOF  
 #undef REDUCED_DOF  
324  /*  /*
325   * $Log$   * $Log$
326   * Revision 1.6  2005/09/15 03:44:21  jgs   * Revision 1.6  2005/09/15 03:44:21  jgs

Legend:
Removed from v.780  
changed lines
  Added in v.781

  ViewVC Help
Powered by ViewVC 1.1.26