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

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

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

branches/trilinos_from_5897/dudley/src/Assemble_setNormal.cpp revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC branches/trilinos_from_5897/dudley/src/Assemble_getNormal.cpp revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 /****************************************************************************  
   
   Assemblage routines: calculates the normal vector at quadrature points on  
   face elements  
   
 *****************************************************************************/  
   
17  #include "Assemble.h"  #include "Assemble.h"
18  #include "ShapeTable.h"  #include "ShapeTable.h"
19  #include "Util.h"  #include "Util.h"
20    
21  namespace dudley {  namespace dudley {
22    
23  void Assemble_setNormal(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, escript::Data* normal)  void Assemble_getNormal(const NodeFile* nodes, const ElementFile* elements,
24                            escript::Data& normal)
25  {  {
26      if (!nodes || !elements)      if (!nodes || !elements)
27          return;          return;
28    
29      const int NN = elements->numNodes;      const int NN = elements->numNodes;
30      const int numDim = nodes->numDim;      const int numDim = nodes->numDim;
31      bool reduced_integration = Assemble_reducedIntegrationOrder(normal);      const int numQuad = (hasReducedIntegrationOrder(normal) ? 1 : NN);
     const int numQuad = (!reduced_integration) ? (elements->numDim + 1) : 1;  
32      const int numDim_local = elements->numLocalDim;      const int numDim_local = elements->numLocalDim;
33      const int NS = elements->numDim + 1;      const int NS = elements->numDim + 1;
34    
# Line 55  void Assemble_setNormal(Dudley_NodeFile* Line 48  void Assemble_setNormal(Dudley_NodeFile*
48      // check the dimensions of normal      // check the dimensions of normal
49      if (!(numDim == numDim_local || numDim - 1 == numDim_local)) {      if (!(numDim == numDim_local || numDim - 1 == numDim_local)) {
50          throw DudleyException("Assemble_setNormal: Cannot calculate normal vector");          throw DudleyException("Assemble_setNormal: Cannot calculate normal vector");
51      } else if (!normal->isDataPointShapeEqual(1, &numDim)) {      } else if (!normal.isDataPointShapeEqual(1, &numDim)) {
52          throw DudleyException("Assemble_setNormal: illegal point data shape of normal Data object");          throw DudleyException("Assemble_setNormal: illegal point data shape of normal Data object");
53      } else if (!normal->numSamplesEqual(numQuad, elements->numElements)) {      } else if (!normal.numSamplesEqual(numQuad, elements->numElements)) {
54          throw DudleyException("Assemble_setNormal: illegal number of samples of normal Data object");          throw DudleyException("Assemble_setNormal: illegal number of samples of normal Data object");
55      } else if (!normal->actsExpanded()) {      } else if (!normal.actsExpanded()) {
56          throw DudleyException("Assemble_setNormal: expanded Data object is expected for normal.");          throw DudleyException("Assemble_setNormal: expanded Data object is expected for normal.");
57      }      }
58    
59      normal->requireWrite();      normal.requireWrite();
60  #pragma omp parallel  #pragma omp parallel
61      {      {
62          std::vector<double> local_X(NS * numDim);          std::vector<double> local_X(NS * numDim);
# Line 71  void Assemble_setNormal(Dudley_NodeFile* Line 64  void Assemble_setNormal(Dudley_NodeFile*
64  #pragma omp for  #pragma omp for
65          for (index_t e = 0; e < elements->numElements; e++) {          for (index_t e = 0; e < elements->numElements; e++) {
66              // gather local coordinates of nodes into local_X              // gather local coordinates of nodes into local_X
67              Dudley_Util_Gather_double(NS,              util::gather(NS, &elements->Nodes[INDEX2(0, e, NN)], numDim,
68                      &elements->Nodes[INDEX2(0, e, NN)], numDim,                           nodes->Coordinates, &local_X[0]);
                     nodes->Coordinates, &local_X[0]);  
69    
70              // calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q)              // calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q)
71              Dudley_Util_SmallMatMult(numDim, numDim_local * numQuad,              util::smallMatMult(numDim, numDim_local * numQuad,
72                                       &dVdv[0], NS, &local_X[0], dSdv);                                       &dVdv[0], NS, &local_X[0], dSdv);
73              // get normalized vector              // get normalized vector
74              double* normal_array = normal->getSampleDataRW(e);              double* normal_array = normal.getSampleDataRW(e);
75              Dudley_NormalVector(numQuad, numDim, numDim_local, &dVdv[0], normal_array);              util::normalVector(numQuad, numDim, numDim_local, &dVdv[0], normal_array);
76          }          }
77      }      }
78  }  }
Legend:
Removed from v.6078  
changed lines
  Added in v.6079
  ViewVC Help
Powered by ViewVC 1.1.26