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

Contents of /branches/domexper/dudley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_findMatchingFaces.c
File MIME type: text/plain
File size: 8191 byte(s)
*** empty log message ***

1 /**************************************************************/
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=TMPMEMALLOC(NN*numDim*faces->numElements,double);
63 center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
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.2 2004/12/14 05:39:30 jgs
175 * *** empty log message ***
176 *
177 * Revision 1.1.1.1.2.1 2004/11/24 01:37:14 gross
178 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
179 *
180 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
181 * initial import of project esys2
182 *
183 * Revision 1.2 2004/07/02 04:21:13 gross
184 * Finley C code has been included
185 *
186 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
187 * Initial version of eys using boost-python.
188 *
189 *
190 */
191

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26