#include "Utils_CorbaException.hxx"
#include "utilities.h"
+#include "MEDMEM_EnsightMeshDriver.hxx"
#include <string>
#include <deque>
#include <map>
#include <TCollection_AsciiString.hxx>
#include <TColStd_SequenceOfAsciiString.hxx>
#include <HDFascii.hxx>
-#include "SALOMEDS_Tool.hxx"
using namespace std;
- Filter_Gen_i* Filter_Gen_i::_FILTERGen = NULL;
+using namespace SALOME_FILTER;
+Filter_Gen_i* Filter_Gen_i::_FILTERGen = NULL;
//=============================================================================
/*!
PortableServer::ObjectId * contId,
const char *instanceName,
const char *interfaceName) :
- Engines_Component_i(orb, poa, contId, instanceName, interfaceName)
+ Engines_Component_i(orb, poa, contId, instanceName, interfaceName),_med(NULL),_newMed(NULL),_mesh(NULL),_newMesh(NULL),_myGradient(NULL),_myDField(NULL),_myIField(NULL),_criteria(NULL)
{
MESSAGE("activate object");
_thisObj = this ;
ASSERT(SINGLETON_<SALOME_NamingService>::IsAlreadyExisting()) ;
_NS->init_orb( _orb ) ;
- //_myFilterI = 0;
_FILTERGen = this;
}
Filter_Gen_i::~Filter_Gen_i()
{
MESSAGE("Filter_Gen_i::~Filter_Gen_i");
+
+ // destruction of gradient field
+ if(_myGradient)
+ delete _myGradient;
+
+ // destruction of support of reference field: reference field is destroyed
+ // by destruction of med object
+ if(_myIField)
+ delete _myIField->getSupport();
+ if(_myDField)
+ delete _myDField->getSupport();
+
+ // destruction of criteria: support and field
+ if(_criteria){
+ delete _criteria->getSupport();
+ delete _criteria;
+ }
+
+ if(_med)
+ delete _med;
+}
+
+void Filter_Gen_i::loadMED(const char* inMedFile)
+{
+ SCRUTE(inMedFile);
+ _file = inMedFile;
+ _med = new ::MED(MED_DRIVER,_file);
+}
+
+void Filter_Gen_i::unloadMED()
+{
+ MESSAGE("unloadMED called");
+ // destruction of gradient field
+ if(_myGradient){
+ delete _myGradient;
+ _myGradient=NULL;
+ }
+
+ // destruction of support of reference field: reference field is destroyed
+ // by destruction of med object
+ if(_myIField){
+ delete _myIField->getSupport();
+ _myIField=NULL;
+ }
+ if(_myDField){
+ delete _myDField->getSupport();
+ _myDField=NULL;
+ }
+
+ // destruction of criteria: support and field
+ if(_criteria){
+ delete _criteria->getSupport();
+ delete _criteria;
+ _criteria=NULL;
+ }
+
+ if(_med){
+ delete _med;
+ _med=NULL;
+ }
+}
+
+StrSeq* Filter_Gen_i::getMeshNames()
+{
+ StrSeq *seq = new StrSeq();
+ deque<string> deq = _med->getMeshNames();
+ seq->length(deq.size());
+ for(int i=0;i<deq.size();i++)
+ (*seq)[i] = deq[i].c_str();
+ return seq;
+}
+
+StrSeq * Filter_Gen_i::getFieldNames()
+{
+ StrSeq *seq = new StrSeq();
+ deque<string> deq = _med->getFieldNames();
+ seq->length(deq.size());
+ for(int i=0;i<deq.size();i++)
+ (*seq)[i] = deq[i].c_str();
+ return seq;
+}
+
+CORBA::Long Filter_Gen_i::getMeshDimension(const char* meshName)
+{
+ return _med->getMesh(meshName)->getMeshDimension();
+}
+
+DTITSeq* Filter_Gen_i::getFieldIteration(const char* fieldName)
+{
+ DTITSeq *seq = new DTITSeq();
+ deque<DT_IT_> deq = _med->getFieldIteration(fieldName);
+ seq->length(deq.size());
+ for(int i=0;i<deq.size();i++){
+ (*seq)[i].dt = deq[i].dt;
+ (*seq)[i].it = deq[i].it;
+ }
+ return seq;
+}
+
+char* Filter_Gen_i::getMeshName(const char* fieldName,CORBA::Long dt,CORBA::Long it)
+{
+ return (char*)(_med->getField(fieldName,dt,it)->getSupport()->getMesh()->getName().c_str());
+}
+
+void Filter_Gen_i::readReferenceField(const char* meshName, const char* fieldName, CORBA::Long ts)
+{
+ // read of input mesh
+ _mesh = _med->getMesh(meshName);
+ _mesh->read();
+
+ // read of input field
+ deque<DT_IT_> myIteration = _med->getFieldIteration (fieldName);
+ MEDMEM::FIELD_* field = _med->getField(fieldName,myIteration[ts].dt,myIteration[ts].it);
+ if (dynamic_cast<MEDMEM::FIELD<double>*>(field)){
+ _myDField = (MEDMEM::FIELD<double>*)field;
+ _myDField->read();
+ _myIField = NULL;
+ }
+ else{
+ _myIField = (MEDMEM::FIELD<int>*)field;
+ _myIField->read();
+ _myDField = NULL;
+ }
+}
+
+void Filter_Gen_i::buildGradient() throw(SALOME_FILTER::FILTER_Gen::FilterError)
+{
+ if(!_myGradient){
+ FIELD<double> * gradient;
+ try{
+ if(_myDField)
+ gradient = _myDField->buildGradient();
+ else
+ gradient = _myIField->buildGradient();
+ _myGradient = gradient->buildNorm2Field();
+ delete gradient;
+ }
+ catch(MEDEXCEPTION& Mex){
+ MESSAGE("SALOME_Exception: Can't calculate gradient");
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Can't calculate gradient");
+ }
+ }
+}
+
+void Filter_Gen_i::getMinMax(CORBA::Double& imin, CORBA::Double& imax,ref_func rf)
+{
+ double min, max;
+
+ switch(rf){
+ case F_FIELD:
+ if (_myDField)
+ _myDField->getMinMax(min,max);
+ else{
+ int xmin, xmax;
+ _myIField->getMinMax(xmin,xmax);
+ min = (double)xmin;
+ max = (double)xmax;
+ }
+ break;
+ case F_GRAD:
+ _myGradient->getMinMax(min,max);
+ break;
+ }
+ imin = min;
+ imax = max;
+}
+
+LongSeq* Filter_Gen_i::getHistogram(CORBA::Long size,ref_func rf)
+{
+ int mysize = size;
+ vector<int> myh;
+
+ switch(rf){
+ case F_FIELD:
+ if (_myDField)
+ myh = _myDField->getHistogram(mysize);
+ else
+ myh = _myIField->getHistogram(mysize);
+ break;
+ case F_GRAD:
+ myh = _myGradient->getHistogram(mysize);
+ break;
+ }
+
+ LongSeq *seq = new LongSeq();
+ seq->length(myh.size());
+ for(int i=0;i<myh.size();i++)
+ (*seq)[i] = myh[i];
+ return seq;
+}
+
+void Filter_Gen_i::generateCriteria(CORBA::Long nbthresh,CORBA::Double fthresh,CORBA::Double sthresh,CORBA::Long areaFlag) throw(SALOME_FILTER::FILTER_Gen::FilterError)
+{
+ double val, min, max;
+ double sigmaV;
+ bool isGVal;
+ MED_EN::medEntityMesh typ;
+ const int *revC;
+ const int *indC;
+ FIELD<double> *volume;
+ set <int> listElements;
+ set <int>::iterator elemIt ;
+
+ if(_myDField)
+ typ = _myDField->getSupport()->getEntity();
+ else
+ typ = _myIField->getSupport()->getEntity();
+
+ // create support on nodes
+ SUPPORT *sup = new SUPPORT(_mesh,"Support",MED_NODE);
+
+ // create integer field on nodes
+ _criteria = new FIELD<int>(sup,1);
+
+ _criteria->setName("Criteria");
+
+ // read number of nodes
+ int NumberOf = sup->getNumberOfElements(MED_ALL_ELEMENTS);
+
+ // if reference field is on elements get reference field on nodes
+ switch(typ){
+ case MED_CELL:
+ if(_myDField){
+ // calculate reverse connectivity to have the list of elements which contains node i
+ revC = _myDField->getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
+ indC = _myDField->getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
+ // calculate volume field on mesh
+ volume = _myDField->getSupport()->getMesh()->getVolume(_myDField->getSupport());
+ }
+ else{
+ // calculate reverse connectivity to have the list of elements which contains node i
+ revC = _myIField->getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
+ indC = _myIField->getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
+ // calculate volume field on mesh
+ volume = _myIField->getSupport()->getMesh()->getVolume(_myIField->getSupport());
+ }
+ break;
+ default:
+ break;
+ }
+
+ for (int i=1; i<NumberOf+1; i++){
+
+ // if reference field is on elements get reference field on nodes
+ switch(typ){
+ case MED_CELL:
+ // listElements contains elements which contains a node of element i
+ listElements.clear();
+
+ for(int j=indC[i-1];j<indC[i];j++){
+ // c element contains node i
+ int c=revC[j-1];
+ listElements.insert(c);
+ }
+
+ // calculate field value on node i
+ sigmaV = 0.;
+ val = 0.;
+ for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
+ int elem = *elemIt;
+ double vol = volume->getValueIJ(elem,1);
+ if( vol != 0. ){
+ sigmaV += 1./vol;
+ if(_myDField)
+ val += _myDField->getValueIJ(elem,1)/vol;
+ else
+ val += ((double)_myIField->getValueIJ(elem,1))/vol;
+ }
+ }
+ val /= sigmaV;
+ break;
+ case MED_FACE:
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on faces");
+ break;
+ case MED_EDGE:
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on edges");
+ break;
+ case MED_NODE:
+ // read reference field value
+ if(_myDField)
+ val = _myDField->getValueIJ(i,1);
+ else
+ val = (double)_myIField->getValueIJ(i,1);
+ break;
+ case MED_ALL_ENTITIES:
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Filter doesn't run on reference field on all entities");
+ break;
+ }
+
+ // set criteria field value
+ if( nbthresh == 1 ){
+ if( areaFlag )
+ if( val >= fthresh ) isGVal = true;
+ else isGVal = false;
+ else
+ if( val <= fthresh ) isGVal = true;
+ else isGVal = false;
+ }
+ else{
+ min = fthresh;
+ max = sthresh;
+ if(sthresh < fthresh){
+ min = sthresh;
+ max = fthresh;
+ }
+ if( areaFlag )
+ if( (val <= min) || (val >= max) ) isGVal = true;
+ else isGVal = false;
+ else
+ if( (val >= min) && (val <= max) ) isGVal = true;
+ else isGVal = false;
+ }
+ if( isGVal )
+ _criteria->setValueIJ(i,1,1);
+ else
+ _criteria->setValueIJ(i,1,0);
+ }
+}
+
+void Filter_Gen_i::createEnsightInputFiles()
+{
+ int id;
+
+ MESSAGE("Calculate boundary");
+ // generate MED boundary of geometry
+ SUPPORT *supB = _mesh->getBoundaryElements(MED_FACE);
+
+ // call ensight driver MED to generate input ensight mesh ,
+ // input ensight boundary mesh and input criteria ensight field for filtoo
+ MESSAGE("Create ensight mesh");
+ ENSIGHT_MESH_DRIVER myMeshDriver("/tmp/input.geom",_mesh);
+ myMeshDriver.addSupport(supB);
+ id=_mesh->addDriver(myMeshDriver);
+ _mesh->write(id);
+
+ MESSAGE("Create ensight field");
+ ENSIGHT_FIELD_DRIVER<int> myFieldDriver("/tmp/input.data",_criteria);
+ id=_criteria->addDriver(myFieldDriver);
+ _criteria->write(id);
+}
+
+void Filter_Gen_i::filtering()
+{
+ string command;
+
+ MESSAGE("call filtoo");
+ // send filtoo command
+ command = "cd /tmp;filtoo -f input -o output > /tmp/filter.log";
+ MESSAGE(command);
+ system(command.c_str());
+
+ // destroy filtoo input files
+// command = "cd /tmp;rm -f input.*";
+// MESSAGE(command);
+// system(command.c_str());
+
+}
+
+void Filter_Gen_i::projectFieldsOnDecimateMesh() throw(SALOME_FILTER::FILTER_Gen::FilterError)
+{
+ string command;
+
+ // have to call ensight driver MED to generate output MED file from filtoo output
+ MED _newMed(ENSIGHT_DRIVER,"/tmp/output.case");
+ _newMed.read();
+
+ // destroy filtoo output files
+// command = "cd /tmp;rm -f output.*";
+// MESSAGE(command);
+// system(command.c_str());
+
+ // get new mesh name
+ deque<string> meshesNames = _newMed.getMeshNames();
+ int numberOfMeshes = meshesNames.size();
+ if( numberOfMeshes != 1)
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid number of meshes in filtoo output");
+
+ // new mesh generated by filtoo
+ _newMesh = _newMed.getMesh(meshesNames[0]);
+
+ // create support on nodes on all new mesh
+ SUPPORT *newSup = new SUPPORT(_newMesh,"Support",MED_NODE);
+
+ // read the id of nodes of output mesh, in input mesh
+ readMapping();
+
+ // read connectivity of new mesh to get neighbour node of created node
+ _connL = _newMesh->getConnectivityLength(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
+ _conn = _newMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
+ _connI = _newMesh->getConnectivityIndex(MED_NODAL,MED_CELL);
+
+ // read number of nodes on new mesh
+ int numberOfNodes = newSup->getNumberOfElements(MED_ALL_ELEMENTS);
+ int numberOfComponents;
+
+ deque<string> fieldsNames = _med->getFieldNames();
+ int numberOfFields = fieldsNames.size();
+
+ try{
+
+ // loop on fields
+ for (int i=0; i<numberOfFields; i++){
+
+ // is the input field the reference field?
+ bool isReferenceField= false;
+ if(_myDField){
+ if( strcmp(_myDField->getName().c_str(),fieldsNames[i].c_str()) == 0)
+ isReferenceField = true;
+ }
+ else
+ if( strcmp(_myIField->getName().c_str(),fieldsNames[i].c_str()) == 0)
+ isReferenceField = true;
+
+ deque<DT_IT_> myIteration = _med->getFieldIteration (fieldsNames[i]);
+ string meshName = _med->getField(fieldsNames[i],myIteration[0].dt,myIteration[0].it)->getSupport()->getMesh()->getName();
+
+ // we process only fields on input mesh
+ if( strcmp(meshName.c_str(),_mesh->getName().c_str()) == 0){
+
+ // loop on time steps
+ int numberOfIteration = myIteration.size();
+ for(int j=0;j<numberOfIteration;j++){
+
+ // select input field
+ MEDMEM::FIELD_* field = _med->getField(fieldsNames[i],myIteration[j].dt,myIteration[j].it);
+ FIELD<double> *myDField = NULL;
+ FIELD<double> *newDField = NULL;
+ FIELD<int> *myIField = NULL;
+ FIELD<int> *newIField = NULL;
+
+ if (dynamic_cast<MEDMEM::FIELD<double>*>(field)){
+ if(!isReferenceField)
+ // read input field on input file
+ myDField = new FIELD<double>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
+ else
+ myDField = _myDField;
+ // create new output field
+ newDField = new FIELD<double>(newSup,field->getNumberOfComponents());
+ newDField->setName(myDField->getName());
+ }
+ else{
+ if(!isReferenceField)
+ // read input field on input file
+ myIField = new FIELD<int>(MEDMEM::MED_DRIVER,_file,fieldsNames[i],myIteration[j].dt,myIteration[j].it);
+ else
+ myIField = _myIField;
+ // create new output field
+ newIField = new FIELD<int>(newSup,field->getNumberOfComponents());
+ newIField->setName(myIField->getName());
+ }
+
+ // loop on nodes on new field
+ for (int k=1; k<numberOfNodes+1; k++){
+ // read number of nodes on input field
+ int l = getNodeNumber(k);
+ double dval;
+ int ival;
+
+ for(int c=1;c<=numberOfComponents;c++){
+ // read value on input field
+ if(myDField)
+ dval = myDField->getValueIJ(l,c);
+ else
+ ival = myIField->getValueIJ(l,c);
+
+ // write value on new field
+ if(newDField)
+ newDField->setValueIJ(k,c,dval);
+ else
+ newIField->setValueIJ(k,c,ival);
+ }
+
+ }
+ if(newDField)
+ _newMed.addField(newDField);
+ else
+ _newMed.addField(newIField);
+
+ // Destroy input field if not referecne field
+ if(!isReferenceField)
+ if(myDField){
+ delete myDField->getSupport();
+ delete myDField;
+ }
+ else{
+ delete myIField->getSupport();
+ delete myIField;
+ }
+
+ }
+ }
+ }
+ }
+ catch(SALOME_Exception){
+ throw SALOME_FILTER::FILTER_Gen::FilterError("Unvalid decimate mlesh created by filtoo");
+ }
+}
+
+void Filter_Gen_i::readMapping()
+{
+ ifstream mapFile("/tmp/output.renum");
+
+ // read number of vertices to map
+ mapFile >> _nbvmap;
+ _map = new int[_nbvmap];
+
+ for(int i=0;i<_nbvmap;i++)
+ mapFile >> _map[i];
+
+}
+
+int Filter_Gen_i::getNodeNumber(int num)
+{
+ int oldnum = _map[num];
+
+ // if new node get neighbour node
+ if(oldnum == 0)
+ oldnum = _map[getNeighbourVertex(num)];
+
+ return oldnum;
}
+int Filter_Gen_i::getNeighbourVertex(int num) throw(SALOME_Exception)
+{
+ int nnum;
+ int numberOfElements = _newMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
+
+ int index;
+ // get index of created node in connectivity array
+ for(index=0;index<_connL;index++){
+ if(_conn[index] == num){
+
+ // get index of element which contain created node
+ int i;
+ for(i=1;i<=numberOfElements;i++)
+ if(_connI[i] > index)
+ break;
+
+ // search neighbour node which are in old mesh
+ for(int j=_connI[i-1];j<_connI[i];j++){
+ nnum = _conn[j-1];
+ if( _map[nnum] != 0)
+ break;
+ }
+
+ // if neighbour node in old mesh: go out loop, else continue
+ if(_map[nnum]!=0)
+ break;
+ }
+ }
+
+ // if no neighbour node in old mesh: throw exception
+ if(_map[nnum]==0)
+ throw SALOME_Exception("None of the neighbour node are in old mesh!!");
+
+ return nnum;
+}
+
+void Filter_Gen_i::createMedOutputFile(const char *outMedFile)
+{
+ int id;
+
+ MESSAGE("Create ensight mesh");
+ id = _newMed->addDriver(MED_DRIVER,outMedFile);
+ _newMed->write(id);
+}
+
+
+//=============================================================================
+/*!
+ * C factory, accessible with dlsym, after dlopen
+ */
+//=============================================================================
+
+extern "C"
+{
+ PortableServer::ObjectId * FILTEREngine_factory(
+ CORBA::ORB_ptr orb,
+ PortableServer::POA_ptr poa,
+ PortableServer::ObjectId * contId,
+ const char *instanceName,
+ const char *interfaceName)
+ {
+ MESSAGE("PortableServer::ObjectId * FilterEngine_factory()");
+ SCRUTE(interfaceName);
+ Filter_Gen_i * myFilter_Gen
+ = new Filter_Gen_i(orb, poa, contId, instanceName, interfaceName);
+ return myFilter_Gen->getId() ;
+ }
+}
