Salome HOME
update for take accounbt of field on partial support
[modules/filter.git] / src / FILTER / field2nodes.cxx
index 06394789d156ee0452f4cb8066166fa8cb064880..b6c275812d9394694ca172bf12484dcc666e69e7 100644 (file)
@@ -43,8 +43,8 @@ int main (int argc, char ** argv) {
 
   string filenameIN ;
   string filenameOUT;
-  set <int> listElements;
-  set <int>::iterator elemIt ;
+  set <int> listElements, supset;
+  set <int>::iterator elemIt, ssIt;
   
   if ( argc == 3 ) {
     filenameIN  = argv[1] ;
@@ -56,8 +56,10 @@ int main (int argc, char ** argv) {
     cout << "-> Read all meshes "  ;
     int NumberOfMeshes = myMed.getNumberOfMeshes() ;
     cout << "( "<<NumberOfMeshes << " ) :" << endl ;
+    FIELD<double> **volume = new FIELD<double>*[NumberOfMeshes];
     deque<string> MeshName = myMed.getMeshNames() ;
     for (int i=0; i<NumberOfMeshes; i++) {
+      volume[i] = NULL;
       myMed.getMesh(MeshName[i])->read() ;
       cout << "-> Mesh "<<i+1<<", named "<<MeshName[i]<<" is read !" << endl;
       int id = myMed.getMesh(MeshName[i])->addDriver(MED_DRIVER,filenameOUT,MeshName[i],MED_EN::MED_ECRI);
@@ -86,83 +88,181 @@ int main (int argc, char ** argv) {
        switch(myField->getSupport()->getEntity()){
        case MED_CELL:
          cout << "*************************** CHAMP AUX CELLULES" << endl;
-         MESH *mesh = myMed.getField(FieldName[i],FieldIteration[j].dt,FieldIteration[j].it)->getSupport()->getMesh();
-         SUPPORT *sup = new SUPPORT(mesh,"Support",MED_NODE);
+         MESH *mesh = myField->getSupport()->getMesh();
+
+         // create new support for new field
+         SUPPORT *newSup;
+         if( myField->getSupport()->isOnAllElements() )
+           newSup = new SUPPORT(mesh,"Support",MED_NODE);
+         else{
+           int nbe = myField->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
+           const int *numb = myField->getSupport()->getNumber(MED_ALL_ELEMENTS);
+           list<int> myList;
+           for(int k=0;k<nbe;k++){
+             myList.push_back(numb[k]);
+             supset.insert(numb[k]);
+           }
+           newSup = mesh->buildSupportOnNodeFromElementList(myList,MED_CELL);
+         }
+
          // read number of nodes
-         int NumberOfNodes = sup->getNumberOfElements(MED_ALL_ELEMENTS);
+         int NumberOfNodes = newSup->getNumberOfElements(MED_ALL_ELEMENTS);
          // calculate reverse connectivity to have the list of elements which contains node i
          const int *revC = myField->getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
          const int *indC = myField->getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
+         int ivol;
          // calculate volume field on mesh
-         FIELD<double> *volume = myField->getSupport()->getMesh()->getVolume(myField->getSupport());
+         for(int k=0;k<NumberOfMeshes;k++)
+           if(strcmp(mesh->getName().c_str(),MeshName[k].c_str())==0)
+             ivol = k;
+         if( volume[ivol] == NULL )
+           volume[ivol] = myField->getSupport()->getMesh()->getVolume(myField->getSupport());
          if (dynamic_cast<MEDMEM::FIELD<double>*>(myField)){
            FIELD<double> *myDField = (MEDMEM::FIELD<double>*)myField;
-           FIELD<double> *newDField = new FIELD<double>(sup,NumberOfComponents);
+           FIELD<double> *newDField = new FIELD<double>(newSup,NumberOfComponents);
            newDField->setName(myField->getName());
+           newDField->setIterationNumber(FieldIteration[j].dt);
+           newDField->setOrderNumber(FieldIteration[j].it);
+           newDField->setTime(myDField->getTime());
            double *val = new double[NumberOfComponents];
-           for (int k=1; k<NumberOfNodes+1; k++){
-             // listElements contains elements which contains a node of element i
-             listElements.clear();
-             for(int j=indC[k-1];j<indC[k];j++){
-               // c element contains node i
-               int c=revC[j-1];
-               listElements.insert(c);
+           if( newSup->isOnAllElements() ){
+             for (int k=1; k<NumberOfNodes+1; k++){
+               // listElements contains elements which contains a node of element i
+               listElements.clear();
+               for(int l=indC[k-1];l<indC[k];l++){
+                 // c element contains node k
+                 int c=revC[l-1];
+                 listElements.insert(c);
+               }
+               
+               // calculate field value on node k 
+               double sigmaV = 0.;
+               for(int l=0;l<NumberOfComponents;l++)
+                 val[l] = 0.;
+               for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
+                 int elem = *elemIt;
+                 double vol = volume[ivol]->getValueIJ(elem,1);
+                 if( vol != 0. ){
+                   sigmaV += 1./vol;
+                   for(int l=1;l<=NumberOfComponents;l++)
+                     val[l-1] += myDField->getValueIJ(elem,l)/vol;
+                 }
+               }
+               for(int l=1;l<=NumberOfComponents;l++)
+                 newDField->setValueIJ(k,l,val[l-1]/sigmaV);
              }
-
-             // calculate field value on node i 
-             double sigmaV = 0.;
-             for(int j=0;j<NumberOfComponents;j++)
-               val[j] = 0.;
-             for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
-               int elem = *elemIt;
-               double vol = volume->getValueIJ(elem,1);
-               if( vol != 0. ){
-                 sigmaV += 1./vol;
-                 for(int j=1;j<=NumberOfComponents;j++)
-                   val[j] += myDField->getValueIJ(elem,j)/vol;
+           }
+           else{
+             const int *numb = newSup->getNumber(MED_ALL_ELEMENTS);
+             for (int k=1; k<=NumberOfNodes; k++){
+               int kk = numb[k-1];
+               // listElements contains elements which contains a node of element i
+               listElements.clear();
+               for(int l=indC[kk-1];l<indC[kk];l++){
+                 // c element contains node k
+                 int c=revC[l-1];
+                 // add c element only if it is in partial support
+                 ssIt = supset.find(c);
+                 if(ssIt!=supset.end())
+                   listElements.insert(c);
+               }
+               
+               // calculate field value on node k 
+               double sigmaV = 0.;
+               for(int l=0;l<NumberOfComponents;l++)
+                 val[l] = 0.;
+               for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
+                 int elem = *elemIt;
+                 double vol = volume[ivol]->getValueIJ(elem,1);
+                 if( vol != 0. ){
+                   sigmaV += 1./vol;
+                   for(int l=1;l<=NumberOfComponents;l++)
+                     val[l-1] += myDField->getValueIJ(elem,l)/vol;
+                 }
                }
+               for(int l=1;l<=NumberOfComponents;l++)
+                 newDField->setValueIJ(k,l,val[l-1]/sigmaV);
              }
-             for(int j=1;j<=NumberOfComponents;j++)
-               newDField->setValueIJ(k,j,val[j]/sigmaV);
            }
            delete [] val;
            int id = newDField->addDriver(MED_DRIVER,filenameOUT,FieldName[i],MED_EN::MED_ECRI);
            newDField->write(id);
+           delete newSup;
+           delete newDField;
            cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is written !" << endl;
          }
          else{
            FIELD<int> *myIField = (MEDMEM::FIELD<int>*)myField;
-           FIELD<int> *newIField = new FIELD<int>(sup,NumberOfComponents);
+           FIELD<int> *newIField = new FIELD<int>(newSup,NumberOfComponents);
            newIField->setName(myField->getName());
+           newIField->setIterationNumber(FieldIteration[j].dt);
+           newIField->setOrderNumber(FieldIteration[j].it);
+           newIField->setTime(myIField->getTime());
            double *val = new double[NumberOfComponents];
-           for (int k=1; k<NumberOfNodes+1; k++){
-             // listElements contains elements which contains a node of element i
-             listElements.clear();
-             for(int j=indC[k-1];j<indC[k];j++){
-               // c element contains node i
-               int c=revC[j-1];
-               listElements.insert(c);
-             }
+           if( newSup->isOnAllElements() ){
+             for (int k=1; k<NumberOfNodes+1; k++){
+               // listElements contains elements which contains a node of element i
+               listElements.clear();
+               for(int l=indC[k-1];l<indC[k];l++){
+                 // c element contains node i
+                 int c=revC[l-1];
+                 listElements.insert(c);
+               }
 
-             // calculate field value on node i 
-             double sigmaV = 0.;
-             for(int j=0;j<NumberOfComponents;j++)
-               val[j] = 0.;
-             for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
-               int elem = *elemIt;
-               double vol = volume->getValueIJ(elem,1);
-               if( vol != 0. ){
-                 sigmaV += 1./vol;
-                 for(int j=1;j<=NumberOfComponents;j++)
-                   val[j] += ((double)myIField->getValueIJ(elem,j))/vol;
+               // calculate field value on node k 
+               double sigmaV = 0.;
+               for(int l=0;l<NumberOfComponents;l++)
+                 val[l] = 0.;
+               for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
+                 int elem = *elemIt;
+                 double vol = volume[ivol]->getValueIJ(elem,1);
+                 if( vol != 0. ){
+                   sigmaV += 1./vol;
+                   for(int l=1;l<=NumberOfComponents;l++)
+                     val[l-1] += ((double)myIField->getValueIJ(elem,l))/vol;
+                 }
+               }
+               for(int l=1;l<=NumberOfComponents;l++)
+                 newIField->setValueIJ(k,l,(int)(val[l-1]/sigmaV));
+             }
+           }
+           else{
+             const int *numb = newSup->getNumber(MED_ALL_ELEMENTS);
+             for (int k=1; k<=NumberOfNodes; k++){
+               int kk = numb[k-1];
+               // listElements contains elements which contains a node of element i
+               listElements.clear();
+               for(int l=indC[kk-1];l<indC[kk];l++){
+                 // c element contains node k
+                 int c=revC[l-1];
+                 // add c element only if it is in partial support
+                 ssIt = supset.find(c);
+                 if(ssIt!=supset.end())
+                   listElements.insert(c);
+               }
+               
+               // calculate field value on node k 
+               double sigmaV = 0.;
+               for(int l=0;l<NumberOfComponents;l++)
+                 val[l] = 0.;
+               for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
+                 int elem = *elemIt;
+                 double vol = volume[ivol]->getValueIJ(elem,1);
+                 if( vol != 0. ){
+                   sigmaV += 1./vol;
+                   for(int l=1;l<=NumberOfComponents;l++)
+                     val[l-1] += myIField->getValueIJ(elem,l)/vol;
+                 }
                }
+               for(int l=1;l<=NumberOfComponents;l++)
+                 newIField->setValueIJ(k,l,(int)(val[l-1]/sigmaV));
              }
-             for(int j=1;j<=NumberOfComponents;j++)
-               newIField->setValueIJ(k,j,(int)(val[j]/sigmaV));
            }
            delete [] val;
            int id = newIField->addDriver(MED_DRIVER,filenameOUT,FieldName[i],MED_EN::MED_ECRI);
            newIField->write(id);
+           delete newSup;
+           delete newIField;
            cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is written !" << endl;
          }
          break;