/[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 6620 - (show annotations)
Mon Aug 14 08:31:59 2017 UTC (5 weeks, 2 days ago) by gross
File size: 8899 byte(s)
gather function fixed.
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2017 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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
20 Finley: Domain
21
22 searches for faces in the mesh which are matching.
23
24 *****************************************************************************/
25
26 #include "FinleyDomain.h"
27 #include "Util.h"
28
29 #include <escript/index.h>
30
31 namespace finley {
32
33 static double lockingGridSize = 0.;
34
35 // this structure is used for matching surface elements
36 struct FaceCenter
37 {
38 int refId;
39 std::vector<double> x;
40 };
41 #define Finley_TRACE
42
43 /// comparison function for findMatchingFaces
44 bool FaceCenterCompare(const FaceCenter& e1, const FaceCenter& e2)
45 {
46 for (int i = 0; i < e1.x.size(); i++) {
47 bool l = (e1.x[i] < e2.x[i]+lockingGridSize);
48 bool g = (e2.x[i] < e1.x[i]+lockingGridSize);
49 if (! (l && g)) {
50 if (l) return true;
51 if (g) return false;
52 }
53 }
54 return e1.refId < e2.refId; // strict order to e1==e2 is false
55 }
56
57 inline double getDist(int e0, int i0, int e1, int i1, int numDim, int NN,
58 const double* X)
59 {
60 double dist = 0.;
61 for (int i = 0; i < numDim; i++) {
62 dist = std::max(dist, std::abs(X[INDEX3(i, i0, e0, numDim, NN)]
63 - X[INDEX3(i, i1, e1, numDim, NN)]));
64 }
65 return dist;
66 }
67
68 void FinleyDomain::findMatchingFaces(double safety_factor, double tolerance,
69 int* numPairs, int* elem0, int* elem1,
70 int* matching_nodes_in_elem1) const
71 {
72 const_ReferenceElement_ptr refElement(m_faceElements->referenceElementSet->
73 borrowReferenceElement(false));
74 const int numDim = m_nodes->numDim;
75 const int NN = m_faceElements->numNodes;
76 const int numNodesOnFace = refElement->Type->numNodesOnFace;
77 const int* faceNodes = refElement->Type->faceNodes;
78 const int* shiftNodes = refElement->Type->shiftNodes;
79 const int* reverseNodes = refElement->Type->reverseNodes;
80
81 if (numNodesOnFace <= 0) {
82 std::stringstream ss;
83 ss << "Mesh::findMatchingFaces: matching faces cannot be applied to "
84 "face elements of type " << refElement->Type->Name;
85 throw escript::ValueError(ss.str());
86 }
87 double* X = new double[NN * numDim * m_faceElements->numElements];
88 std::vector<FaceCenter> center(m_faceElements->numElements);
89 int* a1 = new int[NN];
90 int* a2 = new int[NN];
91 double h = std::numeric_limits<double>::max();
92
93 // TODO: OMP
94 for (index_t e = 0; e < m_faceElements->numElements; e++) {
95 // get the coordinates of the nodes
96 util::gather(NN, &(m_faceElements->Nodes[INDEX2(0,e,NN)]), numDim,
97 m_nodes->Coordinates, &X[INDEX3(0,0,e,numDim,NN)]);
98 // get the element center
99 center[e].refId = e;
100 center[e].x.assign(numDim, 0);
101 for (int i0 = 0; i0 < numNodesOnFace; i0++) {
102 for (int i = 0; i < numDim; i++)
103 center[e].x[i] += X[INDEX3(i,faceNodes[i0],e,numDim,NN)];
104 }
105 for (int i = 0; i < numDim; i++)
106 center[e].x[i] /= numNodesOnFace;
107 // get the minimum distance between nodes in the element
108 for (int i0 = 0; i0 < numNodesOnFace; i0++) {
109 for (int i1 = i0+1; i1 < numNodesOnFace; i1++) {
110 double h_local = getDist(e, faceNodes[i0], e, faceNodes[i1], numDim, NN, X);
111 h = std::min(h, h_local);
112 }
113 }
114 }
115 lockingGridSize = h*std::max(safety_factor, 0.);
116 #ifdef Finley_TRACE
117 printf("locking grid size is %e\n", lockingGridSize);
118 printf("absolute tolerance is %e.\n", h * tolerance);
119 printf("number of face elements is %d\n", m_faceElements->numElements);
120 #endif
121 // sort the elements by center coordinates (lexicographical)
122 std::sort(center.begin(), center.end(), FaceCenterCompare);
123 // find elements with matching center
124 *numPairs = 0;
125
126 // TODO: OMP
127 for (index_t e = 0; e < m_faceElements->numElements-1; e++) {
128 double dist = 0.;
129 for (int i = 0; i < numDim; i++)
130 dist = std::max(dist, std::abs(center[e].x[i]-center[e+1].x[i]));
131 if (dist < h * tolerance) {
132 const int e_0 = center[e].refId;
133 const int e_1 = center[e+1].refId;
134 elem0[*numPairs] = e_0;
135 elem1[*numPairs] = e_1;
136 // now the element e_1 is rotated such that the first node in
137 // element e_0 and e_1 have the same coordinates
138 int* perm = a1;
139 int* perm_tmp = a2;
140 for (int i = 0; i < NN; i++)
141 perm[i] = i;
142 while (1) {
143 // if node 0 and perm[0] are the same we are ready
144 dist = getDist(e_0, 0, e_1, perm[0], numDim, NN, X);
145 if (dist <= h*tolerance)
146 break;
147 if (shiftNodes[0] >= 0) {
148 // rotate the nodes
149 int* itmp_ptr = perm;
150 perm = perm_tmp;
151 perm_tmp = itmp_ptr;
152 #pragma ivdep
153 for (int i = 0; i < NN; i++)
154 perm[i] = perm_tmp[shiftNodes[i]];
155 }
156 // if the permutation is back at the identity, i.e. perm[0]=0,
157 // the faces don't match:
158 if (perm[0] == 0) {
159 std::stringstream ss;
160 ss << "Mesh::findMatchingFaces: couldn't match first node "
161 "of element " << e_0 << " to touching element " << e_1;
162 throw escript::ValueError(ss.str());
163 }
164 }
165 // now we check if the second nodes match
166 if (numNodesOnFace > 1) {
167 dist = getDist(e_0, 1, e_1, perm[faceNodes[1]], numDim, NN, X);
168 // if the second node does not match we reverse the
169 // direction of the nodes
170 if (dist > h*tolerance) {
171 // rotate the nodes
172 if (reverseNodes[0] < 0) {
173 std::stringstream ss;
174 ss << "Mesh::findMatchingFaces: couldn't match the"
175 " second node of element " << e_0
176 << " to touching element " << e_1;
177 throw escript::ValueError(ss.str());
178 } else {
179 int* itmp_ptr = perm;
180 perm = perm_tmp;
181 perm_tmp = itmp_ptr;
182 #pragma ivdep
183 for (int i = 0; i < NN; i++)
184 perm[i] = perm_tmp[reverseNodes[i]];
185 dist = getDist(e_0, 1, e_1, perm[faceNodes[1]], numDim, NN, X);
186 if (dist > h*tolerance) {
187 std::stringstream ss;
188 ss << "Mesh::findMatchingFaces: couldn't match the"
189 " second node of element " << e_0
190 << " to touching element " << e_1;
191 throw escript::ValueError(ss.str());
192 }
193 }
194 }
195 }
196 // we check if the rest of the face nodes match
197 for (int i = 2; i < numNodesOnFace; i++) {
198 const int n = faceNodes[i];
199 dist = getDist(e_0, n, e_1, perm[n], numDim, NN, X);
200 if (dist > h*tolerance) {
201 std::stringstream ss;
202 ss << "Mesh::findMatchingFaces: couldn't match the "
203 << i << "-th node of element " << e_0
204 << " to touching element " << e_1;
205 throw escript::ValueError(ss.str());
206 }
207 }
208 // copy over the permuted nodes of e_1 into matching_nodes_in_elem1
209 for (int i = 0; i < NN; i++)
210 matching_nodes_in_elem1[INDEX2(i,*numPairs,NN)] =
211 m_faceElements->Nodes[INDEX2(perm[i],e_1,NN)];
212 (*numPairs)++;
213 }
214 }
215 #ifdef Finley_TRACE
216 printf("number of pairs of matching faces %d\n", *numPairs);
217 #endif
218
219 delete[] X;
220 delete[] a1;
221 delete[] a2;
222 }
223
224 } // namespace finley
225

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Mesh_findMatchingFaces.cpp:5567-5588 /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 /branches/trilinos_from_5897/finley/src/Mesh_findMatchingFaces.cpp:5898-6118 /release/3.0/finley/src/Mesh_findMatchingFaces.cpp:2591-2601 /release/4.0/finley/src/Mesh_findMatchingFaces.cpp:5380-5406 /trunk/finley/src/Mesh_findMatchingFaces.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26