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

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

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

revision 6008 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  /************************************************************************************/  /****************************************************************************
18    
19  /*    assemblage routines: calculates the normal vector at quadrature points on face elements */    Assemblage routines: calculates the normal vector at quadrature points on
20      face elements
21    
22  /************************************************************************************/  *****************************************************************************/
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
23    
24  #include "Assemble.h"  #include "Assemble.h"
 #include "Util.h"  
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
   
25  #include "ShapeTable.h"  #include "ShapeTable.h"
26    #include "Util.h"
27    
28  /************************************************************************************/  namespace dudley {
29    
30  void Dudley_Assemble_setNormal(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escript::Data* normal)  void Assemble_setNormal(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, escript::Data* normal)
31  {  {
32      double *local_X = NULL, *dVdv = NULL, *normal_array;      if (!nodes || !elements)
33      index_t sign;          return;
     dim_t e, q, NN, NS, numDim, numQuad, numDim_local;  
     bool reduced_integration;  
     const double *dSdv = 0;  
     if (nodes == NULL || elements == NULL)  
     return;  
34    
35      switch (elements->numDim)      const int NN = elements->numNodes;
36      {      const int numDim = nodes->numDim;
37      case 2:      bool reduced_integration = Assemble_reducedIntegrationOrder(normal);
38      dSdv = &(DTDV_2D[0][0]);      const int numQuad = (!reduced_integration) ? (elements->numDim + 1) : 1;
39      break;      const int numDim_local = elements->numLocalDim;
40      case 3:      const int NS = elements->numDim + 1;
41      dSdv = &(DTDV_3D[0][0]);  
42      break;      const double *dSdv = NULL;
43      default:      switch (elements->numDim) {
44      dSdv = &(DTDV_1D[0][0]);          case 2:
45      break;              dSdv = &DTDV_2D[0][0];
46      }          break;
47      Dudley_resetError();          case 3:
48      NN = elements->numNodes;              dSdv = &DTDV_3D[0][0];
49      numDim = nodes->numDim;          break;
50      reduced_integration = Dudley_Assemble_reducedIntegrationOrder(normal);          default:
51      numQuad = (!reduced_integration) ? (elements->numDim + 1) : 1;              dSdv = &DTDV_1D[0][0];
52      numDim_local = elements->numLocalDim;          break;
53      NS = elements->numDim + 1;      }
54    
55      /* set some parameters */      // check the dimensions of normal
56        if (!(numDim == numDim_local || numDim - 1 == numDim_local)) {
57      sign = 1;          throw DudleyException("Assemble_setNormal: Cannot calculate normal vector");
58      /* check the dimensions of normal */      } else if (!normal->isDataPointShapeEqual(1, &numDim)) {
59      if (!(numDim == numDim_local || numDim - 1 == numDim_local))          throw DudleyException("Assemble_setNormal: illegal point data shape of normal Data object");
60      {      } else if (!normal->numSamplesEqual(numQuad, elements->numElements)) {
61      Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: Cannot calculate normal vector");          throw DudleyException("Assemble_setNormal: illegal number of samples of normal Data object");
62      }      } else if (!normal->actsExpanded()) {
63      else if (!isDataPointShapeEqual(normal, 1, &(numDim)))          throw DudleyException("Assemble_setNormal: expanded Data object is expected for normal.");
64      {      }
65      Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal number of samples of normal Data object");  
66      }      normal->requireWrite();
67      else if (!numSamplesEqual(normal, numQuad, elements->numElements))  #pragma omp parallel
68      {      {
69      Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal number of samples of normal Data object");          std::vector<double> local_X(NS * numDim);
70      }          std::vector<double> dVdv(numQuad * numDim * numDim_local);
71      else if (!isDataPointShapeEqual(normal, 1, &(numDim)))  #pragma omp for
72      {          for (index_t e = 0; e < elements->numElements; e++) {
73      Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal point data shape of normal Data object");              // gather local coordinates of nodes into local_X
74      }              Dudley_Util_Gather_double(NS,
75      else if (!isExpanded(normal))                      &elements->Nodes[INDEX2(0, e, NN)], numDim,
76      {                      nodes->Coordinates, &local_X[0]);
77      Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: expanded Data object is expected for normal.");  
78      }              // calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q)
79                Dudley_Util_SmallMatMult(numDim, numDim_local * numQuad,
80      /* now we can start */                                       &dVdv[0], NS, &local_X[0], dSdv);
81      if (Dudley_noError())              // get normalized vector
82      {              double* normal_array = normal->getSampleDataRW(e);
83      requireWrite(normal);              Dudley_NormalVector(numQuad, numDim, numDim_local, &dVdv[0], normal_array);
84  #pragma omp parallel private(local_X,dVdv)          }
     {  
         local_X = dVdv = NULL;  
         /* allocation of work arrays */  
         local_X = new  double[NS * numDim];  
         dVdv = new  double[numQuad * numDim * numDim_local];  
         if (!(Dudley_checkPtr(local_X) || Dudley_checkPtr(dVdv)))  
         {  
         /* open the element loop */  
 #pragma omp for private(e,q,normal_array) schedule(static)  
         for (e = 0; e < elements->numElements; e++)  
         {  
             /* gather local coordinates of nodes into local_X: */  
             Dudley_Util_Gather_double(NS, &(elements->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,  
                           local_X);  
             /*  calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */  
             Dudley_Util_SmallMatMult(numDim, numDim_local * numQuad, dVdv, NS, local_X, dSdv);  
             /* get normalized vector:      */  
             normal_array = getSampleDataRW(normal, e);  
             Dudley_NormalVector(numQuad, numDim, numDim_local, dVdv, normal_array);  
             for (q = 0; q < numQuad * numDim; q++)  
             normal_array[q] *= sign;  
         }  
         }  
         delete[] local_X;  
         delete[] dVdv;  
     }  
85      }      }
86  }  }
87    
88    } // namespace dudley
89    

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

  ViewVC Help
Powered by ViewVC 1.1.26