00001 #include "SundanceUniformRefinementPair.hpp"
00002 using std::cout;
00003 using std::endl;
00004 using std::setw;
00005
00006 namespace Sundance
00007 {
00008 using namespace Teuchos;
00009
00010
00011 UniformRefinementPair::UniformRefinementPair() {}
00012
00013 UniformRefinementPair::UniformRefinementPair(const MeshType& meshType,
00014 const Mesh& coarse)
00015 : meshType_(meshType),
00016 coarse_(coarse),
00017 fine_()
00018 {
00019 if (coarse.spatialDim() == 2)
00020 {
00021 refineTriMesh();
00022 }
00023 TEUCHOS_TEST_FOR_EXCEPT(coarse.spatialDim() != 2);
00024 }
00025
00026 void UniformRefinementPair::refineTriMesh()
00027 {
00028 const MPIComm& comm = coarse_.comm();
00029 fine_ = meshType_.createEmptyMesh(2, comm);
00030
00031 TEUCHOS_TEST_FOR_EXCEPT(fine_.spatialDim() != 2);
00032 TEUCHOS_TEST_FOR_EXCEPT(comm.getNProc() > 1);
00033
00034 int numVerts = coarse_.numCells(0);
00035 int numEdges = coarse_.numCells(1);
00036 int numElems = coarse_.numCells(2);
00037
00038 int numNewVerts = numVerts + numEdges;
00039 int numNewEdges = 2*numEdges + 3*numElems;
00040 int numNewElems = 4*numElems;
00041
00042 Array<int> vertDone(numVerts, 0);
00043 Array<int> oldEdgeDone(numEdges, 0);
00044 Array<int> newEdgeDone(numNewEdges, 0);
00045
00046 oldToNewVertMap_ = Array<int>(numVerts, -1);
00047 newVertToOldLIDMap_ = Array<int>(numVerts+numEdges, -1);
00048 newVertIsOnEdge_ = Array<int>(numNewVerts, 0);
00049 oldEdgeToNewVertMap_ = Array<int>(numEdges, -1);
00050
00051 ArrayOfTuples<int> elemVerts(numElems, 6);
00052 oldToNewElemMap_ = ArrayOfTuples<int>(numElems, 4);
00053 newToOldElemMap_ = Array<int>(numNewElems, -1);
00054 oldEdgeChildren_ = Array<Array<int> >(numEdges);
00055 oldEdgeParallels_ = Array<Array<int> >(numEdges);
00056 newEdgeParents_ = Array<int>(numNewEdges, -1);
00057 newEdgeParallels_ = Array<int>(numNewEdges, -1);
00058
00059 interiorEdges_ = ArrayOfTuples<int>(numElems, 3);
00060
00061 int newVertLID = 0;
00062
00063 for (int c=0; c<numElems; c++)
00064 {
00065
00066 const int* verts = coarse_.elemZeroFacetView(c);
00067 for (int i=0; i<3; i++)
00068 {
00069 int v = verts[i];
00070 if (!vertDone[v])
00071 {
00072 int vertOwner = coarse_.ownerProcID(0, v);
00073 int vertLabel = coarse_.label(0, v);
00074 fine_.addVertex(newVertLID, coarse_.nodePosition(v), vertOwner, vertLabel);
00075
00076 oldToNewVertMap_[v] = newVertLID;
00077 newVertToOldLIDMap_[newVertLID] = v;
00078 elemVerts.value(c, i) = newVertLID;
00079 newVertLID++;
00080 vertDone[v]=true;
00081 }
00082 else
00083 {
00084 elemVerts.value(c, i) = oldToNewVertMap_[v];
00085 }
00086 }
00087
00088
00089 for (int i=0; i<3; i++)
00090 {
00091 int ori;
00092 int e = coarse_.facetLID(2, c, 1, i, ori);
00093 if (!oldEdgeDone[e])
00094 {
00095 int edgeOwner = coarse_.ownerProcID(1, e);
00096 Point centroid = coarse_.centroid(1, e);
00097 fine_.addVertex(newVertLID, centroid, edgeOwner, 0);
00098 elemVerts.value(c, 3+i) = newVertLID;
00099 oldEdgeToNewVertMap_[e] = newVertLID;
00100 newVertIsOnEdge_[newVertLID] = true;
00101 newVertToOldLIDMap_[newVertLID] = e;
00102 oldEdgeDone[e] = true;
00103 newVertLID++;
00104 }
00105 else
00106 {
00107 elemVerts.value(c, 3+i) = oldEdgeToNewVertMap_[e];
00108 }
00109 }
00110 }
00111
00112 Array<Array<int> > pairs
00113 = tuple<Array<int> >(
00114 tuple(0,1),
00115 tuple(1,2),
00116 tuple(2,0)
00117 );
00118
00119 Array<Array<int> > vertPtrs
00120 = tuple<Array<int> >(
00121 tuple(0,4,5),
00122 tuple(1,3,5),
00123 tuple(3,4,5),
00124 tuple(2,3,4)
00125 );
00126
00127 int newElemLID=0;
00128 for (int c=0; c<numElems; c++)
00129 {
00130
00131 Array<int> edgeLIDs(3);
00132 for (int s=0; s<3; s++)
00133 {
00134 int ori = 0;
00135 edgeLIDs[s] = coarse_.facetLID(2, c, 1, s, ori);
00136 }
00137
00138 int interiorEdgeCount = 0;
00139
00140
00141 int elemOwner = coarse_.ownerProcID(2, c);
00142 int elemLabel = coarse_.label(2, c);
00143 for (int p=0; p<vertPtrs.size(); p++)
00144 {
00145 Array<int> verts = tuple(elemVerts.value(c, vertPtrs[p][0]),
00146 elemVerts.value(c, vertPtrs[p][1]),
00147 elemVerts.value(c, vertPtrs[p][2]));
00148 fine_.addElement(newElemLID, verts, elemOwner, elemLabel);
00149
00150 oldToNewElemMap_.value(c, p) = newElemLID;
00151 newToOldElemMap_[newElemLID] = c;
00152
00153
00154 for (int side=0; side<3; side++)
00155 {
00156 int v1 = verts[pairs[side][0]];
00157 int v2 = verts[pairs[side][1]];
00158 int newEdge = lookupEdge(fine_, v1, v2);
00159 if (newEdgeDone[newEdge]) continue;
00160
00161
00162 if (newVertIsOnEdge_[v1] && newVertIsOnEdge_[v2])
00163 {
00164
00165
00166
00167 TEUCHOS_TEST_FOR_EXCEPT(interiorEdgeCount >= 3);
00168 interiorEdges_.value(c, interiorEdgeCount++) = newEdge;
00169
00170
00171 int parEdge = -1;
00172 for (int i=0; i<3; i++)
00173 {
00174 int ev = oldEdgeToNewVertMap_[edgeLIDs[i]];
00175 if (ev==v1 || ev==v2) continue;
00176 parEdge = edgeLIDs[i];
00177 break;
00178 }
00179 TEUCHOS_TEST_FOR_EXCEPT(parEdge==-1);
00180 oldEdgeParallels_[parEdge].append(newEdge);
00181 newEdgeParallels_[newEdge] = parEdge;
00182 }
00183 else
00184 {
00185 for (int n=0; n<2; n++)
00186 {
00187 int v = verts[pairs[side][n]];
00188 if (newVertIsOnEdge_[v])
00189 {
00190 int oldEdge = newVertToOldLIDMap_[v];
00191 oldEdgeChildren_[oldEdge].append(newEdge);
00192 newEdgeParents_[newEdge] = oldEdge;
00193 }
00194 }
00195 }
00196 newEdgeDone[newEdge] = true;
00197 }
00198 newElemLID++;
00199 }
00200 }
00201
00202 TEUCHOS_TEST_FOR_EXCEPT(newElemLID != 4*numElems);
00203
00204
00205 }
00206
00207 int UniformRefinementPair::lookupEdge(const Mesh& mesh,
00208 int v1, int v2) const
00209 {
00210 Array<int> cofacetLIDs;
00211 mesh.getCofacets(0, v1, 1, cofacetLIDs);
00212 for (int c=0; c<cofacetLIDs.size(); c++)
00213 {
00214 int ori;
00215 for (int f=0; f<2; f++)
00216 {
00217 int v = mesh.facetLID(1, cofacetLIDs[c], 0, f, ori);
00218 if (v == v2) return cofacetLIDs[c];
00219 }
00220 }
00221
00222 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00223 "edge (" << v1 << ", " << v2 << ") not found in mesh");
00224 return -1;
00225 }
00226
00227
00228 int UniformRefinementPair::check() const
00229 {
00230 int bad = 0;
00231 double tol = 1.0e-12;
00232 cout << "old vertex to new vertex map" << endl;
00233 cout << "----------------------------------------" << endl;
00234 for (int v=0; v<coarse_.numCells(0); v++)
00235 {
00236 int vNew = oldToNewVertMap()[v];
00237 double d = coarse_.centroid(0, v).distance(fine_.centroid(0, vNew));
00238 string state = "GOOD";
00239 if (d > tol) {state="BAD"; bad++;}
00240 if (fine_.label(0,vNew) != coarse_.label(0,v)) {state="BAD"; bad++;}
00241 cout << setw(4) << v
00242 << setw(16) << coarse_.centroid(0, v)
00243 << setw(4) << coarse_.label(0,v)
00244 << setw(4) << vNew
00245 << setw(16) << fine_.centroid(0, vNew)
00246 << setw(4) << fine_.label(0,vNew)
00247 << setw(6) << state << endl;
00248 }
00249
00250 cout << endl << endl;
00251
00252 cout << "old edge to new vertex map" << endl;
00253 cout << "----------------------------------------" << endl;
00254 for (int e=0; e<coarse_.numCells(1); e++)
00255 {
00256 int vNew = oldEdgeToNewVertMap()[e];
00257 double d = coarse_.centroid(1, e).distance(fine_.centroid(0, vNew));
00258 string state = "GOOD";
00259 if (d > tol) {state="BAD"; bad++;}
00260 cout << setw(4) << e << setw(16)
00261 << coarse_.centroid(1, e) << "\t" << setw(10)
00262 << vNew << setw(16)
00263 << fine_.centroid(0, vNew) << setw(6) << state << endl;
00264 }
00265
00266 cout << endl << endl;
00267 cout << "old vertex to new vertex map inversion" << endl;
00268 cout << "----------------------------------------" << endl;
00269 for (int v=0; v<coarse_.numCells(0); v++)
00270 {
00271 int vNew = oldToNewVertMap()[v];
00272 int vOld = newVertToOldLIDMap()[vNew];
00273
00274 string state = "GOOD";
00275 if (vOld != v) {state="BAD"; bad++;}
00276 cout << setw(4) << v << setw(4) << vNew << setw(4) << vOld
00277 << setw(6) << state << endl;
00278 }
00279
00280 cout << endl << endl;
00281 cout << "old edge to new vertex map inversion" << endl;
00282 cout << "----------------------------------------" << endl;
00283 for (int e=0; e<coarse_.numCells(1); e++)
00284 {
00285 int vNew = oldEdgeToNewVertMap()[e];
00286 int eOld = newVertToOldLIDMap()[vNew];
00287
00288 string state = "GOOD";
00289 if (eOld != e) {state="BAD"; bad++;}
00290 cout << setw(4) << e << setw(4) << vNew << setw(4) << eOld
00291 << setw(6) << state << endl;
00292 }
00293
00294
00295
00296
00297 cout << endl << endl;
00298 cout << "old to new elem map" << endl;
00299 cout << "----------------------------------------" << endl;
00300 for (int c=0; c<coarse_.numCells(2); c++)
00301 {
00302 cout << setw(6) << c;
00303 for (int e=0; e<4; e++)
00304 {
00305 int eNew = oldToNewElemMap().value(c, e);
00306 int eOld = newToOldElemMap()[eNew];
00307 cout << setw(4) << eNew ;
00308 string state = "OK";
00309 if (eOld != c) {state="BAD"; bad++;}
00310 cout << setw(4) << state;
00311 }
00312 cout << endl;
00313
00314 }
00315
00316 cout << endl << endl;
00317 cout << "old edge children and parallels" << endl;
00318 cout << "----------------------------------------" << endl;
00319 for (int e=0; e<coarse_.numCells(1); e++)
00320 {
00321 string state = "GOOD";
00322 const Array<int>& children = oldEdgeChildren()[e];
00323 const Array<int>& parallels = oldEdgeParallels()[e];
00324 Array<int> parents(children.size());
00325 Array<int> coarsePar(parallels.size());
00326 for (int i=0; i<children.size(); i++)
00327 {
00328 parents[i] = newEdgeParents()[children[i]];
00329 if (parents[i] != e) {state="BAD"; bad++;}
00330 }
00331 for (int i=0; i<parallels.size(); i++)
00332 {
00333 coarsePar[i] = newEdgeParallels()[parallels[i]];
00334 if (coarsePar[i] != e) {state="BAD"; bad++;}
00335 }
00336
00337 cout << setw(4) << e << setw(12) << children
00338 << setw(12) << parents
00339 << setw(12) << parallels
00340 << setw(12) << coarsePar
00341 << setw(6) << state << endl;
00342 }
00343
00344 cout << endl << endl;
00345 cout << "interior edges" << endl;
00346 cout << "----------------------------------------" << endl;
00347 for (int c=0; c<coarse_.numCells(2); c++)
00348 {
00349 string state = "GOOD";
00350 Array<int> edges(3);
00351 for (int i=0; i<3; i++)
00352 {
00353 edges[i] = interiorEdges_.value(c,i);
00354 int ori;
00355 int v1 = fine_.facetLID(1, edges[i], 0, 0, ori);
00356 if (!newVertIsOnEdge_[v1])
00357 {bad++; state="BAD";}
00358 int v2 = fine_.facetLID(1, edges[i], 0, 1, ori);
00359 if (!newVertIsOnEdge_[v2])
00360 {bad++; state="BAD";}
00361
00362
00363 Array<int> cofs;
00364 fine_.getCofacets(1, edges[i], 2, cofs);
00365 if (cofs.size() != 2) {bad++; state="BAD";}
00366 for (int cf=0; cf<cofs.size(); cf++)
00367 {
00368 if (newToOldElemMap_[cofs[cf]] != c)
00369 {bad++; state="BAD";}
00370 }
00371 }
00372
00373 cout << setw(4) << c << setw(20) << edges
00374 << setw(6) << state << endl;
00375 }
00376
00377 return bad;
00378 }
00379
00380 }