/[escript]/temp/finley/src/Mesh_findMatchingFaces.c
ViewVC logotype

Annotation of /temp/finley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_findMatchingFaces.c
File MIME type: text/plain
File size: 8013 byte(s)
Initial revision

1 jgs 82 /**************************************************************/
2    
3     /* Finley: Mesh */
4    
5     /* searches for faces in the mesh which are matching */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2003/04 */
10     /* Author: gross@access.edu.au */
11     /* Version: $Id$ */
12    
13     /**************************************************************/
14    
15     #include "Common.h"
16     #include "Finley.h"
17     #include "Util.h"
18     #include "Mesh.h"
19    
20     /**************************************************************/
21    
22     static double Finley_Mesh_lockingGridSize=0;
23    
24     int Finley_Mesh_findMatchingFaces_compar(const void *arg1 , const void *arg2 ) {
25     Finley_Mesh_findMatchingFaces_center *e1,*e2;
26     int i,l,g;
27     e1=(Finley_Mesh_findMatchingFaces_center*) arg1;
28     e2=(Finley_Mesh_findMatchingFaces_center*) arg2;
29     for (i=0;i<MAX_numDim;i++) {
30     l= (e1->x[i]<e2->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
31     g= (e2->x[i]<e1->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
32     if (! (l && g)) {
33     if (l) return -1;
34     if (g) return 1;
35     }
36     }
37     return 0;
38     }
39    
40     void Finley_Mesh_findMatchingFaces(Finley_NodeFile *nodes, Finley_ElementFile *faces, double safety_factor,double tolerance,
41     int* numPairs, int* elem0,int* elem1,int* matching_nodes_in_elem1) {
42     #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
43     {int i; \
44     _dist_=0; \
45     for (i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)])); \
46     }
47    
48     #define SWAP(_i1_,_i2_) \
49     {int* i; \
50     i=(_i2_); \
51     (_i2_)=(_i1_); \
52     (_i1_)=i; \
53     }
54    
55     int e,i,i0,i1,e_0,e_1;
56     double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
57     int NN=faces->ReferenceElement->Type->numNodes;
58     int numDim=nodes->numDim;
59     int a1[NN],a2[NN],*perm,*perm_tmp,n;
60     Finley_Mesh_findMatchingFaces_center *center;
61    
62     X=(double*) TMPMEMALLOC(sizeof(double)*NN*numDim*faces->numElements);
63     center=(Finley_Mesh_findMatchingFaces_center*) TMPMEMALLOC(sizeof(Finley_Mesh_findMatchingFaces_center)*faces->numElements);
64     if (!(Finley_checkPtr(X) || Finley_checkPtr(center)) ) {
65     /* OMP */
66     for (e=0;e<faces->numElements;e++) {
67     /* get the coordinates of the nodes */
68     Finley_Util_Gather_double(NN,&(faces->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,&(X[INDEX3(0,0,e,numDim,NN)]));
69     /* get the element center */
70     center[e].refId=e;
71     for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
72     for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
73     for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faces->ReferenceElement->Type->faceNode[i0],e,numDim,NN)];
74     }
75     for (i=0;i<numDim;i++) center[e].x[i]/=faces->ReferenceElement->Type->numNodesOnFace;
76     /* get the minimum distance between nodes in the element */
77     for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
78     for (i1=i0+1;i1<faces->ReferenceElement->Type->numNodesOnFace;i1++) {
79     getDist(h_local,e,faces->ReferenceElement->Type->faceNode[i0],e,faces->ReferenceElement->Type->faceNode[i1]);
80     h=MIN(h,h_local);
81     }
82     }
83     }
84     /* set the */
85     Finley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
86     #ifdef Finley_TRACE
87     printf("locking grid size is %e\n",Finley_Mesh_lockingGridSize);
88     #endif
89     /* sort the elements by center center coordinates (lexigraphical)*/
90     qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
91     /* find elements with matching center */
92     *numPairs=0;
93     /* OMP */
94     for (e=0;e<faces->numElements-1 && Finley_ErrorCode==NO_ERROR;e++) {
95     if (Finley_Mesh_findMatchingFaces_compar((void*) &(center[e]),(void*) &(center[e+1]))==0) {
96     e_0=center[e].refId;
97     e_1=center[e+1].refId;
98     elem0[*numPairs]=e_0;
99     elem1[*numPairs]=e_1;
100     /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
101     perm=a1;
102     perm_tmp=a2;
103     for (i=0;i<NN;i++) perm[i]=i;
104     while (Finley_ErrorCode==NO_ERROR) {
105     /* if node 0 and perm[0] are the same we are ready */
106     getDist(dist,e_0,0,e_1,perm[0]);
107     if (dist<=h*tolerance) break;
108     if (faces->ReferenceElement->Type->shiftNodes[0]>=0) {
109     /* rotate the nodes */
110     for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->shiftNodes[i]];
111     SWAP(perm,perm_tmp);
112     }
113     /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
114     if (perm[0]==0) {
115     Finley_ErrorCode=VALUE_ERROR;
116     sprintf(Finley_ErrorMsg,"couldn't match first node of element %d to touching element %d",e_0,e_1);
117     }
118     }
119     /* now we check if the second nodes match */
120     if (Finley_ErrorCode==NO_ERROR) {
121     if (faces->ReferenceElement->Type->numNodesOnFace>1) {
122     getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
123     /* if the second node does not match we reverse the direction of the nodes */
124     if (dist>h*tolerance) {
125     /* rotate the nodes */
126     if (faces->ReferenceElement->Type->reverseNodes[0]<0) {
127     Finley_ErrorCode=VALUE_ERROR;
128     sprintf(Finley_ErrorMsg,"couldn't match the second node of element %d to touching element %d",e_0,e_1);
129     } else {
130     for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->reverseNodes[i]];
131     SWAP(perm,perm_tmp);
132     getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
133     if (dist>h*tolerance) {
134     Finley_ErrorCode=VALUE_ERROR;
135     sprintf(Finley_ErrorMsg,"couldn't match the second node of element %d to touching element %d",e_0,e_1);
136     }
137     }
138     }
139     }
140     }
141     /* we check if the rest of the face nodes match: */
142     if (Finley_ErrorCode==NO_ERROR) {
143     for (i=2;i<faces->ReferenceElement->Type->numNodesOnFace;i++) {
144     n=faces->ReferenceElement->Type->faceNode[i];
145     getDist(dist,e_0,n,e_1,perm[n]);
146     if (dist>h*tolerance) {
147     Finley_ErrorCode=VALUE_ERROR;
148     sprintf(Finley_ErrorMsg,"couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
149     break;
150     }
151     }
152     }
153     /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
154     if (Finley_ErrorCode==NO_ERROR) {
155     for (i=0;i<NN;i++) matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
156     }
157     (*numPairs)++;
158     }
159     }
160     #ifdef Finley_TRACE
161     printf("number of pairs of matching faces %d\n",*numPairs);
162     #endif
163     }
164     /* clean up */
165     TMPMEMFREE(X);
166     TMPMEMFREE(center);
167    
168     #undef getDist
169     #undef SWAP
170     }
171    
172     /*
173     * $Log$
174     * Revision 1.1 2004/10/26 06:53:57 jgs
175     * Initial revision
176     *
177     * Revision 1.2 2004/07/02 04:21:13 gross
178     * Finley C code has been included
179     *
180     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
181     * Initial version of eys using boost-python.
182     *
183     *
184     */
185    

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26