/[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 1219 by gross, Sat Jun 23 11:43:12 2007 UTC revision 1220 by btully, Thu Aug 2 04:46:20 2007 UTC
# Line 22  Line 22 
22  /**************************************************************/  /**************************************************************/
23    
24  /*  Author: gross@access.edu.au */  /*  Author: gross@access.edu.au */
25  /*  Version: $Id:$ */  /*  Version: $Id$ */
26    
27  /**************************************************************/  /**************************************************************/
28    
# Line 33  Line 33 
33  #endif  #endif
34    
35    
36    #define NEW_LUMPING /* */
37    
38  /**************************************************************/  /**************************************************************/
39    
40  void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)  void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
# Line 46  void Finley_Assemble_LumpedSystem(Finley Line 48  void Finley_Assemble_LumpedSystem(Finley
48    type_t funcspace;    type_t funcspace;
49    index_t color,*row_index=NULL;    index_t color,*row_index=NULL;
50    double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;    double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
51    register double rtmp;    register double rtmp, m_t, diagS;
52    size_t len_EM_lumpedMat_size;    size_t len_EM_lumpedMat_size;
53    
54    Finley_resetError();    Finley_resetError();
# Line 109  void Finley_Assemble_LumpedSystem(Finley Line 111  void Finley_Assemble_LumpedSystem(Finley
111      expandedD=isExpanded(D);      expandedD=isExpanded(D);
112      S=p.row_jac->ReferenceElement->S;      S=p.row_jac->ReferenceElement->S;
113    
114      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
115      {      {
116         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);         EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
117         row_index=THREAD_MEMALLOC(p.row_NN,index_t);         row_index=THREAD_MEMALLOC(p.row_NN,index_t);
# Line 130  void Finley_Assemble_LumpedSystem(Finley Line 132  void Finley_Assemble_LumpedSystem(Finley
132                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
133                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
134                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
135                              #ifdef NEW_LUMPING /* HRZ lumping */
136                              /*
137                               *           Number of PDEs: 1
138                               *  D_p varies over element: True
139                               */
140                              m_t=0; /* mass of the element: m_t */
141                              for (q=0;q<p.numQuad;q++) {
142                                  m_t+=Vol[q]*D_p[q];
143                              }                          
144                              diagS=0; /* diagonal sum: S */
145                              for (s=0;s<p.row_NS;s++) {
146                                  rtmp=0;
147                                  for (q=0;q<p.numQuad;q++) {
148                                      rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
149                                  }
150                                  EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
151                                  diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
152                              }
153                              /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
154                              for (s=0;s<p.row_NS;s++) {
155                                  EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
156                              }
157                              #else /* row-sum lumping */
158                            for (s=0;s<p.row_NS;s++) {                            for (s=0;s<p.row_NS;s++) {
159                                rtmp=0;                                rtmp=0;
160                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
161                                EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;                                EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
162                            }                            }
163                              #endif
164                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
165                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
166                         } /* end color check */                         } /* end color check */
# Line 152  void Finley_Assemble_LumpedSystem(Finley Line 178  void Finley_Assemble_LumpedSystem(Finley
178                      for(e=0;e<elements->numElements;e++) {                      for(e=0;e<elements->numElements;e++) {
179                         {                         {
180                   #endif                   #endif
181                             Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
182                             memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
183                             D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
184                             for (s=0;s<p.row_NS;s++) {                            #ifdef NEW_LUMPING /* HRZ lumping */
185                                 rtmp=0;                            /*
186                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                             *           Number of PDEs: 1
187                                 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];                             *  D_p varies over element: False
188                             }                             */
189                             for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                            m_t=0; /* mass of the element: m_t */
190                             Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                            for (q=0;q<p.numQuad;q++) {
191                                  m_t+=Vol[q]*D_p[0];
192                              }                          
193                              diagS=0; /* diagonal sum: S */
194                              for (s=0;s<p.row_NS;s++) {
195                                  rtmp=0;
196                                  for (q=0;q<p.numQuad;q++) {
197                                      rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
198                                  }
199                                  EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
200                                  diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
201                              }
202                              /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
203                              for (s=0;s<p.row_NS;s++) {
204                                  EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
205                              }
206                              #else /* row-sum lumping */
207                              for (s=0;s<p.row_NS;s++) {
208                                  rtmp=0;
209                                  for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
210                                  EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
211                              }
212                              #endif
213                              for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
214                              Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
215                        } /* end color check */                        } /* end color check */
216                      } /* end element loop */                      } /* end element loop */
217                   } /* end color loop */                   } /* end color loop */
# Line 182  void Finley_Assemble_LumpedSystem(Finley Line 232  void Finley_Assemble_LumpedSystem(Finley
232                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
233                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
234                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
235                              #ifdef NEW_LUMPING /* HRZ lumping */
236                              /*
237                               *           Number of PDEs: Multiple
238                               *  D_p varies over element: True
239                               */
240                              for (k=0;k<p.numEqu;k++) {
241                                  m_t=0; /* mass of the element: m_t */
242                                  for (q=0;q<p.numQuad;q++) {
243                                      m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
244                                  }                          
245                                  diagS=0; /* diagonal sum: S */
246                                  for (s=0;s<p.row_NS;s++) {
247                                      rtmp=0;
248                                      for (q=0;q<p.numQuad;q++) {
249                                          rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
250                                      }
251                                      EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
252                                      diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
253                                  }
254                                  /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
255                                  for (s=0;s<p.row_NS;s++) {
256                                      EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
257                                  }
258                              }
259                              #else /* row-sum lumping */
260                            for (s=0;s<p.row_NS;s++) {                            for (s=0;s<p.row_NS;s++) {
261                                for (k=0;k<p.numEqu;k++) {                                for (k=0;k<p.numEqu;k++) {
262                                    rtmp=0.;                                    rtmp=0.;
# Line 189  void Finley_Assemble_LumpedSystem(Finley Line 264  void Finley_Assemble_LumpedSystem(Finley
264                                    EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;                                    EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
265                                }                                }
266                            }                            }
267                              #endif
268                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
269                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
270                         } /* end color check */                         } /* end color check */
# Line 209  void Finley_Assemble_LumpedSystem(Finley Line 285  void Finley_Assemble_LumpedSystem(Finley
285                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
286                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
287                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
288                              #ifdef NEW_LUMPING /* HRZ lumping */
289                              /*
290                               *           Number of PDEs: Multiple
291                               *  D_p varies over element: False
292                               */
293                              for (k=0;k<p.numEqu;k++) {
294                                  m_t=0; /* mass of the element: m_t */
295                                  for (q=0;q<p.numQuad;q++) {
296                                      m_t+=Vol[q]*D_p[k];
297                                  }                          
298                                  diagS=0; /* diagonal sum: S */
299                                  for (s=0;s<p.row_NS;s++) {
300                                      rtmp=0;
301                                      for (q=0;q<p.numQuad;q++) {
302                                          rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
303                                      }
304                                      EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
305                                      diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
306                                  }
307                                  /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
308                                  for (s=0;s<p.row_NS;s++) {
309                                      EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
310                                  }
311                              }
312                              #else /* row-sum lumping */
313                            for (s=0;s<p.row_NS;s++) {                            for (s=0;s<p.row_NS;s++) {
314                                rtmp=0;                                rtmp=0;
315                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
316                                for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];                                for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
317                            }                            }
318                              #endif
319                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                            for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
320                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);                            Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
321                         } /* end color check */                         } /* end color check */
Legend:
Removed from v.1219  
changed lines
  Added in v.1220
  ViewVC Help
Powered by ViewVC 1.1.26