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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 4 months ago) by caltinay
File size: 10589 byte(s)
sync and fix.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_findMatchingFaces.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_findMatchingFaces.cpp:2682-2741 /branches/pasowrap/dudley/src/Mesh_findMatchingFaces.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_findMatchingFaces.cpp:3871-3891 /branches/restext/finley/src/Mesh_findMatchingFaces.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Mesh_findMatchingFaces.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_findMatchingFaces.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Mesh_findMatchingFaces.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_findMatchingFaces.cpp:3517-3974 /release/3.0/finley/src/Mesh_findMatchingFaces.cpp:2591-2601 /release/4.0/dudley/src/Mesh_findMatchingFaces.cpp:5380-5406 /trunk/dudley/src/Mesh_findMatchingFaces.cpp:4257-4344,5898-5962 /trunk/ripley/test/python/dudley/src/Mesh_findMatchingFaces.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26