00001 #include "SundanceMeshIOUtils.hpp"
00002 #include "PlayaOut.hpp"
00003 #include "PlayaTabs.hpp"
00004
00005
00006 namespace Sundance
00007 {
00008 using namespace Sundance;
00009 using namespace Sundance;
00010 using namespace Sundance;
00011 using namespace Teuchos;
00012
00013
00014 Expr readNodalFields(const MeshSource& mesher, const Mesh& mesh,
00015 const VectorType<double>& vecType, int verb)
00016 {
00017 Tabs tab0(0);
00018 PLAYA_ROOT_MSG1(verb, tab0 << "begin readNodalFields()");
00019 Tabs tab1;
00020
00021
00022 RCP<Array<Array<double> > > elemAttr;
00023 RCP<Array<Array<double> > > vertAttr;
00024
00025 PLAYA_ROOT_MSG2(verb, tab1 << "reading attributes");
00026 mesher.getAttributes(vertAttr, elemAttr);
00027 int nAttrs = vertAttr->size();
00028 PLAYA_ROOT_MSG2(verb, tab1 << "found " << nAttrs << " attributes");
00029 const Array<Array<double> >& funcVals = *vertAttr;
00030
00031
00032 Array<BasisFamily> bas(nAttrs);
00033 for (int i=0; i<bas.size(); i++) bas[i] = new Lagrange(1);
00034 PLAYA_ROOT_MSG2(verb, tab1 << "forming discrete space");
00035 DiscreteSpace discSpace(mesh, bas, vecType);
00036 PLAYA_ROOT_MSG2(verb, tab1 << "forming discrete function");
00037 Expr u0 = new DiscreteFunction(discSpace, 0.0, "u0");
00038
00039
00040
00041 Vector<double> vec = DiscreteFunction::discFunc(u0)->getVector();
00042 const RCP<DOFMapBase>& dofMap
00043 = DiscreteFunction::discFunc(u0)->map();
00044
00045
00046
00047 PLAYA_ROOT_MSG2(verb, tab1 << "filling discrete function");
00048 Array<int> dofs(1);
00049 for (int i=0; i<mesh.numCells(0); i++)
00050 {
00051 Tabs tab2;
00052 PLAYA_ROOT_MSG3(verb, tab2 << "node i=" << i);
00053 for (int f=0; f<nAttrs; f++)
00054 {
00055 Tabs tab3;
00056
00057 dofMap->getDOFsForCell(0, i, f, dofs);
00058 int dof = dofs[0];
00059 PLAYA_ROOT_MSG3(verb, tab3 << "f=" << f << " dof=" << dof << " val="
00060 << funcVals[f][i]);
00061 loadable(vec)->setElement(dof, funcVals[f][i]);
00062 }
00063 }
00064
00065
00066 DiscreteFunction::discFunc(u0)->setVector(vec);
00067
00068 PLAYA_ROOT_MSG1(verb, tab0 << "done readNodalFields()");
00069 return u0;
00070 }
00071
00072
00073
00074
00075
00076 Expr readSerialGridField(const std::string& gridFile,
00077 double ax, double bx,
00078 double ay, double by,
00079 const VectorType<double>& vecType,
00080 const MeshType& meshType,
00081 Mesh& mesh)
00082 {
00083 ifstream is(gridFile.c_str());
00084 int nNodes ;
00085 int nx;
00086 int ny;
00087 int nAttrs;
00088 is >> nx >> ny >> nAttrs;
00089 nNodes = nx*ny;
00090
00091 Array<Array<double> > funcVals(nAttrs);
00092 for (int i=0; i<nAttrs; i++) funcVals[i].resize(nNodes);
00093 for (int i=0; i<nNodes; i++)
00094 {
00095 for (int f=0; f<nAttrs; f++)
00096 {
00097 is >> funcVals[f][i];
00098 }
00099 }
00100
00101
00102 MeshSource mesher = new PartitionedRectangleMesher(ax, bx, nx-1, 1,
00103 ay, by, ny-1, 1,
00104 meshType);
00105
00106 mesh = mesher.getMesh();
00107
00108
00109
00110
00111 Array<BasisFamily> bas(nAttrs);
00112 for (int i=0; i<bas.size(); i++) bas[i] = new Lagrange(1);
00113 DiscreteSpace discSpace(mesh, bas, vecType);
00114 Expr u0 = new DiscreteFunction(discSpace, 0.0, "u0");
00115
00116
00117
00118
00119 Vector<double> vec = DiscreteFunction::discFunc(u0)->getVector();
00120 const RCP<DOFMapBase>& dofMap
00121 = DiscreteFunction::discFunc(u0)->map();
00122
00123
00124
00125 Array<int> dofs(1);
00126 for (int i=0; i<mesh.numCells(0); i++)
00127 {
00128 for (int f=0; f<nAttrs; f++)
00129 {
00130
00131 dofMap->getDOFsForCell(0, i, f, dofs);
00132 int dof = dofs[0];
00133 loadable(vec)->setElement(dof, funcVals[f][i]);
00134 }
00135 }
00136
00137
00138 DiscreteFunction::discFunc(u0)->setVector(vec);
00139
00140 return u0;
00141 }
00142
00143 double readbackTester(const std::string& infile, const MPIComm& comm, int verb)
00144 {
00145
00146 VectorType<double> vecType = new EpetraVectorType();
00147
00148
00149 Out::root() << "starting to read mesh " << std::endl;
00150 MeshType meshType = new BasicSimplicialMeshType();
00151 MeshSource mesher = new ExodusMeshReader(infile, meshType, verb, comm);
00152 Mesh mesh = mesher.getMesh();
00153
00154 Out::root() << "done reading mesh " << std::endl;
00155 int dim = mesh.spatialDim();
00156
00157 Expr x = new CoordExpr(0);
00158 Expr y = new CoordExpr(1);
00159 Expr z = new CoordExpr(2);
00160
00161
00162
00163 DiscreteSpace discSpace(mesh, new Lagrange(1), vecType);
00164 Out::root() << "done making discfunc " << std::endl;
00165 Expr u;
00166 if (dim==3)
00167 {
00168 L2Projector proj(discSpace, x*sin(x)*y*y+exp(z)*cos(y));
00169 u = proj.project();
00170 }
00171 else if (dim==2)
00172 {
00173 L2Projector proj(discSpace, x*sin(x)+x*x*cos(y));
00174 u = proj.project();
00175 }
00176 else
00177 {
00178 TEUCHOS_TEST_FOR_EXCEPT(dim != 2 && dim != 3);
00179 }
00180 Out::root() << "done projecting function " << std::endl;
00181
00182 CellFilter interior = new MaximalCellFilter();
00183 QuadratureFamily quad = new GaussianQuadrature(2);
00184 Expr F = Integral(interior, u*u, quad);
00185 double fVal = evaluateIntegral(mesh, F);
00186 Out::root() << "done evaluating functional on original mesh" << std::endl;
00187
00188
00189 FieldWriter w = new ExodusWriter("./readbackTest-out");
00190 w.addMesh(mesh);
00191 w.addField("u", new ExprFieldWrapper(u));
00192 w.write();
00193
00194
00195 MeshSource mesher2 = new ExodusMeshReader("./readbackTest-out", meshType, verb, comm);
00196 Mesh mesh2 = mesher2.getMesh();
00197 Expr u2 = readNodalFields(mesher2, mesh2, vecType, verb);
00198 Out::root() << "done readback " << std::endl;
00199
00200
00201 Expr F2 = Integral(interior, u2*u2, quad);
00202 double fVal2 = evaluateIntegral(mesh2, F2);
00203
00204 Out::root() << "functional on original mesh = " << fVal << std::endl;
00205 Out::root() << "functional on readback mesh = " << fVal2 << std::endl;
00206
00207 double diff = std::fabs(fVal - fVal2);
00208 Out::root() << "diff = " << diff << std::endl;
00209
00210 return diff;
00211 }
00212
00213
00214
00215 void invertMap(const Map<int,int>& in, Map<int,int>& out)
00216 {
00217 typedef Map<int,int>::const_iterator iter;
00218
00219 for (iter i=in.begin(); i!=in.end(); i++)
00220 {
00221 out.put(i->second, i->first);
00222 }
00223 }
00224
00225 void serialPartition(
00226 const RCP<SerialPartitionerBase>& part,
00227 int numProc,
00228 const MeshSource& mesher,
00229 const std::string& outfile)
00230 {
00231
00232
00233
00234 TEUCHOS_TEST_FOR_EXCEPTION(mesher.comm().getNProc() > 1, std::runtime_error,
00235 "serialPartition() should only be called from a "
00236 "single-process communicator");
00237
00238
00239 Mesh mesh = mesher.getMesh();
00240 int dim = mesh.spatialDim();
00241
00242
00243 FieldWriter origWriter = new ExodusWriter("orig-writback");
00244 origWriter.addMesh(mesh);
00245 origWriter.write();
00246
00247
00248
00249
00250 Array<Sundance::Map<int, int> > oldElemLIDToNewLIDMaps;
00251 Array<Sundance::Map<int, int> > oldVertLIDToNewLIDMaps;
00252
00253
00254 Array<Mesh> submesh = part->makeMeshParts(mesh, numProc,
00255 oldElemLIDToNewLIDMaps,
00256 oldVertLIDToNewLIDMaps
00257 );
00258
00259
00260 RCP<Array<Array<double> > > oldElemData;
00261 RCP<Array<Array<double> > > oldNodeData;
00262 mesher.getAttributes(oldNodeData, oldElemData);
00263
00264
00265 Array<Set<int> > vertViewers(mesh.numCells(0));
00266 for (int p=0; p<numProc; p++)
00267 {
00268 for (int v=0; v<submesh[p].numCells(0); v++)
00269 {
00270 int gid = submesh[p].mapLIDToGID(0,v);
00271 vertViewers[gid].put(p);
00272 }
00273 }
00274
00275
00276 for (int p=0; p<numProc; p++)
00277 {
00278 Out::os() << "writing part=" << p << " of " << numProc << std::endl;
00279 FieldWriter writer = new ExodusWriter(outfile);
00280
00281 writer.impersonateParallelProc(numProc, p);
00282 writer.addMesh(submesh[p]);
00283
00284
00285 Map<int, int> newElemLIDToOldLIDMap;
00286 Map<int, int> newVertLIDToOldLIDMap;
00287 Out::os() << "preparing maps to submeshes" << std::endl;
00288 invertMap(oldElemLIDToNewLIDMaps[p], newElemLIDToOldLIDMap);
00289 invertMap(oldVertLIDToNewLIDMaps[p], newVertLIDToOldLIDMap);
00290
00291
00292 Out::os() << "mapping element data to submeshes" << std::endl;
00293 RCP<Array<double> > elemData;
00294 int nElemFuncs = oldElemData->size();
00295 int nElems = submesh[p].numCells(dim);
00296 for (int f=0; f<nElemFuncs; f++)
00297 {
00298 elemData = rcp(new Array<double>(nElems));
00299 for (int lid=0; lid<nElems; lid++)
00300 {
00301 int oldLID = newElemLIDToOldLIDMap.get(lid);
00302 (*elemData)[lid] = (*oldElemData)[f][oldLID];
00303 }
00304
00305 writer.addField("uElem["+ Teuchos::toString(f)+"]",
00306 new CellLIDMappedFieldWrapper(dim, 1, elemData));
00307 }
00308
00309
00310 Out::os() << "mapping node data to submeshes" << std::endl;
00311 RCP<Array<double> > nodeData;
00312 int nNodeFuncs = oldNodeData->size();
00313 int nNodes = submesh[p].numCells(0);
00314 for (int f=0; f<nNodeFuncs; f++)
00315 {
00316 nodeData = rcp(new Array<double>(nNodes));
00317 for (int lid=0; lid<nNodes; lid++)
00318 {
00319 int oldLID = newVertLIDToOldLIDMap.get(lid);
00320 (*nodeData)[lid] = (*oldNodeData)[f][oldLID];
00321 }
00322
00323 writer.addField("uNode["+ Teuchos::toString(f)+"]",
00324 new CellLIDMappedFieldWrapper(0, 1, nodeData));
00325 }
00326
00327
00328 Out::os() << "doing write()" << std::endl;
00329 writer.write();
00330 Out::os() << "done part=" << p << " of " << numProc << std::endl;
00331 }
00332
00333
00334 for (int p=0; p<numProc; p++)
00335 {
00336 Out::os() << "writing part=" << p << " of " << numProc << std::endl;
00337
00338 string vvFile = outfile + "-" + Teuchos::toString(numProc)
00339 + "-" + Teuchos::toString(p) + ".pvv";
00340 ofstream of(vvFile.c_str());
00341 of << submesh[p].numCells(0) << endl;
00342 for (int v=0; v<submesh[p].numCells(0); v++)
00343 {
00344 int gid = submesh[p].mapLIDToGID(0, v);
00345 const Set<int>& viewers = vertViewers[gid];
00346 of << v << " " << gid << " " << viewers.size();
00347 for (Set<int>::const_iterator
00348 i=viewers.begin(); i!=viewers.end(); i++)
00349 {
00350 of << " " << *i;
00351 }
00352 of << endl;
00353 }
00354 }
00355 }
00356
00357
00358
00359
00360
00361 }