/[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 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (2 years, 11 months ago) by caltinay
File size: 9108 byte(s)
Much needed sync with trunk...

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 #include "Mesh.h"
26 #include "ShapeTable.h"
27 #include "Util.h"
28
29 namespace dudley {
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 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 < 3; i++) {
41 l = (e1->x[i] < e2->x[i] + Dudley_Mesh_lockingGridSize);
42 g = (e2->x[i] < e1->x[i] + Dudley_Mesh_lockingGridSize);
43 if (!(l && g)) {
44 if (l)
45 return -1;
46 if (g)
47 return 1;
48 }
49 }
50 if (e1->refId < e2->refId)
51 return -1;
52 else if (e1->refId > e2->refId)
53 return 1;
54 else
55 return 0;
56 }
57
58 inline double getDist(int e0, int i0, int e1, int i1, int numDim, int NN,
59 const double* X)
60 {
61 double dist=0.;
62 for (int i=0; i<numDim; i++) {
63 dist=std::max(dist, std::abs(X[INDEX3(i, i0, e0, numDim, NN)]
64 - X[INDEX3(i, i1, e1, numDim, NN)]));
65 }
66 return dist;
67 }
68
69 void Dudley_Mesh_findMatchingFaces(Dudley_NodeFile* nodes,
70 Dudley_ElementFile* faces, double safety_factor, double tolerance,
71 dim_t* numPairs, index_t* elem0, index_t* elem1,
72 index_t* matching_nodes_in_elem1)
73 {
74 double h = HUGE_VAL, dist, *X = NULL;
75 Dudley_Mesh_findMatchingFaces_center *center;
76 index_t e_0, e_1, *a1 = NULL, *a2 = NULL, *perm = NULL, *perm_tmp = NULL, *itmp_ptr = NULL;
77 const index_t *shiftNodes = NULL, *reverseNodes = NULL;
78 dim_t e, i, i0, i1, n, numNodesOnFace;
79
80 dim_t numDim = nodes->numDim;
81 int NN = faces->numNodes;
82 numNodesOnFace = numNodesOnFaceMap[faces->etype];
83 shiftNodes = shiftNodesMap[faces->etype];
84 reverseNodes = reverseNodesMap[faces->etype];
85
86 if (numNodesOnFace <= 0) {
87 std::stringstream ss;
88 ss << "Dudley_Mesh_findMatchingFaces: matching faces cannot be "
89 "applied to face elements of type " << getElementName(faces->etype);
90 const std::string errorMsg(ss.str());
91 throw DudleyException(errorMsg);
92 }
93 X = new double[NN * numDim * faces->numElements];
94 center = new Dudley_Mesh_findMatchingFaces_center[faces->numElements];
95 a1 = new int[NN];
96 a2 = new int[NN];
97
98 /* OMP */
99 for (e = 0; e < faces->numElements; e++) {
100 /* get the coordinates of the nodes */
101 Dudley_Util_Gather_double(NN, &(faces->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
102 &(X[INDEX3(0, 0, e, numDim, NN)]));
103 /* get the element center */
104 center[e].refId = e;
105 for (i = 0; i < 3; i++)
106 center[e].x[i] = 0;
107 for (i0 = 0; i0 < numNodesOnFace; i0++) {
108 for (i = 0; i < numDim; i++)
109 center[e].x[i] += X[INDEX3(i, i0, e, numDim, NN)];
110 }
111 for (i = 0; i < numDim; i++)
112 center[e].x[i] /= numNodesOnFace;
113 /* get the minimum distance between nodes in the element */
114 for (i0 = 0; i0 < numNodesOnFace; i0++) {
115 for (i1 = i0 + 1; i1 < numNodesOnFace; i1++) {
116 const double h_local=getDist(e, i0, e, i1, numDim, NN, X);
117 h=std::min(h, h_local);
118 }
119 }
120 }
121
122 Dudley_Mesh_lockingGridSize = h * std::max(safety_factor, 0.);
123 #ifdef Dudley_TRACE
124 printf("locking grid size is %e\n", Dudley_Mesh_lockingGridSize);
125 printf("absolute tolerance is %e.\n", h * tolerance);
126 #endif
127 /* sort the elements by center coordinates (lexicographical) */
128 qsort(center, faces->numElements, sizeof(Dudley_Mesh_findMatchingFaces_center),
129 Dudley_Mesh_findMatchingFaces_compar);
130 /* find elements with matching center */
131 *numPairs = 0;
132 /* OMP */
133 for (e = 0; e < faces->numElements - 1; e++) {
134 dist = 0;
135 for (i = 0; i < numDim; i++)
136 dist = std::max(dist, std::abs(center[e].x[i] - center[e + 1].x[i]));
137 if (dist < h * tolerance) {
138 e_0 = center[e].refId;
139 e_1 = center[e + 1].refId;
140 elem0[*numPairs] = e_0;
141 elem1[*numPairs] = e_1;
142 /* now the element e_1 is rotated such that the first node in element e_0 and e_1 have the same coordinates */
143 perm = a1;
144 perm_tmp = a2;
145 for (i = 0; i < NN; i++)
146 perm[i] = i;
147 while (1) {
148 /* if node 0 and perm[0] are the same we are ready */
149 dist=getDist(e_0, 0, e_1, perm[0], numDim, NN, X);
150 if (dist <= h * tolerance)
151 break;
152 if (shiftNodes[0] >= 0) {
153 /* rotate the nodes */
154 itmp_ptr = perm;
155 perm = perm_tmp;
156 perm_tmp = itmp_ptr;
157 #pragma ivdep
158 for (i = 0; i < NN; i++)
159 perm[i] = perm_tmp[shiftNodes[i]];
160 }
161 /* if the permutation is back at the identity, ie. perm[0]=0, the faces don't match: */
162 if (perm[0] == 0)
163 {
164 std::stringstream ss;
165 ss << "Mesh_findMatchingFaces: couldn't match first "
166 "node of element " << e_0 << " to touching element "
167 << e_1;
168 std::string errorMsg(ss.str());
169 throw DudleyException(errorMsg);
170 }
171 }
172 /* now we check if the second nodes match */
173 if (numNodesOnFace > 1) {
174 dist=getDist(e_0, 1, e_1, perm[1], numDim, NN, X);
175 /* if the second node does not match we reverse the direction of the nodes */
176 if (dist > h * tolerance) {
177 /* rotate the nodes */
178 if (reverseNodes[0] < 0) {
179 std::stringstream ss;
180 ss << "Mesh_findMatchingFaces: couldn't match "
181 "the second node of element " << e_0
182 << " to touching element " << e_1;
183 const std::string errorMsg(ss.str());
184 throw DudleyException(errorMsg);
185 } else {
186 itmp_ptr = perm;
187 perm = perm_tmp;
188 perm_tmp = itmp_ptr;
189 #pragma ivdep
190 for (i = 0; i < NN; i++)
191 perm[i] = perm_tmp[reverseNodes[i]];
192 dist=getDist(e_0, 1, e_1, perm[1], numDim, NN, X);
193 if (dist > h * tolerance) {
194 std::stringstream ss;
195 ss << "Mesh_findMatchingFaces: couldn't match "
196 "the second node of element " << e_0
197 << " to touching element " << e_1;
198 const std::string errorMsg(ss.str());
199 throw DudleyException(errorMsg);
200 }
201 }
202 }
203 }
204 /* we check if the rest of the face nodes match: */
205 for (i = 2; i < numNodesOnFace; i++) {
206 n = i;
207 dist=getDist(e_0, n, e_1, perm[n], numDim, NN, X);
208 if (dist > h * tolerance) {
209 std::stringstream ss;
210 ss << "Mesh_findMatchingFaces: couldn't match the "
211 << i << "-th node of element " << e_0
212 << " to touching element " << e_1;
213 const std::string errorMsg(ss.str());
214 throw DudleyException(errorMsg);
215 }
216 }
217 /* copy over the permuted nodes of e_1 into matching_nodes_in_elem1 */
218 for (i = 0; i < NN; i++)
219 matching_nodes_in_elem1[INDEX2(i, *numPairs, NN)] = faces->Nodes[INDEX2(perm[i], e_1, NN)];
220 (*numPairs)++;
221 }
222 }
223 #ifdef Dudley_TRACE
224 printf("number of pairs of matching faces %d\n", *numPairs);
225 #endif
226
227 delete[] X;
228 delete[] center;
229 delete[] a1;
230 delete[] a1;
231 }
232
233 } // namespace dudley
234

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-6007 /trunk/ripley/test/python/dudley/src/Mesh_findMatchingFaces.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26