/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: text/plain
File size: 8494 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Finley: Mesh */
18
19 /* searches for faces in the mesh which are matching */
20
21 /**************************************************************/
22
23 #include "Util.h"
24 #include "Mesh.h"
25 /**************************************************************/
26
27 static double Finley_Mesh_lockingGridSize=0;
28
29 int Finley_Mesh_findMatchingFaces_compar(const void *arg1 , const void *arg2 ) {
30 Finley_Mesh_findMatchingFaces_center *e1,*e2;
31 bool_t l,g;
32 dim_t i;
33 e1=(Finley_Mesh_findMatchingFaces_center*) arg1;
34 e2=(Finley_Mesh_findMatchingFaces_center*) arg2;
35 for (i=0;i<MAX_numDim;i++) {
36 l= (e1->x[i]<e2->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
37 g= (e2->x[i]<e1->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
38 if (! (l && g)) {
39 if (l) return -1;
40 if (g) return 1;
41 }
42 }
43 if ( e1->refId < e2->refId ) {
44 return -1;
45 } else if ( e1->refId > e2->refId ) {
46 return 1;
47 } else {
48 return 0;
49 }
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 char error_msg[LenErrorMsg_MAX];
61 double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
62 dim_t NN=faces->ReferenceElement->Type->numNodes;
63 dim_t numDim=nodes->numDim;
64 Finley_Mesh_findMatchingFaces_center *center;
65 index_t e_0,e_1,*a1=NULL,*a2=NULL,*perm=NULL,*perm_tmp=NULL,*itmp_ptr=NULL;
66 dim_t e,i,i0,i1,n;
67
68 X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
69 center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
70 a1=TMPMEMALLOC(NN,int);
71 a2=TMPMEMALLOC(NN,int);
72 if (!(Finley_checkPtr(X) || Finley_checkPtr(center) || Finley_checkPtr(a1) || Finley_checkPtr(a2)) ) {
73 /* OMP */
74 for (e=0;e<faces->numElements;e++) {
75 /* get the coordinates of the nodes */
76 Finley_Util_Gather_double(NN,&(faces->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,&(X[INDEX3(0,0,e,numDim,NN)]));
77 /* get the element center */
78 center[e].refId=e;
79 for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
80 for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
81 for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faces->ReferenceElement->Type->faceNode[i0],e,numDim,NN)];
82 }
83 for (i=0;i<numDim;i++) center[e].x[i]/=faces->ReferenceElement->Type->numNodesOnFace;
84 /* get the minimum distance between nodes in the element */
85 for (i0=0;i0<faces->ReferenceElement->Type->numNodesOnFace;i0++) {
86 for (i1=i0+1;i1<faces->ReferenceElement->Type->numNodesOnFace;i1++) {
87 getDist(h_local,e,faces->ReferenceElement->Type->faceNode[i0],e,faces->ReferenceElement->Type->faceNode[i1]);
88 h=MIN(h,h_local);
89 }
90 }
91 }
92 /* set the */
93 Finley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
94 #ifdef Finley_TRACE
95 printf("locking grid size is %e\n",Finley_Mesh_lockingGridSize);
96 printf("absolute tolerance is %e.\n",h * tolerance);
97 #endif
98 /* sort the elements by center center coordinates (lexigraphical)*/
99 qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
100 /* find elements with matching center */
101 *numPairs=0;
102 /* OMP */
103 for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
104 dist=0;
105 for (i=0;i<numDim;i++) dist=MAX(dist,ABS(center[e].x[i]-center[e+1].x[i]));
106 if (dist < h * tolerance) {
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 itmp_ptr=perm;
122 perm=perm_tmp;
123 perm_tmp=itmp_ptr;
124 #pragma ivdep
125 for (i=0;i<NN;i++) perm[i]=perm_tmp[faces->ReferenceElement->Type->shiftNodes[i]];
126 }
127 /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
128 if (perm[0]==0) {
129 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match first node of element %d to touching element %d",e_0,e_1);
130 Finley_setError(VALUE_ERROR,error_msg);
131 }
132 }
133 /* now we check if the second nodes match */
134 if (Finley_noError()) {
135 if (faces->ReferenceElement->Type->numNodesOnFace>1) {
136 getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
137 /* if the second node does not match we reverse the direction of the nodes */
138 if (dist>h*tolerance) {
139 /* rotate the nodes */
140 if (faces->ReferenceElement->Type->reverseNodes[0]<0) {
141 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
142 Finley_setError(VALUE_ERROR,error_msg);
143 } else {
144 itmp_ptr=perm;
145 perm=perm_tmp;
146 perm_tmp=itmp_ptr;
147 #pragma ivdep
148 for (i=0;i<NN;i++) perm[i]=perm_tmp[faces->ReferenceElement->Type->reverseNodes[i]];
149 getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
150 if (dist>h*tolerance) {
151 sprintf(error_msg,"Mesh_findMatchingFaces: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,"Mesh_findMatchingFaces: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 TMPMEMFREE(a1);
185 TMPMEMFREE(a1);
186
187 #undef getDist
188 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26