00001 #include "SundanceRivaraDriver.hpp"
00002 #include "SundanceMesh.hpp"
00003 #include "SundanceRivaraMesh.hpp"
00004 #include "SundanceExprFieldWrapper.hpp"
00005
00006 using namespace Sundance;
00007 using namespace Sundance;
00008 using namespace Teuchos;
00009 using namespace Sundance::Rivara;
00010 using Sundance::ExprFieldWrapper;
00011 using std::endl;
00012
00013 void sort(const Array<int>& in, Array<int>& rtn)
00014 {
00015 rtn.resize(in.size());
00016 const int* pIn = &(in[0]);
00017 int* pRtn = &(rtn[0]);
00018 for (int i=0; i<in.size(); i++) pRtn[i] = pIn[i];
00019
00020 for (int j=1; j<in.size(); j++)
00021 {
00022 int key = pRtn[j];
00023 int i=j-1;
00024 while (i>=0 && pRtn[i]>key)
00025 {
00026 pRtn[i+1]=pRtn[i];
00027 i--;
00028 }
00029 pRtn[i+1]=key;
00030 }
00031
00032 }
00033
00034 static Time& refTimer()
00035 {
00036 static RCP<Time> rtn
00037 = TimeMonitor::getNewTimer("mesh refinement");
00038 return *rtn;
00039 }
00040
00041 static Time& m2rTimer()
00042 {
00043 static RCP<Time> rtn
00044 = TimeMonitor::getNewTimer("mesh to rivara copy");
00045 return *rtn;
00046 }
00047
00048 static Time& r2mTimer()
00049 {
00050 static RCP<Time> rtn
00051 = TimeMonitor::getNewTimer("rivara to mesh copy");
00052 return *rtn;
00053 }
00054
00055 static Time& volSetTimer()
00056 {
00057 static RCP<Time> rtn
00058 = TimeMonitor::getNewTimer("refinement stack building");
00059 return *rtn;
00060 }
00061
00062
00063 Mesh RefinementTransformation::apply(const Mesh& inputMesh) const
00064 {
00065 TimeMonitor timer(refTimer());
00066
00067 int dim = inputMesh.spatialDim();
00068 MPIComm comm = inputMesh.comm();
00069 int numElems = inputMesh.numCells(dim);
00070
00071 RCP<RivaraMesh> rivMesh = rcp(new RivaraMesh(dim, comm));
00072
00073 Array<int> lidMap;
00074
00075 meshToRivara(inputMesh, lidMap,rivMesh);
00076
00077 ExprFieldWrapper expr(errExpr_);
00078 TEUCHOS_TEST_FOR_EXCEPTION(expr.isPointData(), std::runtime_error,
00079 "Expected cell-based discrete function for area specification");
00080
00081 {
00082 TimeMonitor timer(volSetTimer());
00083 numRefined_ = 0;
00084 for (int c=0; c<numElems; c++)
00085 {
00086 double err = expr.getData(dim,c,0);
00087 int rivLID = lidMap[c];
00088 Element* e = rivMesh->element(rivLID).get();
00089 double vol = e->volume();
00090 double newVol = vol * std::pow(reqErr_/(err+1.0e-12), 0.5*dim);
00091
00092 if (newVol >= vol) continue;
00093 rivMesh->refinementSet().push(e);
00094 rivMesh->refinementAreas().push(newVol);
00095 numRefined_ ++;
00096 }
00097 }
00098 Out::os() << "refining n=" << numRefined_ << " cells" << std::endl;
00099 rivMesh->refine();
00100
00101
00102 Mesh outputMesh = rivaraToMesh(rivMesh, comm);
00103
00104 return outputMesh;
00105 }
00106
00107
00108 void RefinementTransformation::meshToRivara(
00109 const Mesh& mesh,
00110 Array<int>& lidMap,
00111 RCP<RivaraMesh>& rivMesh) const
00112 {
00113 TimeMonitor timer(m2rTimer());
00114 int dim = mesh.spatialDim();
00115 int numNodes = mesh.numCells(0);
00116 int numElems = mesh.numCells(dim);
00117
00118 for (int n=0; n<numNodes; n++)
00119 {
00120 int gid = n;
00121 int label = mesh.label(0,n);
00122 int ownerProcID = mesh.ownerProcID(0,n);
00123 Point x = mesh.nodePosition(n);
00124 rivMesh->addVertex(gid, x, ownerProcID, label);
00125 }
00126
00127 lidMap.resize(numElems);
00128 for (int e=0; e<numElems; e++)
00129 {
00130 int gid = e;
00131 int label = mesh.label(dim,e);
00132 int ownerProcID = mesh.ownerProcID(dim,e);
00133 Array<int> verts;
00134 Array<int> fo;
00135 mesh.getFacetArray(dim, e, 0, verts, fo);
00136 int elemLID = rivMesh->addElement(gid, verts, ownerProcID, label);
00137 lidMap[e] = elemLID;
00138
00139 if (dim==2)
00140 {
00141 Array<int> edgeLIDs;
00142 mesh.getFacetArray(dim, e, 1, edgeLIDs, fo);
00143 for (int f=0; f<3; f++)
00144 {
00145 int edgeLabel = mesh.label(1,edgeLIDs[f]);
00146 rivMesh->element(elemLID)->edge(f)->setLabel(edgeLabel);
00147 }
00148 }
00149 else if (dim==3)
00150 {
00151 Array<int> faceLIDs;
00152 mesh.getFacetArray(dim, e, 2, faceLIDs, fo);
00153 for (int f=0; f<4; f++)
00154 {
00155 int faceLabel = mesh.label(2,faceLIDs[f]);
00156 rivMesh->element(elemLID)->face(f)->setLabel(faceLabel);
00157 }
00158 }
00159 }
00160 }
00161
00162
00163 Mesh RefinementTransformation::rivaraToMesh(const RCP<RivaraMesh>& rivMesh,
00164 const MPIComm& comm) const
00165 {
00166 TimeMonitor timer(r2mTimer());
00167 int dim = rivMesh->spatialDim();
00168 int numNodes = rivMesh->numNodes();
00169
00170 Mesh mesh = meshType_.createEmptyMesh(dim, comm);
00171
00172 for (int n=0; n<numNodes; n++)
00173 {
00174 const RCP<Node>& node = rivMesh->node(n);
00175 const Point& x = node->pt();
00176 int gid = node->globalIndex();
00177 int ownerProcID = node->ownerProc();
00178 int label = node->label();
00179 mesh.addVertex(gid, x, ownerProcID, label);
00180 }
00181
00182
00183 ElementIterator iter(rivMesh.get());
00184
00185 int gid=0;
00186
00187 Array<int> verts(dim+1);
00188 Array<int> fo;
00189 Array<int> edgeLIDs;
00190 Array<int> faceLIDs;
00191
00192 while (iter.hasMoreElements())
00193 {
00194 Tabs tab;
00195 const Element* e = iter.getNextElement();
00196 int ownerProcID = e->ownerProc();
00197 int label = e->label();
00198 const Array<RCP<Node> >& nodes = e->nodes();
00199
00200 for (int i=0; i<nodes.size(); i++)
00201 {
00202 verts[i] = nodes[i]->globalIndex();
00203 }
00204 int lid = mesh.addElement(gid, verts, ownerProcID, label);
00205 gid++;
00206
00207 Array<int> sortedVerts(verts.size());
00208 sort(verts, sortedVerts);
00209
00210 Out::os() << tab << "elem LID=" << lid << " verts=" << sortedVerts << endl;
00211
00212 if (dim==2)
00213 {
00214 if (e->hasNoEdgeLabels()) continue;
00215 mesh.getFacetArray(dim, lid, 1, edgeLIDs, fo);
00216 for (int f=0; f<3; f++)
00217 {
00218 int edgeLabel = e->edge(f)->label();
00219 if (edgeLabel < 0) continue;
00220 mesh.setLabel(1, edgeLIDs[f], edgeLabel);
00221 }
00222 }
00223 else if (dim==3)
00224 {
00225 Tabs tab1;
00226 if (e->hasNoFaceLabels()) continue;
00227 mesh.getFacetArray(dim, lid, 2, faceLIDs, fo);
00228 for (int f=0; f<4; f++)
00229 {
00230
00231 int faceLabel = e->face(f)->label();
00232 if (faceLabel <= 0) continue;
00233 const FaceNodes& nodes = e->face(f)->nodes();
00234 Array<int> nodeLIDs = nodes.nodes().elements();
00235 Out::os() << tab1 << "face nodes = " << nodeLIDs << endl;
00236 int faceNum = -1;
00237 if (nodeLIDs[0]==sortedVerts[0] && nodeLIDs[1]==sortedVerts[1]
00238 && nodeLIDs[2]==sortedVerts[2])
00239 {
00240 faceNum = 3;
00241 }
00242 else if (nodeLIDs[0]==sortedVerts[0] && nodeLIDs[1]==sortedVerts[1]
00243 && nodeLIDs[2]==sortedVerts[3])
00244 {
00245 faceNum = 2;
00246 }
00247 else if (nodeLIDs[0]==sortedVerts[0] && nodeLIDs[1]==sortedVerts[2]
00248 && nodeLIDs[2]==sortedVerts[3])
00249 {
00250 faceNum = 1;
00251 }
00252 else if (nodeLIDs[0]==sortedVerts[1] && nodeLIDs[1]==sortedVerts[2]
00253 && nodeLIDs[2]==sortedVerts[3])
00254 {
00255 faceNum = 0;
00256 }
00257 else
00258 {
00259 TEUCHOS_TEST_FOR_EXCEPT(true);
00260 }
00261 Out::os() << tab1 << "faceLID=" << faceLIDs[faceNum] << " UFC face num="
00262 << faceNum
00263 << " label = " << faceLabel
00264 << " parent = " << lid << std::endl;
00265 mesh.setLabel(2, faceLIDs[faceNum], faceLabel);
00266 }
00267 }
00268 }
00269
00270 return mesh;
00271
00272 }
00273
00274
