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

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)