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

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

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

revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 777 by gross, Wed Jul 12 08:54:45 2006 UTC
# Line 33  Line 33 
33  void Finley_Assemble_integrate(Finley_NodeFile* nodes, Finley_ElementFile* elements,escriptDataC* data,double* out) {  void Finley_Assemble_integrate(Finley_NodeFile* nodes, Finley_ElementFile* elements,escriptDataC* data,double* out) {
34      type_t datacase;      type_t datacase;
35      index_t node_offset;      index_t node_offset;
36        bool_t reducedIntegrationOrder;
37      dim_t q,e,i;      dim_t q,e,i;
     double *local_X, *dVdv, *Vol,*out_local,*data_array,rtmp;  
38      if (nodes==NULL || elements==NULL) return;      if (nodes==NULL || elements==NULL) return;
     dim_t NN=elements->ReferenceElement->Type->numNodes;  
     dim_t NS=elements->ReferenceElement->Type->numShapes;  
39      type_t data_type=getFunctionSpaceType(data);      type_t data_type=getFunctionSpaceType(data);
40      dim_t numComps=getDataPointSize(data);      dim_t numComps=getDataPointSize(data);
     dim_t numQuad=elements->ReferenceElement->numQuadNodes;  
     dim_t numDim=nodes->numDim;  
     dim_t numDim_local=elements->ReferenceElement->Type->numDim;  
41      Finley_resetError();      Finley_resetError();
42                                                                                                                                                                                                                                                                                                
43      /* set some parameter */      /* set some parameter */
44                                                                                                                                                                                                                                                                                                
45      if (data_type==FINLEY_ELEMENTS) {      if (data_type==FINLEY_ELEMENTS) {
46          datacase=0;          reducedIntegrationOrder=FALSE;
         node_offset=0;  
47      } else if (data_type==FINLEY_FACE_ELEMENTS)  {      } else if (data_type==FINLEY_FACE_ELEMENTS)  {
48          datacase=1;          reducedIntegrationOrder=FALSE;
         node_offset=0;  
49      } else if (data_type==FINLEY_CONTACT_ELEMENTS_1)  {      } else if (data_type==FINLEY_CONTACT_ELEMENTS_1)  {
50          datacase=1;          reducedIntegrationOrder=FALSE;
         node_offset=0;  
51      } else if (data_type==FINLEY_CONTACT_ELEMENTS_2)  {      } else if (data_type==FINLEY_CONTACT_ELEMENTS_2)  {
52          datacase=1;          reducedIntegrationOrder=FALSE;
         node_offset=NN-NS;  
53      } else {      } else {
54         Finley_setError(TYPE_ERROR,"__FILE__: integration of data is not possible.");         Finley_setError(TYPE_ERROR,"Finley_Assemble_integrate: integration of data is not possible.");
55      }      }
56    
57      /* check the shape of the data  */      /* get access to jacobean */
58        Finley_ElementFile_Jacobeans* jac=Finley_ElementFile_borrowJacobeans(elements,nodes,FALSE,reducedIntegrationOrder);
59    
60      if (! numSamplesEqual(data,numQuad,elements->numElements)) {      if (Finley_noError()) {
        Finley_setError(TYPE_ERROR,"__FILE__: illegal number of samples of integrant kernel Data object");  
     }  
61    
62      /* now we can start */          /* check the shape of the data  */
63            if (! numSamplesEqual(data,jac->ReferenceElement->numQuadNodes,elements->numElements)) {
64               Finley_setError(TYPE_ERROR,"Finley_Assemble_integrate: illegal number of samples of integrant kernel Data object");
65            }
66    
67      if (Finley_noError()) {  
68          for (q=0;q<numComps;q++) out[q]=0;          /* now we can start */
69          #pragma omp parallel private(local_X,dVdv,Vol,out_local,i)  
70          {          if (Finley_noError()) {
71                for (q=0;q<numComps;q++) out[q]=0;
72              /* allocation of work arrays */              #pragma omp parallel private(q,i)
73                {
74              local_X=dVdv=Vol=out_local=NULL;                  double out_local[numComps], rtmp;
75              out_local=THREAD_MEMALLOC(numComps,double);                  double* data_array;
76              dVdv=THREAD_MEMALLOC(numQuad*numDim_local*numDim,double);  
77              Vol=THREAD_MEMALLOC(numQuad,double);                  /* initialize local result */
78              local_X=THREAD_MEMALLOC(NS*numDim,double);  
79                for (i=0;i<numComps;i++) out_local[i]=0;
80              if (! (Finley_checkPtr(Vol) || Finley_checkPtr(local_X) || Finley_checkPtr(dVdv)  || Finley_checkPtr(out_local))) {  
81                    /* open the element loop */
82                 /* initialize local result */  
83                    if (isExpanded(data)) {
84             for (i=0;i<numComps;i++) out_local[i]=0;                      #pragma omp for private(e) schedule(static)
85                        for(e=0;e<elements->numElements;e++) {
86                 /* open the element loop */                           data_array=getSampleData(data,e);
87                             for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
88                 #pragma omp for private(e,q,data_array,rtmp) schedule(static)                                 for (i=0;i<numComps;i++) out_local[i]+=data_array[INDEX2(i,q,numComps)]*jac->volume[INDEX2(q,e,jac->ReferenceElement->numQuadNodes)];
                for(e=0;e<elements->numElements;e++) {  
                     /* gather local coordinates of nodes into local_X: */  
                     Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);  
                     /*  calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */  
                     Finley_Util_SmallMatMult(numDim,numDim_local*numQuad,dVdv,NS,local_X,elements->ReferenceElement->dSdv);  
                     /* calculate volume /area elements */  
                     switch (datacase) {  
                         case 0:  
                            Finley_Util_DetOfSmallMat(numQuad,numDim,dVdv,Vol);  
                            break;  
                         case 1:  
                            Finley_LengthOfNormalVector(numQuad,numDim,numDim_local,dVdv,Vol);  
                            break;  
                     }  
                     /* out=out+ integrals */  
                     data_array=getSampleData(data,e);  
                     if (isExpanded(data)) {  
                          for (q=0;q<numQuad;q++) {  
                             rtmp=ABS(Vol[q])*elements->ReferenceElement->QuadWeights[q];  
                             for (i=0;i<numComps;i++) out_local[i]+=data_array[INDEX2(i,q,numComps)]*rtmp;  
89                           }                           }
                     } else {  
                          rtmp=0.;  
                          for (q=0;q<numQuad;q++) rtmp+=ABS(Vol[q])*elements->ReferenceElement->QuadWeights[q];  
                          for (i=0;i<numComps;i++) out_local[i]+=data_array[i]*rtmp;  
90                      }                      }
91                 }                  } else {
92                       #pragma omp for private(e) schedule(static)
93                       for(e=0;e<elements->numElements;e++) {
94                            data_array=getSampleData(data,e);
95                            rtmp=0.;
96                            for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) rtmp+=jac->volume[INDEX2(q,e,jac->ReferenceElement->numQuadNodes)];
97                            for (i=0;i<numComps;i++) out_local[i]+=data_array[i]*rtmp;
98                       }
99                    }
100                    /* add local results to global result */
101                    #pragma omp critical
102                    for (i=0;i<numComps;i++) out[i]+=out_local[i];
103    
104                 /* add local results to global result */              }
                #pragma omp critical  
                for (i=0;i<numComps;i++) out[i]+=out_local[i];  
   
            }  
            THREAD_MEMFREE(local_X);  
            THREAD_MEMFREE(dVdv);  
            THREAD_MEMFREE(Vol);  
            THREAD_MEMFREE(out_local);  
105         }         }
106      }     }
107  }  }
108    
109  /*  /*

Legend:
Removed from v.616  
changed lines
  Added in v.777

  ViewVC Help
Powered by ViewVC 1.1.26