1 // FILTER FILTER : implemetation of FILTER idl descriptions
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
24 // File : Filter_Gen_i.cxx
25 // Author : Paul RASCLE, EDF
29 #include "Filter_Gen_i.hxx"
31 #include "Utils_SINGLETON.hxx"
33 #include "Utils_CorbaException.hxx"
34 #include "utilities.h"
36 #include "MEDMEM_EnsightMeshDriver.hxx"
41 #include <TCollection_AsciiString.hxx>
42 #include <TColStd_SequenceOfAsciiString.hxx>
43 #include <HDFascii.hxx>
46 using namespace SALOME_FILTER;
47 Filter_Gen_i* Filter_Gen_i::_FILTERGen = NULL;
49 //=============================================================================
51 * default constructor: not for use
53 //=============================================================================
55 Filter_Gen_i::Filter_Gen_i()
57 MESSAGE("Filter_Gen_i::Filter_Gen_i");
60 //=============================================================================
62 * standard constructor
64 //=============================================================================
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)
73 MESSAGE("activate object");
75 _id = _poa->activate_object(_thisObj);
78 // get an NamingService interface
79 _NS = SINGLETON_<SALOME_NamingService>::Instance() ;
80 ASSERT(SINGLETON_<SALOME_NamingService>::IsAlreadyExisting()) ;
81 _NS->init_orb( _orb ) ;
86 //=============================================================================
88 * default destructor: not for use
90 //=============================================================================
92 Filter_Gen_i::~Filter_Gen_i()
94 MESSAGE("Filter_Gen_i::~Filter_Gen_i");
96 // destruction of gradient field
100 // destruction of support of reference field: reference field is destroyed
101 // by destruction of med object
103 delete _myIField->getSupport();
105 delete _myDField->getSupport();
107 // destruction of criteria: support and field
109 delete _criteria->getSupport();
121 void Filter_Gen_i::loadMED(const char* inMedFile)
125 _med = new ::MED(MED_DRIVER,_file);
128 void Filter_Gen_i::unloadMED()
130 MESSAGE("unloadMED called");
131 // destruction of gradient field
137 // destruction of support of reference field: reference field is destroyed
138 // by destruction of med object
140 delete _myIField->getSupport();
144 delete _myDField->getSupport();
148 // destruction of criteria: support and field
150 delete _criteria->getSupport();
161 StrSeq* Filter_Gen_i::getMeshNames()
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();
171 StrSeq * Filter_Gen_i::getFieldNames()
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();
181 CORBA::Long Filter_Gen_i::getMeshDimension(const char* meshName)
183 return _med->getMesh(meshName)->getMeshDimension();
186 CORBA::Long Filter_Gen_i::getFieldEntity(const char* fieldName,CORBA::Long dt,CORBA::Long it)
188 return _med->getField(fieldName,dt,it)->getSupport()->getEntity();
191 CORBA::Boolean Filter_Gen_i::fieldIsOnAllElements(const char* fieldName,CORBA::Long dt,CORBA::Long it)
193 return _med->getField(fieldName,dt,it)->getSupport()->isOnAllElements();
196 DTITSeq* Filter_Gen_i::getFieldIteration(const char* fieldName)
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;
208 char* Filter_Gen_i::getMeshName(const char* fieldName,CORBA::Long dt,CORBA::Long it)
210 return (char*)(_med->getField(fieldName,dt,it)->getSupport()->getMesh()->getName().c_str());
213 void Filter_Gen_i::readReferenceField(const char* meshName, const char* fieldName, CORBA::Long ts)
215 // read of input mesh
216 _mesh = _med->getMesh(meshName);
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;
228 _myIField = (MEDMEM::FIELD<int>*)field;
234 void Filter_Gen_i::buildGradient() throw(SALOME_FILTER::FILTER_Gen::FilterError)
237 FIELD<double> * gradient;
240 gradient = _myDField->buildGradient();
242 gradient = _myIField->buildGradient();
243 _myGradient = gradient->buildNorm2Field();
246 catch(MEDEXCEPTION& Mex){
247 MESSAGE("SALOME_Exception: Can't calculate gradient");
248 throw SALOME_FILTER::FILTER_Gen::FilterError("Can't calculate gradient");
253 void Filter_Gen_i::getMinMax(CORBA::Double& imin, CORBA::Double& imax,ref_func rf)
260 _myDField->getMinMax(min,max);
263 _myIField->getMinMax(xmin,xmax);
269 _myGradient->getMinMax(min,max);
276 LongSeq* Filter_Gen_i::getHistogram(CORBA::Long size,ref_func rf)
284 myh = _myDField->getHistogram(mysize);
286 myh = _myIField->getHistogram(mysize);
289 myh = _myGradient->getHistogram(mysize);
293 LongSeq *seq = new LongSeq();
294 seq->length(myh.size());
295 for(int i=0;i<myh.size();i++)
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)
302 double val, min, max;
304 MED_EN::medEntityMesh typ;
308 typ = _myDField->getSupport()->getEntity();
310 typ = _myIField->getSupport()->getEntity();
312 // create support on nodes
313 SUPPORT *sup = new SUPPORT(_mesh,"Support",MED_NODE);
315 // create integer field on nodes
316 _criteria = new FIELD<int>(sup,1);
318 _criteria->setName("Criteria");
320 // read number of nodes
321 int NumberOf = sup->getNumberOfElements(MED_ALL_ELEMENTS);
323 for (int i=1; i<NumberOf+1; i++){
325 // if reference field is on elements get reference field on nodes
328 throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on cells");
331 throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on faces");
334 throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on edges");
337 // read reference field value
341 val = _myDField->getValueIJ(i,1);
343 val = (double)_myIField->getValueIJ(i,1);
346 val = _myGradient->getValueIJ(i,1);
350 case MED_ALL_ENTITIES:
351 throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on all entities");
355 // set criteria field value
358 if( val >= fthresh ) isGVal = true;
361 if( val <= fthresh ) isGVal = true;
367 if(sthresh < fthresh){
372 if( (val <= min) || (val >= max) ) isGVal = true;
375 if( (val >= min) && (val <= max) ) isGVal = true;
379 _criteria->setValueIJ(i,1,1);
381 _criteria->setValueIJ(i,1,0);
384 catch(SALOME_Exception& ex){
385 throw SALOME_FILTER::FILTER_Gen::FilterError(ex.what());
389 void Filter_Gen_i::createEnsightInputFiles() throw(SALOME_FILTER::FILTER_Gen::FilterError)
392 // call ensight driver MED to generate input ensight mesh ,
393 // input ensight boundary mesh and input criteria ensight field for filtoo
394 MESSAGE("Create ensight mesh");
395 ENSIGHT_MESH_WRONLY_DRIVER myMeshDriver("/tmp/input.case",_mesh);
397 myMeshDriver.write();
398 myMeshDriver.close();
400 MESSAGE("Create ensight field");
401 ENSIGHT_FIELD_WRONLY_DRIVER<int> myFieldDriver("/tmp/input.case",_criteria);
402 myFieldDriver.open();
403 myFieldDriver.write();
404 myFieldDriver.close();
406 catch(SALOME_Exception& ex){
407 throw SALOME_FILTER::FILTER_Gen::FilterError(ex.what());
411 void Filter_Gen_i::filtering() throw(SALOME_FILTER::FILTER_Gen::FilterError)
415 MESSAGE("call filtoo");
416 // send filtoo command
417 command = "cd /tmp;filtoo -s -m 1 -f input -o output > /tmp/filter.log";
419 int status = system(command.c_str());
421 throw SALOME_FILTER::FILTER_Gen::FilterError("filtoo error");
423 // destroy filtoo input files
424 command = "cd /tmp;rm -f input.*";
426 system(command.c_str());
430 void Filter_Gen_i::projectFieldsOnDecimateMesh() throw(SALOME_FILTER::FILTER_Gen::FilterError)
434 // read of new mesh in ensight file
435 MESH* myMesh = new MESH();
436 ENSIGHT_MESH_RDONLY_DRIVER myEnsightMeshDriver("/tmp/output.case", myMesh) ;
437 myMesh->addDriver(ENSIGHT_DRIVER,"/tmp/output.case","myMesh",MED_EN::MED_LECT);
440 // have to call ensight driver MED to generate output MED file from filtoo output
441 _newMed = new ::MED();
442 _newMed->addMesh(myMesh);
445 deque<string> meshesNames = _newMed->getMeshNames();
446 int numberOfMeshes = meshesNames.size();
447 if( numberOfMeshes != 1)
448 throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid number of meshes in filtoo output");
450 // new mesh generated by filtoo
451 _newMesh = _newMed->getMesh(meshesNames[0]);
453 // create support on nodes on all new mesh
454 SUPPORT *newSup = new SUPPORT(_newMesh,"Support",MED_NODE);
456 // read the id of nodes of output mesh, in input mesh
459 // read connectivity of new mesh to get neighbour node of created node
460 _connL = _newMesh->getConnectivityLength(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
461 _conn = _newMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
462 _connI = _newMesh->getConnectivityIndex(MED_NODAL,MED_CELL);
464 // read number of nodes on new mesh
465 int numberOfNodes = newSup->getNumberOfElements(MED_ALL_ELEMENTS);
466 int numberOfComponents;
468 deque<string> fieldsNames = _med->getFieldNames();
469 int numberOfFields = fieldsNames.size();
474 for (int i=0; i<numberOfFields; i++){
476 // is the input field the reference field?
477 bool isReferenceField= false;
479 if( strcmp(_myDField->getName().c_str(),fieldsNames[i].c_str()) == 0)
480 isReferenceField = true;
483 if( strcmp(_myIField->getName().c_str(),fieldsNames[i].c_str()) == 0)
484 isReferenceField = true;
486 deque<DT_IT_> myIteration = _med->getFieldIteration (fieldsNames[i]);
487 string meshName = _med->getField(fieldsNames[i],myIteration[0].dt,myIteration[0].it)->getSupport()->getMesh()->getName();
489 // we process only fields on input mesh
490 if( strcmp(meshName.c_str(),_mesh->getName().c_str()) == 0){
492 // loop on time steps
493 int numberOfIteration = myIteration.size();
494 for(int j=0;j<numberOfIteration;j++){
496 // select input field
497 MEDMEM::FIELD_* field = _med->getField(fieldsNames[i],myIteration[j].dt,myIteration[j].it);
499 // if field not on nodes, take following field
500 if( field->getSupport()->getEntity() != MED_NODE )
503 FIELD<double> *myDField = NULL;
504 FIELD<double> *newDField = NULL;
505 FIELD<int> *myIField = NULL;
506 FIELD<int> *newIField = NULL;
508 if (dynamic_cast<MEDMEM::FIELD<double>*>(field)){
509 if(!isReferenceField){
510 // read input field on input file
511 myDField = new FIELD<double>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
514 myDField = _myDField;
516 // create new output field
517 newDField = new FIELD<double>(newSup,field->getNumberOfComponents());
518 newDField->setName(myDField->getName());
519 newDField->setIterationNumber(myIteration[j].dt);
520 newDField->setOrderNumber(myIteration[j].it);
521 newDField->setTime(myDField->getTime());
524 if(!isReferenceField)
525 // read input field on input file
526 myIField = new FIELD<int>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
528 myIField = _myIField;
529 // create new output field
530 newIField = new FIELD<int>(newSup,field->getNumberOfComponents());
531 newIField->setName(myIField->getName());
532 newIField->setIterationNumber(myIteration[j].dt);
533 newIField->setOrderNumber(myIteration[j].it);
534 newIField->setTime(myIField->getTime());
536 numberOfComponents = field->getNumberOfComponents();
538 // loop on nodes on new field
539 for (int k=1; k<=numberOfNodes; k++){
540 // read number of nodes on input field
541 int l = getNodeNumber(k);
545 for(int c=1;c<=numberOfComponents;c++){
546 // read value on input field
548 dval = myDField->getValueIJ(l,c);
550 ival = myIField->getValueIJ(l,c);
552 // write value on new field
554 newDField->setValueIJ(k,c,dval);
557 newIField->setValueIJ(k,c,ival);
562 _newMed->addField(newDField);
564 _newMed->addField(newIField);
566 // Destroy input field if not reference field
567 if(!isReferenceField)
569 delete myDField->getSupport();
573 delete myIField->getSupport();
581 catch(SALOME_Exception){
582 throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid decimate mlesh created by filtoo");
585 // destroy filtoo output files
586 command = "cd /tmp;rm -f output.*";
588 system(command.c_str());
592 void Filter_Gen_i::readMapping()
596 ifstream mapFile("/tmp/output.renum");
598 // read first line of file
599 getline(mapFile, firstline);
601 // read number of vertices to map
602 string::size_type pos = firstline.find(":");
603 istringstream iss(firstline.substr(pos+1));
606 // read each vertices in array
607 _map = new int[_nbvmap];
609 for(int i=0;i<_nbvmap;i++){
616 int Filter_Gen_i::getNodeNumber(int num)
618 int oldnum = _map[num-1];
620 // if new node get neighbour node
622 oldnum = _map[getNeighbourVertex(num)-1];
627 int Filter_Gen_i::getNeighbourVertex(int num) throw(SALOME_Exception)
630 int numberOfElements = _newMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
633 // get index of created node in connectivity array
634 for(index=0;index<_connL;index++){
635 if(_conn[index] == num){
637 // get index of element which contain created node
639 for(i=1;i<=numberOfElements;i++)
640 if(_connI[i] > index)
643 // search neighbour node which are in old mesh
644 for(int j=_connI[i-1];j<_connI[i];j++){
646 if( _map[nnum-1] != 0)
650 // if neighbour node in old mesh: go out loop, else continue
656 // if no neighbour node in old mesh: throw exception
658 throw SALOME_Exception("None of the neighbour node are in old mesh!!");
663 void Filter_Gen_i::createMedOutputFile(const char *outMedFile)
667 MESSAGE("Create MED mesh: "<<outMedFile);
668 id = _newMed->addDriver(MED_DRIVER,outMedFile,MED_EN::MED_ECRI);
673 //=============================================================================
675 * C factory, accessible with dlsym, after dlopen
677 //=============================================================================
681 PortableServer::ObjectId * FILTEREngine_factory(
683 PortableServer::POA_ptr poa,
684 PortableServer::ObjectId * contId,
685 const char *instanceName,
686 const char *interfaceName)
688 MESSAGE("PortableServer::ObjectId * FilterEngine_factory()");
689 SCRUTE(interfaceName);
690 Filter_Gen_i * myFilter_Gen
691 = new Filter_Gen_i(orb, poa, contId, instanceName, interfaceName);
692 return myFilter_Gen->getId() ;