/[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 753 - (show annotations)
Mon Jun 26 02:51:25 2006 UTC (14 years, 1 month ago) by woo409
File MIME type: text/plain
File size: 9127 byte(s)
+ Made qsortG usage conditional on USE_QSORTG being defined

Note both altix and win32 fail the test_normal_onFunctionContactOne/Zero tests when using QSORTG.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26