/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp

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

trunk/finley/src/Assemble_PDE_Single2_1D.c revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.c revision 4257 by jfenwick, Wed Feb 27 03:42:40 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  /**************************************************************/  /************************************************************************************/
17    
18  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */
19  /*    the shape functions for test and solution must be identical */  /*    the shape functions for test and solution must be identical */
20    
   
21  /*      -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m  and -(X_,i)_i + Y */  /*      -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m  and -(X_,i)_i + Y */
22    
23  /*    in a 1D domain. The shape functions for test and solution must be identical  */  /*    in a 1D domain. The shape functions for test and solution must be identical  */
# Line 32  Line 32 
32  /*      X = 1  */  /*      X = 1  */
33  /*      Y = scalar   */  /*      Y = scalar   */
34    
35    /************************************************************************************/
 /**************************************************************/  
   
36    
37  #include "Assemble.h"  #include "Assemble.h"
38  #include "Util.h"  #include "Util.h"
# Line 42  Line 40 
40  #include <omp.h>  #include <omp.h>
41  #endif  #endif
42    
43    /************************************************************************************/
44    
45  /**************************************************************/  void Dudley_Assemble_PDE_Single2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46                        Paso_SystemMatrix * Mat, escriptDataC * F,
47  void  Finley_Assemble_PDE_Single2_1D(Assemble_Parameters p, Finley_ElementFile* elements,                      escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48                                       Paso_SystemMatrix* Mat, escriptDataC* F,                      escriptDataC * X, escriptDataC * Y)
49                                       escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {  {
50    
51      #define DIM 1  #define DIM 1
52      index_t color;      index_t color;
53      dim_t e;      dim_t e;
54      double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;      __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
55        double *EM_S, *EM_F, *DSDX;
56      index_t *row_index;      index_t *row_index;
57      register dim_t q, s,r;      register dim_t q, s, r;
58      register double rtmp;      register double rtmp;
59      bool_t add_EM_F, add_EM_S;      bool_t add_EM_F, add_EM_S;
60    
61      bool_t extendedA=isExpanded(A);      bool_t extendedA = isExpanded(A);
62      bool_t extendedB=isExpanded(B);      bool_t extendedB = isExpanded(B);
63      bool_t extendedC=isExpanded(C);      bool_t extendedC = isExpanded(C);
64      bool_t extendedD=isExpanded(D);      bool_t extendedD = isExpanded(D);
65      bool_t extendedX=isExpanded(X);      bool_t extendedX = isExpanded(X);
66      bool_t extendedY=isExpanded(Y);      bool_t extendedY = isExpanded(Y);
67      double *F_p=getSampleData(F,0);      double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68      double *S=p.row_jac->ReferenceElement->S;      double *S = p.row_jac->BasisFunctions->S;
69      dim_t len_EM_S=p.row_NN*p.col_NN;      dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal;
70      dim_t len_EM_F=p.row_NN;      dim_t len_EM_F = p.row_numShapesTotal;
71    
72    #pragma omp parallel private(color, EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q, row_index, q, s,r,rtmp,add_EM_F, add_EM_S)
     #pragma omp parallel private(color, EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, row_index, q, s,r,rtmp,add_EM_F, add_EM_S)  
73      {      {
74         EM_S=THREAD_MEMALLOC(len_EM_S,double);      EM_S = THREAD_MEMALLOC(len_EM_S, double);
75         EM_F=THREAD_MEMALLOC(len_EM_F,double);      EM_F = THREAD_MEMALLOC(len_EM_F, double);
76         row_index=THREAD_MEMALLOC(p.row_NN,index_t);      row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
77    
78         if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {      if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79        {
80            for (color=elements->minColor;color<=elements->maxColor;color++) {  
81               /*  open loop over all elements: */          for (color = elements->minColor; color <= elements->maxColor; color++)
82               #pragma omp for private(e) schedule(static)          {
83               for(e=0;e<elements->numElements;e++){          /*  open loop over all elements: */
84                  if (elements->Color[e]==color) {  #pragma omp for private(e) schedule(static)
85                     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);          for (e = 0; e < elements->numElements; e++)
86                     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);          {
87                     for (q=0;q<len_EM_S;++q) EM_S[q]=0;              if (elements->Color[e] == color)
88                     for (q=0;q<len_EM_F;++q) EM_F[q]=0;              {
89                     add_EM_F=FALSE;  
90                     add_EM_S=FALSE;              A_p = getSampleDataRO(A, e);
91                     /**************************************************************/              C_p = getSampleDataRO(C, e);
92                     /*   process A: */              B_p = getSampleDataRO(B, e);
93                     /**************************************************************/              D_p = getSampleDataRO(D, e);
94                     A_p=getSampleData(A,e);              X_p = getSampleDataRO(X, e);
95                     if (NULL!=A_p) {              Y_p = getSampleDataRO(Y, e);
96                        add_EM_S=TRUE;  
97                        if (extendedA) {              double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
98                           for (s=0;s<p.row_NS;s++) {              DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
99                             for (r=0;r<p.col_NS;r++) {              for (q = 0; q < len_EM_S; ++q)
100                                rtmp=0;                  EM_S[q] = 0;
101                                for (q=0;q<p.numQuad;q++) {              for (q = 0; q < len_EM_F; ++q)
102                                   rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                  EM_F[q] = 0;
103                               }              add_EM_F = FALSE;
104                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;              add_EM_S = FALSE;
105                             }                /************************************************************************************/
106                           }              /*   process A: */
107                        } else {                /************************************************************************************/
108                           for (s=0;s<p.row_NS;s++) {              if (NULL != A_p)
109                             for (r=0;r<p.col_NS;r++) {              {
110                                 rtmp=0;                  add_EM_S = TRUE;
111                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                  if (extendedA)
112                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX2(0,0,DIM)];                  {
113                             }                  A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
114                           }                  for (s = 0; s < p.row_numShapes; s++)
115                         }                  {
116                     }                      for (r = 0; r < p.row_numShapes; r++)
117                     /**************************************************************/                      {
118                     /*   process B: */                      rtmp = 0;
119                     /**************************************************************/                      for (q = 0; q < p.numQuadTotal; q++)
120                     B_p=getSampleData(B,e);                      {
121                     if (NULL!=B_p) {                          rtmp +=
122                        add_EM_S=TRUE;                          vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
123                        if (extendedB) {                          A_q[INDEX3(0, 0, q, DIM, DIM)] *
124                           for (s=0;s<p.row_NS;s++) {                          DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
125                             for (r=0;r<p.col_NS;r++) {                      }
126                               rtmp=0;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
127                               for (q=0;q<p.numQuad;q++) {                      }
128                                  rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX2(0,q,DIM)]*S[INDEX2(r,q,p.row_NS)];                  }
129                               }                  }
130                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                  else
131                             }                  {
132                           }                  for (s = 0; s < p.row_numShapes; s++)
133                        } else {                  {
134                           for (s=0;s<p.row_NS;s++) {                      for (r = 0; r < p.row_numShapes; r++)
135                             for (r=0;r<p.col_NS;r++) {                      {
136                                 rtmp=0;                      rtmp = 0;
137                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*S[INDEX2(r,q,p.row_NS)];                      for (q = 0; q < p.numQuadTotal; q++)
138                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[0];                          rtmp +=
139                             }                          vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
140                           }                          DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
141                        }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
142                     }                          rtmp * A_p[INDEX2(0, 0, DIM)];
143                     /**************************************************************/                      }
144                     /*   process C: */                  }
145                     /**************************************************************/                  }
146                     C_p=getSampleData(C,e);              }
147                     if (NULL!=C_p) {                /************************************************************************************/
148                        add_EM_S=TRUE;              /*   process B: */
149                       if (extendedC) {                /************************************************************************************/
150                           for (s=0;s<p.row_NS;s++) {              if (NULL != B_p)
151                             for (r=0;r<p.col_NS;r++) {              {
152                               rtmp=0;                  add_EM_S = TRUE;
153                               for (q=0;q<p.numQuad;q++) {                  if (extendedB)
154                                  rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                  {
155                               }                  B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
156                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                  for (s = 0; s < p.row_numShapes; s++)
157                             }                  {
158                           }                      for (r = 0; r < p.row_numShapes; r++)
159                       } else {                      {
160                           for (s=0;s<p.row_NS;s++) {                      rtmp = 0;
161                             for (r=0;r<p.col_NS;r++) {                      for (q = 0; q < p.numQuadTotal; q++)
162                                rtmp=0;                      {
163                                for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];                          rtmp +=
164                                EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[0];                          vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
165                             }                          B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
166                           }                      }
167                       }                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
168                     }                      }
169                     /************************************************************* */                  }
170                     /* process D */                  }
171                     /**************************************************************/                  else
172                     D_p=getSampleData(D,e);                  {
173                     if (NULL!=D_p) {                  for (s = 0; s < p.row_numShapes; s++)
174                       add_EM_S=TRUE;                  {
175                       if (extendedD) {                      for (r = 0; r < p.row_numShapes; r++)
176                           for (s=0;s<p.row_NS;s++) {                      {
177                             for (r=0;r<p.col_NS;r++) {                      rtmp = 0;
178                                rtmp=0;                      for (q = 0; q < p.numQuadTotal; q++)
179                                for (q=0;q<p.numQuad;q++) {                          rtmp +=
180                                   rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q]*S[INDEX2(r,q,p.row_NS)];                          vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
181                               }                          S[INDEX2(r, q, p.row_numShapes)];
182                               EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
183                             }                          rtmp * B_p[0];
184                           }                      }
185                       } else {                  }
186                           for (s=0;s<p.row_NS;s++) {                  }
187                             for (r=0;r<p.col_NS;r++) {              }
188                                 rtmp=0;                /************************************************************************************/
189                                 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];              /*   process C: */
190                                 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[0];                /************************************************************************************/
191                             }              if (NULL != C_p)
192                           }              {
193                       }                  add_EM_S = TRUE;
194                     }                  if (extendedC)
195                     /**************************************************************/                  {
196                     /*   process X: */                  C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
197                     /**************************************************************/                  for (s = 0; s < p.row_numShapes; s++)
198                     X_p=getSampleData(X,e);                  {
199                     if (NULL!=X_p) {                      for (r = 0; r < p.row_numShapes; r++)
200                       add_EM_F=TRUE;                      {
201                       if (extendedX) {                      rtmp = 0;
202                          for (s=0;s<p.row_NS;s++) {                      for (q = 0; q < p.numQuadTotal; q++)
203                            rtmp=0;                      {
204                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX2(0,q,DIM)];                          rtmp +=
205                            EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;                          vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
206                          }                          DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
207                       } else {                      }
208                          for (s=0;s<p.row_NS;s++) {                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
209                            rtmp=0;                      }
210                            for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];                  }
211                            EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*X_p[0];                  }
212                          }                  else
213                       }                  {
214                    }                  for (s = 0; s < p.row_numShapes; s++)
215                    /**************************************************************/                  {
216                    /*   process Y: */                      for (r = 0; r < p.row_numShapes; r++)
217                    /**************************************************************/                      {
218                     Y_p=getSampleData(Y,e);                      rtmp = 0;
219                     if (NULL!=Y_p) {                      for (q = 0; q < p.numQuadTotal; q++)
220                       add_EM_F=TRUE;                          rtmp +=
221                       if (extendedY) {                          vol * S[INDEX2(s, q, p.row_numShapes)] *
222                          for (s=0;s<p.row_NS;s++) {                          DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
223                             rtmp=0;                      EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
224                             for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[q];                          rtmp * C_p[0];
225                             EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;                      }
226                          }                  }
227                        } else {                  }
228                          for (s=0;s<p.row_NS;s++) {              }
229                              rtmp=0;                /*********************************************************************************** */
230                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];              /* process D */
231                              EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];                /************************************************************************************/
232                          }              if (NULL != D_p)
233                        }              {
234                      }                  add_EM_S = TRUE;
235                      /***********************************************************************************************/                  if (extendedD)
236                      /* add the element matrices onto the matrix and right hand side                                */                  {
237                      /***********************************************************************************************/                  D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
238                      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 (s = 0; s < p.row_numShapes; s++)
239                      if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);                  {
240                      if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);                      for (r = 0; r < p.row_numShapes; r++)
241                          {
242                  } /* end color check */                      rtmp = 0;
243               } /* end element loop */                      for (q = 0; q < p.numQuadTotal; q++)
244           } /* end color loop */                      {
245                                      rtmp +=
246           THREAD_MEMFREE(EM_S);                          vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
247           THREAD_MEMFREE(EM_F);                          S[INDEX2(r, q, p.row_numShapes)];
248           THREAD_MEMFREE(row_index);                      }
249                        EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
250                        }
251                    }
252                    }
253                    else
254                    {
255                    for (s = 0; s < p.row_numShapes; s++)
256                    {
257                        for (r = 0; r < p.row_numShapes; r++)
258                        {
259                        rtmp = 0;
260                        for (q = 0; q < p.numQuadTotal; q++)
261                            rtmp +=
262                            vol * S[INDEX2(s, q, p.row_numShapes)] *
263                            S[INDEX2(r, q, p.row_numShapes)];
264                        EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
265                            rtmp * D_p[0];
266                        }
267                    }
268                    }
269                }
270                  /************************************************************************************/
271                /*   process X: */
272                  /************************************************************************************/
273                if (NULL != X_p)
274                {
275                    add_EM_F = TRUE;
276                    if (extendedX)
277                    {
278                    X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
279                    for (s = 0; s < p.row_numShapes; s++)
280                    {
281                        rtmp = 0;
282                        for (q = 0; q < p.numQuadTotal; q++)
283                        rtmp +=
284                            vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
285                            X_q[INDEX2(0, q, DIM)];
286                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
287                    }
288                    }
289                    else
290                    {
291                    for (s = 0; s < p.row_numShapes; s++)
292                    {
293                        rtmp = 0;
294                        for (q = 0; q < p.numQuadTotal; q++)
295                        rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
296                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
297                    }
298                    }
299                }
300                 /************************************************************************************/
301                /*   process Y: */
302                 /************************************************************************************/
303                if (NULL != Y_p)
304                {
305                    add_EM_F = TRUE;
306                    if (extendedY)
307                    {
308                    Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
309                    for (s = 0; s < p.row_numShapes; s++)
310                    {
311                        rtmp = 0;
312                        for (q = 0; q < p.numQuadTotal; q++)
313                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
314                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
315                    }
316                    }
317                    else
318                    {
319                    for (s = 0; s < p.row_numShapes; s++)
320                    {
321                        rtmp = 0;
322                        for (q = 0; q < p.numQuadTotal; q++)
323                        rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
324                        EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
325                    }
326                    }
327                }
328                   /*********************************************************************************************************************/
329                /* add the element matrices onto the matrix and right hand side                                */
330                   /*********************************************************************************************************************/
331                for (q = 0; q < p.row_numShapesTotal; q++)
332                    row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
333    
334                if (add_EM_F)
335                    Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
336                               p.row_DOF_UpperBound);
337                if (add_EM_S)
338                    Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
339                                      p.row_numShapesTotal, row_index, p.numComp, EM_S);
340    
341                }       /* end color check */
342            }       /* end element loop */
343            }           /* end color loop */
344    
345            THREAD_MEMFREE(EM_S);   /* these FREEs appear to be inside the if because if any of the allocs */
346            THREAD_MEMFREE(EM_F);   /* failed it means an out of memory (which is not recoverable anyway) */
347            THREAD_MEMFREE(row_index);
348    
349        } /* end of pointer check */      }           /* end of pointer check */
350     } /* end parallel region */      }               /* end parallel region */
351  }  }
352    
353  /*  /*
354   * $Log$   * $Log$
355   */   */

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

  ViewVC Help
Powered by ViewVC 1.1.26