/[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 4496 - (show annotations)
Mon Jul 15 06:53:44 2013 UTC (6 years ago) by caltinay
File size: 8638 byte(s)
finley (WIP):
-moved all of finley into its namespace
-introduced some shared pointers
-Mesh is now a class
-other bits and pieces...

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

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/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 /release/3.0/finley/src/Mesh_findMatchingFaces.cpp:2591-2601 /trunk/finley/src/Mesh_findMatchingFaces.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26