/[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 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 11 months ago) by jgs
Original Path: trunk/finley/src/finleyC/Mesh_findMatchingFaces.c
File MIME type: text/plain
File size: 9250 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 jgs 150 /*
2     ******************************************************************************
3     * *
4     * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5     * *
6     * This software is the property of ACcESS. No part of this code *
7     * may be copied in any form or by any means without the expressed written *
8     * consent of ACcESS. Copying, use or modification of this software *
9     * by any unauthorised person is illegal unless that person has a software *
10     * license agreement with ACcESS. *
11     * *
12     ******************************************************************************
13     */
14    
15 jgs 82 /**************************************************************/
16    
17     /* Finley: Mesh */
18    
19     /* searches for faces in the mesh which are matching */
20    
21     /**************************************************************/
22    
23 jgs 150 /* Author: gross@access.edu.au */
24     /* Version: $Id$
25 jgs 82
26     /**************************************************************/
27    
28     #include "Util.h"
29     #include "Mesh.h"
30    
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 jgs 123 bool_t l,g;
38     dim_t i;
39 jgs 82 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 jgs 123 dim_t* numPairs, index_t* elem0,index_t* elem1,index_t* matching_nodes_in_elem1) {
54 jgs 82 #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
55 jgs 123 {dim_t i; \
56 jgs 82 _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 jgs 123 {index_t* i; \
62 jgs 82 i=(_i2_); \
63     (_i2_)=(_i1_); \
64     (_i1_)=i; \
65     }
66    
67 jgs 150 char error_msg[LenErrorMsg_MAX];
68 jgs 82 double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
69 jgs 123 dim_t NN=faces->ReferenceElement->Type->numNodes;
70     dim_t numDim=nodes->numDim;
71 jgs 82 Finley_Mesh_findMatchingFaces_center *center;
72 jgs 123 index_t e_0,e_1,a1[NN],a2[NN],*perm,*perm_tmp;
73     dim_t e,i,i0,i1,n;
74 jgs 82
75 jgs 102 X=TMPMEMALLOC(NN*numDim*faces->numElements,double);
76     center=TMPMEMALLOC(faces->numElements,Finley_Mesh_findMatchingFaces_center);
77 jgs 82 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     qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
104     /* find elements with matching center */
105     *numPairs=0;
106     /* OMP */
107 jgs 150 for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
108 jgs 82 if (Finley_Mesh_findMatchingFaces_compar((void*) &(center[e]),(void*) &(center[e+1]))==0) {
109     e_0=center[e].refId;
110     e_1=center[e+1].refId;
111     elem0[*numPairs]=e_0;
112     elem1[*numPairs]=e_1;
113     /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
114     perm=a1;
115     perm_tmp=a2;
116     for (i=0;i<NN;i++) perm[i]=i;
117 jgs 150 while (Finley_noError()) {
118 jgs 82 /* if node 0 and perm[0] are the same we are ready */
119     getDist(dist,e_0,0,e_1,perm[0]);
120     if (dist<=h*tolerance) break;
121     if (faces->ReferenceElement->Type->shiftNodes[0]>=0) {
122     /* rotate the nodes */
123     for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->shiftNodes[i]];
124     SWAP(perm,perm_tmp);
125     }
126     /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
127     if (perm[0]==0) {
128 jgs 150 sprintf(error_msg,"__FILE__:couldn't match first node of element %d to touching element %d",e_0,e_1);
129     Finley_setError(VALUE_ERROR,error_msg);
130 jgs 82 }
131     }
132     /* now we check if the second nodes match */
133 jgs 150 if (Finley_noError()) {
134 jgs 82 if (faces->ReferenceElement->Type->numNodesOnFace>1) {
135     getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
136     /* if the second node does not match we reverse the direction of the nodes */
137     if (dist>h*tolerance) {
138     /* rotate the nodes */
139     if (faces->ReferenceElement->Type->reverseNodes[0]<0) {
140 jgs 150 sprintf(error_msg,"__FILE__:couldn't match the second node of element %d to touching element %d",e_0,e_1);
141     Finley_setError(VALUE_ERROR,error_msg);
142 jgs 82 } else {
143     for (i=0;i<NN;i++) perm_tmp[i]=perm[faces->ReferenceElement->Type->reverseNodes[i]];
144     SWAP(perm,perm_tmp);
145     getDist(dist,e_0,1,e_1,perm[faces->ReferenceElement->Type->faceNode[1]]);
146     if (dist>h*tolerance) {
147 jgs 150 sprintf(error_msg,"__FILE__:couldn't match the second node of element %d to touching element %d",e_0,e_1);
148     Finley_setError(VALUE_ERROR,error_msg);
149 jgs 82 }
150     }
151     }
152     }
153     }
154     /* we check if the rest of the face nodes match: */
155 jgs 150 if (Finley_noError()) {
156 jgs 82 for (i=2;i<faces->ReferenceElement->Type->numNodesOnFace;i++) {
157     n=faces->ReferenceElement->Type->faceNode[i];
158     getDist(dist,e_0,n,e_1,perm[n]);
159     if (dist>h*tolerance) {
160 jgs 150 sprintf(error_msg,"__FILE__:couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
161     Finley_setError(VALUE_ERROR,error_msg);
162 jgs 82 break;
163     }
164     }
165     }
166     /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
167 jgs 150 if (Finley_noError()) {
168 jgs 82 for (i=0;i<NN;i++) matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
169     }
170     (*numPairs)++;
171     }
172     }
173     #ifdef Finley_TRACE
174     printf("number of pairs of matching faces %d\n",*numPairs);
175     #endif
176     }
177     /* clean up */
178     TMPMEMFREE(X);
179     TMPMEMFREE(center);
180    
181     #undef getDist
182     #undef SWAP
183     }
184    
185     /*
186     * $Log$
187 jgs 150 * Revision 1.6 2005/09/15 03:44:22 jgs
188     * Merge of development branch dev-02 back to main trunk on 2005-09-15
189     *
190     * Revision 1.5.2.1 2005/09/07 06:26:19 gross
191     * the solver from finley are put into the standalone package paso now
192     *
193 jgs 123 * Revision 1.5 2005/07/08 04:07:51 jgs
194     * Merge of development branch back to main trunk on 2005-07-08
195     *
196 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
197 jgs 97 * *** empty log message ***
198 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:51 gross
199     * some changes towards 64 integers in finley
200 jgs 82 *
201 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:14 gross
202     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
203 jgs 97 *
204 jgs 82 *
205 jgs 123 *
206 jgs 82 */
207    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26