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

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

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

revision 2988 by gross, Tue Mar 16 06:14:30 2010 UTC revision 2989 by gross, Thu Mar 18 01:45:55 2010 UTC
# Line 30  Line 30 
30    
31    
32  /* Disabled until the tests pass */  /* Disabled until the tests pass */
33  #define NEW_LUMPING  /* #define NEW_LUMPING */
34    
35  /**************************************************************/  /**************************************************************/
36    
# Line 50  void Finley_Assemble_LumpedSystem(Finley Line 50  void Finley_Assemble_LumpedSystem(Finley
50  #ifdef NEW_LUMPING  #ifdef NEW_LUMPING
51    register double m_t=0., diagS=0.;    register double m_t=0., diagS=0.;
52  #endif  #endif
53    
54    Finley_resetError();    Finley_resetError();
55    
56    if (nodes==NULL || elements==NULL) return;    if (nodes==NULL || elements==NULL) return;
# Line 77  void Finley_Assemble_LumpedSystem(Finley Line 77  void Finley_Assemble_LumpedSystem(Finley
77    /* set all parameters in p*/    /* set all parameters in p*/
78    Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);    Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79    if (! Finley_noError()) return;    if (! Finley_noError()) return;
80    
81    /* check if all function spaces are the same */    /* check if all function spaces are the same */
82      if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
   if (! numSamplesEqual(D,p.numQuadSub,elements->numElements) ) {  
83          sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);          sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
84          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
85    }    }
86    
87    /*  check the dimensions: */    /*  check the dimensions: */
     
88    if (p.numEqu==1) {    if (p.numEqu==1) {
89      if (!isEmpty(D)) {      if (!isEmpty(D)) {
90         if (!isDataPointShapeEqual(D,0,dimensions)) {         if (!isDataPointShapeEqual(D,0,dimensions)) {
# Line 103  void Finley_Assemble_LumpedSystem(Finley Line 101  void Finley_Assemble_LumpedSystem(Finley
101        }        }
102      }      }
103    }    }
   
104    if (Finley_noError()) {    if (Finley_noError()) {
105         printf("SADSADSADAS 2\n");
106      requireWrite(lumpedMat);      requireWrite(lumpedMat);
107      lumpedMat_p=getSampleDataRW(lumpedMat,0);      lumpedMat_p=getSampleDataRW(lumpedMat,0);
108      len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;      len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
# Line 112  void Finley_Assemble_LumpedSystem(Finley Line 110  void Finley_Assemble_LumpedSystem(Finley
110            
111      expandedD=isExpanded(D);      expandedD=isExpanded(D);
112      S=p.row_jac->BasisFunctions->S;      S=p.row_jac->BasisFunctions->S;
113    
114  #ifdef NEW_LUMPING  #ifdef NEW_LUMPING
115      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
116  #else  #else
# Line 130  void Finley_Assemble_LumpedSystem(Finley Line 128  void Finley_Assemble_LumpedSystem(Finley
128                      for(e=0;e<elements->numElements;e++){                                    for(e=0;e<elements->numElements;e++){              
129              if (elements->Color[e]==color) {              if (elements->Color[e]==color) {
130                              for (isub=0; isub<p.numSub; isub++) {                              for (isub=0; isub<p.numSub; isub++) {
131                                 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);;                                 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
                   
132                                 D_p=getSampleDataRO(D,e);                                                           D_p=getSampleDataRO(D,e);                          
133                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
134    
# Line 153  void Finley_Assemble_LumpedSystem(Finley Line 150  void Finley_Assemble_LumpedSystem(Finley
150                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
151                                         rtmp=0;                                         rtmp=0;
152                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
153                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
154                                     }                                     }
155    
156                                 #endif                                 #endif
157                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
158                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
# Line 163  void Finley_Assemble_LumpedSystem(Finley Line 161  void Finley_Assemble_LumpedSystem(Finley
161                      } /* end element loop */                      } /* end element loop */
162                    } /* end color loop */                    } /* end color loop */
163               } else  {  /* with constant D */                 } else  {  /* with constant D */  
164    printf("SADSADSADAS\n");
165                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
166                      /*  open loop over all elements: */                      /*  open loop over all elements: */
167                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
# Line 174  void Finley_Assemble_LumpedSystem(Finley Line 173  void Finley_Assemble_LumpedSystem(Finley
173                     #ifdef NEW_LUMPING /* HRZ lumping */                     #ifdef NEW_LUMPING /* HRZ lumping */
174                                     m_t=0; /* mass of the element: m_t */                                     m_t=0; /* mass of the element: m_t */
175                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];                                     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
                     
176                                     diagS=0; /* diagonal sum: S */                                     diagS=0; /* diagonal sum: S */
177                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
178                                        rtmp=0;                                        rtmp=0;
# Line 185  void Finley_Assemble_LumpedSystem(Finley Line 183  void Finley_Assemble_LumpedSystem(Finley
183                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */                                     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
184                     rtmp=m_t/diagS*D_p[0];                     rtmp=m_t/diagS*D_p[0];
185                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;                                     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
                             
186                                 #else /* row-sum lumping */                                 #else /* row-sum lumping */
187                                     for (s=0;s<p.row_numShapes;s++) {                                     for (s=0;s<p.row_numShapes;s++) {
188                                         rtmp=0;                                         rtmp=0;
189                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                         for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
190                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];                                         EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
191                                     }                                     }
192    for (s=0;s<p.row_numShapes;s++) printf(" %d %d : %e\n",isub,s,EM_lumpedMat[INDEX2(0,s,p.numEqu)]);
193                                 #endif                                 #endif
194                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];                                 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
195                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                                 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
# Line 251  void Finley_Assemble_LumpedSystem(Finley Line 249  void Finley_Assemble_LumpedSystem(Finley
249                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
250                 for (isub=0; isub<p.numSub; isub++) {                     for (isub=0; isub<p.numSub; isub++) {    
251                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);                                Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
                               memset(EM_lumpedMat,0,len_EM_lumpedMat_size);  
252                                D_p=getSampleDataRO(D,e);                                D_p=getSampleDataRO(D,e);
253                                
254                                #ifdef NEW_LUMPING /* HRZ lumping */                                #ifdef NEW_LUMPING /* HRZ lumping */
# Line 274  void Finley_Assemble_LumpedSystem(Finley Line 271  void Finley_Assemble_LumpedSystem(Finley
271                      for (k=0;k<p.numEqu;k++) {                      for (k=0;k<p.numEqu;k++) {
272                                          rtmp=0.;                                          rtmp=0.;
273                                          for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];                                          for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
274                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];                                          EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
275                                       }                                       }
276                   }                   }
277                    #endif                    #endif

Legend:
Removed from v.2988  
changed lines
  Added in v.2989

  ViewVC Help
Powered by ViewVC 1.1.26