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

Contents of /trunk/finley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1376 - (show annotations)
Wed Jan 9 01:38:18 2008 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 8503 byte(s)
inserted sys.exit(1) into the tests so scons can detect the failure of the test. 
A similar statement has been removed from an earlier as it produces problems on 64bit Linux. Previously exit(0) was called in case of success but now this is not done in order to avoid a fatal end of the program. in the case of an error in the test there could be a fatal error so but I guess that this not really a problem.

PS: the fact that signal 0 was returned even for the case of an error lead to the illusion that all tests have been completed successfully.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26