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

Diff of /branches/trilinos_from_5897/dudley/src/Mesh_createNodeFileMappings.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  /*   Dudley: NodeFile : creates the mappings using the indexReducedNodes */  /*   Dudley: NodeFile : creates the mappings using the indexReducedNodes */
20  /*                 no distribution is happening                          */  /*                 no distribution is happening                          */
21    
22  /************************************************************************************/  /****************************************************************************/
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
23    
24  #include "Mesh.h"  #include "Mesh.h"
 #define UNUSED -1  
25    
26  /************************************************************************************/  namespace dudley {
27    
28    #define UNUSED -1
29    
30  void Dudley_Mesh_createDOFMappingAndCoupling(Dudley_Mesh * in, bool use_reduced_elements)  void Dudley_Mesh_createDOFMappingAndCoupling(Dudley_Mesh* in, bool use_reduced_elements)
31  {  {
32      index_t min_DOF, max_DOF, *shared = NULL, *offsetInShared = NULL, *locDOFMask =      index_t min_DOF, max_DOF, *shared = NULL, *offsetInShared = NULL, *locDOFMask =
33      NULL, i, k, myFirstDOF, myLastDOF, *nodeMask = NULL, firstDOF, lastDOF, *globalDOFIndex, *wanted_DOFs = NULL;          NULL, i, k, myFirstDOF, myLastDOF, *nodeMask = NULL, firstDOF, lastDOF, *globalDOFIndex, *wanted_DOFs = NULL;
34      dim_t mpiSize, len_loc_dof, numNeighbors, n, lastn, numNodes, *rcv_len = NULL, *snd_len = NULL, count;      dim_t mpiSize, len_loc_dof, numNeighbors, n, lastn, numNodes, *rcv_len = NULL, *snd_len = NULL, count;
35      Esys_MPI_rank myRank, p, p_min, p_max, *neighbor = NULL;      int myRank, p, p_min, p_max, *neighbor = NULL;
36      paso::SharedComponents_ptr rcv_shcomp, snd_shcomp;      paso::SharedComponents_ptr rcv_shcomp, snd_shcomp;
37      Dudley_NodeMapping *this_mapping = NULL;      Dudley_NodeMapping *this_mapping = NULL;
38      paso::Connector_ptr this_connector;      paso::Connector_ptr this_connector;
39      paso::Distribution_ptr dof_distribution;      paso::Distribution_ptr dof_distribution;
40      esysUtils::JMPI& mpi_info = in->MPIInfo;      escript::JMPI mpi_info(in->MPIInfo);
41  #ifdef ESYS_MPI  #ifdef ESYS_MPI
42      MPI_Request *mpi_requests = NULL;      MPI_Request *mpi_requests = NULL;
43      MPI_Status *mpi_stati = NULL;      MPI_Status *mpi_stati = NULL;
# Line 50  void Dudley_Mesh_createDOFMappingAndCoup Line 48  void Dudley_Mesh_createDOFMappingAndCoup
48      numNodes = in->Nodes->numNodes;      numNodes = in->Nodes->numNodes;
49      if (use_reduced_elements)      if (use_reduced_elements)
50      {      {
51      dof_distribution = in->Nodes->reducedDegreesOfFreedomDistribution;          dof_distribution = in->Nodes->reducedDegreesOfFreedomDistribution;
52      globalDOFIndex = in->Nodes->globalReducedDOFIndex;          globalDOFIndex = in->Nodes->globalReducedDOFIndex;
53      }      }
54      else      else
55      {      {
56      dof_distribution = in->Nodes->degreesOfFreedomDistribution;          dof_distribution = in->Nodes->degreesOfFreedomDistribution;
57      globalDOFIndex = in->Nodes->globalDegreesOfFreedom;          globalDOFIndex = in->Nodes->globalDegreesOfFreedom;
58      }      }
59      myFirstDOF = dof_distribution->getFirstComponent();      myFirstDOF = dof_distribution->getFirstComponent();
60      myLastDOF = dof_distribution->getLastComponent();      myLastDOF = dof_distribution->getLastComponent();
# Line 69  void Dudley_Mesh_createDOFMappingAndCoup Line 67  void Dudley_Mesh_createDOFMappingAndCoup
67    
68      if (max_DOF < min_DOF)      if (max_DOF < min_DOF)
69      {      {
70      min_DOF = myFirstDOF;          min_DOF = myFirstDOF;
71      max_DOF = myLastDOF - 1;          max_DOF = myLastDOF - 1;
72      }      }
73    
74      p_min = mpiSize;      p_min = mpiSize;
75      p_max = -1;      p_max = -1;
76      if (max_DOF >= min_DOF)      if (max_DOF >= min_DOF)
77      {      {
78      for (p = 0; p < mpiSize; ++p)          for (p = 0; p < mpiSize; ++p)
79      {          {
80          if (dof_distribution->first_component[p] <= min_DOF)              if (dof_distribution->first_component[p] <= min_DOF)
81          p_min = p;                  p_min = p;
82          if (dof_distribution->first_component[p] <= max_DOF)              if (dof_distribution->first_component[p] <= max_DOF)
83          p_max = p;                  p_max = p;
84      }          }
85      }      }
86    
87      len_loc_dof = max_DOF - min_DOF + 1;      len_loc_dof = max_DOF - min_DOF + 1;
88      if (!((min_DOF <= myFirstDOF) && (myLastDOF - 1 <= max_DOF)))      if (!((min_DOF <= myFirstDOF) && (myLastDOF - 1 <= max_DOF)))
89      {      {
90      Dudley_setError(SYSTEM_ERROR, "Local elements do not span local degrees of freedom.");          throw DudleyException("Local elements do not span local degrees of freedom.");
     return;  
91      }      }
92      rcv_len = new  dim_t[mpiSize];      rcv_len = new  dim_t[mpiSize];
93      snd_len = new  dim_t[mpiSize];      snd_len = new  dim_t[mpiSize];
# Line 103  void Dudley_Mesh_createDOFMappingAndCoup Line 100  void Dudley_Mesh_createDOFMappingAndCoup
100  #endif  #endif
101      wanted_DOFs = new  index_t[numNodes];      wanted_DOFs = new  index_t[numNodes];
102      nodeMask = new  index_t[numNodes];      nodeMask = new  index_t[numNodes];
103      neighbor = new  Esys_MPI_rank[mpiSize];      neighbor = new  int[mpiSize];
104      shared = new  index_t[numNodes * (p_max - p_min + 1)];      shared = new  index_t[numNodes * (p_max - p_min + 1)];
105      offsetInShared = new  index_t[mpiSize + 1];      offsetInShared = new  index_t[mpiSize + 1];
106      locDOFMask = new  index_t[len_loc_dof];      locDOFMask = new  index_t[len_loc_dof];
     if (!  
     (Dudley_checkPtr(neighbor) || Dudley_checkPtr(shared) || Dudley_checkPtr(offsetInShared)  
      || Dudley_checkPtr(locDOFMask) || Dudley_checkPtr(nodeMask) || Dudley_checkPtr(rcv_len)  
      || Dudley_checkPtr(snd_len) || Dudley_checkPtr(mpi_requests) || Dudley_checkPtr(mpi_stati)  
      || Dudley_checkPtr(mpi_stati)))  
     {  
107    
108      memset(rcv_len, 0, sizeof(dim_t) * mpiSize);      memset(rcv_len, 0, sizeof(dim_t) * mpiSize);
109  #pragma omp parallel  #pragma omp parallel
110      {      {
111  #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
112          for (i = 0; i < len_loc_dof; ++i)          for (i = 0; i < len_loc_dof; ++i)
113          locDOFMask[i] = UNUSED;              locDOFMask[i] = UNUSED;
114  #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
115          for (i = 0; i < numNodes; ++i)          for (i = 0; i < numNodes; ++i)
116          nodeMask[i] = UNUSED;              nodeMask[i] = UNUSED;
117  #pragma omp for private(i,k) schedule(static)  #pragma omp for private(i,k) schedule(static)
118          for (i = 0; i < numNodes; ++i)          for (i = 0; i < numNodes; ++i)
119          {          {
120          k = globalDOFIndex[i];              k = globalDOFIndex[i];
121          if (k > -1)              if (k > -1)
122          {              {
123              locDOFMask[k - min_DOF] = UNUSED - 1;                  locDOFMask[k - min_DOF] = UNUSED - 1;
124  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
125              if ((k - min_DOF) >= len_loc_dof)                  if ((k - min_DOF) >= len_loc_dof)
126              {                  {
127              printf("BOUNDS_CHECK %s %d i=%d k=%d min_DOF=%d\n", __FILE__, __LINE__, i, k, min_DOF);                      printf("BOUNDS_CHECK %s %d i=%d k=%d min_DOF=%d\n", __FILE__, __LINE__, i, k, min_DOF);
128              exit(1);                      exit(1);
129              }                  }
130  #endif  #endif
131          }              }
132          }          }
133    
134  #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
135          for (i = myFirstDOF - min_DOF; i < myLastDOF - min_DOF; ++i)          for (i = myFirstDOF - min_DOF; i < myLastDOF - min_DOF; ++i)
136          {          {
137          locDOFMask[i] = i - myFirstDOF + min_DOF;              locDOFMask[i] = i - myFirstDOF + min_DOF;
138  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
139          if (i < 0 || i >= len_loc_dof)              if (i < 0 || i >= len_loc_dof)
140          {              {
141              printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);                  printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
142              exit(1);                  exit(1);
143          }              }
144  #endif  #endif
145          }          }
146      }      }
147    
148      numNeighbors = 0;      numNeighbors = 0;
149      n = 0;      n = 0;
150      lastn = n;      lastn = n;
151      for (p = p_min; p <= p_max; ++p)      for (p = p_min; p <= p_max; ++p)
152      {      {
153          firstDOF = MAX(min_DOF, dof_distribution->first_component[p]);          firstDOF = std::max(min_DOF, dof_distribution->first_component[p]);
154          lastDOF = MIN(max_DOF + 1, dof_distribution->first_component[p + 1]);          lastDOF = std::min(max_DOF + 1, dof_distribution->first_component[p + 1]);
155          if (p != myRank)          if (p != myRank)
156          {          {
157          for (i = firstDOF - min_DOF; i < lastDOF - min_DOF; ++i)              for (i = firstDOF - min_DOF; i < lastDOF - min_DOF; ++i)
158          {              {
159  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
160              if (i < 0 || i >= len_loc_dof)                  if (i < 0 || i >= len_loc_dof)
161              {                  {
162              printf("BOUNDS_CHECK %s %d p=%d i=%d\n", __FILE__, __LINE__, p, i);                      printf("BOUNDS_CHECK %s %d p=%d i=%d\n", __FILE__, __LINE__, p, i);
163              exit(1);                      exit(1);
164              }                  }
165  #endif  #endif
166              if (locDOFMask[i] == UNUSED - 1)                  if (locDOFMask[i] == UNUSED - 1)
167              {                  {
168              locDOFMask[i] = myLastDOF - myFirstDOF + n;                      locDOFMask[i] = myLastDOF - myFirstDOF + n;
169              wanted_DOFs[n] = i + min_DOF;                      wanted_DOFs[n] = i + min_DOF;
170              ++n;                      ++n;
171              }                  }
172          }              }
173          if (n > lastn)              if (n > lastn)
174          {              {
175              rcv_len[p] = n - lastn;                  rcv_len[p] = n - lastn;
176              neighbor[numNeighbors] = p;                  neighbor[numNeighbors] = p;
177  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
178              if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)                  if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
179              {                  {
180              printf("BOUNDS_CHECK %s %d p=%d numNeighbors=%d n=%d\n", __FILE__, __LINE__, p, numNeighbors,                      printf("BOUNDS_CHECK %s %d p=%d numNeighbors=%d n=%d\n", __FILE__, __LINE__, p, numNeighbors,
181                     n);                             n);
182              exit(1);                      exit(1);
183              }                  }
184  #endif  #endif
185              offsetInShared[numNeighbors] = lastn;                  offsetInShared[numNeighbors] = lastn;
186              numNeighbors++;                  numNeighbors++;
187              lastn = n;                  lastn = n;
188          }              }
189          }          }
190      }      }
191  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
192      if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)      if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
193      {      {
194          printf("BOUNDS_CHECK %s %d numNeighbors=%d\n", __FILE__, __LINE__, numNeighbors);          printf("BOUNDS_CHECK %s %d numNeighbors=%d\n", __FILE__, __LINE__, numNeighbors);
195          exit(1);          exit(1);
196      }      }
197  #endif  #endif
198      offsetInShared[numNeighbors] = lastn;      offsetInShared[numNeighbors] = lastn;
199    
200      /* assign new DOF labels to nodes */      /* assign new DOF labels to nodes */
201  #pragma omp parallel for private(i,k) schedule(static)  #pragma omp parallel for private(i,k) schedule(static)
202      for (i = 0; i < numNodes; ++i)      for (i = 0; i < numNodes; ++i)
203      {      {
204          k = globalDOFIndex[i];          k = globalDOFIndex[i];
205          if (k > -1)          if (k > -1)
206          nodeMask[i] = locDOFMask[k - min_DOF];              nodeMask[i] = locDOFMask[k - min_DOF];
207      }      }
208    
209      /* now we can set the mapping from nodes to local DOFs */      /* now we can set the mapping from nodes to local DOFs */
210      this_mapping = Dudley_NodeMapping_alloc(numNodes, nodeMask, UNUSED);      this_mapping = Dudley_NodeMapping_alloc(numNodes, nodeMask, UNUSED);
211      /* define how to get DOF values for controlled bu other processors */      /* define how to get DOF values for controlled bu other processors */
212  #ifdef BOUNDS_CHECK  #ifdef BOUNDS_CHECK
213      for (i = 0; i < offsetInShared[numNeighbors]; ++i)      for (i = 0; i < offsetInShared[numNeighbors]; ++i)
214      {      {
215          if (i < 0 || i >= numNodes * (p_max - p_min + 1))          if (i < 0 || i >= numNodes * (p_max - p_min + 1))
216          {          {
217          printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);              printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
218          exit(1);              exit(1);
219          }          }
220      }      }
221  #endif  #endif
222  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
223      for (i = 0; i < offsetInShared[numNeighbors]; ++i)      for (i = 0; i < offsetInShared[numNeighbors]; ++i)
224          shared[i] = myLastDOF - myFirstDOF + i;          shared[i] = myLastDOF - myFirstDOF + i;
225    
226      rcv_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,      rcv_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,
227                  numNeighbors, neighbor, shared, offsetInShared, 1, 0,              numNeighbors, neighbor, shared, offsetInShared, 1, 0,
228                  mpi_info));              mpi_info));
229    
230      /*      /*
231       *    now we build the sender       *    now we build the sender
232       */       */
233  #ifdef ESYS_MPI  #ifdef ESYS_MPI
234      MPI_Alltoall(rcv_len, 1, MPI_INT, snd_len, 1, MPI_INT, mpi_info->comm);      MPI_Alltoall(rcv_len, 1, MPI_INT, snd_len, 1, MPI_INT, mpi_info->comm);
235  #else  #else
236      for (p = 0; p < mpiSize; ++p)      for (p = 0; p < mpiSize; ++p)
237          snd_len[p] = rcv_len[p];          snd_len[p] = rcv_len[p];
238  #endif  #endif
239      count = 0;      count = 0;
240      for (p = 0; p < rcv_shcomp->numNeighbors; p++)      for (p = 0; p < rcv_shcomp->numNeighbors; p++)
241      {      {
242  #ifdef ESYS_MPI  #ifdef ESYS_MPI
243          MPI_Isend(&(wanted_DOFs[rcv_shcomp->offsetInShared[p]]),          MPI_Isend(&(wanted_DOFs[rcv_shcomp->offsetInShared[p]]),
244                rcv_shcomp->offsetInShared[p + 1] - rcv_shcomp->offsetInShared[p], MPI_INT,                    rcv_shcomp->offsetInShared[p + 1] - rcv_shcomp->offsetInShared[p], MPI_INT,
245                rcv_shcomp->neighbor[p], mpi_info->msg_tag_counter + myRank, mpi_info->comm,                    rcv_shcomp->neighbor[p], mpi_info->counter() + myRank, mpi_info->comm,
246                &mpi_requests[count]);                    &mpi_requests[count]);
247  #endif  #endif
248          count++;          count++;
249      }      }
250      n = 0;      n = 0;
251      numNeighbors = 0;      numNeighbors = 0;
252      for (p = 0; p < mpiSize; p++)      for (p = 0; p < mpiSize; p++)
253      {      {
254          if (snd_len[p] > 0)          if (snd_len[p] > 0)
255          {          {
256  #ifdef ESYS_MPI  #ifdef ESYS_MPI
257          MPI_Irecv(&(shared[n]), snd_len[p],              MPI_Irecv(&(shared[n]), snd_len[p],
258                MPI_INT, p, mpi_info->msg_tag_counter + p, mpi_info->comm, &mpi_requests[count]);                        MPI_INT, p, mpi_info->counter() + p, mpi_info->comm, &mpi_requests[count]);
259  #endif  #endif
260          count++;              count++;
261          neighbor[numNeighbors] = p;              neighbor[numNeighbors] = p;
262          offsetInShared[numNeighbors] = n;              offsetInShared[numNeighbors] = n;
263          numNeighbors++;              numNeighbors++;
264          n += snd_len[p];              n += snd_len[p];
265          }          }
266      }      }
267      mpi_info->incCounter(mpi_info->size);      offsetInShared[numNeighbors] = n;
     offsetInShared[numNeighbors] = n;  
268  #ifdef ESYS_MPI  #ifdef ESYS_MPI
269      MPI_Waitall(count, mpi_requests, mpi_stati);      mpi_info->incCounter(mpi_info->size);
270        MPI_Waitall(count, mpi_requests, mpi_stati);
271  #endif  #endif
272      /* map global ids to local id's */      /* map global ids to local id's */
273  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
274      for (i = 0; i < offsetInShared[numNeighbors]; ++i)      for (i = 0; i < offsetInShared[numNeighbors]; ++i)
275      {      {
276          shared[i] = locDOFMask[shared[i] - min_DOF];          shared[i] = locDOFMask[shared[i] - min_DOF];
     }  
   
     snd_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,  
                 numNeighbors, neighbor, shared, offsetInShared, 1, 0,  
                 dof_distribution->mpi_info));  
   
     if (Dudley_noError())  
         this_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));  
     /* assign new DOF labels to nodes */  
277      }      }
278    
279        snd_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,
280                numNeighbors, neighbor, shared, offsetInShared, 1, 0,
281                dof_distribution->mpi_info));
282    
283        this_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
284    
285        /* assign new DOF labels to nodes */
286      delete[] rcv_len;      delete[] rcv_len;
287      delete[] snd_len;      delete[] snd_len;
288      delete[] mpi_requests;      delete[] mpi_requests;
# Line 303  void Dudley_Mesh_createDOFMappingAndCoup Line 293  void Dudley_Mesh_createDOFMappingAndCoup
293      delete[] shared;      delete[] shared;
294      delete[] offsetInShared;      delete[] offsetInShared;
295      delete[] locDOFMask;      delete[] locDOFMask;
296      if (Dudley_noError())      if (use_reduced_elements)
297      {      {
298      if (use_reduced_elements)          in->Nodes->reducedDegreesOfFreedomMapping = this_mapping;
299      {          in->Nodes->reducedDegreesOfFreedomConnector = this_connector;
         in->Nodes->reducedDegreesOfFreedomMapping = this_mapping;  
         in->Nodes->reducedDegreesOfFreedomConnector = this_connector;  
     }  
     else  
     {  
         in->Nodes->degreesOfFreedomMapping = this_mapping;  
         in->Nodes->degreesOfFreedomConnector = this_connector;  
     }  
300      }      }
301      else      else
302      {      {
303      Dudley_NodeMapping_free(this_mapping);          in->Nodes->degreesOfFreedomMapping = this_mapping;
304            in->Nodes->degreesOfFreedomConnector = this_connector;
305      }      }
306  }  }
307    
308  void Dudley_Mesh_createMappings(Dudley_Mesh * mesh, index_t * dof_distribution, index_t * node_distribution)  void Dudley_Mesh_createMappings(Dudley_Mesh* mesh, index_t* dof_distribution, index_t* node_distribution)
309  {  {
310      int i;      int i;
311      index_t *maskReducedNodes = NULL, *indexReducedNodes = NULL;      index_t *maskReducedNodes = NULL, *indexReducedNodes = NULL;
312      dim_t numReducedNodes;      dim_t numReducedNodes;
313    
314      maskReducedNodes = new  index_t[mesh->Nodes->numNodes];      maskReducedNodes = new index_t[mesh->Nodes->numNodes];
315      indexReducedNodes = new  index_t[mesh->Nodes->numNodes];      indexReducedNodes = new index_t[mesh->Nodes->numNodes];
316    
     if (!(Dudley_checkPtr(maskReducedNodes) || Dudley_checkPtr(indexReducedNodes)))  
     {  
317  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
318      for (i = 0; i < mesh->Nodes->numNodes; ++i)      for (i = 0; i < mesh->Nodes->numNodes; ++i)
319          maskReducedNodes[i] = -1;          maskReducedNodes[i] = -1;
320      Dudley_Mesh_markNodes(maskReducedNodes, 0, mesh, TRUE);      Dudley_Mesh_markNodes(maskReducedNodes, 0, mesh, true);
321    
322      numReducedNodes = Dudley_Util_packMask(mesh->Nodes->numNodes, maskReducedNodes, indexReducedNodes);      numReducedNodes = Dudley_Util_packMask(mesh->Nodes->numNodes, maskReducedNodes, indexReducedNodes);
323      if (Dudley_noError())      Dudley_Mesh_createNodeFileMappings(mesh, numReducedNodes,
324          Dudley_Mesh_createNodeFileMappings(mesh, numReducedNodes, indexReducedNodes, dof_distribution,              indexReducedNodes, dof_distribution, node_distribution);
                            node_distribution);  
     }  
   
325      delete[] maskReducedNodes;      delete[] maskReducedNodes;
326      delete[] indexReducedNodes;      delete[] indexReducedNodes;
327  }  }
328    
329  void Dudley_Mesh_createNodeFileMappings(Dudley_Mesh * in, dim_t numReducedNodes, index_t * indexReducedNodes,  void Dudley_Mesh_createNodeFileMappings(Dudley_Mesh* in, dim_t numReducedNodes,
330                      index_t * dof_first_component, index_t * nodes_first_component)          index_t* indexReducedNodes, index_t* dof_first_component, index_t* nodes_first_component)
331  {  {
332    
333      index_t myFirstDOF, myLastDOF, myFirstNode, myLastNode, *reduced_dof_first_component = NULL, *nodeMask = NULL,      index_t myFirstDOF, myLastDOF, myFirstNode, myLastNode, *reduced_dof_first_component = NULL, *nodeMask = NULL,
334      *reduced_nodes_first_component = NULL, k, *maskMyReducedDOF = NULL, *indexMyReducedDOF =          *reduced_nodes_first_component = NULL, k, *maskMyReducedDOF = NULL, *indexMyReducedDOF =
335      NULL, *maskMyReducedNodes = NULL, *indexMyReducedNodes = NULL;          NULL, *maskMyReducedNodes = NULL, *indexMyReducedNodes = NULL;
336      dim_t myNumDOF, myNumNodes, myNumReducedNodes, myNumReducedDOF, globalNumReducedNodes, globalNumReducedDOF, i,      dim_t myNumDOF, myNumNodes, myNumReducedNodes, myNumReducedDOF, globalNumReducedNodes, globalNumReducedDOF, i,
337      mpiSize;          mpiSize;
338      Esys_MPI_rank myRank;      int myRank;
339    
340      mpiSize = in->Nodes->MPIInfo->size;      mpiSize = in->Nodes->MPIInfo->size;
341      myRank = in->Nodes->MPIInfo->rank;      myRank = in->Nodes->MPIInfo->rank;
# Line 368  void Dudley_Mesh_createNodeFileMappings( Line 345  void Dudley_Mesh_createNodeFileMappings(
345      reduced_dof_first_component = new  index_t[mpiSize + 1];      reduced_dof_first_component = new  index_t[mpiSize + 1];
346      reduced_nodes_first_component = new  index_t[mpiSize + 1];      reduced_nodes_first_component = new  index_t[mpiSize + 1];
347    
348      if (!(Dudley_checkPtr(reduced_dof_first_component) || Dudley_checkPtr(reduced_nodes_first_component)))      myFirstDOF = dof_first_component[myRank];
349      {      myLastDOF = dof_first_component[myRank + 1];
350        myNumDOF = myLastDOF - myFirstDOF;
351      myFirstDOF = dof_first_component[myRank];  
352      myLastDOF = dof_first_component[myRank + 1];      myFirstNode = nodes_first_component[myRank];
353      myNumDOF = myLastDOF - myFirstDOF;      myLastNode = nodes_first_component[myRank + 1];
354        myNumNodes = myLastNode - myFirstNode;
355      myFirstNode = nodes_first_component[myRank];  
356      myLastNode = nodes_first_component[myRank + 1];      maskMyReducedDOF = new index_t[myNumDOF];
357      myNumNodes = myLastNode - myFirstNode;      indexMyReducedDOF = new index_t[myNumDOF];
358        maskMyReducedNodes = new index_t[myNumNodes];
359      maskMyReducedDOF = new  index_t[myNumDOF];      indexMyReducedNodes = new index_t[myNumNodes];
     indexMyReducedDOF = new  index_t[myNumDOF];  
     maskMyReducedNodes = new  index_t[myNumNodes];  
     indexMyReducedNodes = new  index_t[myNumNodes];  
   
     if (!  
         (Dudley_checkPtr(maskMyReducedDOF) || Dudley_checkPtr(indexMyReducedDOF)  
          || Dudley_checkPtr(maskMyReducedNodes) || Dudley_checkPtr(indexMyReducedNodes)))  
     {  
360    
361  #pragma omp parallel private(i)  #pragma omp parallel private(i)
362          {      {
363  #pragma omp for schedule(static)  #pragma omp for schedule(static)
364          for (i = 0; i < myNumNodes; ++i)          for (i = 0; i < myNumNodes; ++i)
365              maskMyReducedNodes[i] = -1;              maskMyReducedNodes[i] = -1;
366  #pragma omp for schedule(static)  #pragma omp for schedule(static)
367          for (i = 0; i < myNumDOF; ++i)          for (i = 0; i < myNumDOF; ++i)
368              maskMyReducedDOF[i] = -1;              maskMyReducedDOF[i] = -1;
369  #pragma omp for private(k) schedule(static)  #pragma omp for private(k) schedule(static)
370          for (i = 0; i < numReducedNodes; ++i)          for (i = 0; i < numReducedNodes; ++i)
371          {          {
372              k = in->Nodes->globalNodesIndex[indexReducedNodes[i]];              k = in->Nodes->globalNodesIndex[indexReducedNodes[i]];
373              if ((k >= myFirstNode) && (myLastNode > k))              if ((k >= myFirstNode) && (myLastNode > k))
374              maskMyReducedNodes[k - myFirstNode] = i;                  maskMyReducedNodes[k - myFirstNode] = i;
375              k = in->Nodes->globalDegreesOfFreedom[indexReducedNodes[i]];              k = in->Nodes->globalDegreesOfFreedom[indexReducedNodes[i]];
376              if ((k >= myFirstDOF) && (myLastDOF > k))              if ((k >= myFirstDOF) && (myLastDOF > k))
377              {              {
378              maskMyReducedDOF[k - myFirstDOF] = i;                  maskMyReducedDOF[k - myFirstDOF] = i;
379              }              }
380          }          }
381          }      }
382          myNumReducedNodes = Dudley_Util_packMask(myNumNodes, maskMyReducedNodes, indexMyReducedNodes);      myNumReducedNodes = Dudley_Util_packMask(myNumNodes, maskMyReducedNodes, indexMyReducedNodes);
383          myNumReducedDOF = Dudley_Util_packMask(myNumDOF, maskMyReducedDOF, indexMyReducedDOF);      myNumReducedDOF = Dudley_Util_packMask(myNumDOF, maskMyReducedDOF, indexMyReducedDOF);
384    
385  #ifdef ESYS_MPI  #ifdef ESYS_MPI
386          MPI_Allgather(&myNumReducedNodes, 1, MPI_INT, reduced_nodes_first_component, 1, MPI_INT,      MPI_Allgather(&myNumReducedNodes, 1, MPI_INT, reduced_nodes_first_component, 1, MPI_INT,
387                in->Nodes->MPIInfo->comm);                    in->Nodes->MPIInfo->comm);
388          MPI_Allgather(&myNumReducedDOF, 1, MPI_INT, reduced_dof_first_component, 1, MPI_INT,      MPI_Allgather(&myNumReducedDOF, 1, MPI_INT, reduced_dof_first_component, 1, MPI_INT,
389                in->Nodes->MPIInfo->comm);                    in->Nodes->MPIInfo->comm);
390  #else  #else
391          reduced_nodes_first_component[0] = myNumReducedNodes;      reduced_nodes_first_component[0] = myNumReducedNodes;
392          reduced_dof_first_component[0] = myNumReducedDOF;      reduced_dof_first_component[0] = myNumReducedDOF;
393  #endif  #endif
394          globalNumReducedNodes = 0;      globalNumReducedNodes = 0;
395          globalNumReducedDOF = 0;      globalNumReducedDOF = 0;
396          for (i = 0; i < mpiSize; ++i)      for (i = 0; i < mpiSize; ++i)
397          {      {
398          k = reduced_nodes_first_component[i];          k = reduced_nodes_first_component[i];
399          reduced_nodes_first_component[i] = globalNumReducedNodes;          reduced_nodes_first_component[i] = globalNumReducedNodes;
400          globalNumReducedNodes += k;          globalNumReducedNodes += k;
401    
402          k = reduced_dof_first_component[i];          k = reduced_dof_first_component[i];
403          reduced_dof_first_component[i] = globalNumReducedDOF;          reduced_dof_first_component[i] = globalNumReducedDOF;
404          globalNumReducedDOF += k;          globalNumReducedDOF += k;
405          }      }
406          reduced_nodes_first_component[mpiSize] = globalNumReducedNodes;      reduced_nodes_first_component[mpiSize] = globalNumReducedNodes;
407          reduced_dof_first_component[mpiSize] = globalNumReducedDOF;      reduced_dof_first_component[mpiSize] = globalNumReducedDOF;
408          /* ==== distribution of Nodes =============================== */      /* ==== distribution of Nodes =============================== */
409          in->Nodes->nodesDistribution.reset(new paso::Distribution(in->Nodes->MPIInfo, nodes_first_component, 1, 0));      in->Nodes->nodesDistribution.reset(new paso::Distribution(in->Nodes->MPIInfo, nodes_first_component, 1, 0));
410    
411          /* ==== distribution of DOFs =============================== */      /* ==== distribution of DOFs =============================== */
412          in->Nodes->degreesOfFreedomDistribution.reset(      in->Nodes->degreesOfFreedomDistribution.reset(
413                      new paso::Distribution(in->Nodes->MPIInfo, dof_first_component, 1, 0));              new paso::Distribution(in->Nodes->MPIInfo, dof_first_component, 1, 0));
414    
415          /* ==== distribution of reduced Nodes =============================== */      /* ==== distribution of reduced Nodes =============================== */
416          in->Nodes->reducedNodesDistribution.reset(new paso::Distribution(      in->Nodes->reducedNodesDistribution.reset(new paso::Distribution(
417                      in->Nodes->MPIInfo, reduced_nodes_first_component, 1, 0));              in->Nodes->MPIInfo, reduced_nodes_first_component, 1, 0));
418    
419          /* ==== distribution of reduced DOF =============================== */      /* ==== distribution of reduced DOF =============================== */
420          in->Nodes->reducedDegreesOfFreedomDistribution.reset(      in->Nodes->reducedDegreesOfFreedomDistribution.reset(
421                      new paso::Distribution(in->Nodes->MPIInfo, reduced_dof_first_component, 1, 0));              new paso::Distribution(in->Nodes->MPIInfo, reduced_dof_first_component, 1, 0));
422      }  
423      delete[] maskMyReducedDOF;      delete[] maskMyReducedDOF;
424      delete[] indexMyReducedDOF;      delete[] indexMyReducedDOF;
425      delete[] maskMyReducedNodes;      delete[] maskMyReducedNodes;
426      delete[] indexMyReducedNodes;      delete[] indexMyReducedNodes;
     }  
427      delete[] reduced_dof_first_component;      delete[] reduced_dof_first_component;
428      delete[] reduced_nodes_first_component;      delete[] reduced_nodes_first_component;
429    
430      nodeMask = new  index_t[in->Nodes->numNodes];      nodeMask = new index_t[in->Nodes->numNodes];
431      if (!Dudley_checkPtr(nodeMask) && Dudley_noError())      /* ==== nodes mapping which is a dummy structure ======== */
     {  
   
     /* ==== nodes mapping which is a dummy structure ======== */  
432  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
433      for (i = 0; i < in->Nodes->numNodes; ++i)      for (i = 0; i < in->Nodes->numNodes; ++i)
434          nodeMask[i] = i;          nodeMask[i] = i;
435      in->Nodes->nodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);      in->Nodes->nodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
436    
437      /* ==== mapping between nodes and reduced nodes ========== */      /* ==== mapping between nodes and reduced nodes ========== */
438  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
439      for (i = 0; i < in->Nodes->numNodes; ++i)      for (i = 0; i < in->Nodes->numNodes; ++i)
440          nodeMask[i] = UNUSED;          nodeMask[i] = UNUSED;
441  #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
442      for (i = 0; i < numReducedNodes; ++i)      for (i = 0; i < numReducedNodes; ++i)
443          nodeMask[indexReducedNodes[i]] = i;          nodeMask[indexReducedNodes[i]] = i;
444      in->Nodes->reducedNodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);      in->Nodes->reducedNodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
     }  
445      delete[] nodeMask;      delete[] nodeMask;
446      /* ==== mapping between nodes and DOFs + DOF connector ========== */      /* ==== mapping between nodes and DOFs + DOF connector ========== */
447      if (Dudley_noError())      Dudley_Mesh_createDOFMappingAndCoupling(in, false);
448          Dudley_Mesh_createDOFMappingAndCoupling(in, FALSE);      /* == mapping between nodes and reduced DOFs + reduced DOF connector == */
449      /* ==== mapping between nodes and reduced DOFs + reduced DOF connector ========== */      Dudley_Mesh_createDOFMappingAndCoupling(in, true);
     if (Dudley_noError())  
         Dudley_Mesh_createDOFMappingAndCoupling(in, TRUE);  
450    
451      /* get the Ids for DOFs and reduced nodes */      /* get the Ids for DOFs and reduced nodes */
     if (Dudley_noError())  
     {  
452  #pragma omp parallel private(i)  #pragma omp parallel private(i)
453      {      {
454  #pragma omp for  #pragma omp for
455          for (i = 0; i < in->Nodes->reducedNodesMapping->numTargets; ++i)          for (i = 0; i < in->Nodes->reducedNodesMapping->numTargets; ++i)
456          in->Nodes->reducedNodesId[i] = in->Nodes->Id[in->Nodes->reducedNodesMapping->map[i]];              in->Nodes->reducedNodesId[i] = in->Nodes->Id[in->Nodes->reducedNodesMapping->map[i]];
457  #pragma omp for  #pragma omp for
458          for (i = 0; i < in->Nodes->degreesOfFreedomMapping->numTargets; ++i)          for (i = 0; i < in->Nodes->degreesOfFreedomMapping->numTargets; ++i)
459          in->Nodes->degreesOfFreedomId[i] = in->Nodes->Id[in->Nodes->degreesOfFreedomMapping->map[i]];              in->Nodes->degreesOfFreedomId[i] = in->Nodes->Id[in->Nodes->degreesOfFreedomMapping->map[i]];
460  #pragma omp for  #pragma omp for
461          for (i = 0; i < in->Nodes->reducedDegreesOfFreedomMapping->numTargets; ++i)          for (i = 0; i < in->Nodes->reducedDegreesOfFreedomMapping->numTargets; ++i)
462          in->Nodes->reducedDegreesOfFreedomId[i] =              in->Nodes->reducedDegreesOfFreedomId[i] =
463              in->Nodes->Id[in->Nodes->reducedDegreesOfFreedomMapping->map[i]];                  in->Nodes->Id[in->Nodes->reducedDegreesOfFreedomMapping->map[i]];
     }  
     }  
     else  
     {  
     Dudley_NodeMapping_free(in->Nodes->nodesMapping);  
     Dudley_NodeMapping_free(in->Nodes->reducedNodesMapping);  
     Dudley_NodeMapping_free(in->Nodes->degreesOfFreedomMapping);  
     Dudley_NodeMapping_free(in->Nodes->reducedDegreesOfFreedomMapping);  
     in->Nodes->nodesDistribution.reset();  
     in->Nodes->reducedNodesDistribution.reset();  
     in->Nodes->degreesOfFreedomDistribution.reset();  
     in->Nodes->reducedDegreesOfFreedomDistribution.reset();  
     in->Nodes->degreesOfFreedomConnector.reset();  
     in->Nodes->reducedDegreesOfFreedomConnector.reset();  
     in->Nodes->nodesMapping = NULL;  
     in->Nodes->reducedNodesMapping = NULL;  
     in->Nodes->degreesOfFreedomMapping = NULL;  
     in->Nodes->reducedDegreesOfFreedomMapping = NULL;  
464      }      }
465  }  }
466    
467    } // namespace dudley
468    

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

  ViewVC Help
Powered by ViewVC 1.1.26