/[escript]/branches/domexper/dudley/src/Mesh_findMatchingFaces.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3086 - (hide annotations)
Thu Aug 5 05:07:58 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 8606 byte(s)
Another pass at removing finley

1 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 jfenwick 3086 /* Dudley: Mesh */
18 jgs 82
19     /* searches for faces in the mesh which are matching */
20    
21     /**************************************************************/
22    
23     #include "Util.h"
24     #include "Mesh.h"
25 gross 2748
26 jgs 82 /**************************************************************/
27    
28 jfenwick 3086 static double Dudley_Mesh_lockingGridSize=0;
29 jgs 82
30 jfenwick 3086 int Dudley_Mesh_findMatchingFaces_compar(const void *arg1 , const void *arg2 ) {
31     Dudley_Mesh_findMatchingFaces_center *e1,*e2;
32 jgs 123 bool_t l,g;
33     dim_t i;
34 jfenwick 3086 e1=(Dudley_Mesh_findMatchingFaces_center*) arg1;
35     e2=(Dudley_Mesh_findMatchingFaces_center*) arg2;
36 jgs 82 for (i=0;i<MAX_numDim;i++) {
37 jfenwick 3086 l= (e1->x[i]<e2->x[i]+Dudley_Mesh_lockingGridSize) ? TRUE : FALSE;
38     g= (e2->x[i]<e1->x[i]+Dudley_Mesh_lockingGridSize) ? TRUE : FALSE;
39 jgs 82 if (! (l && g)) {
40     if (l) return -1;
41     if (g) return 1;
42     }
43     }
44 gross 888 if ( e1->refId < e2->refId ) {
45     return -1;
46     } else if ( e1->refId > e2->refId ) {
47 gross 1134 return 1;
48 gross 888 } else {
49     return 0;
50     }
51 jgs 82 }
52    
53 jfenwick 3086 void Dudley_Mesh_findMatchingFaces(Dudley_NodeFile *nodes, Dudley_ElementFile *faces, double safety_factor,double tolerance,
54 jgs 123 dim_t* numPairs, index_t* elem0,index_t* elem1,index_t* matching_nodes_in_elem1) {
55 jgs 82 #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
56 jgs 123 {dim_t i; \
57 jgs 82 _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 jfenwick 3086 Dudley_ReferenceElement* refElement=NULL;
61 jgs 150 char error_msg[LenErrorMsg_MAX];
62 gross 2748 double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
63 jfenwick 3086 Dudley_Mesh_findMatchingFaces_center *center;
64 gross 2748 index_t e_0,e_1,*a1=NULL,*a2=NULL,*perm=NULL,*perm_tmp=NULL,*itmp_ptr=NULL, *faceNodes=NULL, *shiftNodes=NULL, *reverseNodes=NULL ;
65     dim_t e,i,i0,i1,n,NN,numNodesOnFace;
66    
67     dim_t numDim=nodes->numDim;
68 jgs 82
69 jfenwick 3086 refElement= Dudley_ReferenceElementSet_borrowReferenceElement(faces->referenceElementSet, FALSE);
70 gross 2748 NN=faces->numNodes;
71    
72     numNodesOnFace=refElement->Type->numNodesOnFace;
73     faceNodes=refElement->Type->faceNodes;
74     shiftNodes=refElement->Type->shiftNodes;
75     reverseNodes=refElement->Type->reverseNodes;
76    
77     if (numNodesOnFace<=0) {
78 jfenwick 3086 sprintf(error_msg,"Dudley_Mesh_findMatchingFaces: matching faces cannot be applied to face elements of type %s",refElement->Type->Name);
79     Dudley_setError(TYPE_ERROR,error_msg);
80 gross 2748 return;
81     }
82 jgs 102 X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
83 jfenwick 3086 center=TMPMEMALLOC(faces->numElements,Dudley_Mesh_findMatchingFaces_center);
84 gross 1028 a1=TMPMEMALLOC(NN,int);
85     a2=TMPMEMALLOC(NN,int);
86 jfenwick 3086 if (!(Dudley_checkPtr(X) || Dudley_checkPtr(center) || Dudley_checkPtr(a1) || Dudley_checkPtr(a2)) ) {
87 jgs 82 /* OMP */
88     for (e=0;e<faces->numElements;e++) {
89     /* get the coordinates of the nodes */
90 jfenwick 3086 Dudley_Util_Gather_double(NN,&(faces->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,&(X[INDEX3(0,0,e,numDim,NN)]));
91 jgs 82 /* get the element center */
92     center[e].refId=e;
93     for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
94 gross 2748 for (i0=0;i0<numNodesOnFace;i0++) {
95     for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
96 jgs 82 }
97 gross 2748 for (i=0;i<numDim;i++) center[e].x[i]/=numNodesOnFace;
98 jgs 82 /* get the minimum distance between nodes in the element */
99 gross 2748 for (i0=0;i0<numNodesOnFace;i0++) {
100     for (i1=i0+1;i1<numNodesOnFace;i1++) {
101     getDist(h_local,e,faceNodes[i0],e,faceNodes[i1]);
102 jgs 82 h=MIN(h,h_local);
103     }
104     }
105     }
106     /* set the */
107 jfenwick 3086 Dudley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
108     #ifdef Dudley_TRACE
109     printf("locking grid size is %e\n",Dudley_Mesh_lockingGridSize);
110 gross 888 printf("absolute tolerance is %e.\n",h * tolerance);
111 jgs 82 #endif
112     /* sort the elements by center center coordinates (lexigraphical)*/
113 jfenwick 3086 qsort(center,faces->numElements,sizeof(Dudley_Mesh_findMatchingFaces_center),Dudley_Mesh_findMatchingFaces_compar);
114 jgs 82 /* find elements with matching center */
115     *numPairs=0;
116     /* OMP */
117 jfenwick 3086 for (e=0;e<faces->numElements-1 && Dudley_noError();e++) {
118 gross 888 dist=0;
119     for (i=0;i<numDim;i++) dist=MAX(dist,ABS(center[e].x[i]-center[e+1].x[i]));
120     if (dist < h * tolerance) {
121 jgs 82 e_0=center[e].refId;
122     e_1=center[e+1].refId;
123     elem0[*numPairs]=e_0;
124     elem1[*numPairs]=e_1;
125     /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
126     perm=a1;
127     perm_tmp=a2;
128     for (i=0;i<NN;i++) perm[i]=i;
129 jfenwick 3086 while (Dudley_noError()) {
130 jgs 82 /* if node 0 and perm[0] are the same we are ready */
131     getDist(dist,e_0,0,e_1,perm[0]);
132     if (dist<=h*tolerance) break;
133 gross 2748 if (shiftNodes[0]>=0) {
134 jgs 82 /* rotate the nodes */
135 gross 1377 itmp_ptr=perm;
136     perm=perm_tmp;
137     perm_tmp=itmp_ptr;
138     #pragma ivdep
139 gross 2748 for (i=0;i<NN;i++) perm[i]=perm_tmp[shiftNodes[i]];
140 jgs 82 }
141     /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
142     if (perm[0]==0) {
143 gross 888 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match first node of element %d to touching element %d",e_0,e_1);
144 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
145 jgs 82 }
146     }
147     /* now we check if the second nodes match */
148 jfenwick 3086 if (Dudley_noError()) {
149 gross 2748 if (numNodesOnFace>1) {
150     getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
151 jgs 82 /* if the second node does not match we reverse the direction of the nodes */
152     if (dist>h*tolerance) {
153     /* rotate the nodes */
154 gross 2748 if (reverseNodes[0]<0) {
155 gross 888 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
156 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
157 jgs 82 } else {
158 gross 1377 itmp_ptr=perm;
159     perm=perm_tmp;
160     perm_tmp=itmp_ptr;
161     #pragma ivdep
162 gross 2748 for (i=0;i<NN;i++) perm[i]=perm_tmp[reverseNodes[i]];
163     getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
164 jgs 82 if (dist>h*tolerance) {
165 gross 888 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",e_0,e_1);
166 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
167 jgs 82 }
168     }
169     }
170     }
171     }
172     /* we check if the rest of the face nodes match: */
173 jfenwick 3086 if (Dudley_noError()) {
174 gross 2748 for (i=2;i<numNodesOnFace;i++) {
175     n=faceNodes[i];
176 jgs 82 getDist(dist,e_0,n,e_1,perm[n]);
177     if (dist>h*tolerance) {
178 gross 888 sprintf(error_msg,"Mesh_findMatchingFaces:couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
179 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
180 jgs 82 break;
181     }
182     }
183     }
184     /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
185 jfenwick 3086 if (Dudley_noError()) {
186 jgs 82 for (i=0;i<NN;i++) matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
187     }
188     (*numPairs)++;
189     }
190     }
191 jfenwick 3086 #ifdef Dudley_TRACE
192 jgs 82 printf("number of pairs of matching faces %d\n",*numPairs);
193     #endif
194     }
195     /* clean up */
196     TMPMEMFREE(X);
197     TMPMEMFREE(center);
198 gross 1028 TMPMEMFREE(a1);
199     TMPMEMFREE(a1);
200 jgs 82
201     #undef getDist
202     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26