src/FILTER/Filter_Gen_i.cxx

   1 //  FILTER FILTER : implemetation of FILTER idl descriptions
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //
  23 //
  24 //  File   : Filter_Gen_i.cxx
  25 //  Author : Paul RASCLE, EDF
  26 //  Module : FILTER
  27 //  $Header$
  28
  29 #include "Filter_Gen_i.hxx"
  30
  31 #include "Utils_SINGLETON.hxx"
  32 #include "OpUtil.hxx"
  33 #include "Utils_CorbaException.hxx"
  34 #include "utilities.h"
  35
  36 #include "MEDMEM_EnsightMeshDriver.hxx"
  37 #include <string>
  38 #include <deque>
  39 #include <map>
  40
  41 #include <TCollection_AsciiString.hxx>
  42 #include <TColStd_SequenceOfAsciiString.hxx>
  43 #include <HDFascii.hxx>
  44
  45 using namespace std;
  46 using namespace SALOME_FILTER;
  47 Filter_Gen_i* Filter_Gen_i::_FILTERGen = NULL;
  48
  49 //=============================================================================
  50 /*!
  51  *  default constructor: not for use
  52  */
  53 //=============================================================================
  54
  55 Filter_Gen_i::Filter_Gen_i()
  56 {
  57   MESSAGE("Filter_Gen_i::Filter_Gen_i");
  58 }
  59
  60 //=============================================================================
  61 /*!
  62  *  standard constructor
  63  */
  64 //=============================================================================
  65
  66 Filter_Gen_i:: Filter_Gen_i(CORBA::ORB_ptr orb,
  67                                 PortableServer::POA_ptr poa,
  68                                 PortableServer::ObjectId * contId,
  69                                 const char *instanceName,
  70                                 const char *interfaceName) :
  71   Engines_Component_i(orb, poa, contId, instanceName, interfaceName),_med(NULL),_newMed(NULL),_mesh(NULL),_newMesh(NULL),_myGradient(NULL),_myDField(NULL),_myIField(NULL),_criteria(NULL),_map(NULL)
  72 {
  73   MESSAGE("activate object");
  74   _thisObj = this ;
  75   _id = _poa->activate_object(_thisObj);
  76
  77   _duringLoad=false;
  78   // get an NamingService interface
  79   _NS = SINGLETON_<SALOME_NamingService>::Instance() ;
  80   ASSERT(SINGLETON_<SALOME_NamingService>::IsAlreadyExisting()) ;
  81   _NS->init_orb( _orb ) ;
  82
  83   _FILTERGen = this;
  84 }
  85
  86 //=============================================================================
  87 /*!
  88  *  default destructor: not for use
  89  */
  90 //=============================================================================
  91
  92 Filter_Gen_i::~Filter_Gen_i()
  93 {
  94   MESSAGE("Filter_Gen_i::~Filter_Gen_i");
  95
  96   // destruction of gradient field
  97   if(_myGradient)
  98     delete _myGradient;
  99
 100   // destruction of support of reference field: reference field is destroyed
 101   // by destruction of med object
 102   if(_myIField)
 103     delete _myIField->getSupport();
 104   if(_myDField)
 105     delete _myDField->getSupport();
 106
 107  // destruction of criteria: support and field
 108   if(_criteria){
 109     delete _criteria->getSupport();
 110     delete _criteria;
 111   }
 112
 113   if(_map)
 114     delete [] _map;
 115   if(_med)
 116     delete _med;
 117   if(_newMed)
 118     delete _newMed;
 119 }
 120
 121 void Filter_Gen_i::loadMED(const char* inMedFile)
 122 {
 123   SCRUTE(inMedFile);
 124   _file = inMedFile;
 125   _med = new ::MED(MED_DRIVER,_file);
 126 }
 127
 128 void Filter_Gen_i::unloadMED()
 129 {
 130   MESSAGE("unloadMED called");
 131   // destruction of gradient field
 132   if(_myGradient){
 133     delete _myGradient;
 134     _myGradient=NULL;
 135   }
 136
 137   // destruction of support of reference field: reference field is destroyed
 138   // by destruction of med object
 139   if(_myIField){
 140     delete _myIField->getSupport();
 141     _myIField=NULL;
 142   }
 143   if(_myDField){
 144     delete _myDField->getSupport();
 145     _myDField=NULL;
 146   }
 147
 148  // destruction of criteria: support and field
 149   if(_criteria){
 150     delete _criteria->getSupport();
 151     delete _criteria;
 152     _criteria=NULL;
 153   }
 154
 155   if(_med){
 156     delete _med;
 157     _med=NULL;
 158   }
 159 }
 160
 161 StrSeq* Filter_Gen_i::getMeshNames()
 162 {
 163   StrSeq *seq = new StrSeq();
 164   deque<string> deq = _med->getMeshNames();
 165   seq->length(deq.size());
 166   for(int i=0;i<deq.size();i++)
 167     (*seq)[i] = deq[i].c_str();
 168   return seq;
 169 }
 170
 171 StrSeq * Filter_Gen_i::getFieldNames()
 172 {
 173   StrSeq *seq = new StrSeq();
 174   deque<string> deq = _med->getFieldNames();
 175   seq->length(deq.size());
 176   for(int i=0;i<deq.size();i++)
 177      (*seq)[i] = deq[i].c_str();
 178   return seq;
 179 }
 180
 181 CORBA::Long  Filter_Gen_i::getMeshDimension(const char* meshName)
 182 {
 183   return _med->getMesh(meshName)->getMeshDimension();
 184 }
 185
 186 CORBA::Long  Filter_Gen_i::getFieldEntity(const char* fieldName,CORBA::Long dt,CORBA::Long it)
 187 {
 188   return _med->getField(fieldName,dt,it)->getSupport()->getEntity();
 189 }
 190
 191 CORBA::Boolean Filter_Gen_i::fieldIsOnAllElements(const char* fieldName,CORBA::Long dt,CORBA::Long it)
 192 {
 193   return _med->getField(fieldName,dt,it)->getSupport()->isOnAllElements();
 194 }
 195
 196 DTITSeq* Filter_Gen_i::getFieldIteration(const char* fieldName)
 197 {
 198   DTITSeq *seq = new DTITSeq();
 199   deque<DT_IT_> deq = _med->getFieldIteration(fieldName);
 200   seq->length(deq.size());
 201   for(int i=0;i<deq.size();i++){
 202     (*seq)[i].dt = deq[i].dt;
 203     (*seq)[i].it = deq[i].it;
 204   }
 205   return seq;
 206 }
 207
 208 char* Filter_Gen_i::getMeshName(const char* fieldName,CORBA::Long dt,CORBA::Long it)
 209 {
 210   return (char*)(_med->getField(fieldName,dt,it)->getSupport()->getMesh()->getName().c_str());
 211 }
 212
 213 void Filter_Gen_i::readReferenceField(const char* meshName, const char* fieldName, CORBA::Long ts)
 214 {
 215   // read of input mesh
 216   _mesh = _med->getMesh(meshName);
 217   _mesh->read();
 218
 219   // read of input field
 220   deque<DT_IT_> myIteration = _med->getFieldIteration (fieldName);
 221   MEDMEM::FIELD_* field = _med->getField(fieldName,myIteration[ts].dt,myIteration[ts].it);
 222   if (dynamic_cast<MEDMEM::FIELD<double>*>(field)){
 223     _myDField = (MEDMEM::FIELD<double>*)field;
 224     _myDField->read();
 225     _myIField = NULL;
 226   }
 227   else{
 228     _myIField = (MEDMEM::FIELD<int>*)field;
 229     _myIField->read();
 230     _myDField = NULL;
 231   }
 232 }
 233
 234 void Filter_Gen_i::buildGradient() throw(SALOME_FILTER::FILTER_Gen::FilterError)
 235 {
 236   if(!_myGradient){
 237     FIELD<double> * gradient;
 238     try{
 239       if(_myDField)
 240         gradient = _myDField->buildGradient();
 241       else
 242         gradient = _myIField->buildGradient();
 243       _myGradient = gradient->buildNorm2Field();
 244       delete gradient;
 245     }
 246     catch(MEDEXCEPTION& Mex){
 247       MESSAGE("SALOME_Exception: Can't calculate gradient");
 248       throw SALOME_FILTER::FILTER_Gen::FilterError("Can't calculate gradient");
 249     }
 250   }
 251 }
 252
 253 void Filter_Gen_i::getMinMax(CORBA::Double& imin, CORBA::Double& imax,ref_func rf)
 254 {
 255   double min, max;
 256
 257   switch(rf){
 258   case F_FIELD:
 259     if (_myDField)
 260       _myDField->getMinMax(min,max);
 261     else{
 262       int xmin, xmax;
 263       _myIField->getMinMax(xmin,xmax);
 264       min = (double)xmin;
 265       max = (double)xmax;
 266     }
 267     break;
 268   case F_GRAD:
 269     _myGradient->getMinMax(min,max);
 270     break;
 271   }
 272   imin = min;
 273   imax = max;
 274 }
 275
 276 LongSeq* Filter_Gen_i::getHistogram(CORBA::Long size,ref_func rf)
 277 {
 278   int mysize = size;
 279   vector<int> myh;
 280
 281   switch(rf){
 282   case F_FIELD:
 283     if (_myDField)
 284       myh = _myDField->getHistogram(mysize);
 285     else
 286       myh = _myIField->getHistogram(mysize);
 287     break;
 288   case F_GRAD:
 289     myh = _myGradient->getHistogram(mysize);
 290     break;
 291   }
 292
 293   LongSeq *seq = new LongSeq();
 294   seq->length(myh.size());
 295   for(int i=0;i<myh.size();i++)
 296     (*seq)[i] = myh[i];
 297   return seq;
 298 }
 299
 300 void Filter_Gen_i::generateCriteria(CORBA::Long nbthresh,CORBA::Double fthresh,CORBA::Double sthresh,CORBA::Boolean areaFlag,ref_func rf) throw(SALOME_FILTER::FILTER_Gen::FilterError)
 301 {
 302   double val, min, max;
 303   bool isGVal;
 304   MED_EN::medEntityMesh typ;
 305
 306   if(_myDField)
 307     typ = _myDField->getSupport()->getEntity();
 308   else
 309     typ = _myIField->getSupport()->getEntity();
 310
 311   // create support on nodes
 312   SUPPORT *sup = new SUPPORT(_mesh,"Support",MED_NODE);
 313
 314   // create integer field on nodes
 315   _criteria = new FIELD<int>(sup,1);
 316
 317   _criteria->setName("Criteria");
 318
 319   // read number of nodes
 320   int NumberOf = sup->getNumberOfElements(MED_ALL_ELEMENTS);
 321
 322   for (int i=1; i<NumberOf+1; i++){
 323
 324     // if reference field is on elements get reference field on nodes
 325     switch(typ){
 326     case MED_CELL:
 327       throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on cells");
 328       break;
 329     case MED_FACE:
 330       throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on faces");
 331       break;
 332     case MED_EDGE:
 333       throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on edges");
 334       break;
 335     case MED_NODE:
 336       // read reference field value
 337       switch(rf){
 338       case F_FIELD:
 339         if(_myDField)
 340           val = _myDField->getValueIJ(i,1);
 341         else
 342           val = (double)_myIField->getValueIJ(i,1);
 343         break;
 344       case F_GRAD:
 345         val = _myGradient->getValueIJ(i,1);
 346         break;
 347       }
 348       break;
 349     case MED_ALL_ENTITIES:
 350       throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on all entities");
 351       break;
 352     }
 353
 354     // set criteria field value
 355     if( nbthresh == 1 ){
 356       if( areaFlag )
 357         if( val >= fthresh ) isGVal = true;
 358         else isGVal = false;
 359       else
 360         if( val <= fthresh ) isGVal = true;
 361         else isGVal = false;
 362     }
 363     else{
 364       min = fthresh;
 365       max = sthresh;
 366       if(sthresh < fthresh){
 367         min = sthresh;
 368         max = fthresh;
 369       }
 370       if( areaFlag )
 371         if( (val <= min) || (val >= max) ) isGVal = true;
 372         else isGVal = false;
 373       else
 374         if( (val >= min) && (val <= max) ) isGVal = true;
 375         else isGVal = false;
 376     }
 377     if( isGVal )
 378       _criteria->setValueIJ(i,1,1);
 379     else
 380       _criteria->setValueIJ(i,1,0);
 381   }
 382 }
 383
 384 void Filter_Gen_i::createEnsightInputFiles()
 385 {
 386   MESSAGE("Calculate boundary");
 387   // generate MED boundary of geometry
 388   SUPPORT *supB = _mesh->getBoundaryElements(MED_FACE);
 389
 390   // call ensight driver MED to generate input ensight mesh ,
 391   // input ensight boundary mesh and input criteria ensight field for filtoo
 392   MESSAGE("Create ensight mesh");
 393   ENSIGHT_MESH_WRONLY_DRIVER myMeshDriver("/tmp/input.case",_mesh);
 394   myMeshDriver.addSupport(supB);
 395   myMeshDriver.open();
 396   myMeshDriver.write();
 397   myMeshDriver.close();
 398
 399   MESSAGE("Create ensight field");
 400   ENSIGHT_FIELD_WRONLY_DRIVER<int> myFieldDriver("/tmp/input.case",_criteria);
 401   myFieldDriver.open();
 402   myFieldDriver.write();
 403   myFieldDriver.close();
 404 }
 405
 406 void Filter_Gen_i::filtering()
 407 {
 408   string command;
 409
 410   MESSAGE("call filtoo");
 411   // send filtoo command
 412   command = "cd /tmp;filtoo -f input -o output > /tmp/filter.log";
 413   MESSAGE(command);
 414   system(command.c_str());
 415
 416   // destroy filtoo input files
 417   command = "cd /tmp;rm -f input.*";
 418   MESSAGE(command);
 419   system(command.c_str());
 420
 421 }
 422
 423 void Filter_Gen_i::projectFieldsOnDecimateMesh() throw(SALOME_FILTER::FILTER_Gen::FilterError)
 424 {
 425   string command;
 426
 427   // read of new mesh in ensight file
 428   MESH* myMesh = new MESH();
 429   ENSIGHT_MESH_RDONLY_DRIVER myEnsightMeshDriver("/tmp/output.case", myMesh) ;
 430   myMesh->addDriver(ENSIGHT_DRIVER,"/tmp/output.case","myMesh",MED_EN::MED_LECT);
 431   myMesh->read() ;
 432
 433   // have to call ensight driver MED to generate output MED file from filtoo output
 434   _newMed = new ::MED();
 435   _newMed->addMesh(myMesh);
 436
 437   // get new mesh name
 438   deque<string> meshesNames = _newMed->getMeshNames();
 439   int numberOfMeshes = meshesNames.size();
 440   if( numberOfMeshes != 1)
 441     throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid number of meshes in filtoo output");
 442
 443   // new mesh generated by filtoo
 444   _newMesh = _newMed->getMesh(meshesNames[0]);
 445
 446   // create support on nodes on all new mesh
 447   SUPPORT *newSup = new SUPPORT(_newMesh,"Support",MED_NODE);
 448
 449   // read the id of nodes of output mesh, in input mesh
 450   readMapping();
 451
 452   // read connectivity of new mesh to get neighbour node of created node
 453   _connL = _newMesh->getConnectivityLength(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
 454   _conn = _newMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
 455   _connI = _newMesh->getConnectivityIndex(MED_NODAL,MED_CELL);
 456
 457   // read number of nodes on new mesh
 458   int numberOfNodes = newSup->getNumberOfElements(MED_ALL_ELEMENTS);
 459   int numberOfComponents;
 460
 461   deque<string> fieldsNames = _med->getFieldNames();
 462   int numberOfFields = fieldsNames.size();
 463
 464   try{
 465
 466     // loop on fields
 467     for (int i=0; i<numberOfFields; i++){
 468
 469       // is the input field the reference field?
 470       bool isReferenceField= false;
 471       if(_myDField){
 472         if( strcmp(_myDField->getName().c_str(),fieldsNames[i].c_str()) == 0)
 473           isReferenceField = true;
 474       }
 475       else
 476         if( strcmp(_myIField->getName().c_str(),fieldsNames[i].c_str()) == 0)
 477           isReferenceField = true;
 478
 479       deque<DT_IT_> myIteration = _med->getFieldIteration (fieldsNames[i]);
 480       string meshName = _med->getField(fieldsNames[i],myIteration[0].dt,myIteration[0].it)->getSupport()->getMesh()->getName();
 481
 482       // we process only fields on input mesh
 483       if( strcmp(meshName.c_str(),_mesh->getName().c_str()) == 0){
 484
 485         // loop on time steps
 486         int numberOfIteration = myIteration.size();
 487         for(int j=0;j<numberOfIteration;j++){
 488
 489           // select input field
 490           MEDMEM::FIELD_* field = _med->getField(fieldsNames[i],myIteration[j].dt,myIteration[j].it);
 491
 492           // if field not on nodes, take following field
 493           if( field->getSupport()->getEntity() != MED_NODE )
 494             break;
 495
 496           FIELD<double> *myDField = NULL;
 497           FIELD<double> *newDField = NULL;
 498           FIELD<int> *myIField = NULL;
 499           FIELD<int> *newIField = NULL;
 500
 501           if (dynamic_cast<MEDMEM::FIELD<double>*>(field)){
 502             if(!isReferenceField){
 503               // read input field on input file
 504               myDField = new FIELD<double>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
 505             }
 506             else{
 507               myDField = _myDField;
 508             }
 509             // create new output field
 510             newDField = new FIELD<double>(newSup,field->getNumberOfComponents());
 511             newDField->setName(myDField->getName());
 512           }
 513           else{
 514             if(!isReferenceField)
 515               // read input field on input file
 516               myIField = new FIELD<int>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
 517             else
 518               myIField = _myIField;
 519             // create new output field
 520             newIField = new FIELD<int>(newSup,field->getNumberOfComponents());
 521             newIField->setName(myIField->getName());
 522           }
 523           numberOfComponents = field->getNumberOfComponents();
 524
 525           // loop on nodes on new field
 526           for (int k=1; k<=numberOfNodes; k++){
 527             // read number of nodes on input field
 528             int l = getNodeNumber(k);
 529             double dval;
 530             int ival;
 531
 532             for(int c=1;c<=numberOfComponents;c++){
 533               // read value on input field
 534               if(myDField)
 535                 dval = myDField->getValueIJ(l,c);
 536               else
 537                 ival = myIField->getValueIJ(l,c);
 538
 539               // write value on new field
 540               if(newDField){
 541                 newDField->setValueIJ(k,c,dval);
 542               }
 543               else
 544                 newIField->setValueIJ(k,c,ival);
 545             }
 546
 547           }
 548           if(newDField)
 549             _newMed->addField(newDField);
 550           else
 551             _newMed->addField(newIField);
 552
 553           // Destroy input field if not reference field
 554           if(!isReferenceField)
 555             if(myDField){
 556               delete myDField->getSupport();
 557               delete myDField;
 558             }
 559             else{
 560               delete myIField->getSupport();
 561               delete myIField;
 562             }
 563
 564         }
 565       }
 566     }
 567   }
 568   catch(SALOME_Exception){
 569     throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid decimate mlesh created by filtoo");
 570   }
 571
 572   // destroy filtoo output files
 573   command = "cd /tmp;rm -f output.*";
 574   MESSAGE(command);
 575   system(command.c_str());
 576
 577 }
 578
 579 void Filter_Gen_i::readMapping()
 580 {
 581   ifstream mapFile("/tmp/output.renum");
 582
 583   // read number of vertices to map
 584   mapFile >> _nbvmap;
 585   _map = new int[_nbvmap];
 586
 587   for(int i=0;i<_nbvmap;i++)
 588     mapFile >> _map[i];
 589
 590 }
 591
 592 int Filter_Gen_i::getNodeNumber(int num)
 593 {
 594   int oldnum = _map[num-1];
 595
 596   // if new node get neighbour node
 597   if(oldnum == 0)
 598     oldnum = _map[getNeighbourVertex(num)-1];
 599
 600   return oldnum;
 601 }
 602
 603 int Filter_Gen_i::getNeighbourVertex(int num) throw(SALOME_Exception)
 604 {
 605   int nnum;
 606   int numberOfElements = _newMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
 607
 608   int index;
 609   // get index of created node in connectivity array
 610   for(index=0;index<_connL;index++){
 611     if(_conn[index] == num){
 612
 613       // get index of element which contain created node
 614       int i;
 615       for(i=1;i<=numberOfElements;i++)
 616         if(_connI[i] > index)
 617           break;
 618
 619       // search neighbour node which are in old mesh
 620       for(int j=_connI[i-1];j<_connI[i];j++){
 621         nnum = _conn[j-1];
 622         if( _map[nnum-1] != 0)
 623           break;
 624       }
 625
 626       // if neighbour node in old mesh: go out loop, else continue
 627       if(_map[nnum-1]!=0)
 628         break;
 629     }
 630   }
 631
 632   // if no neighbour node in old mesh: throw exception
 633   if(_map[nnum-1]==0)
 634     throw SALOME_Exception("None of the neighbour node are in old mesh!!");
 635
 636   return nnum;
 637 }
 638
 639 void Filter_Gen_i::createMedOutputFile(const char *outMedFile)
 640 {
 641   int id;
 642
 643   MESSAGE("Create MED mesh: "<<outMedFile);
 644   id = _newMed->addDriver(MED_DRIVER,outMedFile,MED_EN::MED_ECRI);
 645   _newMed->write(id);
 646 }
 647
 648
 649 //=============================================================================
 650 /*!
 651  * C factory, accessible with dlsym, after dlopen
 652  */
 653 //=============================================================================
 654
 655 extern "C"
 656 {
 657   PortableServer::ObjectId * FILTEREngine_factory(
 658                                CORBA::ORB_ptr orb,
 659                                PortableServer::POA_ptr poa,
 660                                PortableServer::ObjectId * contId,
 661                                const char *instanceName,
 662                                const char *interfaceName)
 663   {
 664     MESSAGE("PortableServer::ObjectId * FilterEngine_factory()");
 665     SCRUTE(interfaceName);
 666     Filter_Gen_i * myFilter_Gen
 667       = new Filter_Gen_i(orb, poa, contId, instanceName, interfaceName);
 668     return myFilter_Gen->getId() ;
 669   }
 670 }