/[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 1220 by btully, Thu Aug 2 04:46:20 2007 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006,2007 by ACcENULLNULL MNRF              *  *
4   *                                                          *  * Copyright (c) 2003-2008 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open NULLoftware License version 3.0    *  *
8   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Primary Business: Queensland, Australia
9   *                                                          *  * Licensed under the Open Software License version 3.0
10   ************************************************************  * http://www.opensource.org/licenses/osl-3.0.php
11  */  *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
# Line 21  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$ */  
   
 /**************************************************************/  
   
25  #include "Assemble.h"  #include "Assemble.h"
26  #include "Util.h"  #include "Util.h"
27  #ifdef _OPENMP  #ifdef _OPENMP
# Line 33  Line 29 
29  #endif  #endif
30    
31    
32  #define NEW_LUMPING /* */  /* Disabled until the tests pass */
33    /* #define NEW_LUMPING */ /* */
34    
35  /**************************************************************/  /**************************************************************/
36    
# Line 43  void Finley_Assemble_LumpedSystem(Finley Line 40  void Finley_Assemble_LumpedSystem(Finley
40    bool_t reducedIntegrationOrder=FALSE, expandedD;    bool_t reducedIntegrationOrder=FALSE, expandedD;
41    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
42    Assemble_Parameters p;    Assemble_Parameters p;
   double time0;  
43    dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;    dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
44    type_t funcspace;    type_t funcspace;
45    index_t color,*row_index=NULL;    index_t color,*row_index=NULL;
46    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;
47    register double rtmp, m_t, diagS;    register double rtmp;
48    size_t len_EM_lumpedMat_size;    size_t len_EM_lumpedMat_size;
49    
50    #ifdef NEW_LUMPING
51      register double m_t, diagS;
52    #endif
53    
54    Finley_resetError();    Finley_resetError();
55    
56    if (nodes==NULL || elements==NULL) return;    if (nodes==NULL || elements==NULL) return;
# Line 111  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, m_t, diagS)      #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
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);
118         if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {         if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
119            if (p.numEqu == 1) {            if (p.numEqu == 1) {
120               if (expandedD) {               if (expandedD) {
                  #ifndef PASO_MPI  
121                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
122                      /*  open loop over all elements: */                      /*  open loop over all elements: */
123                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
124                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
125                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
                  #else  
                  {  
                     for(e=0;e<elements->numElements;e++) {  
                        {  
                  #endif  
126                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
127                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
128                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
# Line 137  void Finley_Assemble_LumpedSystem(Finley Line 131  void Finley_Assemble_LumpedSystem(Finley
131                             *           Number of PDEs: 1                             *           Number of PDEs: 1
132                             *  D_p varies over element: True                             *  D_p varies over element: True
133                             */                             */
134                              #pragma omp parallel private(m_t)
135                            m_t=0; /* mass of the element: m_t */                            m_t=0; /* mass of the element: m_t */
136                            for (q=0;q<p.numQuad;q++) {                            for (q=0;q<p.numQuad;q++) {
137                                m_t+=Vol[q]*D_p[q];                                m_t+=Vol[q]*D_p[q];
138                            }                                                      }                          
139                              #pragma omp parallel private(diagS)
140                            diagS=0; /* diagonal sum: S */                            diagS=0; /* diagonal sum: S */
141                            for (s=0;s<p.row_NS;s++) {                            for (s=0;s<p.row_NS;s++) {
142                                rtmp=0;                                rtmp=0;
# Line 167  void Finley_Assemble_LumpedSystem(Finley Line 163  void Finley_Assemble_LumpedSystem(Finley
163                      } /* end element loop */                      } /* end element loop */
164                    } /* end color loop */                    } /* end color loop */
165               } else  {               } else  {
                  #ifndef PASO_MPI  
166                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
167                   /*  open loop over all elements: */                   /*  open loop over all elements: */
168                   #pragma omp for private(e) schedule(static)                   #pragma omp for private(e) schedule(static)
169                   for(e=0;e<elements->numElements;e++){                   for(e=0;e<elements->numElements;e++){
170                      if (elements->Color[e]==color) {                      if (elements->Color[e]==color) {
                  #else  
                  {  
                     for(e=0;e<elements->numElements;e++) {  
                        {  
                  #endif  
171                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
172                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
173                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
# Line 218  void Finley_Assemble_LumpedSystem(Finley Line 208  void Finley_Assemble_LumpedSystem(Finley
208               }               }
209            } else {            } else {
210               if (expandedD) {               if (expandedD) {
                  #ifndef PASO_MPI  
211                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
212                      /*  open loop over all elements: */                      /*  open loop over all elements: */
213                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
214                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
215                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
                  #else  
                  {  
                     for(e=0;e<elements->numElements;e++) {  
                        {  
                  #endif  
216                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
217                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
218                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
# Line 271  void Finley_Assemble_LumpedSystem(Finley Line 255  void Finley_Assemble_LumpedSystem(Finley
255                      } /* end element loop */                      } /* end element loop */
256                  } /* end color loop */                  } /* end color loop */
257               } else {               } else {
258                   #ifndef PASO_MPI                   /*  open loop over all elements: */
259                   for (color=elements->minColor;color<=elements->maxColor;color++) {                   for (color=elements->minColor;color<=elements->maxColor;color++) {
                     /*  open loop over all elements: */  
260                      #pragma omp for private(e) schedule(static)                      #pragma omp for private(e) schedule(static)
261                      for(e=0;e<elements->numElements;e++){                      for(e=0;e<elements->numElements;e++){
262                         if (elements->Color[e]==color) {                         if (elements->Color[e]==color) {
                  #else  
                  {  
                     for(e=0;e<elements->numElements;e++) {  
                        {  
                  #endif  
263                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);                            Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
264                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);                            memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
265                            D_p=getSampleData(D,e);                            D_p=getSampleData(D,e);
# Line 328  void Finley_Assemble_LumpedSystem(Finley Line 306  void Finley_Assemble_LumpedSystem(Finley
306         THREAD_MEMFREE(row_index);         THREAD_MEMFREE(row_index);
307      } /* end parallel region */      } /* end parallel region */
308    }    }
   #ifdef Finley_TRACE  
   printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);  
   #endif  
309  }  }

Legend:
Removed from v.1220  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26