/[escript]/trunk/finley/src/Mesh_findMatchingFaces.cpp
ViewVC logotype

Contents of /trunk/finley/src/Mesh_findMatchingFaces.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 8497 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Finley: Mesh */
20
21 /* searches for faces in the mesh which are matching */
22
23 /************************************************************************************/
24
25 #include "Util.h"
26 #include "Mesh.h"
27
28 /************************************************************************************/
29
30 using namespace finley;
31
32 static double Finley_Mesh_lockingGridSize=0;
33
34 int Finley_Mesh_findMatchingFaces_compar(const void *arg1 , const void *arg2 ) {
35 Finley_Mesh_findMatchingFaces_center *e1,*e2;
36 bool_t l,g;
37 dim_t i;
38 e1=(Finley_Mesh_findMatchingFaces_center*) arg1;
39 e2=(Finley_Mesh_findMatchingFaces_center*) arg2;
40 for (i=0;i<MAX_numDim;i++) {
41 l= (e1->x[i]<e2->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
42 g= (e2->x[i]<e1->x[i]+Finley_Mesh_lockingGridSize) ? TRUE : FALSE;
43 if (! (l && g)) {
44 if (l) return -1;
45 if (g) return 1;
46 }
47 }
48 if ( e1->refId < e2->refId ) {
49 return -1;
50 } else if ( e1->refId > e2->refId ) {
51 return 1;
52 } else {
53 return 0;
54 }
55 }
56
57 void Finley_Mesh_findMatchingFaces(NodeFile *nodes, ElementFile *faces,
58 double safety_factor,double tolerance,
59 int* numPairs, int* elem0, int* elem1,
60 int* matching_nodes_in_elem1)
61 {
62 #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
63 {\
64 _dist_=0;\
65 for (int i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)]));\
66 }
67 ReferenceElement* refElement=NULL;
68 char error_msg[LenErrorMsg_MAX];
69 double h=DBLE(HUGE_VAL),h_local,dist,*X=NULL;
70 Finley_Mesh_findMatchingFaces_center *center;
71 index_t e_0,e_1,*a1=NULL,*a2=NULL,*perm=NULL,*perm_tmp=NULL,*itmp_ptr=NULL, *faceNodes=NULL, *shiftNodes=NULL, *reverseNodes=NULL ;
72 dim_t e,i,i0,i1,n,NN,numNodesOnFace;
73
74 dim_t numDim=nodes->numDim;
75
76 refElement= ReferenceElementSet_borrowReferenceElement(faces->referenceElementSet, FALSE);
77 NN=faces->numNodes;
78
79 numNodesOnFace=refElement->Type->numNodesOnFace;
80 faceNodes=refElement->Type->faceNodes;
81 shiftNodes=refElement->Type->shiftNodes;
82 reverseNodes=refElement->Type->reverseNodes;
83
84 if (numNodesOnFace<=0) {
85 sprintf(error_msg,"Finley_Mesh_findMatchingFaces: matching faces cannot be applied to face elements of type %s",refElement->Type->Name);
86 Finley_setError(TYPE_ERROR,error_msg);
87 return;
88 }
89 X=new double[NN*numDim*faces->numElements];
90 center=new Finley_Mesh_findMatchingFaces_center[faces->numElements];
91 a1=new int[NN];
92 a2=new int[NN];
93 /* TODO: OMP */
94 for (e=0;e<faces->numElements;e++) {
95 /* get the coordinates of the nodes */
96 util::gather(NN, &(faces->Nodes[INDEX2(0,e,NN)]), numDim,
97 nodes->Coordinates, &(X[INDEX3(0,0,e,numDim,NN)]));
98 /* get the element center */
99 center[e].refId=e;
100 for (i=0;i<MAX_numDim;i++) center[e].x[i]=0;
101 for (i0=0;i0<numNodesOnFace;i0++) {
102 for (i=0;i<numDim;i++) center[e].x[i]+=X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
103 }
104 for (i=0;i<numDim;i++) center[e].x[i]/=numNodesOnFace;
105 /* get the minimum distance between nodes in the element */
106 for (i0=0;i0<numNodesOnFace;i0++) {
107 for (i1=i0+1;i1<numNodesOnFace;i1++) {
108 getDist(h_local,e,faceNodes[i0],e,faceNodes[i1]);
109 h=MIN(h,h_local);
110 }
111 }
112 }
113 Finley_Mesh_lockingGridSize=h*MAX(safety_factor,0);
114 #ifdef Finley_TRACE
115 printf("locking grid size is %e\n",Finley_Mesh_lockingGridSize);
116 printf("absolute tolerance is %e.\n",h * tolerance);
117 #endif
118 /* sort the elements by center coordinates (lexicographical)*/
119 qsort(center,faces->numElements,sizeof(Finley_Mesh_findMatchingFaces_center),Finley_Mesh_findMatchingFaces_compar);
120 /* find elements with matching center */
121 *numPairs=0;
122 /* TODO: OMP */
123 for (e=0;e<faces->numElements-1 && Finley_noError();e++) {
124 dist=0;
125 for (i=0;i<numDim;i++) dist=MAX(dist,ABS(center[e].x[i]-center[e+1].x[i]));
126 if (dist < h * tolerance) {
127 e_0=center[e].refId;
128 e_1=center[e+1].refId;
129 elem0[*numPairs]=e_0;
130 elem1[*numPairs]=e_1;
131 /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
132 perm=a1;
133 perm_tmp=a2;
134 for (i=0;i<NN;i++) perm[i]=i;
135 while (Finley_noError()) {
136 /* if node 0 and perm[0] are the same we are ready */
137 getDist(dist,e_0,0,e_1,perm[0]);
138 if (dist<=h*tolerance) break;
139 if (shiftNodes[0]>=0) {
140 /* rotate the nodes */
141 itmp_ptr=perm;
142 perm=perm_tmp;
143 perm_tmp=itmp_ptr;
144 #pragma ivdep
145 for (i=0;i<NN;i++) perm[i]=perm_tmp[shiftNodes[i]];
146 }
147 /* if the permutation is back at the identity, i.e. perm[0]=0, the faces don't match: */
148 if (perm[0]==0) {
149 sprintf(error_msg,"Mesh_findMatchingFaces: couldn't match first node of element %d to touching element %d",e_0,e_1);
150 Finley_setError(VALUE_ERROR,error_msg);
151 }
152 }
153 /* now we check if the second nodes match */
154 if (Finley_noError()) {
155 if (numNodesOnFace>1) {
156 getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
157 /* if the second node does not match we reverse the direction of the nodes */
158 if (dist>h*tolerance) {
159 /* rotate the nodes */
160 if (reverseNodes[0]<0) {
161 sprintf(error_msg,"Mesh_findMatchingFaces: couldn't match the second node of element %d to touching element %d",e_0,e_1);
162 Finley_setError(VALUE_ERROR,error_msg);
163 } else {
164 itmp_ptr=perm;
165 perm=perm_tmp;
166 perm_tmp=itmp_ptr;
167 #pragma ivdep
168 for (i=0;i<NN;i++) perm[i]=perm_tmp[reverseNodes[i]];
169 getDist(dist,e_0,1,e_1,perm[faceNodes[1]]);
170 if (dist>h*tolerance) {
171 sprintf(error_msg,"Mesh_findMatchingFaces: couldn't match the second node of element %d to touching element %d",e_0,e_1);
172 Finley_setError(VALUE_ERROR,error_msg);
173 }
174 }
175 }
176 }
177 }
178 /* we check if the rest of the face nodes match: */
179 if (Finley_noError()) {
180 for (i=2;i<numNodesOnFace;i++) {
181 n=faceNodes[i];
182 getDist(dist,e_0,n,e_1,perm[n]);
183 if (dist>h*tolerance) {
184 sprintf(error_msg,"Mesh_findMatchingFaces: couldn't match the %d-th node of element %d to touching element %d",i,e_0,e_1);
185 Finley_setError(VALUE_ERROR,error_msg);
186 break;
187 }
188 }
189 }
190 // copy over the permuted nodes of e_1 into matching_nodes_in_elem1
191 if (Finley_noError()) {
192 for (i=0;i<NN;i++)
193 matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)]=faces->Nodes[INDEX2(perm[i],e_1,NN)];
194 }
195 (*numPairs)++;
196 }
197 }
198 #ifdef Finley_TRACE
199 printf("number of pairs of matching faces %d\n",*numPairs);
200 #endif
201
202 /* clean up */
203 delete[] X;
204 delete[] center;
205 delete[] a1;
206 delete[] a2;
207 #undef getDist
208 }
209

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/Mesh_findMatchingFaces.cpp:2682-2741 /branches/pasowrap/finley/src/Mesh_findMatchingFaces.cpp:3661-3674 /branches/py3_attempt2/finley/src/Mesh_findMatchingFaces.cpp:3871-3891 /branches/restext/finley/src/Mesh_findMatchingFaces.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Mesh_findMatchingFaces.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_findMatchingFaces.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Mesh_findMatchingFaces.cpp:3471-3974 /release/3.0/finley/src/Mesh_findMatchingFaces.cpp:2591-2601 /trunk/finley/src/Mesh_findMatchingFaces.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26