/[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 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 9040 byte(s)
Updated the copyright header.


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

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