00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundancePeriodicMesh1D.hpp"
00032
00033 #include "SundanceCellJacobianBatch.hpp"
00034 #include "SundanceMaximalCofacetBatch.hpp"
00035 #include "SundanceDebug.hpp"
00036 #include "SundanceOut.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 #include "Teuchos_Time.hpp"
00039 #include "Teuchos_TimeMonitor.hpp"
00040 #include "SundanceObjectWithVerbosity.hpp"
00041
00042 using namespace Sundance;
00043 using namespace Teuchos;
00044 using namespace std;
00045 using Playa::MPIComm;
00046 using Playa::MPIContainerComm;
00047
00048
00049 PeriodicMesh1D::PeriodicMesh1D(double xMin, double xMax, int numElems)
00050 : MeshBase(1, MPIComm::self(), ExodusMeshOrder),
00051 numElems_(numElems),
00052 xMin_(xMin),
00053 xMax_(xMax),
00054 x_(numElems+1),
00055 verts_(numElems),
00056 labels_(2)
00057 {
00058 labels_[0].resize(numElems_);
00059 labels_[1].resize(numElems_);
00060
00061 for (int i=0; i<numElems_; i++)
00062 {
00063 labels_[0][i] = 0;
00064 labels_[1][i] = 1;
00065 verts_[i].resize(2);
00066 verts_[i][0] = i;
00067 verts_[i][1] = (i+1) % numElems;
00068 }
00069
00070 for (int i=0; i<=numElems_; i++)
00071 {
00072 x_[i] = Point(xMin_ + ((double) i)/((double) numElems_)*(xMax_-xMin_));
00073 }
00074 }
00075
00076 int PeriodicMesh1D::numCells(int cellDim) const
00077 {
00078 switch(cellDim)
00079 {
00080 case 0 :
00081 return numElems_;
00082 case 1:
00083 return numElems_;
00084 default:
00085 TEUCHOS_TEST_FOR_EXCEPT(true);
00086 }
00087 return -1;
00088 }
00089
00090
00091 Point PeriodicMesh1D::nodePosition(int i) const
00092 {
00093 return x_[i];
00094 }
00095
00096
00097 const double* PeriodicMesh1D::nodePositionView(int i) const
00098 {
00099 return &(x_[i][0]);
00100 }
00101
00102 void PeriodicMesh1D::getJacobians(int cellDim, const Array<int>& cellLID,
00103 CellJacobianBatch& jBatch) const
00104 {
00105 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00106 "cellDim=" << cellDim
00107 << " is not in expected range [0, " << spatialDim()
00108 << "]");
00109
00110 jBatch.resize(cellLID.size(), spatialDim(), cellDim);
00111
00112 int nCells = cellLID.size();
00113
00114 if (cellDim==0)
00115 {
00116 for (int i=0; i<nCells; i++)
00117 {
00118 double* detJ = jBatch.detJ(i);
00119 *detJ = 1.0;
00120 }
00121 }
00122 else
00123 {
00124 for (int i=0; i<nCells; i++)
00125 {
00126 int lid = cellLID[i];
00127 double* J = jBatch.jVals(i);
00128 double x0 = x_[lid][0];
00129 double x1 = x_[lid+1][0];
00130 J[0] = fabs(x1-x0);
00131 }
00132 }
00133 }
00134
00135
00136 void PeriodicMesh1D::getCellDiameters(int cellDim, const Array<int>& cellLID,
00137 Array<double>& cellDiameters) const
00138 {
00139 cellDiameters.resize(cellLID.size());
00140
00141 TEUCHOS_TEST_FOR_EXCEPT(cellDim < 0 || cellDim > 1);
00142
00143 if (cellDim==1)
00144 {
00145 for (int i=0; i<cellLID.size(); i++)
00146 {
00147 int c = cellLID[i];
00148 double h = x_[c+1][0]-x_[c][0];
00149 cellDiameters[i] = h;
00150 }
00151 }
00152 else
00153 {
00154 for (int i=0; i<cellLID.size(); i++)
00155 {
00156 cellDiameters[i] = 1.0;
00157 }
00158 }
00159 }
00160
00161 void PeriodicMesh1D::pushForward(int cellDim, const Array<int>& cellLID,
00162 const Array<Point>& refQuadPts,
00163 Array<Point>& physQuadPts) const
00164 {
00165 TEUCHOS_TEST_FOR_EXCEPT(cellDim < 0 || cellDim > 1);
00166
00167 if (cellDim==1)
00168 {
00169 if (physQuadPts.size() > 0) physQuadPts.resize(0);
00170 physQuadPts.reserve(refQuadPts.size() * cellLID.size());
00171
00172 for (int i=0; i<cellLID.size(); i++)
00173 {
00174 int c = cellLID[i];
00175 Point h = x_[c+1]-x_[c];
00176 for (int q=0; q<refQuadPts.size(); q++)
00177 {
00178 physQuadPts.append(Point(x_[c] + refQuadPts[q][0] * h));
00179 }
00180 }
00181 }
00182 else
00183 {
00184 for (int i=0; i<cellLID.size(); i++)
00185 {
00186 physQuadPts.append(x_[cellLID[i]]);
00187 }
00188 }
00189 }
00190
00191 int PeriodicMesh1D::numFacets(int cellDim, int cellLID,
00192 int facetDim) const
00193 {
00194 TEUCHOS_TEST_FOR_EXCEPT(cellLID < 0 || cellLID >= numElems_);
00195 if (cellDim == 1 && facetDim==0) return 2;
00196 return 0;
00197 }
00198
00199
00200
00201 int PeriodicMesh1D::facetLID(int cellDim, int cellLID,
00202 int facetDim, int facetIndex,
00203 int& facetOrientation) const
00204 {
00205 TEUCHOS_TEST_FOR_EXCEPT(cellLID < 0 || cellLID >= numElems_);
00206
00207 TEUCHOS_TEST_FOR_EXCEPT(cellDim != 1);
00208 TEUCHOS_TEST_FOR_EXCEPT(facetDim != 0);
00209 TEUCHOS_TEST_FOR_EXCEPT(facetIndex < 0);
00210 TEUCHOS_TEST_FOR_EXCEPT(facetIndex > 1);
00211
00212 return verts_[cellLID][facetIndex];
00213 }
00214
00215 void PeriodicMesh1D::getFacetLIDs(int cellDim,
00216 const Array<int>& cellLID,
00217 int facetDim,
00218 Array<int>& facetLID,
00219 Array<int>& facetSign) const
00220 {
00221 facetLID.resize(2*cellLID.size());
00222 facetSign.resize(2*cellLID.size());
00223
00224 for (int i=0; i<cellLID.size(); i++)
00225 {
00226 facetLID[2*i] = this->facetLID(cellDim, cellLID[i], facetDim, 0, facetSign[2*i]);
00227 facetLID[2*i+1] = this->facetLID(cellDim, cellLID[i], facetDim, 1, facetSign[2*i]);
00228 }
00229 }
00230
00231
00232 const int* PeriodicMesh1D::elemZeroFacetView(int cellLID) const
00233 {
00234 return &(verts_[cellLID][0]);
00235 }
00236
00237
00238 int PeriodicMesh1D::numMaxCofacets(int cellDim, int cellLID) const
00239 {
00240 TEUCHOS_TEST_FOR_EXCEPT(cellDim != 0);
00241 return 2;
00242 }
00243
00244 int PeriodicMesh1D::maxCofacetLID(int cellDim, int cellLID,
00245 int cofacetIndex,
00246 int& facetIndex) const
00247 {
00248 TEUCHOS_TEST_FOR_EXCEPT(cellDim != 0);
00249 if (cofacetIndex==0) return (cellLID-1) % numElems_;
00250 return cellLID;
00251 }
00252
00253 void PeriodicMesh1D::getMaxCofacetLIDs(const Array<int>& cellLIDs,
00254 MaximalCofacetBatch& cofacets) const
00255 {
00256 TEUCHOS_TEST_FOR_EXCEPT(true);
00257 }
00258
00259 void PeriodicMesh1D::getCofacets(int cellDim, int cellLID,
00260 int cofacetDim, Array<int>& cofacetLIDs) const
00261 {
00262 TEUCHOS_TEST_FOR_EXCEPT(cellDim != 0);
00263 TEUCHOS_TEST_FOR_EXCEPT(cofacetDim != 1);
00264
00265 cofacetLIDs.resize(2);
00266 cofacetLIDs[0] = (cellLID-1) % numElems_;
00267 cofacetLIDs[1] = cellLID;
00268 }
00269
00270 int PeriodicMesh1D::mapGIDToLID(int cellDim, int globalIndex) const
00271 {
00272 return globalIndex;
00273 }
00274
00275 bool PeriodicMesh1D::hasGID(int cellDim, int globalIndex) const
00276 {
00277 return true;
00278 }
00279
00280 int PeriodicMesh1D::mapLIDToGID(int cellDim, int localIndex) const
00281 {
00282 return localIndex;
00283 }
00284
00285 CellType PeriodicMesh1D::cellType(int cellDim) const
00286 {
00287 if (cellDim==0) return PointCell;
00288 else if (cellDim==1) return LineCell;
00289 else return NullCell;
00290 }
00291
00292 int PeriodicMesh1D::label(int cellDim, int cellLID) const
00293 {
00294 return labels_[cellDim][cellLID];
00295 }
00296
00297 void PeriodicMesh1D::getLabels(int cellDim, const Array<int>& cellLID,
00298 Array<int>& labels) const
00299 {
00300 labels.resize(cellLID.size());
00301 const Array<int>& ld = labels_[cellDim];
00302 for (int i=0; i<cellLID.size(); i++)
00303 {
00304 labels[i] = ld[cellLID[i]];
00305 }
00306 }
00307
00308 Set<int> PeriodicMesh1D::getAllLabelsForDimension(int cellDim) const
00309 {
00310 Set<int> rtn;
00311
00312 const Array<int>& ld = labels_[cellDim];
00313 for (int i=0; i<ld.size(); i++)
00314 {
00315 rtn.put(ld[i]);
00316 }
00317
00318 return rtn;
00319 }
00320
00321 void PeriodicMesh1D::setLabel(int cellDim, int cellLID, int label)
00322 {
00323 labels_[cellDim][cellLID] = label;
00324 }
00325
00326
00327 void PeriodicMesh1D::getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const
00328 {
00329 cellLIDs.resize(0);
00330 const Array<int>& ld = labels_[cellDim];
00331 for (int i=0; i<ld.size(); i++)
00332 {
00333 if (ld[i] == label) cellLIDs.append(i);
00334 }
00335 }