CompuCell/DOXYGEN/core/CompuCell3D/steppables/PDESolvers/FlexibleDiffusionSolverFE.cpp

Go to the documentation of this file.

#include <CompuCell3D/Simulator.h>
#include <CompuCell3D/Automaton/Automaton.h>
#include <CompuCell3D/Potts3D/Potts3D.h>
#include <CompuCell3D/Potts3D/CellInventory.h>
#include <CompuCell3D/Field3D/WatchableField3D.h>
#include <CompuCell3D/Field3D/Field3DImpl.h>
#include <CompuCell3D/Field3D/Field3D.h>
#include <CompuCell3D/Field3D/Field3DIO.h>
#include <BasicUtils/BasicClassGroup.h>
#include <CompuCell3D/steppables/BoxWatcher/BoxWatcher.h>

#include <BasicUtils/BasicString.h>
#include <BasicUtils/BasicException.h>
#include <BasicUtils/BasicRandomNumberGenerator.h>
#include <PublicUtilities/StringUtils.h>
#include <string>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>

// std::ostream & operator<<(std::ostream & out,CompuCell3D::DiffusionData & diffData){
//
//
// }


using namespace CompuCell3D;
using namespace std;

#define EXP_STL
#include "FlexibleDiffusionSolverFE.h"

void FlexibleDiffusionSolverSerializer::serialize(){

        for(int i = 0 ; i < solverPtr->diffSecrFieldTuppleVec.size() ; ++i){
                ostringstream outName;

                outName<<solverPtr->diffSecrFieldTuppleVec[i].diffData.fieldName<<"_"<<currentStep<<"."<<serializedFileExtension;
                ofstream outStream(outName.str().c_str());
                solverPtr->outputField( outStream,solverPtr->concentrationFieldVector[i]);
        }

}

void FlexibleDiffusionSolverSerializer::readFromFile(){
        try{
                for(int i = 0 ; i < solverPtr->diffSecrFieldTuppleVec.size() ; ++i){
                        ostringstream inName;
                        inName<<solverPtr->diffSecrFieldTuppleVec[i].diffData.fieldName<<"."<<serializedFileExtension;

                        solverPtr->readConcentrationField(inName.str().c_str(),solverPtr->concentrationFieldVector[i]);;
                }

        } catch (BasicException &e) {
                cerr<<"COULD NOT FIND ONE OF THE FILES"<<endl;
                throw BasicException("Error in reading diffusion fields from file",e);
        }


}



FlexibleDiffusionSolverFE::FlexibleDiffusionSolverFE()
: DiffusableVector<float>(),deltaX(1.0),deltaT(1.0)
{
        serializerPtr=0;
        serializeFlag=false;
        readFromFileFlag=false;
        haveCouplingTerms=false;
        serializeFrequency=0;
        boxWatcherSteppable=0;
}
FlexibleDiffusionSolverFE::~FlexibleDiffusionSolverFE()
{
        if(serializerPtr)
                delete serializerPtr ; serializerPtr=0;
}
void FlexibleDiffusionSolverFE::init(Simulator *_simulator, CC3DXMLElement *_xmlData) {


        simPtr=_simulator;
        simulator=_simulator;
        potts = _simulator->getPotts();
        automaton=potts->getAutomaton();

        cellInventoryPtr=& potts->getCellInventory();

        //   cellFieldG=potts->getCellFieldG();
        cellFieldG = (WatchableField3D<CellG *> *)potts->getCellFieldG();
        fieldDim=cellFieldG->getDim();

        cerr<<"INSIDE INIT"<<endl;



        diffusePtr=&FlexibleDiffusionSolverFE::diffuse;
        secretePtr=&FlexibleDiffusionSolverFE::secrete;


        update(_xmlData,true);

        numberOfFields=diffSecrFieldTuppleVec.size();



        vector<string> concentrationFieldNameVectorTmp; //temporary vector for field names
        concentrationFieldNameVectorTmp.assign(diffSecrFieldTuppleVec.size(),string(""));

        cerr<<"diffSecrFieldTuppleVec.size()="<<diffSecrFieldTuppleVec.size()<<endl;

        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
                //       cerr<<" concentrationFieldNameVector[i]="<<diffDataVec[i].fieldName<<endl;
                //       concentrationFieldNameVector.push_back(diffDataVec[i].fieldName);
                concentrationFieldNameVectorTmp[i] = diffSecrFieldTuppleVec[i].diffData.fieldName;
                cerr<<" concentrationFieldNameVector[i]="<<concentrationFieldNameVectorTmp[i]<<endl;
        }

        //setting up couplingData - field-field interaction terms
        vector<CouplingData>::iterator pos;

        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
                pos=diffSecrFieldTuppleVec[i].diffData.couplingDataVec.begin();
                for(int j = 0 ; j < diffSecrFieldTuppleVec[i].diffData.couplingDataVec.size() ; ++j){

                        for(int idx=0; idx<concentrationFieldNameVectorTmp.size() ; ++idx){
                                if( concentrationFieldNameVectorTmp[idx] == diffSecrFieldTuppleVec[i].diffData.couplingDataVec[j].intrFieldName ){
                                        diffSecrFieldTuppleVec[i].diffData.couplingDataVec[j].fieldIdx=idx;
                                        haveCouplingTerms=true; //if this is called at list once we have already coupling terms and need to proceed differently with scratch field initialization
                                        break;
                                }
                                //this means that required interacting field name has not been found
                                if( idx == concentrationFieldNameVectorTmp.size()-1 ){
                                        //remove this interacting term
                                        //                pos=&(diffDataVec[i].degradationDataVec[j]);
                                        diffSecrFieldTuppleVec[i].diffData.couplingDataVec.erase(pos);
                                }
                        }
                        ++pos;
                }
        }


        cerr<<"FIELDS THAT I HAVE"<<endl;
        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
                cerr<<"Field "<<i<<" name: "<<concentrationFieldNameVectorTmp[i]<<endl;
        }




        cerr<<"FlexibleDiffusionSolverFE: extra Init in read XML"<<endl;

        workFieldDim=Dim3D(fieldDim.x+2,fieldDim.y+2,fieldDim.z+2);
        if(!haveCouplingTerms){
                allocateDiffusableFieldVector(diffSecrFieldTuppleVec.size()+1,workFieldDim); //+1 is for additional scratch field
        }else{
                allocateDiffusableFieldVector(2*diffSecrFieldTuppleVec.size(),workFieldDim); //with coupling terms every field need to have its own scratch field
        }

        //here I need to copy field names from concentrationFieldNameVectorTmp to concentrationFieldNameVector
        //because concentrationFieldNameVector is reallocated with default values once I call allocateDiffusableFieldVector


        for(unsigned int i=0 ; i < concentrationFieldNameVectorTmp.size() ; ++i){
                concentrationFieldNameVector[i]=concentrationFieldNameVectorTmp[i];
        }


        //register fields once they have been allocated
        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
                simPtr->registerConcentrationField(concentrationFieldNameVector[i] , concentrationFieldVector[i]);
                cerr<<"registring field: "<<concentrationFieldNameVector[i]<<" field address="<<concentrationFieldVector[i]<<endl;
        }





        //    exit(0);

        periodicBoundaryCheckVector.assign(3,false);
        string boundaryName;
        boundaryName=potts->getBoundaryXName();
        changeToLower(boundaryName);
        if(boundaryName=="periodic")  {
                periodicBoundaryCheckVector[0]=true;
        }
        boundaryName=potts->getBoundaryYName();
        changeToLower(boundaryName);
        if(boundaryName=="periodic")  {
                periodicBoundaryCheckVector[1]=true;
        }

        boundaryName=potts->getBoundaryZName();
        changeToLower(boundaryName);
        if(boundaryName=="periodic")  {
                periodicBoundaryCheckVector[2]=true;
        }

        simulator->registerSteerableObject(this);

}
void FlexibleDiffusionSolverFE::extraInit(Simulator *simulator){

        if((serializeFlag || readFromFileFlag) && !serializerPtr){
                serializerPtr=new FlexibleDiffusionSolverSerializer();
                serializerPtr->solverPtr=this;
        }

        if(serializeFlag){
                simulator->registerSerializer(serializerPtr);
        }

        bool useBoxWatcher=false;
        for (int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i){
                if(diffSecrFieldTuppleVec[i].diffData.useBoxWatcher){
                        useBoxWatcher=true;
                        break;
                }
        }
        bool steppableAlreadyRegisteredFlag;
        if(useBoxWatcher){
                boxWatcherSteppable=(BoxWatcher*)Simulator::steppableManager.get("BoxWatcher",&steppableAlreadyRegisteredFlag);
                if(!steppableAlreadyRegisteredFlag)
                        boxWatcherSteppable->init(simulator);
        }
}

void FlexibleDiffusionSolverFE::start() {
        //     if(diffConst> (1.0/6.0-0.05) ){ //hard coded condtion for stability of the solutions - assume dt=1 dx=dy=dz=1
        //
        //       cerr<<"CANNOT SOLVE DIFFUSION EQUATION: STABILITY PROBLEM - DIFFUSION CONSTANT TOO LARGE. EXITING..."<<endl;
        //       exit(0);
        //
        //    }


        dt_dx2=deltaT/(deltaX*deltaX);
        if(readFromFileFlag){
                try{

                        serializerPtr->readFromFile();

                } catch (BasicException &e){
                        cerr<<"Going to fail-safe initialization"<<endl;
                        initializeConcentration(); //if there was error, initialize using failsafe defaults
                }

        }else{
                initializeConcentration();//Normal reading from User specified files
        }

}


void FlexibleDiffusionSolverFE::initializeConcentration()
{
        for(unsigned int i = 0 ; i <diffSecrFieldTuppleVec.size() ; ++i)
        {
                if(diffSecrFieldTuppleVec[i].diffData.concentrationFileName.empty()) continue;
                cerr << "fail-safe initialization " << diffSecrFieldTuppleVec[i].diffData.concentrationFileName << endl;
                readConcentrationField(diffSecrFieldTuppleVec[i].diffData.concentrationFileName,concentrationFieldVector[i]);
        }


}

void FlexibleDiffusionSolverFE::step(const unsigned int _currentStep) {

        currentStep=_currentStep;

        (this->*secretePtr)();

        (this->*diffusePtr)();


        if(serializeFrequency>0 && serializeFlag && !(_currentStep % serializeFrequency)){
                serializerPtr->setCurrentStep(currentStep);
                serializerPtr->serialize();
        }

}

void FlexibleDiffusionSolverFE::secreteOnContactSingleField(unsigned int idx){


        SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;


        std::map<unsigned char,SecretionOnContactData>::iterator mitr;
        std::map<unsigned char,SecretionOnContactData>::iterator end_mitr=secrData.typeIdSecrOnContactDataMap.end();

        CellG *currentCellPtr;

        ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
        Array3D_t & concentrationArray = concentrationFieldPtr->getContainer();

        float currentConcentration;
        float secrConst;
        float secrConstMedium=0.0;
        std::map<unsigned char, float> * contactCellMapMediumPtr;
        std::map<unsigned char, float> * contactCellMapPtr;
        std::map<unsigned char, float>::iterator mitrTypeConst;

        bool secreteInMedium=false;
        //the assumption is that medium has type ID 0
        mitr=secrData.typeIdSecrOnContactDataMap.find(automaton->getTypeId("Medium"));

        if( mitr != end_mitr ){
                secreteInMedium=true;
                contactCellMapMediumPtr = &(mitr->second.contactCellMap);
        }



        Point3D pt;
        Neighbor n;
        CellG *nCell=0;
        WatchableField3D<CellG *> *fieldG = (WatchableField3D<CellG *> *)potts->getCellFieldG();
        unsigned char type;



        for (int z = 1; z < workFieldDim.z-1; z++)
                for (int y = 1; y < workFieldDim.y-1; y++)
                        for (int x = 1; x < workFieldDim.x-1; x++){
                                pt=Point3D(x-1,y-1,z-1);
                                currentCellPtr=cellFieldG->get(pt);
                                //             currentCellPtr=cellFieldG->get(pt);
                                currentConcentration = concentrationArray[x][y][z];


                                if(secreteInMedium && ! currentCellPtr){

                                        for (int i = 0  ; i<=maxNeighborIndex/*offsetVec.size()*/ ; ++i ){
                                                n=boundaryStrategy->getNeighborDirect(pt,i);
                                                if(!n.distance)//not a valid neighbor
                                                        continue;
                                                nCell = fieldG->get(n.pt);
                                                //                      nCell = fieldG->get(n.pt);
                                                if(nCell)
                                                        type=nCell->type;
                                                else
                                                        type=0;


                                                mitrTypeConst=contactCellMapMediumPtr->find(type);

                                                if(mitrTypeConst != contactCellMapMediumPtr->end()){//OK to secrete, contact detected
                                                        secrConstMedium = mitrTypeConst->second;

                                                        concentrationArray[x][y][z]=currentConcentration+secrConstMedium;
                                                }



                                        }

                                        continue;
                                }

                                if(currentCellPtr){
                                        mitr=secrData.typeIdSecrOnContactDataMap.find(currentCellPtr->type);
                                        if(mitr!=end_mitr){

                                                contactCellMapPtr = &(mitr->second.contactCellMap);

                                                for (int i = 0  ; i<=maxNeighborIndex/*offsetVec.size() */; ++i ){

                                                        n=boundaryStrategy->getNeighborDirect(pt,i);
                                                        if(!n.distance)//not a valid neighbor
                                                                continue;
                                                        nCell = fieldG->get(n.pt);
                                                        //                      nCell = fieldG->get(n.pt);
                                                        if(nCell)
                                                                type=nCell->type;
                                                        else
                                                                type=0;

                                                        if (currentCellPtr==nCell) continue; //skip secretion in pixels belongin to the same cell

                                                        mitrTypeConst=contactCellMapPtr->find(type);
                                                        if(mitrTypeConst != contactCellMapPtr->end()){//OK to secrete, contact detected
                                                                secrConst=mitrTypeConst->second;
                                                                //                         concentrationField->set(pt,currentConcentration+secrConst);
                                                                concentrationArray[x][y][z]=currentConcentration+secrConst;
                                                        }

                                                }


                                        }
                                }
                        }

}

void FlexibleDiffusionSolverFE::secreteSingleField(unsigned int idx){


        SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;


        std::map<unsigned char,float>::iterator mitr;   
        std::map<unsigned char,float>::iterator end_mitr=secrData.typeIdSecrConstMap.end();
        std::map<unsigned char,UptakeData>::iterator mitrUptake;
        std::map<unsigned char,UptakeData>::iterator end_mitrUptake=secrData.typeIdUptakeDataMap.end();



        CellG *currentCellPtr;
        Field3DImpl<float> * concentrationField=concentrationFieldVector[idx];
        float currentConcentration;
        float secrConst;
        float secrConstMedium=0.0;
        float maxUptakeInMedium=0.0;
        float relativeUptakeRateInMedium=0.0;

        ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
        Array3D_t & concentrationArray = concentrationFieldPtr->getContainer();

        bool doUptakeFlag=false;
        bool uptakeInMediumFlag=false;
        bool secreteInMedium=false;
        //the assumption is that medium has type ID 0
        mitr=secrData.typeIdSecrConstMap.find(automaton->getTypeId("Medium"));

        if( mitr != end_mitr){
                secreteInMedium=true;
                secrConstMedium=mitr->second;
        }

        //uptake for medium setup
        if(secrData.typeIdUptakeDataMap.size()){
                doUptakeFlag=true;
        }
        //uptake for medium setup
        if(doUptakeFlag){
                mitrUptake=secrData.typeIdUptakeDataMap.find(automaton->getTypeId("Medium"));
                if(mitrUptake != end_mitrUptake){
                        maxUptakeInMedium=mitrUptake->second.maxUptake;
                        relativeUptakeRateInMedium=mitrUptake->second.relativeUptakeRate;
                        uptakeInMediumFlag=true;

                }
        }

        
        Point3D pt;
        for (int z = 1; z < workFieldDim.z-1; z++)
                for (int y = 1; y < workFieldDim.y-1; y++)
                        for (int x = 1; x < workFieldDim.x-1; x++){
                                pt=Point3D(x-1,y-1,z-1);
                                //             cerr<<"pt="<<pt<<" is valid "<<cellFieldG->isValid(pt)<<endl;
                                currentCellPtr=cellFieldG->get(pt);
                                //             currentCellPtr=cellFieldG->get(pt);
                                //             cerr<<"THIS IS PTR="<<currentCellPtr<<endl;

                                //             if(currentCellPtr)
                                //                cerr<<"This is id="<<currentCellPtr->id<<endl;
                                currentConcentration = concentrationArray[x][y][z];

                                if(secreteInMedium && ! currentCellPtr){
                                        concentrationArray[x][y][z]=currentConcentration+secrConstMedium;
                                }

                                if(currentCellPtr){
                                        mitr=secrData.typeIdSecrConstMap.find(currentCellPtr->type);
                                        if(mitr!=end_mitr){
                                                secrConst=mitr->second;
                                                concentrationArray[x][y][z]=currentConcentration+secrConst;
                                        }
                                }

                                if(doUptakeFlag){
                                        if(uptakeInMediumFlag && ! currentCellPtr){
                                                if(currentConcentration>maxUptakeInMedium){
                                                        concentrationArray[x][y][z]-=maxUptakeInMedium; 
                                                }else{
                                                        concentrationArray[x][y][z]-= currentConcentration*relativeUptakeRateInMedium;  
                                                }
                                        }
                                        if(currentCellPtr){
                                                
                                                mitrUptake=secrData.typeIdUptakeDataMap.find(currentCellPtr->type);
                                                if(mitrUptake!=end_mitrUptake){
                                                        if(currentConcentration > mitrUptake->second.maxUptake){
                                                                concentrationArray[x][y][z]-=mitrUptake->second.maxUptake;
                                                                //cerr<<" uptake concentration="<< currentConcentration<<" relativeUptakeRate="<<mitrUptake->second.relativeUptakeRate<<" subtract="<<mitrUptake->second.maxUptake<<endl;
                                                        }else{
                                                                //cerr<<"concentration="<< currentConcentration<<" relativeUptakeRate="<<mitrUptake->second.relativeUptakeRate<<" subtract="<<currentConcentration*mitrUptake->second.relativeUptakeRate<<endl;
                                                                concentrationArray[x][y][z]-= currentConcentration*mitrUptake->second.relativeUptakeRate;       
                                                        }
                                                }
                                        }
                                }
                        }


}


void FlexibleDiffusionSolverFE::secreteConstantConcentrationSingleField(unsigned int idx){


        SecretionData & secrData=diffSecrFieldTuppleVec[idx].secrData;


        std::map<unsigned char,float>::iterator mitr;
        std::map<unsigned char,float>::iterator end_mitr=secrData.typeIdSecrConstMap.end();

        CellG *currentCellPtr;
        Field3DImpl<float> * concentrationField=concentrationFieldVector[idx];
        float currentConcentration;
        float secrConst;
        float secrConstMedium=0.0;

        ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
        Array3D_t & concentrationArray = concentrationFieldPtr->getContainer();


        bool secreteInMedium=false;
        //the assumption is that medium has type ID 0
        mitr=secrData.typeIdSecrConstMap.find(automaton->getTypeId("Medium"));

        if( mitr != end_mitr){
                secreteInMedium=true;
                secrConstMedium=mitr->second;
        }

        Point3D pt;
        for (int z = 1; z < workFieldDim.z-1; z++)
                for (int y = 1; y < workFieldDim.y-1; y++)
                        for (int x = 1; x < workFieldDim.x-1; x++){
                                pt=Point3D(x-1,y-1,z-1);
                                //             cerr<<"pt="<<pt<<" is valid "<<cellFieldG->isValid(pt)<<endl;
                                currentCellPtr=cellFieldG->get(pt);
                                //             currentCellPtr=cellFieldG->get(pt);
                                //             cerr<<"THIS IS PTR="<<currentCellPtr<<endl;

                                //             if(currentCellPtr)
                                //                cerr<<"This is id="<<currentCellPtr->id<<endl;
                                //currentConcentration = concentrationArray[x][y][z];

                                if(secreteInMedium && ! currentCellPtr){
                                        concentrationArray[x][y][z]=secrConstMedium;
                                }

                                if(currentCellPtr){
                                        mitr=secrData.typeIdSecrConstMap.find(currentCellPtr->type);
                                        if(mitr!=end_mitr){
                                                secrConst=mitr->second;
                                                concentrationArray[x][y][z]=secrConst;
                                        }
                                }
                        }


}



void FlexibleDiffusionSolverFE::secrete() {

        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){

                for(unsigned int j = 0 ; j <diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec.size() ; ++j){
                        (this->*diffSecrFieldTuppleVec[i].secrData.secretionFcnPtrVec[j])(i);

                        //          (this->*secrDataVec[i].secretionFcnPtrVec[j])(i);
                }


        }



}



float FlexibleDiffusionSolverFE::couplingTerm(Point3D & _pt,std::vector<CouplingData> & _couplDataVec,float _currentConcentration){

        float couplingTerm=0.0;
        float coupledConcentration;
        for(int i =  0 ; i < _couplDataVec.size() ; ++i){
                coupledConcentration=concentrationFieldVector[_couplDataVec[i].fieldIdx]->get(_pt);
                couplingTerm+=_couplDataVec[i].couplingCoef*_currentConcentration*coupledConcentration;
        }

        return couplingTerm;



}
void FlexibleDiffusionSolverFE::boundaryConditionInit(ConcentrationField_t *concentrationField){
        Array3D_t & concentrationArray=concentrationField->getContainer();

        //dealing with periodic boundary condition in x direction
        if(periodicBoundaryCheckVector[0]){//periodic boundary conditions were set in x direction
                int x=0;
                for(int y=0 ; y<workFieldDim.y ; ++y)
                        for(int z=0 ; z<workFieldDim.z ; ++z){
                                concentrationArray[x][y][z]=concentrationArray[workFieldDim.x-2][y][z];
                        }
                        x=workFieldDim.x-1;
                        for(int y=0 ; y<workFieldDim.y ; ++y)
                                for(int z=0 ; z<workFieldDim.z ; ++z){
                                        concentrationArray[x][y][z]=concentrationArray[1][y][z];
                                }
        }else{//noFlux BC
                int x=0;
                for(int y=0 ; y<workFieldDim.y ; ++y)
                        for(int z=0 ; z<workFieldDim.z ; ++z){
                                concentrationArray[x][y][z]=concentrationArray[x+1][y][z];
                        }
                        x=workFieldDim.x-1;
                        for(int y=0 ; y<workFieldDim.y ; ++y)
                                for(int z=0 ; z<workFieldDim.z ; ++z){
                                        concentrationArray[x][y][z]=concentrationArray[x-1][y][z];
                                }

        }

        //dealing with periodic boundary condition in y direction
        if(periodicBoundaryCheckVector[1]){//periodic boundary conditions were set in x direction
                int y=0;
                for(int x=0 ; x<workFieldDim.x ; ++x)
                        for(int z=0 ; z<workFieldDim.z ; ++z){
                                concentrationArray[x][y][z]=concentrationArray[x][workFieldDim.y-2][z];
                        }
                        y=workFieldDim.y-1;
                        for(int x=0 ; x<workFieldDim.x ; ++x)
                                for(int z=0 ; z<workFieldDim.z ; ++z){
                                        concentrationArray[x][y][z]=concentrationArray[x][1][z];
                                }
        }else{//NoFlux BC
                int y=0;
                for(int x=0 ; x<workFieldDim.x ; ++x)
                        for(int z=0 ; z<workFieldDim.z ; ++z){
                                concentrationArray[x][y][z]=concentrationArray[x][y+1][z];
                        }
                        y=workFieldDim.y-1;
                        for(int x=0 ; x<workFieldDim.x ; ++x)
                                for(int z=0 ; z<workFieldDim.z ; ++z){
                                        concentrationArray[x][y][z]=concentrationArray[x][y-1][z];
                                }
        }

        //dealing with periodic boundary condition in z direction
        if(periodicBoundaryCheckVector[2]){//periodic boundary conditions were set in x direction
                int z=0;
                for(int x=0 ; x<workFieldDim.x ; ++x)
                        for(int y=0 ; y<workFieldDim.y ; ++y){
                                concentrationArray[x][y][z]=concentrationArray[x][y][workFieldDim.z-2];
                        }
                        z=workFieldDim.z-1;
                        for(int x=0 ; x<workFieldDim.x ; ++x)
                                for(int y=0 ; y<workFieldDim.y ; ++y){
                                        concentrationArray[x][y][z]=concentrationArray[x][y][1];
                                }
        }else{//Noflux BC
                int z=0;
                for(int x=0 ; x<workFieldDim.x ; ++x)
                        for(int y=0 ; y<workFieldDim.y ; ++y){
                                concentrationArray[x][y][z]=concentrationArray[x][y][z+1];
                        }
                        z=workFieldDim.z-1;
                        for(int x=0 ; x<workFieldDim.x ; ++x)
                                for(int y=0 ; y<workFieldDim.y ; ++y){
                                        concentrationArray[x][y][z]=concentrationArray[x][y][z-1];
                                }

        }




}


void FlexibleDiffusionSolverFE::diffuseSingleField(unsigned int idx){





        Point3D pt, n;
        unsigned int token = 0;
        double distance;
        CellG * currentCellPtr=0,*nCell=0;

        short currentCellType=0;
        float concentrationSum=0.0;
        float updatedConcentration=0.0;


        float currentConcentration=0.0;
        short neighborCounter=0;

        DiffusionData & diffData = diffSecrFieldTuppleVec[idx].diffData;
        float diffConst=diffData.diffConst;
        float decayConst=diffData.decayConst;
        float deltaT=diffData.deltaT;
        float deltaX=diffData.deltaX;
        float dt_dx2=deltaT/(deltaX*deltaX);



        std::set<unsigned char>::iterator sitr;
        std::set<unsigned char>::iterator end_sitr=diffData.avoidTypeIdSet.end();
        std::set<unsigned char>::iterator end_sitr_decay=diffData.avoidDecayInIdSet.end();

        Automaton *automaton=potts->getAutomaton();

        ConcentrationField_t * concentrationFieldPtr=concentrationFieldVector[idx];
        ConcentrationField_t * scratchFieldPtr;

        if(!haveCouplingTerms)
                scratchFieldPtr=concentrationFieldVector[diffSecrFieldTuppleVec.size()];
        else
                scratchFieldPtr=concentrationFieldVector[diffSecrFieldTuppleVec.size()+idx];

        Array3D_t & concentrationArray = concentrationFieldPtr->getContainer();
        Array3D_t & scratchArray = scratchFieldPtr->getContainer();

        boundaryConditionInit(concentrationFieldPtr);//initializing boundary conditions

        bool avoidMedium=false;
        bool avoidDecayInMedium=false;
        //the assumption is that medium has type ID 0
        if(diffData.avoidTypeIdSet.find(automaton->getTypeId("Medium")) != end_sitr){
                avoidMedium=true;
        }

        if(diffData.avoidDecayInIdSet.find(automaton->getTypeId("Medium")) != end_sitr_decay){
                avoidDecayInMedium=true;
        }

        unsigned x_min=1,x_max=workFieldDim.x-1;
        unsigned y_min=1,y_max=workFieldDim.y-1;
        unsigned z_min=1,z_max=workFieldDim.z-1;

        if(diffData.useBoxWatcher){
                Point3D minCoordinates=*(boxWatcherSteppable->getMinCoordinatesPtr());
                Point3D maxCoordinates=*(boxWatcherSteppable->getMaxCoordinatesPtr());
                cerr<<"FLEXIBLE DIFF SOLVER maxCoordinates="<<maxCoordinates<<" minCoordinates="<<minCoordinates<<endl;
                x_min=minCoordinates.x+1;
                x_max=maxCoordinates.x+1;
                y_min=minCoordinates.y+1;
                y_max=maxCoordinates.y+1;
                z_min=minCoordinates.z+1;
                z_max=maxCoordinates.z+1;
        }
        cerr<<"x_min="<<x_min<<" x_max="<<x_max<<"y_min="<<y_min<<" y_max="<<y_max<<"z_min="<<z_min<<" z_max="<<z_max<<endl;

        for (int z = z_min; z < z_max; z++)
                for (int y = y_min; y < y_max; y++)
                        for (int x = x_min; x < x_max; x++){
                                //for (int z = 1; z < workFieldDim.z-1; z++)
                                // for (int y = 1; y < workFieldDim.y-1; y++)
                                //   for (int x = 1; x < workFieldDim.x-1; x++){
                                pt=Point3D(x-1,y-1,z-1);
                                currentCellPtr=cellFieldG->get(pt);
                                //          currentCellPtr=cellFieldG->get(pt);

                                currentConcentration = concentrationArray[x][y][z];


                                if(avoidMedium && !currentCellPtr){//if medium is to be avoided
                                        if(avoidDecayInMedium){
                                                scratchArray[x][y][z]=currentConcentration;
                                        }
                                        else{
                                                scratchArray[x][y][z]=currentConcentration-deltaT*(decayConst*currentConcentration);
                                        }
                                        continue;
                                }

                                if(currentCellPtr && diffData.avoidTypeIdSet.find(currentCellPtr->type)!=end_sitr){
                                        if(diffData.avoidDecayInIdSet.find(currentCellPtr->type)!=end_sitr_decay){
                                                scratchArray[x][y][z]=currentConcentration;
                                        }else{
                                                scratchArray[x][y][z]=currentConcentration-deltaT*(decayConst*currentConcentration);
                                        }
                                        continue; // avoid user defined types
                                }

                                updatedConcentration=0.0;
                                concentrationSum=0.0;
                                neighborCounter=0;

                                //loop over nearest neighbors
                                CellG *neighborCellPtr=0;
                                const std::vector<Point3D> & offsetVecRef=boundaryStrategy->getOffsetVec(pt);
                                //          cerr<<"maxNeighborIndex="<<maxNeighborIndex<<endl;
                                for (int i = 0  ; i<=maxNeighborIndex /*offsetVec.size()*/ ; ++i ){ 
                                        const Point3D & offset = offsetVecRef[i];


                                        if(diffData.avoidTypeIdSet.size()||avoidMedium){ //means user defined types to avoid in terms of the diffusion

                                                n=pt+offsetVecRef[i];
                                                neighborCellPtr=cellFieldG->get(n);
                                                if(avoidMedium && !neighborCellPtr) continue; // avoid medium if specified by the user
                                                if(neighborCellPtr && diffData.avoidTypeIdSet.find(neighborCellPtr->type)!=end_sitr) continue;//avoid user specified types
                                        }
                                        concentrationSum += concentrationArray[x+offset.x][y+offset.y][z+offset.z];

                                        ++neighborCounter;

                                }


                                updatedConcentration =  dt_dx2*diffConst*(concentrationSum - neighborCounter*currentConcentration)+currentConcentration;





                                //processing decay depandent on type of the current cell
                                if(currentCellPtr){
                                        if(diffData.avoidDecayInIdSet.find(currentCellPtr->type)!=end_sitr_decay){
                                                ;//decay in this type is forbidden
                                        }else{
                                                updatedConcentration-=deltaT*(decayConst*currentConcentration);//decay in this type is allowed
                                        }
                                }else{
                                        if(avoidDecayInMedium){
                                                ;//decay in Medium is forbidden
                                        }else{
                                                updatedConcentration-=deltaT*(decayConst*currentConcentration); //decay in Medium is allowed
                                        }
                                }
                                if(haveCouplingTerms){
                                        updatedConcentration+=couplingTerm(pt,diffData.couplingDataVec,currentConcentration);
                                }

                                //imposing artificial limits on allowed concentration
                                if(diffData.useThresholds){
                                        if(updatedConcentration>diffData.maxConcentration){
                                                updatedConcentration=diffData.maxConcentration;
                                        }
                                        if(updatedConcentration<diffData.minConcentration){
                                                updatedConcentration=diffData.minConcentration;
                                        }
                                }

                                scratchArray[x][y][z]=updatedConcentration;//updating scratch
                        }

}

bool FlexibleDiffusionSolverFE::isBoudaryRegion(int x, int y, int z, Dim3D dim)
{
        if (x < 2 || x > dim.x - 3 || y < 2 || y > dim.y - 3 || z < 2 || z > dim.z - 3 )
                return true;
        else
                return false;
}

void FlexibleDiffusionSolverFE::diffuse() {

        for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){
                diffuseSingleField(i);
                if(!haveCouplingTerms){ //without coupling terms field swap takes place immediately aftera given field has been diffused
                        ConcentrationField_t * concentrationField=concentrationFieldVector[i];
                        ConcentrationField_t * scratchField=concentrationFieldVector[diffSecrFieldTuppleVec.size()];
                        //copy updated values from scratch to concentration fiel
                        scrarch2Concentration(scratchField,concentrationField);
                }

        }

        if(haveCouplingTerms){  //with coupling terms we swap scratch and concentration field after all fields have been diffused
                for(unsigned int i = 0 ; i < diffSecrFieldTuppleVec.size() ; ++i ){
                        ConcentrationField_t * concentrationField=concentrationFieldVector[i];
                        ConcentrationField_t * scratchField=concentrationFieldVector[diffSecrFieldTuppleVec.size()+i];
                        scrarch2Concentration(scratchField,concentrationField);

                }

        }


}
void FlexibleDiffusionSolverFE::scrarch2Concentration( ConcentrationField_t *scratchField, ConcentrationField_t *concentrationField){
        scratchField->switchContainersQuick(*(concentrationField));

}

void FlexibleDiffusionSolverFE::outputField( std::ostream & _out, ConcentrationField_t *_concentrationField){
        Point3D pt;
        float tempValue;

        for (pt.z = 0; pt.z <