CompuCell/DOXYGEN/core/CompuCell3D/plugins/ExternalPotential/ExternalPotentialPlugin.cpp

Go to the documentation of this file.

#include <CompuCell3D/Simulator.h>
#include <CompuCell3D/Potts3D/Potts3D.h>
#include <CompuCell3D/Field3D/Field3D.h>
#include <CompuCell3D/Field3D/AdjacentNeighbor.h>
#include <CompuCell3D/Automaton/Automaton.h>
#include <BasicUtils/BasicClassGroup.h>


using namespace CompuCell3D;
using namespace std;

#define EXP_STL
#include "ExternalPotentialPlugin.h"

namespace CompuCell3D {

        ExternalPotentialPlugin::ExternalPotentialPlugin():lambdaVec(Coordinates3D<float>(0.,0.,0.)),xmlData(0)
{}
ExternalPotentialPlugin::~ExternalPotentialPlugin()
{
}

void ExternalPotentialPlugin::init(Simulator *_simulator, CC3DXMLElement *_xmlData){
  
  xmlData=_xmlData;
  potts = _simulator->getPotts();
  cellFieldG=(WatchableField3D<CellG *> * )potts->getCellFieldG();

  
  potts->registerEnergyFunctionWithName(this,"ExternalPotential");
  
   Dim3D fieldDim=potts->getCellFieldG()->getDim();
   
   
   

  adjNeighbor.initialize(fieldDim);
  adjNeighbor_ptr=&adjNeighbor;

  if(potts->getBoundaryXName()=="Periodic"){
      adjNeighbor.setPeriodicX();
  }
  if(potts->getBoundaryYName()=="Periodic"){
      adjNeighbor.setPeriodicY();
  }
  if(potts->getBoundaryZName()=="Periodic"){
      adjNeighbor.setPeriodicZ();
  }

  _simulator->registerSteerableObject(this);



}

void ExternalPotentialPlugin::extraInit(Simulator *_simulator){
   // check why registering in init gives segfault
//    _simulator->registerSteerableObject(this);

                update(xmlData,true);
}




void ExternalPotentialPlugin::update(CC3DXMLElement *_xmlData, bool _fullInitFlag){


        if(!_xmlData->getNumberOfChildren()){ 
                functionType=BYCELLID;
        }else{
                if(_xmlData->findElement("ExternalPotentialParameters"))
                        functionType=BYCELLTYPE;
                else if (_xmlData->findElement("Lambda"))
                        functionType=GLOBAL;
                else //in case users put garbage xml use changeEnergyByCellId
                        functionType=BYCELLID;
        }

        Automaton *automaton=potts->getAutomaton();
        cerr<<"automaton="<<automaton<<endl;

        switch(functionType){
                case BYCELLID:
                        //set fcn ptr
                        changeEnergyFcnPtr=&ExternalPotentialPlugin::changeEnergyByCellId;
                        break;

                case BYCELLTYPE:
                        {
                                externalPotentialParamVector.clear();
                                vector<int> typeIdVec;
                                vector<ExternalPotentialParam> externalPotentialParamVectorTmp;

                                CC3DXMLElementList energyVec=_xmlData->getElements("ExternalPotentialParameters");

                                for (int i = 0 ; i<energyVec.size(); ++i){
                                        ExternalPotentialParam extPotentialParam;

                                        extPotentialParam.lambdaVec.x=energyVec[i]->getAttributeAsDouble("x");
                                        extPotentialParam.lambdaVec.y=energyVec[i]->getAttributeAsDouble("y");
                                        extPotentialParam.lambdaVec.z=energyVec[i]->getAttributeAsDouble("z");

                                        extPotentialParam.typeName=energyVec[i]->getAttribute("CellType");
                                        cerr<<"automaton="<<automaton<<endl;
                                        typeIdVec.push_back(automaton->getTypeId(extPotentialParam.typeName));

                                        externalPotentialParamVectorTmp.push_back(extPotentialParam);                           
                                }
                                vector<int>::iterator pos=max_element(typeIdVec.begin(),typeIdVec.end());
                                int maxTypeId=*pos;
                                externalPotentialParamVector.assign(maxTypeId+1,ExternalPotentialParam());
                                for (int i = 0 ; i < externalPotentialParamVectorTmp.size() ; ++i){
                                        externalPotentialParamVector[typeIdVec[i]]=externalPotentialParamVectorTmp[i];
                                }

                                //set fcn ptr
                                changeEnergyFcnPtr=&ExternalPotentialPlugin::changeEnergyByCellType;
                        }
                        break;

                case GLOBAL:
                        //using Global Volume Energy Parameters
                        lambdaVec=Coordinates3D<float>(_xmlData->getFirstElement("Lambda")->getAttributeAsDouble("x"),
                                                                                                         _xmlData->getFirstElement("Lambda")->getAttributeAsDouble("y"),
                                                                                                         _xmlData->getFirstElement("Lambda")->getAttributeAsDouble("z")
                                );
                        //set fcn ptr
                        changeEnergyFcnPtr=&ExternalPotentialPlugin::changeEnergyGlobal;
                        break;

                default:
                        //set fcn ptr
                        changeEnergyFcnPtr=&ExternalPotentialPlugin::changeEnergyByCellId;
        }

}


double ExternalPotentialPlugin::changeEnergyGlobal(const Point3D &pt,  const CellG *newCell,
                                                const CellG *oldCell)
{



   double deltaEnergyOld=0.0;
   double deltaEnergyNew=0.0;
   CellG *neighborPtr;

   Dim3D fieldDim=cellFieldG->getDim();


   const vector<Point3D>  & neighborsOffsetVec = adjNeighbor_ptr->getAdjFace2FaceNeighborOffsetVec(pt);
   unsigned int neighborSize=neighborsOffsetVec.size();

   Point3D ptAdj;
   

   //x, y,or z depending in which direction we want potential gradient to act
   //I add 1 to avoid zero change of energy when z is at the boundary
   short prefferedCoordinate=pt.z+1;
   //short deltaCoordinate;
   Coordinates3D<short> deltaCoordinate(0,0,0);
   
   // have to do fieldDim.z-2 because the neighbor of pt with max_z can be a pt with z=0 but they are spaced by delta z=1
   // thus you need to do deltaCoordinate % (fieldDim.z-2). Otherwise if you do not do fieldDim.z-2
   // then you may get the following (max_z-0) % max_z =0, whereas it should be 1 !

//    cerr<<"******************CHANGE BEGIN**************************"<<endl;
   int counter=0;
   Point3D ptFlipNeighbor=potts->getFlipNeighbor();
   Point3D deltaFlip;
   
   deltaFlip=pt;
   deltaFlip-=ptFlipNeighbor;
      
   for(unsigned int i = 0 ; i < neighborSize ; ++i){
      ptAdj=pt;
      ptAdj+=neighborsOffsetVec[i];
      
      if( cellFieldG->isValid(ptAdj) ){
         ++counter;
         neighborPtr=cellFieldG->get(ptAdj);
         
         if(/*neighborPtr &&*/ oldCell && neighborPtr!=oldCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);

            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){

               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );


             }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }


            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
            deltaEnergyOld+=  deltaCoordinate.X()*lambdaVec.X()+
                              deltaCoordinate.Y()*lambdaVec.Y()+
                              deltaCoordinate.Z()*lambdaVec.Z();

         }

         if(/*neighborPtr && */newCell && neighborPtr!=newCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);
            
            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){
               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );
               
            }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }
                        

            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
            
            deltaEnergyNew+=  deltaCoordinate.X()*lambdaVec.X()+
                              deltaCoordinate.Y()*lambdaVec.Y()+
                              deltaCoordinate.Z()*lambdaVec.Z();
            
         }


      }
       
   }


   return deltaEnergyNew-deltaEnergyOld;
   
}      


double ExternalPotentialPlugin::changeEnergyByCellType(const Point3D &pt,  const CellG *newCell,
                                                const CellG *oldCell)
{

        

   double deltaEnergyOld=0.0;
   double deltaEnergyNew=0.0;
   CellG *neighborPtr;

   Dim3D fieldDim=cellFieldG->getDim();


   const vector<Point3D>  & neighborsOffsetVec = adjNeighbor_ptr->getAdjFace2FaceNeighborOffsetVec(pt);
   unsigned int neighborSize=neighborsOffsetVec.size();

   Point3D ptAdj;
   

   //x, y,or z depending in which direction we want potential gradient to act
   //I add 1 to avoid zero change of energy when z is at the boundary
   short prefferedCoordinate=pt.z+1;
   //short deltaCoordinate;
   Coordinates3D<short> deltaCoordinate(0,0,0);
   
   // have to do fieldDim.z-2 because the neighbor of pt with max_z can be a pt with z=0 but they are spaced by delta z=1
   // thus you need to do deltaCoordinate % (fieldDim.z-2). Otherwise if you do not do fieldDim.z-2
   // then you may get the following (max_z-0) % max_z =0, whereas it should be 1 !

//    cerr<<"******************CHANGE BEGIN**************************"<<endl;
   int counter=0;
   Point3D ptFlipNeighbor=potts->getFlipNeighbor();
   Point3D deltaFlip;
   
   deltaFlip=pt;
   deltaFlip-=ptFlipNeighbor;
      
   for(unsigned int i = 0 ; i < neighborSize ; ++i){
      ptAdj=pt;
      ptAdj+=neighborsOffsetVec[i];
      
      if( cellFieldG->isValid(ptAdj) ){
         ++counter;
         neighborPtr=cellFieldG->get(ptAdj);
         
         if(/*neighborPtr &&*/ oldCell && neighborPtr!=oldCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);

            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){

               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );


             }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }


            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
                                
            deltaEnergyOld+=  deltaCoordinate.X()*externalPotentialParamVector[oldCell->type].lambdaVec.X()+
                              deltaCoordinate.Y()*externalPotentialParamVector[oldCell->type].lambdaVec.Y()+
                              deltaCoordinate.Z()*externalPotentialParamVector[oldCell->type].lambdaVec.Z();

         }

         if(/*neighborPtr && */newCell && neighborPtr!=newCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);
            
            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){
               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );
               
            }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }
                        

            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
            
            deltaEnergyNew+=  deltaCoordinate.X()*externalPotentialParamVector[newCell->type].lambdaVec.X()+
                              deltaCoordinate.Y()*externalPotentialParamVector[newCell->type].lambdaVec.Y()+
                              deltaCoordinate.Z()*externalPotentialParamVector[newCell->type].lambdaVec.Z();
            
         }


      }
       
   }


   return deltaEnergyNew-deltaEnergyOld;
   
}      

double ExternalPotentialPlugin::changeEnergyByCellId(const Point3D &pt,  const CellG *newCell,
                                                const CellG *oldCell)
{

        

   double deltaEnergyOld=0.0;
   double deltaEnergyNew=0.0;
   CellG *neighborPtr;

   Dim3D fieldDim=cellFieldG->getDim();


   const vector<Point3D>  & neighborsOffsetVec = adjNeighbor_ptr->getAdjFace2FaceNeighborOffsetVec(pt);
   unsigned int neighborSize=neighborsOffsetVec.size();

   Point3D ptAdj;
   

   //x, y,or z depending in which direction we want potential gradient to act
   //I add 1 to avoid zero change of energy when z is at the boundary
   short prefferedCoordinate=pt.z+1;
   //short deltaCoordinate;
   Coordinates3D<short> deltaCoordinate(0,0,0);
   
   // have to do fieldDim.z-2 because the neighbor of pt with max_z can be a pt with z=0 but they are spaced by delta z=1
   // thus you need to do deltaCoordinate % (fieldDim.z-2). Otherwise if you do not do fieldDim.z-2
   // then you may get the following (max_z-0) % max_z =0, whereas it should be 1 !

//    cerr<<"******************CHANGE BEGIN**************************"<<endl;
   int counter=0;
   Point3D ptFlipNeighbor=potts->getFlipNeighbor();
   Point3D deltaFlip;
   
   deltaFlip=pt;
   deltaFlip-=ptFlipNeighbor;
      
   for(unsigned int i = 0 ; i < neighborSize ; ++i){
      ptAdj=pt;
      ptAdj+=neighborsOffsetVec[i];
      
      if( cellFieldG->isValid(ptAdj) ){
         ++counter;
         neighborPtr=cellFieldG->get(ptAdj);
         
         if(/*neighborPtr &&*/ oldCell && neighborPtr!=oldCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);

            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){

               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );


             }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }


            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
                                
            deltaEnergyOld+=  deltaCoordinate.X()*oldCell->lambdaVecX+
                              deltaCoordinate.Y()*oldCell->lambdaVecY+
                              deltaCoordinate.Z()*oldCell->lambdaVecZ;

         }

         if(/*neighborPtr && */newCell && neighborPtr!=newCell){

            deltaCoordinate.XRef()=(ptAdj.x-pt.x);
            deltaCoordinate.YRef()=(ptAdj.y-pt.y);
            deltaCoordinate.ZRef()=(ptAdj.z-pt.z);
            
            if(fabs(static_cast<double>(deltaCoordinate.X()))>1){
               deltaCoordinate.XRef()=
               ( deltaCoordinate.X()>0 ? -(deltaCoordinate.X()+1) % (fieldDim.x-1) :-(deltaCoordinate.X()-1) % (fieldDim.x-1) );
               
            }

            if(fabs(static_cast<double>(deltaCoordinate.Y()))>1){
               deltaCoordinate.YRef()=
               ( deltaCoordinate.Y()>0 ? -(deltaCoordinate.Y()+1) % (fieldDim.y-1) :-(deltaCoordinate.Y()-1) % (fieldDim.y-1) );
            }
                        

            if(fabs(static_cast<double>(deltaCoordinate.Z()))>1){
               deltaCoordinate.ZRef()=
               ( deltaCoordinate.Z()>0 ? -(deltaCoordinate.Z()+1) % (fieldDim.z-1) :-(deltaCoordinate.Z()-1) % (fieldDim.z-1) );
            }
            
            deltaEnergyNew+=  deltaCoordinate.X()*newCell->lambdaVecX+
                              deltaCoordinate.Y()*newCell->lambdaVecY+
                              deltaCoordinate.Z()*newCell->lambdaVecZ;
            
         }


      }
       
   }


   return deltaEnergyNew-deltaEnergyOld;
   
} 



double ExternalPotentialPlugin::changeEnergy(const Point3D &pt,  const CellG *newCell,
                                                const CellG *oldCell)
{

        return (this->*changeEnergyFcnPtr)(pt,newCell,oldCell);

}      






std::string ExternalPotentialPlugin::toString(){
      return "ExternalPotential";
}


std::string ExternalPotentialPlugin::steerableName(){
   
   return toString();
}



};

---