X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDLoader%2FMEDFileField.cxx;h=dc81062983ab4c36394b02db050c007211927003;hb=662a2a2393a25baef77e42f74204b11b70a9646c;hp=f4c906f8b77c57c8c5d7377555d8c214258e6809;hpb=68291f2f5498143761e7d2b37e4113fc976517d3;p=tools%2Fmedcoupling.git diff --git a/src/MEDLoader/MEDFileField.cxx b/src/MEDLoader/MEDFileField.cxx index f4c906f8b..dc8106298 100644 --- a/src/MEDLoader/MEDFileField.cxx +++ b/src/MEDLoader/MEDFileField.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D +// Copyright (C) 2007-2023 CEA, EDF // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -35,9 +35,10 @@ #include #include +// From MEDLoader.cxx TU: extern INTERP_KERNEL::NormalizedCellType typmai2[MED_N_CELL_FIXED_GEO]; extern med_geometry_type typmainoeud[1]; -extern med_geometry_type typmai3[34]; +extern med_geometry_type typmai3[INTERP_KERNEL::NORM_MAXTYPE]; using namespace MEDCoupling; @@ -169,7 +170,7 @@ std::vector< std::pair > MEDFileFields::getCommonIterations(bool& areTh int MEDFileFields::getNumberOfFields() const { - return _fields.size(); + return (int)_fields.size(); } std::vector MEDFileFields::getFieldsNames() const @@ -234,7 +235,7 @@ void MEDFileFields::simpleRepr(int bkOffset, std::ostream& oss) const for(std::vector< MCAuto >::const_iterator it=_fields.begin();it!=_fields.end();it++,i++) { const MEDFileAnyTypeFieldMultiTSWithoutSDA *cur=(*it); - std::string chapter(17,'0'+i); + std::string chapter(17,(char)('0'+i)); oss << startLine << chapter << std::endl; if(cur) { @@ -256,7 +257,7 @@ MEDFileFields::MEDFileFields() MEDFileFields::MEDFileFields(med_idt fid, bool loadAll, const MEDFileMeshes *ms, const MEDFileEntities *entities) try:MEDFileFieldGlobsReal(fid) { - int nbFields(MEDnField(fid)); + med_int nbFields(MEDnField(fid)); _fields.resize(nbFields); med_field_type typcha; for(int i=0;i b(_fields.size(),true); int k=bg; for(int i=0;i& oldCode, const std::vector& newCode, const DataArrayInt *renumO2N) +bool MEDFileFields::renumberEntitiesLyingOnMesh(const std::string& meshName, const std::vector& oldCode, const std::vector& newCode, const DataArrayIdType *renumO2N) { bool ret(false); for(std::vector< MCAuto >::iterator it=_fields.begin();it!=_fields.end();it++) @@ -576,7 +582,7 @@ bool MEDFileFields::renumberEntitiesLyingOnMesh(const std::string& meshName, con * * \return A new object that the caller is responsible to deallocate. */ -MEDFileFields *MEDFileFields::extractPart(const std::map >& extractDef, MEDFileMesh *mm) const +MEDFileFields *MEDFileFields::extractPart(const std::map >& extractDef, MEDFileMesh *mm) const { if(!mm) throw INTERP_KERNEL::Exception("MEDFileFields::extractPart : input mesh is NULL !"); @@ -607,6 +613,235 @@ void MEDFileFields::accept(MEDFileFieldVisitor& visitor) const } } +class PFLData +{ +public: + PFLData():_add_pts_in_pfl(0) { } + PFLData(const MCAuto& mat, const MCAuto& pfl, mcIdType nbOfNewPts):_matrix(mat),_pfl(pfl),_add_pts_in_pfl(nbOfNewPts) { } + std::string getPflName() const { if(_pfl.isNull()) { return std::string(); } else { return _pfl->getName(); } } + mcIdType getNbOfAddPtsInPfl() const { return _add_pts_in_pfl; } + MCAuto getProfile() const { return _pfl; } + MCAuto getMatrix() const { return _matrix; } +private: + MCAuto _matrix; + MCAuto _pfl; + mcIdType _add_pts_in_pfl; +}; + +class MEDFileFieldLin2QuadVisitor : public MEDFileFieldVisitor +{ +public: + MEDFileFieldLin2QuadVisitor(const MEDFileUMesh *lin, const MEDFileUMesh *quad, const MEDFileFieldGlobsReal *linGlobs, MEDFileFields* outFs):_lin(lin),_quad(quad),_lin_globs(linGlobs),_out_fs(outFs),_gt(INTERP_KERNEL::NORM_ERROR),_1ts_update_requested(false) { } + void newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field) { if(field->getMeshName()!=_lin->getName()) return; _cur_fmts=MEDFileFieldMultiTS::New(); } + void endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field) { if(_cur_fmts.isNotNull()) { if(_cur_fmts->getNumberOfTS()>0) _out_fs->pushField(_cur_fmts); } } + // + void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts); + void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts); + // + void newMeshEntry(const MEDFileFieldPerMesh *fpm); + void endMeshEntry(const MEDFileFieldPerMesh *fpm) { } + // + void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt); + void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt) { } + // + void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd); +private: + void updateData(MEDFileFieldPerMeshPerTypePerDisc *pmtd, mcIdType deltaNbNodes); +private: + const MEDFileUMesh *_lin; + const MEDFileUMesh *_quad; + const MEDFileFieldGlobsReal *_lin_globs; + MEDFileFields *_out_fs; + MCAuto _cur_fmts; + MCAuto _cur_f1ts; + INTERP_KERNEL::NormalizedCellType _gt; + // Info on 1TS modification + bool _1ts_update_requested; + // Cache of matrix to compute faster the values on newly created points + std::map< std::string, PFLData > _cache; + std::vector _pfls_to_be_updated; +}; + +void MEDFileFieldLin2QuadVisitor::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd) +{ + if(_cur_f1ts.isNull()) + return; + if(pmptpd->getType()!=ON_NODES) + throw INTERP_KERNEL::Exception("Not managed yet for ON_CELLS ON_GAUSS_NE and ON_GAUSS_PT"); + _1ts_update_requested=true; + MEDFileAnyTypeField1TSWithoutSDA *ct(_cur_f1ts->contentNotNullBase()); + int locId(pmptpd->getFather()->locIdOfLeaf(pmptpd)); + MEDFileFieldPerMeshPerTypePerDisc *pmtdToModify(ct->getLeafGivenMeshAndTypeAndLocId(_lin->getName(),_gt,locId)); + std::string pflName(pmptpd->getProfile()); + _pfls_to_be_updated.push_back(pflName); + std::map< std::string, PFLData >::iterator itCache(_cache.find(pflName)); + if(itCache!=_cache.end()) + { + updateData(pmtdToModify,(*itCache).second.getNbOfAddPtsInPfl()); + return ; + } + MCAuto pfl; + if(pflName.empty()) + pfl=DataArrayIdType::Range(0,pmptpd->getNumberOfVals(),1); + else + pfl=_lin_globs->getProfile(pflName)->deepCopy(); + // + MCAuto mesh3D(_lin->getMeshAtLevel(0)),mesh3DQuadratic(_quad->getMeshAtLevel(0)); + MCAuto cellIds(mesh3D->getCellIdsLyingOnNodes(pfl->begin(),pfl->end(),true)); + MCAuto mesh3DQuadraticRestricted(mesh3DQuadratic->buildPartOfMySelf(cellIds->begin(),cellIds->end(),true)); + MCAuto mesh3DQuadraticRestrictedNodeIds(mesh3DQuadraticRestricted->computeFetchedNodeIds()); + mesh3DQuadraticRestrictedNodeIds->checkMonotonic(true); + MCAuto newPtsIds(mesh3DQuadraticRestrictedNodeIds->buildSubstraction(pfl)); + MCAuto allSeg3; + { + MCAuto a,b,c,d; + MCAuto seg3Tmp(mesh3DQuadraticRestricted->explodeIntoEdges(a,b,c,d)); + allSeg3=MEDCoupling1SGTUMesh::New(seg3Tmp); + } + if(allSeg3->getCellModelEnum()!=INTERP_KERNEL::NORM_SEG3) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypePerDisc : invalid situation where SEG3 expected !"); + MCAuto midPts,cellSeg3Ids,matrix; + { + DataArrayIdType *nodeConn(allSeg3->getNodalConnectivity()); + nodeConn->rearrange(3); + { + std::vector v(1,2); + midPts=nodeConn->keepSelectedComponents(v); + } + cellSeg3Ids=DataArrayIdType::FindPermutationFromFirstToSecond(midPts,newPtsIds); + { + std::vector v(2); v[0]=0; v[1]=1; + MCAuto tmp(nodeConn->keepSelectedComponents(v)); + matrix=tmp->selectByTupleId(cellSeg3Ids->begin(),cellSeg3Ids->end()); + } + nodeConn->rearrange(1); + } + MCAuto pflq; + if(!pflName.empty()) + { + std::vector vs(2); + vs[0]=pfl; vs[1]=newPtsIds; + pflq=DataArrayIdType::Aggregate(vs); + pflq->setName(pflName); + } + PFLData pdata(matrix,pflq,newPtsIds->getNumberOfTuples()); + _cache[pflName]=pdata; + updateData(pmtdToModify,pdata.getNbOfAddPtsInPfl()); +} + +void MEDFileFieldLin2QuadVisitor::updateData(MEDFileFieldPerMeshPerTypePerDisc *pmtd, mcIdType deltaNbNodes) +{ + pmtd->incrementNbOfVals(deltaNbNodes); +} + +void MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt) +{ + const MEDFileFieldPerMeshPerType *pmpt2(dynamic_cast(pmpt)); + if(!pmpt2) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry : not managed for structure elements !"); + if(pmpt2->getNumberOfLoc()!=1) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newPerMeshPerTypeEntry : not managed for multi discr per timestep !"); + _gt=pmpt->getGeoType(); +} + +void MEDFileFieldLin2QuadVisitor::newMeshEntry(const MEDFileFieldPerMesh *fpm) +{ + if(fpm->getMeshName()!=_lin->getName()) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::newMeshEntry : mismatch into meshName !"); +} + +void MEDFileFieldLin2QuadVisitor::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts) +{ + _1ts_update_requested=false; _pfls_to_be_updated.clear(); + if(!ts) + return ; + const MEDFileField1TSWithoutSDA *tsd(dynamic_cast(ts)); + if(!tsd) + return ; + MCAuto contentCpy(ts->deepCopy()); + MCAuto contentCpy2(DynamicCastSafe(contentCpy)); + if(contentCpy2.isNull()) + return; + _cur_f1ts=MEDFileField1TS::New(*contentCpy2,true); + _cur_f1ts->deepCpyGlobs(*_lin_globs); +} + +void MEDFileFieldLin2QuadVisitor::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts) +{ + if(_cur_f1ts.isNull()) + return ; + if(_1ts_update_requested) + { + MCAuto matrix,oldPfl; + for(std::vector::const_iterator it=_pfls_to_be_updated.begin();it!=_pfls_to_be_updated.end();it++) + { + std::map< std::string, PFLData >::const_iterator it2(_cache.find(*it)); + if(it2==_cache.end()) + throw INTERP_KERNEL::Exception("MEDFileFieldLin2QuadVisitor::endTimeStepEntry : invalid situation !"); + matrix=(*it2).second.getMatrix(); + if((*it).empty()) + continue; + int locId(_cur_f1ts->getProfileId(*it)); + oldPfl.takeRef(_cur_f1ts->getProfile(*it)); + { + std::vector locToKill(1,locId); + _cur_f1ts->killProfileIds(locToKill); + } + _cur_f1ts->appendProfile((*it2).second.getProfile()); + } + DataArrayDouble *arr(_cur_f1ts->getUndergroundDataArray()); + MCAuto res; + { + std::vector v(1,0),v2(1,1); + MCAuto pts0(matrix->keepSelectedComponents(v)); + MCAuto pts1(matrix->keepSelectedComponents(v2)); + if(oldPfl.isNotNull()) + { + pts0=oldPfl->findIdForEach(pts0->begin(),pts0->end()); + pts1=oldPfl->findIdForEach(pts1->begin(),pts1->end()); + } + MCAuto part0(arr->selectByTupleId(*pts0)); + MCAuto part1(arr->selectByTupleId(*pts1)); + res=DataArrayDouble::Add(part0,part1); + res->applyLin(0.5,0.); + } + res=DataArrayDouble::Aggregate(arr,res); + _cur_f1ts->setArray(res); + } + if(_cur_fmts.isNotNull()) + { _cur_fmts->pushBackTimeStep(_cur_f1ts); } + _1ts_update_requested=false; +} + +/*! + * \a newQuad is expected to be the result of MEDFileUMesh::linearToQuadratic of \a oldLin + */ +MCAuto MEDFileFields::linearToQuadratic(const MEDFileMeshes *oldLin, const MEDFileMeshes *newQuad) const +{ + if(!oldLin || !newQuad) + throw INTERP_KERNEL::Exception("MEDFileFields::linearToQuadratic : input meshes must be non NULL !"); + MCAuto ret(MEDFileFields::New()); + for(int i=0;igetNumberOfMeshes();i++) + { + MEDFileMesh *mm(oldLin->getMeshAtPos(i)); + if(!mm) + continue; + MEDFileUMesh *mmu(dynamic_cast(mm)); + if(!mmu) + continue; + MEDFileMesh *mmq(newQuad->getMeshWithName(mmu->getName())); + MEDFileUMesh *mmqu(dynamic_cast(mmq)); + if(!mmqu) + { + std::ostringstream oss; oss << "MEDFileFields::linearToQuadratic : mismatch of name between input meshes for name \"" << mmu->getName() << "\""; + throw INTERP_KERNEL::Exception(oss.str()); + } + MEDFileFieldLin2QuadVisitor vis(mmu,mmqu,this,ret); + accept(vis); + } + return ret; +} + MEDFileAnyTypeFieldMultiTS *MEDFileFields::getFieldAtPos(int i) const { if(i<0 || i>=(int)_fields.size()) @@ -619,17 +854,20 @@ MEDFileAnyTypeFieldMultiTS *MEDFileFields::getFieldAtPos(int i) const return 0; MCAuto ret; const MEDFileFieldMultiTSWithoutSDA *fmtsC(dynamic_cast(fmts)); - const MEDFileIntFieldMultiTSWithoutSDA *fmtsC2(dynamic_cast(fmts)); + const MEDFileInt32FieldMultiTSWithoutSDA *fmtsC2(dynamic_cast(fmts)); + const MEDFileInt64FieldMultiTSWithoutSDA *fmtsC4(dynamic_cast(fmts)); const MEDFileFloatFieldMultiTSWithoutSDA *fmtsC3(dynamic_cast(fmts)); if(fmtsC) ret=MEDFileFieldMultiTS::New(*fmtsC,false); else if(fmtsC2) - ret=MEDFileIntFieldMultiTS::New(*fmtsC2,false); + ret=MEDFileInt32FieldMultiTS::New(*fmtsC2,false); + else if(fmtsC4) + ret=MEDFileInt64FieldMultiTS::New(*fmtsC4,false); else if(fmtsC3) ret=MEDFileFloatFieldMultiTS::New(*fmtsC3,false); else { - std::ostringstream oss; oss << "MEDFileFields::getFieldAtPos : At pos #" << i << " field is neither double (FLOAT64) nor float (FLOAT32) nor integer (INT32) !"; + std::ostringstream oss; oss << "MEDFileFields::getFieldAtPos : At pos #" << i << " field is neither double (FLOAT64) nor float (FLOAT32) nor integer (INT32) nor integer (INT64) !"; throw INTERP_KERNEL::Exception(oss.str()); } ret->shallowCpyGlobs(*this); @@ -839,6 +1077,125 @@ void MEDFileFields::getMeshSENames(std::vector< std::pairaggregateFieldsOnSameMeshes(ms); +} + +/*! + * This method is dedicated to explosed Structured Elements that can lead to exotic situation. + * Especially when there are several structured elements for a same field. + * + * This method looks into meshes into \a ms if there is presence of multiple mesh having same name. + * If so, these meshes are aggregated in the same order than \a ms. + * The fields in \a this lying on the same meshName are also aggregated in the same order than \a this. + */ +void MEDFileFields::aggregateFieldsOnSameMeshes(MEDFileMeshes *ms) +{ + if(!ms) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : ms is nullptr !"); + // + std::vector msNames(ms->getMeshesNames()); + std::set msNamesSet(msNames.begin(),msNames.end()); + if(msNames.size() == msNamesSet.size()) + return ; + // + std::map> > fsByName; + for(auto fmts : _fields) + { + fsByName[fmts->getMeshName()].push_back(fmts); + } + std::vector fieldsNamesToBeAggregated; + std::vector< MCAuto > otherFields; + std::set expectedMeshNamesToMerge; + for(auto fieldsWithSame : fsByName) + { + if(fieldsWithSame.second.size() > 1) + { + fieldsNamesToBeAggregated.push_back(fieldsWithSame.first); + std::set< std::string > zeMeshNames; + for(auto fmtsWithSameName : fieldsWithSame.second) + zeMeshNames.insert(fmtsWithSameName->getMeshName()); + if(zeMeshNames.size()!=1) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : Presence of multiple MultiTS instances with same name but lying on same meshName. Looks bad !"); + std::string meshNameToMerge = *zeMeshNames.begin(); + if(expectedMeshNamesToMerge.find(meshNameToMerge) != expectedMeshNamesToMerge.end()) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : unexpected situation ! Error in implementation !"); + expectedMeshNamesToMerge.insert(*zeMeshNames.begin()); + } + else + { + otherFields.push_back(fieldsWithSame.second.front()); + } + } + for(auto fieldNameToBeAggregated : fieldsNamesToBeAggregated) + { + auto base_fs = fsByName[fieldNameToBeAggregated].front(); + auto fieldsToBeAggregated = fsByName[fieldNameToBeAggregated]; + std::vector< std::vector< std::pair > > dtsToAggregate; + std::vector< MCAuto > eltsToAggregate; + for(auto fieldToBeAggregated : fieldsToBeAggregated) + { + std::vector< std::pair,std::pair > > entries; + int iteration,order; + { + auto dtits = fieldToBeAggregated->getIterations(); + iteration = dtits.front().first; + order = dtits.front().second; + } + fieldToBeAggregated->getUndergroundDataArrayExt(iteration,order,entries); + std::vector< std::pair > dtsToPush; + for(auto entry : entries) + dtsToPush.push_back({entry.first.first,entry.second.second-entry.second.first}); + dtsToAggregate.push_back(dtsToPush); + MCAuto eltToAggregate = MEDFileAnyTypeFieldMultiTS::BuildNewInstanceFromContent(fieldToBeAggregated); + eltsToAggregate.push_back(eltToAggregate); + } + MCAuto gg = MEDFileAnyTypeFieldMultiTS::Aggregate(FromVecAutoToVecOfConst(eltsToAggregate),dtsToAggregate); + gg->setMeshName(base_fs->getMeshName()); + otherFields.push_back(gg->getContent()); + } + // now deal with meshes + std::map > msByName; + for(auto iMesh = 0 ; iMesh < ms->getNumberOfMeshes() ; ++iMesh) + { + auto curMesh = ms->getMeshAtPos(iMesh); + msByName[curMesh->getName()].push_back(curMesh); + } + std::set meshesNamesToBeAggregated; + std::vector< MCAuto > otherMeshes; + for(auto msWithSameName : msByName) + { + if(msWithSameName.second.size()>1) + meshesNamesToBeAggregated.insert(msWithSameName.first); + else + { + otherMeshes.push_back( MCAuto::TakeRef(msWithSameName.second.front()) ); + } + } + if(meshesNamesToBeAggregated != expectedMeshNamesToMerge) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : mismatch between meshes to be aggregated and meshnames into fields to be aggregated"); + std::vector coos; + for(auto meshNameToBeAggregated : meshesNamesToBeAggregated) + { + for(auto curMesh : msByName[meshNameToBeAggregated]) + { + if(!curMesh->getNonEmptyLevels().empty()) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : only meshes without cells supported."); + MEDFileUMesh *curMeshU(dynamic_cast(curMesh)); + if(!curMeshU) + THROW_IK_EXCEPTION("MEDFileFields::aggregateFieldsOnSameMeshes : only unstructured mesh supported."); + coos.push_back(curMeshU->getCoords()); + } + MCAuto coo=DataArrayDouble::Aggregate(coos); + MCAuto gg = MEDFileUMesh::New(); + gg->setName(meshNameToBeAggregated); + gg->setCoords(coo); + otherMeshes.push_back(DynamicCast(gg)); + } + // until this point nothing has changed in \a this nor in \a ms as if a const method. + ms->resize(0); + for(auto mesh : otherMeshes) + ms->pushMesh(mesh); + _fields = otherFields; } MCAuto MEDFileFields::partOfThisOnStructureElements() const @@ -881,7 +1238,7 @@ int MEDFileFields::getPosFromFieldName(const std::string& fieldName) const { std::string tmp(fieldName); std::vector poss; - for(std::size_t i=0;i<_fields.size();i++) + for(unsigned int i=0;i<_fields.size();i++) { const MEDFileAnyTypeFieldMultiTSWithoutSDA *f(_fields[i]); if(f)