/[escript]/branches/intelc_win32/finley/src/Mesh_findMatchingFaces.c
ViewVC logotype

Contents of /branches/intelc_win32/finley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 752 - (show annotations)
Mon Jun 26 02:25:41 2006 UTC (13 years, 2 months ago) by woo409
File MIME type: text/plain
File size: 8948 byte(s)
+ Added a qsort.c file which contains a drop in replacement for qsort (call it as qsortG). This one appears to be a stable implementation and the test .msh files on windows have been set up to be the same as unix again except for the exponent digits (3 instead of 2).
With ALL the qsorts replaced with qsortG only two tests fail now on win32:
test_normal_onFunctionOnContactOne
test_normal_onFunctionOnContactZero

Both give wrong result errors.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26