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