#include "domain_algorithm_extract.hpp" #include "extract_domain.hpp" #include "domain.hpp" #include "grid.hpp" #include "grid_transformation_factory_impl.hpp" #include "attribute_template.hpp" namespace xios { shared_ptr CDomainAlgorithmExtract::create(bool isSource, CGrid* gridDst, CGrid* gridSrc, CTransformation* transformation, int elementPositionInGrid, std::map& elementPositionInGridSrc2ScalarPosition, std::map& elementPositionInGridSrc2AxisPosition, std::map& elementPositionInGridSrc2DomainPosition, std::map& elementPositionInGridDst2ScalarPosition, std::map& elementPositionInGridDst2AxisPosition, std::map& elementPositionInGridDst2DomainPosition) TRY { std::vector domainListDestP = gridDst->getDomains(); std::vector domainListSrcP = gridSrc->getDomains(); CExtractDomain* extractDomain = dynamic_cast (transformation); int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid]; int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid]; return make_shared(isSource, domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], extractDomain); } CATCH bool CDomainAlgorithmExtract::dummyRegistered_ = CDomainAlgorithmExtract::registerTrans(); bool CDomainAlgorithmExtract::registerTrans() TRY { return CGridTransformationFactory::registerTransformation(TRANS_EXTRACT_DOMAIN, create); } CATCH CDomainAlgorithmExtract::CDomainAlgorithmExtract(bool isSource, CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain) : CAlgorithmTransformationTransfer(isSource), domainSrc_(domainSource), domainDest_(domainDestination) TRY { // Reset geometrical attributes to avoid incompatible (user/domainSource) attributs // attributs will be defined using domainSource and/or transformation attributs domainDestination->type.reset(); domainDestination->ni_glo.reset(); domainDestination->nj_glo.reset(); domainDestination->i_index.reset(); // defined using domainSource->getLocalElement() domainDestination->j_index.reset(); // " domainDestination->ibegin.reset(); // will be computed in domainDestination->checkDomain() (from checkAttributes()) domainDestination->ni.reset(); // " domainDestination->jbegin.reset(); // " domainDestination->nj.reset(); // " domainDestination->mask_1d.reset(); // defined scanning domainSource->getFullView() & domainSource->getWorkflowView() differencies domainDestination->mask_2d.reset(); // in all case domainDestination->mask_1d used as reference // domainDestination->data_* attributes will be computed in : domainDestination->data_dim.reset(); // domainDestination->checkDomainData() (from checkAttributes()) domainDestination->data_ni.reset(); domainDestination->data_nj.reset(); domainDestination->data_ibegin.reset(); domainDestination->data_jbegin.reset(); domainDestination->data_i_index.reset(); // domainDestination->checkCompression() (from checkAttributes()) domainDestination->data_j_index.reset(); // Next attributes will be set using domainSource->attributes domainDestination->lonvalue_1d.reset(); domainDestination->latvalue_1d.reset(); domainDestination->lonvalue_2d.reset(); domainDestination->latvalue_2d.reset(); domainDestination->nvertex.reset(); domainDestination->bounds_lon_1d.reset(); domainDestination->bounds_lat_1d.reset(); domainDestination->bounds_lon_2d.reset(); domainDestination->bounds_lat_2d.reset(); domainDestination->area.reset(); domainDestination->radius.reset(); extractDomain->checkValid(domainSource); extractIBegin_ = extractDomain->ibegin.getValue(); extractJBegin_ = extractDomain->jbegin.getValue(); extractNi_ = extractDomain->ni.getValue(); extractNj_ = extractDomain->nj.getValue(); extractIEnd_ = extractIBegin_ + extractNi_ - 1; extractJEnd_ = extractJBegin_ + extractNj_ - 1; if (extractNi_ > domainSource->ni_glo.getValue()) { ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)", << "Extract size is greater than size of domain source" << "Size ni_glo of domain source " <getId() << " is " << domainSource->ni_glo.getValue() << std::endl << "Extract size is " << extractNi_ ); } if (extractNj_ > domainSource->nj_glo.getValue()) { ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)", << "Extract size is greater than size of domain source" << "Size nj_glo of domain source " <getId() << " is " << domainSource->nj_glo.getValue() << std::endl << "Extract size is " << extractNj_ ); } // Calculate the size of local domain int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0 ; int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc; for (int j = 0; j < domainSrc_->nj.getValue(); j++) { for (int i = 0; i < domainSrc_->ni.getValue(); i++) { ind = j*domainSrc_->ni + i; iIdxSrc = domainSrc_->i_index(ind); if ((iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_)) { jIdxSrc = domainSrc_->j_index(ind); if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_)) { if ((niDest == 0) && (njDest == 0)) { destIBegin = i; destJBegin = j; } if (i == destIBegin) ++njDest; } if (j == destJBegin) ++niDest; } } } // Set attributes for this transformation domainDest_->type = domainSrc_ -> type ; domainDest_->ni_glo.setValue(extractNi_); domainDest_->nj_glo.setValue(extractNj_); domainDest_->i_index.resize(niDest*njDest); domainDest_->j_index.resize(niDest*njDest); // Resize lon/lat, bounds, area arrays to local domain dimensions if (!domainSrc_->lonvalue_1d.isEmpty()) { if (domainDest_->type == CDomain::type_attr::rectilinear) { domainDest_->lonvalue_1d.resize(niDest); domainDest_->latvalue_1d.resize(njDest); } else if (domainDest_->type == CDomain::type_attr::unstructured) { domainDest_->lonvalue_1d.resize(niDest); domainDest_->latvalue_1d.resize(niDest); } } else if (!domainSrc_->lonvalue_2d.isEmpty()) { domainDest_->lonvalue_2d.resize(niDest,njDest); domainDest_->latvalue_2d.resize(niDest,njDest); } if (domainSrc_->hasBounds) { if (!domainSrc_->bounds_lon_2d.isEmpty()) { domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest); domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest); } else if (!domainSrc_->bounds_lon_1d.isEmpty()) { domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest); domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest); } } if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest); // Set attributes required to define domainDestination->localElement_ and associated views, full and workflow) CArray sourceGlobalIdx = domainSource->getLocalElement()->getGlobalIndex(); int indexSize = sourceGlobalIdx.numElements(); domainDest_->data_i_index.resize(niDest*njDest); domainDestination->data_i_index = -1; domainDest_->data_j_index.resize(niDest*njDest); domainDestination->data_j_index = 0; CArray sourceWorkflowIdx = domainSource->getLocalView(CElementView::WORKFLOW)->getIndex(); int srcWorkflowSize = sourceWorkflowIdx.numElements(); int iIdxSrcMin = INT_MAX; int jIdxSrcMin = INT_MAX; int IdxMin = INT_MAX; for (int countSrc = 0; countSrc < indexSize ; ++countSrc) { if ( sourceGlobalIdx(countSrc)%domainSource->ni_glo < iIdxSrcMin ) iIdxSrcMin = sourceGlobalIdx(countSrc)%domainSource->ni_glo; if ( sourceGlobalIdx(countSrc)/domainSource->ni_glo < jIdxSrcMin ) jIdxSrcMin = sourceGlobalIdx(countSrc)/domainSource->ni_glo; if ( sourceGlobalIdx(countSrc) < IdxMin ) IdxMin = sourceGlobalIdx(countSrc); } int countDest(0); // increment of the position in destination domain for (int countSrc = 0; countSrc < indexSize ; ++countSrc) { int iIdxSrc = sourceGlobalIdx(countSrc)%domainSource->ni_glo; int jIdxSrc = sourceGlobalIdx(countSrc)/domainSource->ni_glo; // check that point countSrc concerned by extract if ( (iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_) && (jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_) ) { // if concerned, convert source the global indexation in the extracted frame domainDest_->i_index(countDest) = iIdxSrc-extractIBegin_; domainDest_->j_index(countDest) = jIdxSrc-extractJBegin_; // ------------------ define transformation only if in the WF ------------------ // countSrc+IdxMin is the global position (not index) considering the distributed memory // - can be compared to the workflow view int iIdxSrc2 = (countSrc+IdxMin)%domainSource->ni_glo; int jIdxSrc2 = (countSrc+IdxMin)/domainSource->ni_glo; int convert_locally_global_idx = (jIdxSrc2-jIdxSrcMin)*domainSource->ni + (iIdxSrc2-iIdxSrcMin) ; bool concerned_by_WF(false); for ( int i = 0 ; idata_i_index( countDest ) = countDest; } // ----------------------------------------------------------------------------- int iIdxDestLocal = countDest%niDest; int jIdxDestLocal = countDest/niDest; int iIdxSrcLocal = countSrc%domainSource->ni; int jIdxSrcLocal = countSrc/domainSource->ni; // area if (!domainSrc_->area.isEmpty()) { domainDest_->area(iIdxDestLocal,jIdxDestLocal) = domainSrc_->area(iIdxSrcLocal,jIdxSrcLocal); } // bounds if (!domainDest_->bounds_lon_1d.isEmpty()) { for (int n = 0; n < domainSrc_->nvertex; ++n) { domainDest_->bounds_lon_1d(n, countDest) = domainSrc_->bounds_lon_1d(n,countSrc); domainDest_->bounds_lat_1d(n, countDest) = domainSrc_->bounds_lat_1d(n,countSrc); } } else if (!domainDest_->bounds_lon_2d.isEmpty()) { for (int n = 0; n < domainSrc_->nvertex; ++n) { domainDest_->bounds_lon_2d(n, iIdxDestLocal, jIdxDestLocal) = domainSrc_->bounds_lon_2d(n, iIdxSrcLocal, jIdxSrcLocal); domainDest_->bounds_lat_2d(n, iIdxDestLocal, jIdxDestLocal) = domainSrc_->bounds_lat_2d(n, iIdxSrcLocal, jIdxSrcLocal); } } // lon/lat if (!domainDest_->lonvalue_1d.isEmpty()) { if (domainDest_->type == CDomain::type_attr::rectilinear) { // i : scan nbr of points in src domainDest_->lonvalue_1d(iIdxDestLocal) = domainSrc_->lonvalue_1d(iIdxSrcLocal); domainDest_->latvalue_1d(jIdxDestLocal) = domainSrc_->latvalue_1d(jIdxSrcLocal); } else if (domainDest_->type == CDomain::type_attr::unstructured) { domainDest_->lonvalue_1d(countDest) = domainSrc_->lonvalue_1d(countSrc); domainDest_->latvalue_1d(countDest) = domainSrc_->latvalue_1d(countSrc); } } else if (!domainDest_->lonvalue_2d.isEmpty()) { if (domainDest_->type == CDomain::type_attr::curvilinear) { domainDest_->lonvalue_2d(iIdxDestLocal, jIdxDestLocal) = domainSrc_->lonvalue_2d(iIdxSrcLocal,jIdxSrcLocal); domainDest_->latvalue_2d(iIdxDestLocal, jIdxDestLocal) = domainSrc_->latvalue_2d(iIdxSrcLocal,jIdxSrcLocal); } } // if point i has been identified as extracted, increment position in destination domain for the next point countDest++; } } domainDestination->checkAttributes() ; this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ; } CATCH }