00001 #include "SundancePartitionedLineMesher.hpp"
00002 #include "SundanceOut.hpp"
00003 #include "SundanceCollectiveExceptionCheck.hpp"
00004
00005 using namespace Sundance;
00006 using namespace Sundance;
00007
00008 using namespace Teuchos;
00009 using namespace Sundance;
00010
00011
00012 PartitionedLineMesher::PartitionedLineMesher(const ParameterList& params)
00013 : MeshSourceBase(params),
00014 ax_(params.get<double>("ax")),
00015 bx_(params.get<double>("bx")),
00016 nx_(params.get<int>("nx"))
00017 {;}
00018
00019 Mesh PartitionedLineMesher::fillMesh() const
00020 {
00021 Mesh mesh;
00022
00023 try
00024 {
00025 SUNDANCE_OUT(this->verb() > 0,
00026 "PartitionedLineMesher::fillLocalMesh() is meshing "
00027 "interval [" << ax_ << ", " << bx_ << "]");
00028
00029 mesh = createMesh(1);
00030
00031
00032
00033
00034 int np = nProc();
00035 int nppx = nx_/np;
00036
00037 SUNDANCE_OUT(this->verb() > 0,
00038 "PartitionedLineMesher::fillLocalMesh() has " << nppx
00039 << " points per proc");
00040
00041 int px = myRank();
00042
00043 int lowestVisiblePtX = px*nppx-1;
00044 if (lowestVisiblePtX < 0) lowestVisiblePtX = 0;
00045
00046 int highestVisiblePtX = lowestVisiblePtX + nppx + 1;
00047 if (highestVisiblePtX > nx_) highestVisiblePtX = nx_;
00048
00049 SUNDANCE_OUT(this->verb() > 0,
00050 "index range is [" << lowestVisiblePtX << ", " <<
00051 highestVisiblePtX << "]");
00052
00053 Array<int> pts(highestVisiblePtX-lowestVisiblePtX+1);
00054 int globalIndex = 0;
00055
00056
00057 for (int i=0; i<=nx_; i++, globalIndex++)
00058 {
00059 if (i < lowestVisiblePtX || i > highestVisiblePtX) continue;
00060 int pointOwner = i/nppx;
00061 if (i==nx_) pointOwner--;
00062 Point x( ax_ + ((double) i)*(bx_-ax_)/((double) nx_));
00063
00064 SUNDANCE_OUT(this->verb() > 1, "adding point GID="
00065 << globalIndex << " x=" << x << " owner=" << pointOwner);
00066 int lid = mesh.addVertex(globalIndex, x, pointOwner, 0);
00067 pts[i-lowestVisiblePtX] = globalIndex;
00068 SUNDANCE_OUT(this->verb() > 3,
00069 "point " << x << " registered with LID=" << lid);
00070 }
00071
00072
00073 globalIndex = 0 ;
00074
00075 for (int i=0; i<nx_; i++, globalIndex++)
00076 {
00077 if (i < lowestVisiblePtX || i >= highestVisiblePtX) continue;
00078 int a = pts[i-lowestVisiblePtX];
00079 int b = pts[i-lowestVisiblePtX+1];
00080 int cellOwner = i/nppx;
00081 SUNDANCE_OUT(this->verb() > 1, "adding elem GID="
00082 << globalIndex << " nodes=" << tuple(a,b)
00083 << " owner=" << cellOwner);
00084
00085 int lid = mesh.addElement(globalIndex, tuple(a,b), cellOwner, 0);
00086 SUNDANCE_OUT(this->verb() > 3,
00087 "elem " << tuple(a,b) << " registered with LID=" << lid);
00088 }
00089
00090 mesh.freezeTopology();
00091
00092 if (px==0) mesh.setLabel(0, 0, 1);
00093 if (px==np-1) mesh.setLabel(0, mesh.mapGIDToLID(0, nx_), 2);
00094
00095
00096
00097 }
00098 catch(std::exception& e0)
00099 {
00100 reportFailure(comm());
00101 SUNDANCE_TRACE_MSG(e0, "while meshing a line");
00102 }
00103 TEUCHOS_TEST_FOR_EXCEPTION(checkForFailures(comm()), std::runtime_error,
00104 "off-proc error detected on proc=" << myRank()
00105 << " while meshing line");
00106 return mesh;
00107
00108 }