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

Annotation of /trunk/dudley/src/Mesh_findMatchingFaces.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4348 - (hide annotations)
Tue Apr 2 07:15:43 2013 UTC (6 years, 1 month ago) by caltinay
File size: 7651 byte(s)
Reapply some of the spelling fixes that got lost....

1 jgs 150
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 ksteube 1312
16 jfenwick 3981 /************************************************************************************/
17 jgs 82
18 jfenwick 3086 /* Dudley: Mesh */
19 jgs 82
20     /* searches for faces in the mesh which are matching */
21    
22 jfenwick 3981 /************************************************************************************/
23 jgs 82
24     #include "Util.h"
25     #include "Mesh.h"
26 gross 2748
27 jfenwick 3210 #include "ShapeTable.h"
28    
29 jfenwick 3981 /************************************************************************************/
30 jgs 82
31 jfenwick 3224 static double Dudley_Mesh_lockingGridSize = 0;
32 jgs 82
33 jfenwick 3224 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 jgs 82 }
65    
66 jfenwick 3224 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 jgs 82 #define getDist(_dist_,_e0_,_i0_,_e1_,_i1_) \
71 jgs 123 {dim_t i; \
72 jgs 82 _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 jfenwick 3224 }
75 jgs 150 char error_msg[LenErrorMsg_MAX];
76 jfenwick 3224 double h = DBLE(HUGE_VAL), h_local, dist, *X = NULL;
77 jfenwick 3086 Dudley_Mesh_findMatchingFaces_center *center;
78 jfenwick 3224 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 jgs 82
82 jfenwick 3224 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 gross 2748 }
98 jfenwick 4332 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 jfenwick 3224 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 caltinay 4348 /* sort the elements by center coordinates (lexicographical) */
138 jfenwick 3224 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 jgs 82 }
251     /* clean up */
252 jfenwick 4332 delete[] X;
253     delete[] center;
254     delete[] a1;
255     delete[] a1;
256 jgs 82
257     #undef getDist
258     }

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

  ViewVC Help
Powered by ViewVC 1.1.26