00001
00002
00003
00004
00005 #include "SundanceVertexSort.hpp"
00006 #include "PlayaExceptions.hpp"
00007 #include "SundanceOut.hpp"
00008
00009 namespace Sundance
00010 {
00011 using Teuchos::Array;
00012 using Teuchos::tuple;
00013 using std::endl;
00014
00015
00016 int mapExSideToMissingVertex(int dim, int exFaceID)
00017 {
00018 static Array<int> map3D = tuple(2, 0, 1, 3);
00019 static Array<int> map2D = tuple(2, 0, 1);
00020 switch(dim)
00021 {
00022 case 3:
00023 return map3D[exFaceID];
00024 case 2:
00025 return map2D[exFaceID];
00026 default:
00027 TEUCHOS_TEST_FOR_EXCEPTION(dim != 3, std::runtime_error,
00028 "dimension " << dim << " not supported "
00029 "in mapExSideToMissingVertex()");
00030 }
00031 return -1;
00032 }
00033
00034 Array<int> exSideVertPos(int dim, int sideIndex)
00035 {
00036 static Array<Array<int> > faceVerts = tuple<Array<int> >
00037 (
00038 tuple(0,1,3),
00039 tuple(1,2,3),
00040 tuple(0,2,3),
00041 tuple(0,1,2)
00042 );
00043 static Array<Array<int> > edgeVerts = tuple<Array<int> >
00044 (
00045 tuple(0,1),
00046 tuple(1,2),
00047 tuple(0,2)
00048 );
00049
00050 if (dim==2) return edgeVerts[sideIndex];
00051 else if (dim==3) return faceVerts[sideIndex];
00052 else
00053 {
00054 TEUCHOS_TEST_FOR_EXCEPT(true);
00055 }
00056 return faceVerts[0];
00057 }
00058
00059 Array<int> ufcSideVertPos(int dim, int f)
00060 {
00061 static Array<Array<int> > faceVerts = tuple<Array<int> >
00062 (
00063 tuple(1,2,3),
00064 tuple(0,2,3),
00065 tuple(0,1,3),
00066 tuple(0,1,2)
00067 );
00068 static Array<Array<int> > edgeVerts = tuple<Array<int> >
00069 (
00070 tuple(1,2),
00071 tuple(0,2),
00072 tuple(0,1)
00073 );
00074
00075 if (dim==2) return edgeVerts[f];
00076 else if (dim==3) return faceVerts[f];
00077 else
00078 {
00079 TEUCHOS_TEST_FOR_EXCEPT(true);
00080 }
00081 return faceVerts[0];
00082
00083 }
00084
00085 Array<int> ufcSide(int f, const Array<int>& verts)
00086 {
00087 Array<int> rtn;
00088 for (int i=0; i<verts.size(); i++)
00089 {
00090 if (f==i) continue;
00091 rtn.append(verts[i]);
00092 }
00093 return rtn;
00094 }
00095
00096 Array<int> exSide(int f, const Array<int>& verts)
00097 {
00098 Array<int> rtn;
00099 Array<int> exPos = exSideVertPos(verts.size()-1, f);
00100 for (int i=0; i<exPos.size(); i++)
00101 {
00102 rtn.append(verts[exPos[i]]);
00103 }
00104 return rtn;
00105 }
00106
00107 void vertexSort(Array<int>& verts, int* key)
00108 {
00109 Array<std::pair<int,int> > A(verts.size());
00110 for (int i=0; i<verts.size(); i++)
00111 {
00112 A[i].first = verts[i];
00113 A[i].second = i;
00114 }
00115
00116 insertionSort(A);
00117
00118 *key = 0;
00119 int base = verts.size();
00120 int p = 1;
00121 for (int i=verts.size()-1; i>=0; i--)
00122 {
00123 verts[i] = A[i].first;
00124 *key += p * A[i].second;
00125 p *= base;
00126 }
00127 }
00128
00129 void getKeyedPerm(int key, Array<int>& digits)
00130 {
00131 int r = key;
00132 int B = digits.size();
00133 for (int b=B-1; b>=0; b--)
00134 {
00135 int B_toThe_b = iPow(B,b);
00136 digits[B-b-1] = r/B_toThe_b;
00137 r = r - digits[B-b-1]*B_toThe_b;
00138 }
00139 }
00140
00141
00142 int iPow(int base, int n)
00143 {
00144 int rtn = 1;
00145 for (int i=0; i<n; i++) rtn = rtn*base;
00146 return rtn;
00147 }
00148
00149
00150 int exVertPosToUFCVertPos(int meshDim, int permKey, int exVertPos)
00151 {
00152 Array<int> digits(meshDim+1);
00153 getKeyedPerm(permKey, digits);
00154 for (int i=0; i<digits.size(); i++)
00155 {
00156 if (digits[i] == exVertPos) return i;
00157 }
00158
00159 Out::os() << "vert=" << exVertPos << ", key=" << permKey
00160 << ", digits = " << digits << endl;
00161 return -1;
00162 }
00163
00164 int exFacetIndexToUFCFacetIndex(int meshDim, int permKey,
00165 int exFacetID)
00166 {
00167 int missingExVert = mapExSideToMissingVertex(meshDim, exFacetID);
00168 int missingUFCVert = exVertPosToUFCVertPos(meshDim, permKey, missingExVert);
00169 return missingUFCVert;
00170 }
00171
00172
00173 int ufcFacetIndexToExFacetIndex(int meshDim, int ufcFacetIndex)
00174 {
00175 static Array<int> map2D = tuple(1,2,0);
00176 static Array<int> map3D = tuple(1,2,0,3);
00177
00178 if (meshDim==2)
00179 {
00180 return map2D[ufcFacetIndex];
00181 }
00182 else if (meshDim==3)
00183 {
00184 return map3D[ufcFacetIndex];
00185 }
00186 else
00187 {
00188 TEUCHOS_TEST_FOR_EXCEPT( meshDim!=3 && meshDim!=2 );
00189 }
00190 TEUCHOS_TEST_FOR_EXCEPT(true);
00191 return -1;
00192 }
00193
00194 }
00195
00196
00197