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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26