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

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)
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);
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