/[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

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

Legend:
Removed from v.3144  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26