/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_integrate.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_integrate.cpp

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

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 /****************************************************************************  
   
   Assemblage routines: integrates data on quadrature points  
   
 *****************************************************************************/  
   
17  #include "Assemble.h"  #include "Assemble.h"
18  #include "Util.h"  #include "Util.h"
19    
20  namespace dudley {  namespace dudley {
21    
22  void Assemble_integrate(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, const escript::Data* data, double* out)  void Assemble_integrate(const NodeFile* nodes, const ElementFile* elements,
23                            const escript::Data& data, std::vector<double>& out)
24  {  {
25      if (!nodes || !elements)      if (!nodes || !elements)
26          return;          return;
27    
28      const int my_mpi_rank = nodes->MPIInfo->rank;      const int my_mpi_rank = nodes->MPIInfo->rank;
29      Dudley_ElementFile_Jacobians* jac = Dudley_ElementFile_borrowJacobians(      const ElementFile_Jacobians* jac = elements->borrowJacobians(nodes,
30              elements, nodes, Assemble_reducedIntegrationOrder(data));                                           hasReducedIntegrationOrder(data));
31    
32      const dim_t numQuadTotal = jac->numQuad;      const dim_t numQuadTotal = jac->numQuad;
33      // check the shape of the data      // check the shape of the data
34      if (!data->numSamplesEqual(numQuadTotal, elements->numElements)) {      if (!data.numSamplesEqual(numQuadTotal, elements->numElements)) {
35          throw DudleyException("Assemble_integrate: illegal number of samples of integrant kernel Data object");          throw DudleyException("Assemble_integrate: illegal number of samples of integrant kernel Data object");
36      }      }
37    
38      const int numComps = data->getDataPointSize();      const int numComps = data.getDataPointSize();
39    
40      for (int q = 0; q < numComps; q++)      for (int q = 0; q < numComps; q++)
41          out[q] = 0;          out[q] = 0;
# Line 49  void Assemble_integrate(Dudley_NodeFile* Line 44  void Assemble_integrate(Dudley_NodeFile*
44      {      {
45          std::vector<double> out_local(numComps);          std::vector<double> out_local(numComps);
46    
47          if (data->actsExpanded()) {          if (data.actsExpanded()) {
48  #pragma omp for  #pragma omp for
49              for (index_t e = 0; e < elements->numElements; e++) {              for (index_t e = 0; e < elements->numElements; e++) {
50                  if (elements->Owner[e] == my_mpi_rank) {                  if (elements->Owner[e] == my_mpi_rank) {
51                      const double vol = jac->absD[e] * jac->quadweight;                      const double vol = jac->absD[e] * jac->quadweight;
52                      const double* data_array = data->getSampleDataRO(e);                      const double* data_array = data.getSampleDataRO(e);
53                      for (int q = 0; q < numQuadTotal; q++) {                      for (int q = 0; q < numQuadTotal; q++) {
54                          for (int i = 0; i < numComps; i++)                          for (int i = 0; i < numComps; i++)
55                              out_local[i] += data_array[INDEX2(i, q, numComps)] * vol;                              out_local[i] += data_array[INDEX2(i, q, numComps)] * vol;
# Line 66  void Assemble_integrate(Dudley_NodeFile* Line 61  void Assemble_integrate(Dudley_NodeFile*
61              for (index_t e = 0; e < elements->numElements; e++) {              for (index_t e = 0; e < elements->numElements; e++) {
62                  if (elements->Owner[e] == my_mpi_rank) {                  if (elements->Owner[e] == my_mpi_rank) {
63                      const double vol = jac->absD[e] * jac->quadweight;                      const double vol = jac->absD[e] * jac->quadweight;
64                      const double* data_array = data->getSampleDataRO(e);                      const double* data_array = data.getSampleDataRO(e);
65                      double rtmp = 0.;                      double rtmp = 0.;
66                      for (int q = 0; q < numQuadTotal; q++)                      for (int q = 0; q < numQuadTotal; q++)
67                          rtmp += vol;                          rtmp += vol;

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26