/[escript]/branches/doubleplusgood/dudley/src/Mesh_findMatchingFaces.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Mesh_findMatchingFaces.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 7656 byte(s)
like sand though the hourglass
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 /* Dudley: 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 #include "ShapeTable.h"
28
29 /************************************************************************************/
30
31 static double Dudley_Mesh_lockingGridSize = 0;
32
33 int Dudley_Mesh_findMatchingFaces_compar(const void *arg1, const void *arg2)
34 {
35 Dudley_Mesh_findMatchingFaces_center *e1, *e2;
36 bool_t l, g;
37 dim_t i;
38 e1 = (Dudley_Mesh_findMatchingFaces_center *) arg1;
39 e2 = (Dudley_Mesh_findMatchingFaces_center *) arg2;
40 for (i = 0; i < MAX_numDim; i++)
41 {
42 l = (e1->x[i] < e2->x[i] + Dudley_Mesh_lockingGridSize) ? TRUE : FALSE;
43 g = (e2->x[i] < e1->x[i] + Dudley_Mesh_lockingGridSize) ? TRUE : FALSE;
44 if (!(l && g))
45 {
46 if (l)
47 return -1;
48 if (g)
49 return 1;
50 }
51 }
52 if (e1->refId < e2->refId)
53 {
54 return -1;
55 }
56 else if (e1->refId > e2->refId)
57 {
58 return 1;
59 }
60 else
61 {
62 return 0;
63 }
64 }
65
66 void Dudley_Mesh_findMatchingFaces(Dudley_NodeFile * nodes, Dudley_ElementFile * faces, double safety_factor,
67 double tolerance, dim_t * numPairs, index_t * elem0, index_t * elem1,
68 index_t * matching_nodes_in_elem1)
69 {
70 #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
71 {dim_t i; \
72 _dist_=0; \
73 for (i=0;i<numDim;i++) _dist_=MAX(_dist_,ABS(X[INDEX3(i,_i0_,_e0_,numDim,NN)]-X[INDEX3(i,_i1_,_e1_,numDim,NN)])); \
74 }
75 char error_msg[LenErrorMsg_MAX];
76 double h = DBLE(HUGE_VAL), h_local, dist, *X = NULL;
77 Dudley_Mesh_findMatchingFaces_center *center;
78 index_t e_0, e_1, *a1 = NULL, *a2 = NULL, *perm = NULL, *perm_tmp = NULL, *itmp_ptr = NULL;
79 const index_t *shiftNodes = NULL, *reverseNodes = NULL;
80 dim_t e, i, i0, i1, n, NN, numNodesOnFace;
81
82 dim_t numDim = nodes->numDim;
83
84 NN = faces->numNodes;
85
86 numNodesOnFace = numNodesOnFaceMap[faces->etype];
87 shiftNodes = shiftNodesMap[faces->etype];
88 reverseNodes = reverseNodesMap[faces->etype];
89
90 if (numNodesOnFace <= 0)
91 {
92 sprintf(error_msg,
93 "Dudley_Mesh_findMatchingFaces: matching faces cannot be applied to face elements of type %s",
94 getElementName(faces->etype));
95 Dudley_setError(TYPE_ERROR, error_msg);
96 return;
97 }
98 X = new double[NN * numDim * faces->numElements];
99 center = new Dudley_Mesh_findMatchingFaces_center[faces->numElements];
100 a1 = new int[NN];
101 a2 = new int[NN];
102 if (!(Dudley_checkPtr(X) || Dudley_checkPtr(center) || Dudley_checkPtr(a1) || Dudley_checkPtr(a2)))
103 {
104 /* OMP */
105 for (e = 0; e < faces->numElements; e++)
106 {
107 /* get the coordinates of the nodes */
108 Dudley_Util_Gather_double(NN, &(faces->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
109 &(X[INDEX3(0, 0, e, numDim, NN)]));
110 /* get the element center */
111 center[e].refId = e;
112 for (i = 0; i < MAX_numDim; i++)
113 center[e].x[i] = 0;
114 for (i0 = 0; i0 < numNodesOnFace; i0++)
115 {
116 for (i = 0; i < numDim; i++)
117 center[e].x[i] += X[INDEX3(i, i0, e, numDim, NN)];
118 }
119 for (i = 0; i < numDim; i++)
120 center[e].x[i] /= numNodesOnFace;
121 /* get the minimum distance between nodes in the element */
122 for (i0 = 0; i0 < numNodesOnFace; i0++)
123 {
124 for (i1 = i0 + 1; i1 < numNodesOnFace; i1++)
125 {
126 getDist(h_local, e, i0, e, i1);
127 h = MIN(h, h_local);
128 }
129 }
130 }
131 /* set the */
132 Dudley_Mesh_lockingGridSize = h * MAX(safety_factor, 0);
133 #ifdef Dudley_TRACE
134 printf("locking grid size is %e\n", Dudley_Mesh_lockingGridSize);
135 printf("absolute tolerance is %e.\n", h * tolerance);
136 #endif
137 /* sort the elements by center center coordinates (lexigraphical) */
138 qsort(center, faces->numElements, sizeof(Dudley_Mesh_findMatchingFaces_center),
139 Dudley_Mesh_findMatchingFaces_compar);
140 /* find elements with matching center */
141 *numPairs = 0;
142 /* OMP */
143 for (e = 0; e < faces->numElements - 1 && Dudley_noError(); e++)
144 {
145 dist = 0;
146 for (i = 0; i < numDim; i++)
147 dist = MAX(dist, ABS(center[e].x[i] - center[e + 1].x[i]));
148 if (dist < h * tolerance)
149 {
150 e_0 = center[e].refId;
151 e_1 = center[e + 1].refId;
152 elem0[*numPairs] = e_0;
153 elem1[*numPairs] = e_1;
154 /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
155 perm = a1;
156 perm_tmp = a2;
157 for (i = 0; i < NN; i++)
158 perm[i] = i;
159 while (Dudley_noError())
160 {
161 /* if node 0 and perm[0] are the same we are ready */
162 getDist(dist, e_0, 0, e_1, perm[0]);
163 if (dist <= h * tolerance)
164 break;
165 if (shiftNodes[0] >= 0)
166 {
167 /* rotate the nodes */
168 itmp_ptr = perm;
169 perm = perm_tmp;
170 perm_tmp = itmp_ptr;
171 #pragma ivdep
172 for (i = 0; i < NN; i++)
173 perm[i] = perm_tmp[shiftNodes[i]];
174 }
175 /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
176 if (perm[0] == 0)
177 {
178 sprintf(error_msg,
179 "Mesh_findMatchingFaces:couldn't match first node of element %d to touching element %d",
180 e_0, e_1);
181 Dudley_setError(VALUE_ERROR, error_msg);
182 }
183 }
184 /* now we check if the second nodes match */
185 if (Dudley_noError())
186 {
187 if (numNodesOnFace > 1)
188 {
189 getDist(dist, e_0, 1, e_1, perm[1]);
190 /* if the second node does not match we reverse the direction of the nodes */
191 if (dist > h * tolerance)
192 {
193 /* rotate the nodes */
194 if (reverseNodes[0] < 0)
195 {
196 sprintf(error_msg,
197 "Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",
198 e_0, e_1);
199 Dudley_setError(VALUE_ERROR, error_msg);
200 }
201 else
202 {
203 itmp_ptr = perm;
204 perm = perm_tmp;
205 perm_tmp = itmp_ptr;
206 #pragma ivdep
207 for (i = 0; i < NN; i++)
208 perm[i] = perm_tmp[reverseNodes[i]];
209 getDist(dist, e_0, 1, e_1, perm[1]);
210 if (dist > h * tolerance)
211 {
212 sprintf(error_msg,
213 "Mesh_findMatchingFaces:couldn't match the second node of element %d to touching element %d",
214 e_0, e_1);
215 Dudley_setError(VALUE_ERROR, error_msg);
216 }
217 }
218 }
219 }
220 }
221 /* we check if the rest of the face nodes match: */
222 if (Dudley_noError())
223 {
224 for (i = 2; i < numNodesOnFace; i++)
225 {
226 n = i;
227 getDist(dist, e_0, n, e_1, perm[n]);
228 if (dist > h * tolerance)
229 {
230 sprintf(error_msg,
231 "Mesh_findMatchingFaces:couldn't match the %d-th node of element %d to touching element %d",
232 i, e_0, e_1);
233 Dudley_setError(VALUE_ERROR, error_msg);
234 break;
235 }
236 }
237 }
238 /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
239 if (Dudley_noError())
240 {
241 for (i = 0; i < NN; i++)
242 matching_nodes_in_elem1[INDEX2(i, *numPairs, NN)] = faces->Nodes[INDEX2(perm[i], e_1, NN)];
243 }
244 (*numPairs)++;
245 }
246 }
247 #ifdef Dudley_TRACE
248 printf("number of pairs of matching faces %d\n", *numPairs);
249 #endif
250 }
251 /* clean up */
252 delete[] X;
253 delete[] center;
254 delete[] a1;
255 delete[] a1;
256
257 #undef getDist
258 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26