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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3210 - (show annotations)
Mon Sep 27 04:13:19 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 8716 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26