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
00032
00033
00034
00035
00036
00037 #include "SundancePeanoMesh2D.hpp"
00038
00039 #include "SundanceMeshType.hpp"
00040 #include "SundanceCellJacobianBatch.hpp"
00041 #include "SundanceMaximalCofacetBatch.hpp"
00042 #include "SundanceMeshSource.hpp"
00043 #include "SundanceDebug.hpp"
00044 #include "SundanceOut.hpp"
00045 #include "PlayaMPIContainerComm.hpp"
00046 #include "Teuchos_Time.hpp"
00047 #include "Teuchos_TimeMonitor.hpp"
00048 #include "SundanceObjectWithVerbosity.hpp"
00049 #include "SundanceCollectiveExceptionCheck.hpp"
00050
00051 #ifdef HAVE_SUNDANCE_PEANO
00052
00053 using namespace Sundance;
00054 using namespace Teuchos;
00055 using Playa::MPIComm;
00056 using Playa::MPIContainerComm;
00057
00058
00059
00060
00061 Point PeanoMesh2D::returnPoint(0.0 , 0.0);
00062
00063 PeanoMesh2D::PeanoMesh2D(int dim, const MPIComm& comm ,
00064 const MeshEntityOrder& order)
00065 : MeshBase(dim, comm , order),_dimension(dim), _comm(comm)
00066 ,_peanoMesh(NULL)
00067 {
00068 _uniqueResolution = 1.0;
00069 }
00070
00071 void PeanoMesh2D::createMesh(
00072 double position_x,
00073 double position_y,
00074 double offset_x,
00075 double offset_y,
00076 double resolution
00077 ){
00078 double position[2];
00079 double offset[2];
00080 double res[2];
00081
00082
00083 position[0] = position_x; position[1] = position_y;
00084 offset[0] = offset_x; offset[1] = offset_y;
00085 res[0] = resolution; res[1] = resolution;
00086
00087
00088
00089
00090 SUNDANCE_VERB_LOW(" create Peano Mesh ... ");
00091 _dimension = 2;
00092
00093 _peanoMesh = new SundancePeanoInterface2D( position , offset , res );
00094 _uniqueResolution = _peanoMesh->returnResolution(0);
00095
00096 SUNDANCE_VERB_LOW(" Peano Mesh created ... \n");
00097
00098 SUNDANCE_VERB_LOW(" After Plot ... \n");
00099 }
00100
00101
00102 PeanoMesh2D::~PeanoMesh2D() {
00103
00104 }
00105
00106
00107 int PeanoMesh2D::numCells(int dim) const {
00108
00109 return _peanoMesh->numCells(dim);
00110 }
00111
00112 Point PeanoMesh2D::nodePosition(int i) const {
00113
00114
00115 double* coords;
00116 coords = _peanoMesh->nodePositionView(i);
00117
00118 PeanoMesh2D::returnPoint[0] = coords[0];
00119 PeanoMesh2D::returnPoint[1] = coords[1];
00120 return PeanoMesh2D::returnPoint;
00121 }
00122
00123 const double* PeanoMesh2D::nodePositionView(int i) const {
00124
00125
00126 nodePosition(i);
00127 return &(PeanoMesh2D::returnPoint[0]);
00128 }
00129
00130 void PeanoMesh2D::getJacobians(int cellDim, const Array<int>& cellLID,
00131 CellJacobianBatch& jBatch) const
00132 {
00133
00134 SUNDANCE_VERB_HIGH("getJacobians()");
00135 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00136 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00137 int nCells = cellLID.size();
00138 int tmp , tmp_index , tmp_index1;
00139 Point pnt(0.0,0.0);
00140 Point pnt1(0.0,0.0);
00141
00142 jBatch.resize(cellLID.size(), spatialDim(), cellDim);
00143 if (cellDim < spatialDim())
00144 {
00145
00146 for (int i=0; i<nCells; i++)
00147 {
00148
00149 double* detJ = jBatch.detJ(i);
00150 switch(cellDim)
00151 {
00152 case 0: *detJ = 1.0;
00153 break;
00154 case 1:
00155 tmp_index = this->facetLID(cellDim, cellLID[i] , 0 , 0 , tmp );
00156 tmp_index1 = this->facetLID(cellDim, cellLID[i] , 0 , 1 , tmp );
00157 pnt = nodePosition(tmp_index);
00158 pnt1 = nodePosition(tmp_index1);
00159 pnt = pnt1 - pnt;
00160 *detJ = sqrt(pnt * pnt);
00161 break;
00162 default:
00163 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00164 "cellDim=" << cellDim << " in PeanoMesh2D::getJacobians()");
00165 }
00166 }
00167 }else{
00168
00169 SUNDANCE_VERB_HIGH("cellDim == spatialDim()");
00170 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00171 {
00172
00173 double* J = jBatch.jVals(i);
00174 switch(cellDim)
00175 {
00176 case 2:
00177 J[0] = _peanoMesh->returnResolution(0);
00178 J[1] = 0.0; J[2] = 0.0;
00179 J[3] = _peanoMesh->returnResolution(1);
00180 break;
00181 default:
00182 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00183 "cellDim=" << cellDim
00184 << " in SundancePeano2D::getJacobians()");
00185 }
00186 }
00187 }
00188 }
00189
00190 void PeanoMesh2D::getCellDiameters(int cellDim, const Array<int>& cellLID,
00191 Array<double>& cellDiameters) const {
00192 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00193 "cellDim=" << cellDim << " is not in expected range [0, " << spatialDim() << "]");
00194 SUNDANCE_VERB_HIGH("getCellDiameters()");
00195 cellDiameters.resize(cellLID.size());
00196
00197 int tmp , tmp_index , tmp_index1;
00198 Point pnt(0.0,0.0);
00199 Point pnt1(0.0,0.0);
00200
00201 if (cellDim < spatialDim())
00202 {
00203
00204 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00205 {
00206 switch(cellDim)
00207 {
00208 case 0:
00209 cellDiameters[i] = 1.0;
00210 break;
00211 case 1:
00212 tmp_index = this->facetLID(cellDim, cellLID[i] , 0 , 0 , tmp );
00213 tmp_index1= this->facetLID(cellDim, cellLID[i] , 0 , 1 , tmp );
00214 pnt = nodePosition(tmp_index);
00215 pnt1 = nodePosition(tmp_index1);
00216 pnt = pnt1 - pnt;
00217 cellDiameters[i] = sqrt(pnt * pnt);
00218 break;
00219 default:
00220 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00221 "cellDim=" << cellDim << " in PeanoMesh2D::getCellDiameters()");
00222 }
00223 }
00224 }
00225 else
00226 {
00227
00228 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00229 {
00230 switch(cellDim)
00231 {
00232 case 2:
00233 cellDiameters[i] = (_peanoMesh->returnResolution(0)+_peanoMesh->returnResolution(1))/2.0;
00234 break;
00235 default:
00236 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00237 "cellDim=" << cellDim
00238 << " in PeanoMesh2D::getCellDiameters()");
00239 }
00240 }
00241 }
00242 }
00243
00244 void PeanoMesh2D::pushForward(int cellDim, const Array<int>& cellLID,
00245 const Array<Point>& refQuadPts,
00246 Array<Point>& physQuadPts) const {
00247
00248
00249 TEUCHOS_TEST_FOR_EXCEPTION(cellDim < 0 || cellDim > spatialDim(), std::logic_error,
00250 "cellDim=" << cellDim
00251 << " is not in expected range [0, " << spatialDim()
00252 << "]");
00253
00254 int nQuad = refQuadPts.size();
00255 double start_point[2] , end_point[2];
00256
00257 Array<double> J(cellDim*cellDim);
00258
00259 if (physQuadPts.size() > 0) physQuadPts.resize(0);
00260 physQuadPts.reserve(cellLID.size() * refQuadPts.size());
00261 for (unsigned int i=0; i<(unsigned int)cellLID.size(); i++)
00262 {
00263 int lid = cellLID[i];
00264 _peanoMesh->pushForward( cellDim, lid ,start_point , end_point );
00265 Point pnt( start_point[0] , start_point[1] );
00266 Point pnt1( end_point[0] , end_point[1] );
00267 switch(cellDim)
00268 {
00269 case 0:
00270 physQuadPts.append(pnt);
00271 break;
00272 case 1:{
00273 for (int q=0; q<nQuad; q++) {
00274 physQuadPts.append(pnt + (pnt1-pnt)*refQuadPts[q][0]);
00275 }
00276 break;}
00277 case 2:{
00278 for (int q=0; q<nQuad; q++) {
00279 physQuadPts.append( pnt
00280 + Point(refQuadPts[q][0]*_peanoMesh->returnResolution(0),
00281 refQuadPts[q][1]*_peanoMesh->returnResolution(1)));
00282 }
00283 break;}
00284 default:
00285 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "impossible switch value "
00286 "in PeanoMesh2D::getJacobians()");
00287 }
00288 }
00289 }
00290
00291 int PeanoMesh2D::ownerProcID(int cellDim, int cellLID) const {
00292 SUNDANCE_VERB_HIGH("ownerProcID()"); return 0; }
00293
00294
00295 int PeanoMesh2D::numFacets(int cellDim, int cellLID,
00296 int facetDim) const {
00297 SUNDANCE_VERB_HIGH("numFacets()");
00298 return _peanoMesh->numFacets(cellDim, cellLID, facetDim);
00299 }
00300
00301
00302 int PeanoMesh2D::facetLID(int cellDim, int cellLID,
00303 int facetDim, int facetIndex,
00304 int& facetOrientation) const {
00305 int LID;
00306 LID = _peanoMesh->facetLID( cellDim,cellLID, facetDim, facetIndex, facetOrientation);
00307
00308 return LID;
00309 }
00310
00311 void PeanoMesh2D::getFacetLIDs(int cellDim,
00312 const Array<int>& cellLID,
00313 int facetDim,
00314 Array<int>& facetLID,
00315 Array<int>& facetSign) const {
00316 SUNDANCE_VERB_HIGH("getFacetLIDs()");
00317
00318 int LID = 0 , cLID , facetOrientation ;
00319 int ptr = 0;
00320
00321 int nf = numFacets(cellDim, cellLID[0], facetDim);
00322 facetLID.resize(cellLID.size() * nf);
00323 facetSign.resize(cellLID.size() * nf);
00324
00325 for (unsigned int i = 0 ; i < (unsigned int)cellLID.size() ; i++){
00326 cLID = cellLID[i];
00327 for (int f=0; f<nf; f++, ptr++) {
00328
00329 LID = this->facetLID( cellDim, cLID, facetDim, f , facetOrientation);
00330
00331
00332 facetLID[ptr] = LID;
00333 facetSign[ptr] = facetOrientation;
00334 }
00335 }
00336 }
00337
00338 const int* PeanoMesh2D::elemZeroFacetView(int cellLID) const {
00339 return _peanoMesh->elemZeroFacetView(cellLID);
00340 }
00341
00342 int PeanoMesh2D::numMaxCofacets(int cellDim, int cellLID) const {
00343
00344 int coFacetCounter;
00345 coFacetCounter = _peanoMesh->numMaxCofacets( cellDim, cellLID);
00346
00347 return coFacetCounter;
00348 }
00349
00350 int PeanoMesh2D::maxCofacetLID(int cellDim, int cellLID,
00351 int cofacetIndex,
00352 int& facetIndex) const {
00353 int rtn;
00354 rtn = _peanoMesh->maxCofacetLID(cellDim, cellLID, cofacetIndex, facetIndex);
00355
00356
00357 return rtn;
00358 }
00359
00360 void PeanoMesh2D::getMaxCofacetLIDs(const Array<int>& cellLIDs,
00361 MaximalCofacetBatch& cofacets) const {
00362 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error," PeanoMesh2D::getMaxCofacetLIDs() not implemented yet");
00363
00364 }
00365
00366
00367 void PeanoMesh2D::getCofacets(int cellDim, int cellLID,
00368 int cofacetDim, Array<int>& cofacetLIDs) const {
00369 int tmpVect[12] , nrCofacets;
00370 _peanoMesh->getCofacets( cellDim, cellLID, cofacetDim, &tmpVect[0], nrCofacets);
00371 cofacetLIDs.resize(nrCofacets);
00372 for (int ii = 0 ; ii < nrCofacets ; ii++ ) cofacetLIDs[ii] = tmpVect[ii];
00373 }
00374
00375
00376 int PeanoMesh2D::mapGIDToLID(int cellDim, int globalIndex) const {
00377 SUNDANCE_VERB_HIGH("mapGIDToLID()");
00378
00379
00380 return globalIndex;
00381 }
00382
00383 bool PeanoMesh2D::hasGID(int cellDim, int globalIndex) const {
00384 SUNDANCE_VERB_HIGH("hasGID()");
00385
00386
00387 return true;
00388 }
00389
00390 int PeanoMesh2D::mapLIDToGID(int cellDim, int localIndex) const {
00391 SUNDANCE_VERB_HIGH("mapLIDToGID()");
00392
00393
00394 return localIndex;
00395 }
00396
00397 CellType PeanoMesh2D::cellType(int cellDim) const {
00398
00399 switch(cellDim)
00400 {
00401 case 0: return PointCell;
00402 case 1: return LineCell;
00403 case 2: return QuadCell;
00404 case 3: return BrickCell;
00405 default:
00406 return NullCell;
00407 }
00408 }
00409
00410 int PeanoMesh2D::label(int cellDim, int cellLID) const {
00411 return _peanoMesh->label( cellDim, cellLID);
00412 }
00413
00414 void PeanoMesh2D::getLabels(int cellDim, const Array<int>& cellLID,
00415 Array<int>& labels) const {
00416 int tmpIndex;
00417 SUNDANCE_VERB_HIGH("getLabels()");
00418
00419 labels.resize(cellLID.size());
00420
00421 for (tmpIndex = 0 ; tmpIndex < (int)cellLID.size() ; tmpIndex++){
00422 labels[tmpIndex] = _peanoMesh->label( cellDim, cellLID[tmpIndex]);
00423 }
00424 }
00425
00426 Set<int> PeanoMesh2D::getAllLabelsForDimension(int cellDim) const {
00427 Set<int> rtn;
00428 int tmpIndex;
00429 SUNDANCE_VERB_HIGH("getAllLabelsForDimension()");
00430
00431 for (tmpIndex = 0 ; tmpIndex < _peanoMesh->numCells(cellDim) ; tmpIndex++){
00432 rtn.put( _peanoMesh->label( cellDim, tmpIndex) );
00433 }
00434 return rtn;
00435 }
00436
00437 void PeanoMesh2D::getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const {
00438 int tmpIndex , tmpLabel;
00439 SUNDANCE_VERB_HIGH("getLIDsForLabel()");
00440 for (tmpIndex = 0 ; tmpIndex < _peanoMesh->numCells(cellDim) ; tmpIndex++){
00441 tmpLabel = this->label( cellDim , tmpIndex);
00442 if (tmpLabel == label) cellLIDs.append( tmpIndex );
00443 }
00444 }
00445
00446 void PeanoMesh2D::setLabel(int cellDim, int cellLID, int label) {
00447 _peanoMesh->setLabel(cellDim, cellLID, label);
00448 }
00449
00450
00451 void PeanoMesh2D::assignIntermediateCellGIDs(int cellDim) {
00452 SUNDANCE_VERB_HIGH("assignIntermediateCellGIDs()");
00453
00454 }
00455
00456
00457 bool PeanoMesh2D::hasIntermediateGIDs(int dim) const {
00458 SUNDANCE_VERB_HIGH("hasIntermediateGIDs()");
00459 return true;
00460 }
00461
00462
00463 #endif