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

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

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

revision 6008 by caltinay, Mon Feb 22 06:59:27 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  /************************************************************************************/  /****************************************************************************/
18    
19  /*   Dudley: Mesh */  /*   Dudley: Mesh */
20    
21  /* searches for faces in the mesh which are matching */  /* searches for faces in the mesh which are matching */
22    
23  /************************************************************************************/  /****************************************************************************/
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
24    
25  #include "Mesh.h"  #include "Mesh.h"
26  #include "ShapeTable.h"  #include "ShapeTable.h"
27  #include "Util.h"  #include "Util.h"
28    
29    namespace dudley {
 /************************************************************************************/  
30    
31  static double Dudley_Mesh_lockingGridSize = 0;  static double Dudley_Mesh_lockingGridSize = 0;
32    
# Line 41  int Dudley_Mesh_findMatchingFaces_compar Line 37  int Dudley_Mesh_findMatchingFaces_compar
37      dim_t i;      dim_t i;
38      e1 = (Dudley_Mesh_findMatchingFaces_center *) arg1;      e1 = (Dudley_Mesh_findMatchingFaces_center *) arg1;
39      e2 = (Dudley_Mesh_findMatchingFaces_center *) arg2;      e2 = (Dudley_Mesh_findMatchingFaces_center *) arg2;
40      for (i = 0; i < MAX_numDim; i++)      for (i = 0; i < 3; i++) {
     {  
41          l = (e1->x[i] < e2->x[i] + Dudley_Mesh_lockingGridSize);          l = (e1->x[i] < e2->x[i] + Dudley_Mesh_lockingGridSize);
42          g = (e2->x[i] < e1->x[i] + Dudley_Mesh_lockingGridSize);          g = (e2->x[i] < e1->x[i] + Dudley_Mesh_lockingGridSize);
43          if (!(l && g))          if (!(l && g)) {
         {  
44              if (l)              if (l)
45                  return -1;                  return -1;
46              if (g)              if (g)
# Line 54  int Dudley_Mesh_findMatchingFaces_compar Line 48  int Dudley_Mesh_findMatchingFaces_compar
48          }          }
49      }      }
50      if (e1->refId < e2->refId)      if (e1->refId < e2->refId)
     {  
51          return -1;          return -1;
     }  
52      else if (e1->refId > e2->refId)      else if (e1->refId > e2->refId)
     {  
53          return 1;          return 1;
     }  
54      else      else
     {  
55          return 0;          return 0;
56    }
57    
58    inline double getDist(int e0, int i0, int e1, int i1, int numDim, int NN,
59                          const double* X)
60    {
61        double dist=0.;
62        for (int i=0; i<numDim; i++) {
63            dist=std::max(dist, std::abs(X[INDEX3(i, i0, e0, numDim, NN)]
64                        - X[INDEX3(i, i1, e1, numDim, NN)]));
65      }      }
66        return dist;
67  }  }
68    
69  void Dudley_Mesh_findMatchingFaces(Dudley_NodeFile* nodes,  void Dudley_Mesh_findMatchingFaces(Dudley_NodeFile* nodes,
# Line 72  void Dudley_Mesh_findMatchingFaces(Dudle Line 71  void Dudley_Mesh_findMatchingFaces(Dudle
71          dim_t* numPairs, index_t* elem0, index_t* elem1,          dim_t* numPairs, index_t* elem0, index_t* elem1,
72          index_t* matching_nodes_in_elem1)          index_t* matching_nodes_in_elem1)
73  {  {
74  #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \      double h = HUGE_VAL, dist, *X = NULL;
       {dim_t i;   \  
       _dist_=0; \  
       for (i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)])); \  
       }  
     double h = HUGE_VAL, h_local, dist, *X = NULL;  
75      Dudley_Mesh_findMatchingFaces_center *center;      Dudley_Mesh_findMatchingFaces_center *center;
76      index_t e_0, e_1, *a1 = NULL, *a2 = NULL, *perm = NULL, *perm_tmp = NULL, *itmp_ptr = NULL;      index_t e_0, e_1, *a1 = NULL, *a2 = NULL, *perm = NULL, *perm_tmp = NULL, *itmp_ptr = NULL;
77      const index_t *shiftNodes = NULL, *reverseNodes = NULL;      const index_t *shiftNodes = NULL, *reverseNodes = NULL;
78      dim_t e, i, i0, i1, n, NN, numNodesOnFace;      dim_t e, i, i0, i1, n, numNodesOnFace;
79    
80      dim_t numDim = nodes->numDim;      dim_t numDim = nodes->numDim;
81        int NN = faces->numNodes;
     NN = faces->numNodes;  
   
82      numNodesOnFace = numNodesOnFaceMap[faces->etype];      numNodesOnFace = numNodesOnFaceMap[faces->etype];
83      shiftNodes = shiftNodesMap[faces->etype];      shiftNodes = shiftNodesMap[faces->etype];
84      reverseNodes = reverseNodesMap[faces->etype];      reverseNodes = reverseNodesMap[faces->etype];
85    
86      if (numNodesOnFace <= 0)      if (numNodesOnFace <= 0) {
     {  
87          std::stringstream ss;          std::stringstream ss;
88          ss << "Dudley_Mesh_findMatchingFaces: matching faces cannot be "          ss << "Dudley_Mesh_findMatchingFaces: matching faces cannot be "
89              "applied to face elements of type " << getElementName(faces->etype);              "applied to face elements of type " << getElementName(faces->etype);
90          std::string errorMsg(ss.str());          const std::string errorMsg(ss.str());
91          Dudley_setError(TYPE_ERROR, errorMsg.c_str());          throw DudleyException(errorMsg);
         return;  
92      }      }
93      X = new double[NN * numDim * faces->numElements];      X = new double[NN * numDim * faces->numElements];
94      center = new  Dudley_Mesh_findMatchingFaces_center[faces->numElements];      center = new Dudley_Mesh_findMatchingFaces_center[faces->numElements];
95      a1 = new  int[NN];      a1 = new int[NN];
96      a2 = new  int[NN];      a2 = new int[NN];
97      if (!(Dudley_checkPtr(X) || Dudley_checkPtr(center) || Dudley_checkPtr(a1) || Dudley_checkPtr(a2)))  
98      {      /* OMP */
99          /* OMP */      for (e = 0; e < faces->numElements; e++) {
100          for (e = 0; e < faces->numElements; e++)          /* get the coordinates of the nodes */
101          {          Dudley_Util_Gather_double(NN, &(faces->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
102              /* get the coordinates of the nodes */                                    &(X[INDEX3(0, 0, e, numDim, NN)]));
103              Dudley_Util_Gather_double(NN, &(faces->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,          /* get the element center */
104                                        &(X[INDEX3(0, 0, e, numDim, NN)]));          center[e].refId = e;
105              /* get the element center */          for (i = 0; i < 3; i++)
106              center[e].refId = e;              center[e].x[i] = 0;
107              for (i = 0; i < MAX_numDim; i++)          for (i0 = 0; i0 < numNodesOnFace; i0++) {
                 center[e].x[i] = 0;  
             for (i0 = 0; i0 < numNodesOnFace; i0++)  
             {  
                 for (i = 0; i < numDim; i++)  
                     center[e].x[i] += X[INDEX3(i, i0, e, numDim, NN)];  
             }  
108              for (i = 0; i < numDim; i++)              for (i = 0; i < numDim; i++)
109                  center[e].x[i] /= numNodesOnFace;                  center[e].x[i] += X[INDEX3(i, i0, e, numDim, NN)];
110              /* get the minimum distance between nodes in the element */          }
111              for (i0 = 0; i0 < numNodesOnFace; i0++)          for (i = 0; i < numDim; i++)
112              {              center[e].x[i] /= numNodesOnFace;
113                  for (i1 = i0 + 1; i1 < numNodesOnFace; i1++)          /* get the minimum distance between nodes in the element */
114                  {          for (i0 = 0; i0 < numNodesOnFace; i0++) {
115                      getDist(h_local, e, i0, e, i1);              for (i1 = i0 + 1; i1 < numNodesOnFace; i1++) {
116                      h = MIN(h, h_local);                  const double h_local=getDist(e, i0, e, i1, numDim, NN, X);
117                  }                  h=std::min(h, h_local);
118              }              }
119          }          }
120          /* set the */      }
121          Dudley_Mesh_lockingGridSize = h * MAX(safety_factor, 0);  
122        Dudley_Mesh_lockingGridSize = h * std::max(safety_factor, 0.);
123  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
124          printf("locking grid size is %e\n", Dudley_Mesh_lockingGridSize);      printf("locking grid size is %e\n", Dudley_Mesh_lockingGridSize);
125          printf("absolute tolerance is %e.\n", h * tolerance);      printf("absolute tolerance is %e.\n", h * tolerance);
126  #endif  #endif
127          /* sort the elements by center coordinates (lexicographical) */      /* sort the elements by center coordinates (lexicographical) */
128          qsort(center, faces->numElements, sizeof(Dudley_Mesh_findMatchingFaces_center),      qsort(center, faces->numElements, sizeof(Dudley_Mesh_findMatchingFaces_center),
129                Dudley_Mesh_findMatchingFaces_compar);            Dudley_Mesh_findMatchingFaces_compar);
130          /* find elements with matching center */      /* find elements with matching center */
131          *numPairs = 0;      *numPairs = 0;
132          /* OMP */      /* OMP */
133          for (e = 0; e < faces->numElements - 1 && Dudley_noError(); e++)      for (e = 0; e < faces->numElements - 1; e++) {
134          {          dist = 0;
135              dist = 0;          for (i = 0; i < numDim; i++)
136              for (i = 0; i < numDim; i++)              dist = std::max(dist, std::abs(center[e].x[i] - center[e + 1].x[i]));
137                  dist = MAX(dist, ABS(center[e].x[i] - center[e + 1].x[i]));          if (dist < h * tolerance) {
138              if (dist < h * tolerance)              e_0 = center[e].refId;
139              {              e_1 = center[e + 1].refId;
140                  e_0 = center[e].refId;              elem0[*numPairs] = e_0;
141                  e_1 = center[e + 1].refId;              elem1[*numPairs] = e_1;
142                  elem0[*numPairs] = e_0;              /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
143                  elem1[*numPairs] = e_1;              perm = a1;
144                  /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */              perm_tmp = a2;
145                  perm = a1;              for (i = 0; i < NN; i++)
146                  perm_tmp = a2;                  perm[i] = i;
147                  for (i = 0; i < NN; i++)              while (1) {
148                      perm[i] = i;                  /* if node 0 and perm[0] are the same we are ready */
149                  while (Dudley_noError())                  dist=getDist(e_0, 0, e_1, perm[0], numDim, NN, X);
150                    if (dist <= h * tolerance)
151                        break;
152                    if (shiftNodes[0] >= 0) {
153                        /* rotate the nodes */
154                        itmp_ptr = perm;
155                        perm = perm_tmp;
156                        perm_tmp = itmp_ptr;
157    #pragma ivdep
158                        for (i = 0; i < NN; i++)
159                            perm[i] = perm_tmp[shiftNodes[i]];
160                    }
161                    /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
162                    if (perm[0] == 0)
163                  {                  {
164                      /* if node 0 and perm[0] are the same we are ready */                      std::stringstream ss;
165                      getDist(dist, e_0, 0, e_1, perm[0]);                      ss << "Mesh_findMatchingFaces: couldn't match first "
166                      if (dist <= h * tolerance)                          "node of element " << e_0 << " to touching element "
167                          break;                          << e_1;
168                      if (shiftNodes[0] >= 0)                      std::string errorMsg(ss.str());
169                      {                      throw DudleyException(errorMsg);
170                          /* rotate the nodes */                  }
171                }
172                /* now we check if the second nodes match */
173                if (numNodesOnFace > 1) {
174                    dist=getDist(e_0, 1, e_1, perm[1], numDim, NN, X);
175                    /* if the second node does not match we reverse the direction of the nodes */
176                    if (dist > h * tolerance) {
177                        /* rotate the nodes */
178                        if (reverseNodes[0] < 0) {
179                            std::stringstream ss;
180                            ss << "Mesh_findMatchingFaces: couldn't match "
181                                "the second node of element " << e_0
182                                << " to touching element " << e_1;
183                            const std::string errorMsg(ss.str());
184                            throw DudleyException(errorMsg);
185                        } else {
186                          itmp_ptr = perm;                          itmp_ptr = perm;
187                          perm = perm_tmp;                          perm = perm_tmp;
188                          perm_tmp = itmp_ptr;                          perm_tmp = itmp_ptr;
189  #pragma ivdep  #pragma ivdep
190                          for (i = 0; i < NN; i++)                          for (i = 0; i < NN; i++)
191                              perm[i] = perm_tmp[shiftNodes[i]];                              perm[i] = perm_tmp[reverseNodes[i]];
192                      }                          dist=getDist(e_0, 1, e_1, perm[1], numDim, NN, X);
193                      /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */                          if (dist > h * tolerance) {
                     if (perm[0] == 0)  
                     {  
                         std::stringstream ss;  
                         ss << "Mesh_findMatchingFaces: couldn't match first "  
                             "node of element " << e_0 << " to touching element "  
                             << e_1;  
                         std::string errorMsg(ss.str());  
                         Dudley_setError(VALUE_ERROR, errorMsg.c_str());  
                     }  
                 }  
                 /* now we check if the second nodes match */  
                 if (Dudley_noError())  
                 {  
                     if (numNodesOnFace > 1)  
                     {  
                         getDist(dist, e_0, 1, e_1, perm[1]);  
                         /* if the second node does not match we reverse the direction of the nodes */  
                         if (dist > h * tolerance)  
                         {  
                             /* rotate the nodes */  
                             if (reverseNodes[0] < 0)  
                             {  
                                 std::stringstream ss;  
                                 ss << "Mesh_findMatchingFaces: couldn't match "  
                                     "the second node of element " << e_0  
                                     << " to touching element " << e_1;  
                                 std::string errorMsg(ss.str());  
                                 Dudley_setError(VALUE_ERROR, errorMsg.c_str());  
                             }  
                             else  
                             {  
                                 itmp_ptr = perm;  
                                 perm = perm_tmp;  
                                 perm_tmp = itmp_ptr;  
 #pragma ivdep  
                                 for (i = 0; i < NN; i++)  
                                     perm[i] = perm_tmp[reverseNodes[i]];  
                                 getDist(dist, e_0, 1, e_1, perm[1]);  
                                 if (dist > h * tolerance)  
                                 {  
                                     std::stringstream ss;  
                                     ss << "Mesh_findMatchingFaces: couldn't match "  
                                         "the second node of element " << e_0  
                                         << " to touching element " << e_1;  
                                     std::string errorMsg(ss.str());  
                                     Dudley_setError(VALUE_ERROR, errorMsg.c_str());  
                                 }  
                             }  
                         }  
                     }  
                 }  
                 /* we check if the rest of the face nodes match: */  
                 if (Dudley_noError())  
                 {  
                     for (i = 2; i < numNodesOnFace; i++)  
                     {  
                         n = i;  
                         getDist(dist, e_0, n, e_1, perm[n]);  
                         if (dist > h * tolerance)  
                         {  
194                              std::stringstream ss;                              std::stringstream ss;
195                              ss << "Mesh_findMatchingFaces: couldn't match the "                              ss << "Mesh_findMatchingFaces: couldn't match "
196                                  << i << "-th node of element " << e_0                                  "the second node of element " << e_0
197                                  << " to touching element " << e_1;                                  << " to touching element " << e_1;
198                              std::string errorMsg(ss.str());                              const std::string errorMsg(ss.str());
199                              Dudley_setError(VALUE_ERROR, errorMsg.c_str());                              throw DudleyException(errorMsg);
                             break;  
200                          }                          }
201                      }                      }
202                  }                  }
203                  /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */              }
204                  if (Dudley_noError())              /* we check if the rest of the face nodes match: */
205                  {              for (i = 2; i < numNodesOnFace; i++) {
206                      for (i = 0; i < NN; i++)                  n = i;
207                          matching_nodes_in_elem1[INDEX2(i, *numPairs, NN)] = faces->Nodes[INDEX2(perm[i], e_1, NN)];                  dist=getDist(e_0, n, e_1, perm[n], numDim, NN, X);
208                    if (dist > h * tolerance) {
209                        std::stringstream ss;
210                        ss << "Mesh_findMatchingFaces: couldn't match the "
211                            << i << "-th node of element " << e_0
212                            << " to touching element " << e_1;
213                        const std::string errorMsg(ss.str());
214                        throw DudleyException(errorMsg);
215                  }                  }
                 (*numPairs)++;  
216              }              }
217                /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
218                for (i = 0; i < NN; i++)
219                    matching_nodes_in_elem1[INDEX2(i, *numPairs, NN)] = faces->Nodes[INDEX2(perm[i], e_1, NN)];
220                (*numPairs)++;
221          }          }
222        }
223  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
224          printf("number of pairs of matching faces %d\n", *numPairs);      printf("number of pairs of matching faces %d\n", *numPairs);
225  #endif  #endif
226      }  
     /* clean up */  
227      delete[] X;      delete[] X;
228      delete[] center;      delete[] center;
229      delete[] a1;      delete[] a1;
230      delete[] a1;      delete[] a1;
   
 #undef getDist  
231  }  }
232    
233    } // namespace dudley
234    

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

  ViewVC Help
Powered by ViewVC 1.1.26