/[escript]/temp/finley/src/Mesh_findMatchingFaces.c
ViewVC logotype

Annotation of /temp/finley/src/Mesh_findMatchingFaces.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years ago) by ksteube
Original Path: trunk/finley/src/Mesh_findMatchingFaces.c
File MIME type: text/plain
File size: 9127 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26