00001 #include "SundancePartitionedRectangleMesher.hpp"
00002 #include "SundanceOut.hpp"
00003
00004 using namespace Sundance;
00005 using namespace Sundance;
00006
00007 using namespace Teuchos;
00008 using namespace Sundance;
00009
00010 PartitionedRectangleMesher::PartitionedRectangleMesher(const ParameterList& params)
00011 : MeshSourceBase(params),
00012 ax_(params.get<double>("ax")),
00013 bx_(params.get<double>("bx")),
00014 nx_(params.get<int>("nx")),
00015 npx_(1),
00016 ay_(params.get<double>("ay")),
00017 by_(params.get<double>("by")),
00018 ny_(params.get<int>("ny")),
00019 npy_(1)
00020 {
00021 if (params.isParameter("npx"))
00022 {
00023 npx_ = params.get<int>("npx");
00024 }
00025 if (params.isParameter("npy"))
00026 {
00027 npy_ = params.get<int>("npy");
00028 }
00029 }
00030
00031
00032 void PartitionedRectangleMesher::balanceXY(int n, int* npx, int* npy)
00033 {
00034 int m = (int) floor(sqrt((double)n));
00035 for (int i=m; i>=1; i--)
00036 {
00037 if (n % i == 0)
00038 {
00039 *npx = i;
00040 *npy = n/i;
00041 return ;
00042 }
00043 }
00044
00045 *npx = n;
00046 *npy = 1;
00047 }
00048
00049
00050 Mesh PartitionedRectangleMesher::fillMesh() const
00051 {
00052 SUNDANCE_OUT(this->verb() > 0,
00053 "PartitionedRectangleMesher::fillLocalMesh() is meshing "
00054 "rectangle [" << ax_ << ", " << bx_ << "] by ["
00055 << ay_ << ", " << by_ << "]");
00056
00057 SUNDANCE_OUT(this->verb() >= 3,
00058 "PartitionedRectangleMesher::fillLocalMesh() starting creation "
00059 "of empty mesh");
00060
00061 Mesh mesh = createMesh(2);
00062
00063 SUNDANCE_OUT(this->verb() >= 3,
00064 "PartitionedRectangleMesher::fillLocalMesh() done creation of "
00065 "empty mesh");
00066
00067
00068
00069 int np = nProc();
00070 int rank = myRank();
00071
00072 TEUCHOS_TEST_FOR_EXCEPTION(npx_ * npy_ != np, std::runtime_error,
00073 "PartitionedRectangleMesher::fillLocalMesh(): product "
00074 "of npx=" << npx_ << " and npy=" << npy_
00075 << " is not equal to np=" << np);
00076
00077
00078
00079 int nppx = nx_;
00080 int nppy = ny_;
00081 int nxTot = nx_*npx_;
00082 int nyTot = ny_*npy_;
00083
00084 int px = rank/npy_;
00085 int py = rank % npy_;
00086
00087 int lowestVisiblePtX = px*nppx-1;
00088 int lowestVisiblePtY = py*nppy-1;
00089 if (lowestVisiblePtX < 0) lowestVisiblePtX = 0;
00090 if (lowestVisiblePtY < 0) lowestVisiblePtY = 0;
00091
00092
00093 int highestVisiblePtX = lowestVisiblePtX + nppx + 1;
00094 int highestVisiblePtY = lowestVisiblePtY + nppy + 1;
00095 if (highestVisiblePtX > nxTot) highestVisiblePtX = nxTot;
00096 if (highestVisiblePtY > nyTot) highestVisiblePtY = nyTot;
00097
00098
00099 Array<Array<int> > pts(highestVisiblePtX-lowestVisiblePtX+1);
00100 for (int i=0; i<pts.size(); i++) pts[i].resize(highestVisiblePtY-lowestVisiblePtY+1);
00101
00102 int globalIndex = 0;
00103
00104
00105 for (int i=0; i<=nxTot; i++)
00106 {
00107 for (int j=0; j<=nyTot; j++, globalIndex++)
00108 {
00109 if (i < lowestVisiblePtX || i > highestVisiblePtX) continue;
00110 if (j < lowestVisiblePtY || j > highestVisiblePtY) continue;
00111
00112 int ip = i/nppx;
00113 if (i==nxTot) ip--;
00114 int jp = j/nppy;
00115 if (j==nyTot) jp--;
00116 int pointOwner = ip*npy_ + jp;
00117
00118 Point x( ax_ + ((double) i)*(bx_-ax_)/((double) nxTot) ,
00119 ay_ + ((double) j)*(by_-ay_)/((double) nyTot));
00120 SUNDANCE_OUT(this->verb() > 1, "adding point GID="
00121 << globalIndex << " x=" << x << " owner="
00122 << pointOwner);
00123 int lid = mesh.addVertex(globalIndex, x, pointOwner, 0);
00124 pts[i-lowestVisiblePtX][j-lowestVisiblePtY] = globalIndex;
00125
00126 SUNDANCE_OUT(this->verb() >= 3,
00127 "point " << x << " registered with LID=" << lid);
00128 }
00129 }
00130
00131
00132 globalIndex = 0 ;
00133
00134 for (int i=0; i<nxTot; i++)
00135 {
00136 for (int j=0; j<nyTot; j++, globalIndex+=2)
00137 {
00138 if (i < lowestVisiblePtX || i >= highestVisiblePtX) continue;
00139 if (j < lowestVisiblePtY || j >= highestVisiblePtY) continue;
00140
00141 int a = pts[i-lowestVisiblePtX][j-lowestVisiblePtY];
00142 int b = pts[i+1-lowestVisiblePtX][j-lowestVisiblePtY];
00143 int c = pts[i+1-lowestVisiblePtX][j+1-lowestVisiblePtY];
00144 int d = pts[i-lowestVisiblePtX][j+1-lowestVisiblePtY];
00145
00146 int ip = i/nppx;
00147 int jp = j/nppy;
00148 int cellOwner = ip*npy_ + jp;
00149 Array<int> tri1;
00150 Array<int> tri2;
00151 if (i%2 == j%2)
00152 {
00153 tri1 = tuple(a,b,c);
00154 tri2 = tuple(a,c,d);
00155 }
00156 else
00157 {
00158 tri1 = tuple(a,b,d);
00159 tri2 = tuple(b,c,d);
00160 }
00161 int lid1 = mesh.addElement(globalIndex, tri1, cellOwner, 0);
00162 SUNDANCE_OUT(this->verb() >= 3,
00163 "elem " << tri1
00164 << " registered with LID=" << lid1);
00165 int lid2 = mesh.addElement(globalIndex+1, tri2, cellOwner, 0);
00166 SUNDANCE_OUT(this->verb() >= 3,
00167 "elem " << tri2
00168 << " registered with LID=" << lid2);
00169
00170 }
00171 }
00172
00173 mesh.freezeTopology();
00174
00175 return mesh;
00176 }