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

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_getSize.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: calculates the minimum distance between two vertices  
   of elements and assigns the value to each quadrature point in out  
   
 *****************************************************************************/  
   
17  #include "Assemble.h"  #include "Assemble.h"
18  #include "Util.h"  #include "Util.h"
19    
20  namespace dudley {  namespace dudley {
21    
22  void Assemble_getSize(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, escript::Data* out)  void Assemble_getSize(const NodeFile* nodes, const ElementFile* elements,
23                          escript::Data& out)
24  {  {
25      if (!nodes || !elements)      if (!nodes || !elements)
26          return;          return;
# Line 36  void Assemble_getSize(Dudley_NodeFile* n Line 30  void Assemble_getSize(Dudley_NodeFile* n
30      // now we look up what type of elements we need based on the function space      // now we look up what type of elements we need based on the function space
31      // of out. If it is DUDLEY_REDUCED_ELEMENTS or      // of out. If it is DUDLEY_REDUCED_ELEMENTS or
32      // DUDLEY_REDUCED_FACE_ELEMENTS then we have single quad point      // DUDLEY_REDUCED_FACE_ELEMENTS then we have single quad point
33      int numQuad;      int numQuad = (hasReducedIntegrationOrder(out) ? 1 : elements->numNodes);
     if (Assemble_reducedIntegrationOrder(out)) {  
         numQuad = 1;  
     } else {  
         numQuad = elements->numDim + 1;  
     }  
34      const int NN = elements->numNodes;      const int NN = elements->numNodes;
35      const int NS = elements->numDim + 1;      const int NS = elements->numDim + 1;
36      const int NVertices = elements->numDim + 1;      const int NVertices = elements->numDim + 1;
37    
38      // check the dimensions of out      // check the dimensions of out
39      if (!out->numSamplesEqual(numQuad, elements->numElements)) {      if (!out.numSamplesEqual(numQuad, elements->numElements)) {
40          throw DudleyException("Assemble_getSize: illegal number of samples of element size Data object");          throw DudleyException("Assemble_getSize: illegal number of samples of element size Data object");
41      } else if (!out->isDataPointShapeEqual(0, &numDim)) {      } else if (!out.isDataPointShapeEqual(0, &numDim)) {
42          throw DudleyException("Assemble_getSize: illegal data point shape of element size Data object");          throw DudleyException("Assemble_getSize: illegal data point shape of element size Data object");
43      } else if (!out->actsExpanded()) {      } else if (!out.actsExpanded()) {
44          throw DudleyException("Assemble_getSize: expanded Data object is expected for element size.");          throw DudleyException("Assemble_getSize: expanded Data object is expected for element size.");
45      }      }
46    
47      // now we can start      // now we can start
48      out->requireWrite();      out.requireWrite();
49  #pragma omp parallel  #pragma omp parallel
50      {      {
51          std::vector<double> local_X(NN * numDim);          std::vector<double> local_X(NN * numDim);
52  #pragma omp for  #pragma omp for
53          for (index_t e = 0; e < elements->numElements; e++) {          for (index_t e = 0; e < elements->numElements; e++) {
54              // gather local coordinates of nodes into local_X(numDim,NN)              // gather local coordinates of nodes into local_X(numDim,NN)
55              Dudley_Util_Gather_double(NS, &elements->Nodes[INDEX2(0, e, NN)],              util::gather(NS, &elements->Nodes[INDEX2(0, e, NN)], numDim,
56                                        numDim, nodes->Coordinates, &local_X[0]);                           nodes->Coordinates, &local_X[0]);
57              // calculate minimal differences              // calculate minimal differences
58              double max_diff = 0;              double max_diff = 0;
59              for (int n0 = 0; n0 < NVertices; n0++) {              for (int n0 = 0; n0 < NVertices; n0++) {
# Line 80  void Assemble_getSize(Dudley_NodeFile* n Line 69  void Assemble_getSize(Dudley_NodeFile* n
69              }              }
70              max_diff = sqrt(max_diff);              max_diff = sqrt(max_diff);
71              // set all values to max_diff              // set all values to max_diff
72              double* out_array = out->getSampleDataRW(e);              double* out_array = out.getSampleDataRW(e);
73              for (int q = 0; q < numQuad; q++)              for (int q = 0; q < numQuad; q++)
74                  out_array[q] = max_diff;                  out_array[q] = max_diff;
75          }          }
Legend:
Removed from v.6078  
changed lines
  Added in v.6079
  ViewVC Help
Powered by ViewVC 1.1.26