# Diff of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_Points.cpp

revision 6008 by caltinay, Wed Feb 17 23:53:30 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 15  Line 15
15  *****************************************************************************/  *****************************************************************************/
16
17
18  /************************************************************************************/  /****************************************************************************
19
20  /*    assembles the system of numEq PDEs into the stiffness matrix S right hand side F  */    Assembles the system of numEqu PDEs into the stiffness matrix S right hand
21  /*    the shape functions for test and solution must be identical */    side F
22
23          d_dirac_{k,m} u_m and y_dirac_k
24
25  /*      d_dirac_{k,m} u_m yand _dirac_k */    u has p.numComp components in a 3D domain.
26      The shape functions for test and solution must be identical and
27      row_NS == row_NN.
28
29  /*    u has p.numComp components in a 3D domain. The shape functions for test and solution must be identical  */    Shape of the coefficients:
/*    and row_NS == row_NN                                                                                  */
30
31  /*    Shape of the coefficients: */        d_dirac = p.numEqu x p.numComp
32          y_dirac = p.numEqu
33
/*      d_dirac = p.numEqu x p.numComp  */
/*      y_dirac = p.numEqu   */
34
35    *****************************************************************************/
/************************************************************************************/

#define ESNEEDPYTHON
#include "esysUtils/first.h"
36
37  #include "Assemble.h"  #include "Assemble.h"
38  #include "Util.h"  #include "Util.h"
#ifdef _OPENMP
#include <omp.h>
#endif
39
40    namespace dudley {
41
42  /************************************************************************************/  void Assemble_PDE_Points(const Assemble_Parameters& p,
43                             const Dudley_ElementFile* elements,
44                             escript::ASM_ptr mat, escript::Data& F,
45                             const escript::Data& d_dirac,
46                             const escript::Data& y_dirac)
47    {
48        double* F_p = NULL;
49        if (!F.isEmpty()) {
50            F.requireWrite();
51            F_p = F.getSampleDataRW(0);
52        }
53
54  void  Dudley_Assemble_PDE_Points(Dudley_Assemble_Parameters p,  #pragma omp parallel
Dudley_ElementFile* elements,
escript::ASM_ptr mat, escript::Data* F,
const escript::Data* d_dirac, const escript::Data* y_dirac) {

index_t color, e, row_index;
__const double  *d_dirac_p, *y_dirac_p;

double *F_p=(requireWrite(F), getSampleDataRW(F,0));    /* use comma, to get around the mixed code and declarations thing */

#pragma omp parallel private(color, d_dirac_p, y_dirac_p)
55      {      {
56            for (color=elements->minColor;color<=elements->maxColor;color++) {          for (int color=elements->minColor;color<=elements->maxColor;color++) {
57               /*  open loop over all elements: */              // loop over all elements
58               #pragma omp for private(e) schedule(static)  #pragma omp for
59               for(e=0;e<elements->numElements;e++){              for(index_t e=0; e<elements->numElements; e++) {
60                  if (elements->Color[e]==color) {                  if (elements->Color[e]==color) {
61                        const index_t row_index=p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]];
62                        if (!y_dirac.isEmpty()) {
63                            const double* y_dirac_p=y_dirac.getSampleDataRO(e);
65                                            y_dirac_p, F_p, p.row_DOF_UpperBound);
66                        }
67
68             d_dirac_p=getSampleDataRO(d_dirac, e);                      if (!d_dirac.isEmpty()) {
69                     y_dirac_p=getSampleDataRO(y_dirac, e);                          const double* d_dirac_p=d_dirac.getSampleDataRO(e);
71                     row_index=p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]];                                 p.numEqu, 1, &row_index, p.numComp, d_dirac_p);
72                                  }
73             if (NULL!=y_dirac_p)  Dudley_Util_AddScatter(1,                  } // end color check
74                                                          &row_index,              } // end element loop
75                                                          p.numEqu,          } // end color loop
76                                                          y_dirac_p,      } // end parallel region
F_p,
p.row_DOF_UpperBound);

1,
&row_index,
p.numEqu,
1,
&row_index,
p.numComp,
d_dirac_p);
} /* end color check */
} /* end element loop */
} /* end color loop */
} /* end parallel region */
77  }  }
78
79    } // namespace dudley
80

Legend:
 Removed from v.6008 changed lines Added in v.6009