/[escript]/branches/domexper/dudley/src/Mesh_resolveNodeIds.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Mesh_resolveNodeIds.c

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

trunk/esys2/finley/src/finleyC/Mesh_resolveNodeIds.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC branches/domexper/dudley/src/Mesh_resolveNodeIds.c revision 3114 by jfenwick, Fri Aug 27 05:26:25 2010 UTC
# Line 1  Line 1 
 /**************************************************************/  
1    
2  /*   Finley: Mesh */  /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
 /*   at input the element nodes refers to the numbering defined the Id assigned to the nodes in the */  
 /*   NodeFile. At the output, the numbering of the element nodes is between 0 and numNodes */  
 /*   degreesOfFreedom are not neccessarily referening to a dense numbering */  
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Copyrights by ACcESS Australia 2003 */  /*   Dudley: Mesh */
18  /*   Author: gross@access.edu.au */  
19  /*   Version: $Id$ */  /*   at input the element nodes refers to the numbering defined the global Id assigned to the nodes in the */
20    /*   NodeFile. It is also not ensured that all nodes refered by an element is actually available */
21    /*   on the process.  At the output, a local node labeling is used and all nodes are available */
22    /*   In particular the numbering of the element nodes is between 0 and in->NodeFile->numNodes */
23    /*   The function does not create a distribution of the degrees of freedom. */
24    
25  /**************************************************************/  /**************************************************************/
26    
 #include "Finley.h"  
27  #include "Mesh.h"  #include "Mesh.h"
28    #include "Util.h"
29    
30  /**************************************************************/  /**************************************************************/
31    
32  void  Finley_Mesh_resolveNodeIds(Finley_Mesh* in) {  void  Dudley_Mesh_resolveNodeIds(Dudley_Mesh* in) {
   maybelong k,min_id,max_id,min_id2,max_id2,len,numDim,newNumNodes,n;  
   maybelong *maskNodes=NULL,*maskElements=NULL,*index=NULL;  
   Finley_NodeFile *newNodeFile=NULL;  
   Finley_ErrorCode=NO_ERROR;  
   numDim=Finley_Mesh_getDim(in);  
     
   /*   find the minimum and maximum id used: */  
     
   min_id=MAYBELONG_MAX;  
   max_id=-MAYBELONG_MAX;  
   Finley_NodeFile_setIdRange(&min_id2,&max_id2,in->Nodes);  
   if (min_id2==MAYBELONG_MAX || max_id2==-MAYBELONG_MAX) {  
     Finley_ErrorCode=VALUE_ERROR;  
     sprintf(Finley_ErrorMsg,"Mesh has not been defined completely.");  
     goto clean;  
   }  
33    
34      index_t min_id, max_id,  min_id2, max_id2, global_min_id, global_max_id,
35              *globalToNewLocalNodeLabels=NULL, *newLocalToGlobalNodeLabels=NULL;
36      dim_t len, n, newNumNodes, numDim;
37      Dudley_NodeFile *newNodeFile=NULL;
38      #ifdef PASO_MPI
39      index_t id_range[2], global_id_range[2];
40      #endif
41      numDim=Dudley_Mesh_getDim(in);
42      /*  find the minimum and maximum id used by elements: */
43      min_id=INDEX_T_MAX;
44      max_id=-INDEX_T_MAX;
45    //printf("Trying Elements:\n");
46      Dudley_ElementFile_setNodeRange(&min_id2,&max_id2,in->Elements);
47    max_id=MAX(max_id,max_id2);    max_id=MAX(max_id,max_id2);
48    min_id=MIN(min_id,min_id2);    min_id=MIN(min_id,min_id2);
49    Finley_ElementFile_setNodeRange(&min_id2,&max_id2,in->Elements);  //printf("Trying FaceElements:\n");
50      Dudley_ElementFile_setNodeRange(&min_id2,&max_id2,in->FaceElements);
51    max_id=MAX(max_id,max_id2);    max_id=MAX(max_id,max_id2);
52    min_id=MIN(min_id,min_id2);    min_id=MIN(min_id,min_id2);
53    Finley_ElementFile_setNodeRange(&min_id2,&max_id2,in->FaceElements);  //printf("Trying ContactElements:\n");
54      Dudley_ElementFile_setNodeRange(&min_id2,&max_id2,in->ContactElements);
55    max_id=MAX(max_id,max_id2);    max_id=MAX(max_id,max_id2);
56    min_id=MIN(min_id,min_id2);    min_id=MIN(min_id,min_id2);
57    Finley_ElementFile_setNodeRange(&min_id2,&max_id2,in->ContactElements);    Dudley_ElementFile_setNodeRange(&min_id2,&max_id2,in->Points);
58    max_id=MAX(max_id,max_id2);    max_id=MAX(max_id,max_id2);
59    min_id=MIN(min_id,min_id2);    min_id=MIN(min_id,min_id2);
60    Finley_ElementFile_setNodeRange(&min_id2,&max_id2,in->Points);    #ifdef PASO_MPI
61    max_id=MAX(max_id,max_id2);       id_range[0]=-min_id;
62    min_id=MIN(min_id,min_id2);       id_range[1]=max_id;
63    #ifdef Finley_TRACE       MPI_Allreduce( id_range, global_id_range, 2, MPI_INT, MPI_MAX, in->MPIInfo->comm );
64    printf("Node id range is %d:%d\n",min_id,max_id);       global_min_id=-global_id_range[0];
65         global_max_id=global_id_range[1];
66      #else
67         global_min_id=min_id;
68         global_max_id=max_id;
69    #endif    #endif
70        #ifdef Dudley_TRACE
71    /*  allocate a new node file used to gather existing node file: */    printf("Node id range used by elements is %d:%d\n",global_min_id,global_max_id);
72        #endif
73    len=max_id-min_id+1;    if (min_id>max_id) {
74    newNodeFile=Finley_NodeFile_alloc(numDim);       max_id=-1;
75    if (Finley_ErrorCode!=NO_ERROR) goto clean;       min_id=0;
   
   maskNodes=(maybelong*)TMPMEMALLOC(len*sizeof(maybelong));  
   if (Finley_checkPtr(maskNodes)) goto clean;  
   
   maskElements=(maybelong*)TMPMEMALLOC(len*sizeof(maybelong));  
   if (Finley_checkPtr(maskElements)) goto clean;  
   
   index=(maybelong*)TMPMEMALLOC(in->Nodes->numNodes*sizeof(maybelong));  
   if (Finley_checkPtr(maskElements)) goto clean;  
   
   Finley_NodeFile_allocTable(newNodeFile,len);  
   if (Finley_ErrorCode!=NO_ERROR) goto clean;  
   
   #pragma omp parallel for private(n) schedule(static)  
   for (n=0;n<in->Nodes->numNodes;n++) index[n]=-1;  
   #pragma omp parallel for private(n) schedule(static)  
   for (n=0;n<len;n++) {  
          maskNodes[n]=-1;  
          maskElements[n]=-1;  
76    }    }
   /*  mark the nodes referred by elements in maskElements: */  
77        
78    Finley_Mesh_markNodes(maskElements,min_id,in,FALSE);    /* allocate mappings for new local node labeling to global node labeling (newLocalToGlobalNodeLabels)
79         and global node labeling to the new local node labeling (globalToNewLocalNodeLabels[i-min_id] is the
80    /*  mark defined nodes */       new local id of global node i) */
81      len=(max_id>=min_id) ? max_id-min_id+1 : 0 ;
82    #pragma omp parallel for private(k) schedule(static)    globalToNewLocalNodeLabels=TMPMEMALLOC(len,index_t); /* local mask for used nodes */
83    for (k=0;k<in->Nodes->numNodes;k++) maskNodes[in->Nodes->Id[k]-min_id]=1;    newLocalToGlobalNodeLabels=TMPMEMALLOC(len,index_t);
84      if (! ( (Dudley_checkPtr(globalToNewLocalNodeLabels) && Dudley_checkPtr(newLocalToGlobalNodeLabels) ) ) ) {
85    /*  check if all referenced nodes are actually defined: */  
86           #pragma omp parallel
87    #pragma omp parallel for private(k) schedule(static)         {
88    for (k=0;k<len;k++) {             #pragma omp for private(n) schedule(static)
89       /* if a node is refered by an element is there a node defined ?*/             for (n=0;n<len;n++) newLocalToGlobalNodeLabels[n]=-1;
90       if (maskElements[k]>=0 && maskNodes[k]<0) {             #pragma omp for private(n) schedule(static)
91         Finley_ErrorCode=VALUE_ERROR;             for (n=0;n<len;n++) globalToNewLocalNodeLabels[n]=-1;
92         sprintf(Finley_ErrorMsg,"Node id %d is referenced by element but is not defined.",k+min_id);         }
93       }  
94           /*  mark the nodes referred by elements in globalToNewLocalNodeLabels which is currently used as a mask: */
95           Dudley_Mesh_markNodes(globalToNewLocalNodeLabels,min_id,in,FALSE);
96    
97           /* create a local labeling newLocalToGlobalNodeLabels of the local nodes by packing the mask globalToNewLocalNodeLabels*/
98    
99           newNumNodes=Dudley_Util_packMask(len,globalToNewLocalNodeLabels,newLocalToGlobalNodeLabels);
100    
101           /* invert the new labeling and shift the index newLocalToGlobalNodeLabels to global node ids */
102           #pragma omp parallel for private(n) schedule(static)
103           for (n=0;n<newNumNodes;n++) {
104                  #ifdef BOUNDS_CHECK
105                         if (n >= len || n < 0) { printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n); exit(1); }
106                         if (newLocalToGlobalNodeLabels[n] >= len || newLocalToGlobalNodeLabels[n] < 0) { printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n); exit(1); }
107                  #endif
108                  globalToNewLocalNodeLabels[newLocalToGlobalNodeLabels[n]]=n;
109                  newLocalToGlobalNodeLabels[n]+=min_id;
110            }
111            /* create a new table */
112            newNodeFile=Dudley_NodeFile_alloc(numDim,in->MPIInfo);
113            if (Dudley_noError()) {
114               Dudley_NodeFile_allocTable(newNodeFile,newNumNodes);
115            }
116            if (Dudley_noError()) {
117                Dudley_NodeFile_gather_global(newLocalToGlobalNodeLabels,in->Nodes, newNodeFile);
118            }
119            if (Dudley_noError()) {
120               Dudley_NodeFile_free(in->Nodes);
121               in->Nodes=newNodeFile;
122               /*  relable nodes of the elements: */
123               Dudley_Mesh_relableElementNodes(globalToNewLocalNodeLabels,min_id,in);
124            }
125    }    }
126      TMPMEMFREE(globalToNewLocalNodeLabels);
127    if (Finley_ErrorCode==NO_ERROR) {    TMPMEMFREE(newLocalToGlobalNodeLabels);
128        if (! Dudley_noError()) {
129        /*  scatter the nodefile in->nodes into newNodeFile using index; */         Dudley_NodeFile_free(newNodeFile);
       #pragma omp parallel for private(k) schedule(static)  
       for (k=0;k<in->Nodes->numNodes;k++) index[k]=in->Nodes->Id[k]-min_id;  
       Finley_NodeFile_scatter(index,in->Nodes,newNodeFile);  
     
       /*  relable used nodes: */  
       /* index maps the new node labeling onto the old one */  
       newNumNodes=Finley_Util_packMask(len,maskElements,index);  
       #pragma omp parallel for private(k) schedule(static)  
       for (k=0;k<newNumNodes;k++) maskElements[index[k]]=k;  
   
       /*  create a new table of nodes: */  
       Finley_NodeFile_deallocTable(in->Nodes);  
       Finley_NodeFile_allocTable(in->Nodes,newNumNodes);  
       if (Finley_ErrorCode!=NO_ERROR) goto clean;  
   
       /* gather the new nodefile into in->Nodes */  
       Finley_NodeFile_gather(index,newNodeFile,in->Nodes);  
   
       /*  relable nodes of the elements: */  
       Finley_Mesh_relableElementNodes(maskElements,min_id,in);  
   
130    }    }
   
   /*  clean-up: */  
     
   clean: TMPMEMFREE(maskNodes);  
          TMPMEMFREE(maskElements);  
          TMPMEMFREE(index);  
          Finley_NodeFile_deallocTable(newNodeFile);  
          Finley_NodeFile_dealloc(newNodeFile);  
131  }  }
   
 /*  
 * $Log$  
 * Revision 1.1  2004/10/26 06:53:57  jgs  
 * Initial revision  
 *  
 * Revision 1.1.1.1  2004/06/24 04:00:40  johng  
 * Initial version of eys using boost-python.  
 *  
 *  
 */  
   

Legend:
Removed from v.82  
changed lines
  Added in v.3114

  ViewVC Help
Powered by ViewVC 1.1.26